% Driver for Chapter 3: the digital waveguide and the Karplus-Strong string.
%
% The editable script supplies
%
%     y = ks(x0, D, tuning, M)
%
% a delay loop of length D samples, started from x0.
%
% mode 'equiv'  the lossless loop against the leapfrog scheme of Chapter 2 at
%               lam = 1, on grids from 25 to 3200 intervals.  The leapfrog
%               scheme is run by the driver's own copy; its first 2N output
%               samples seed the loop, and the two are then compared for the
%               rest of the run, sample for sample, and timed.
% mode 'keys'   one note on each of 25 keys from C2 to C8, with each tuning.
%               The fundamental's frequency is measured by wv_peakfreq and
%               compared with the equal-tempered target; its decay time is
%               measured from the envelope.  A two-octave scale is also played
%               with each tuning, and once by exact modal synthesis at the
%               target pitches, for the ear to compare.
p = wv_params();
fs = p.fs;
out = struct();

if strcmp(p.mode, 'equiv')
    Ns = p.Ns(:)';
    M = round(p.T * fs);
    TL = NaN(1, numel(Ns));
    TW = NaN(1, numel(Ns));
    DIFF = NaN(1, numel(Ns));
    for k = 1:numel(Ns)
        N = Ns(k);
        x = (0:N)' / N;
        u0 = min(x / p.beta, (1 - x) / (1 - p.beta));
        um = 0.5 * ([u0(2:end); 0] + [0; u0(1:end-1)]);   % exact at lam = 1
        um(1) = 0;
        um(end) = 0;
        drv_lf(u0, um, 50);
        t0 = tic;
        yl = drv_lf(u0, um, M);
        TL(k) = toc(t0) * 1000;
        seed = yl(1:2 * N);
        ks(seed, 2 * N, 'lossless', 4 * N);
        t0 = tic;
        yw = ks(seed, 2 * N, 'lossless', M);
        TW(k) = toc(t0) * 1000;
        DIFF(k) = max(abs(yw - yl)) / max(abs(yl));
    end
    out.Ns = Ns;
    out.tl = TL;
    out.tw = TW;
    out.diff = DIFF;
    out.M = M;
else
    notes = p.notes(:)';
    tunings = {'integer', 'linear', 'allpass'};
    nt = numel(tunings);
    Mn = round(p.Tn * fs);
    CE = NaN(nt, numel(notes));
    T60 = NaN(nt, numel(notes));
    for i = 1:numel(notes)
        f = wv_note(notes(i));
        D = fs / f;
        x0 = drv_seed(f, fs, p.beta);
        for k = 1:nt
            y = ks(x0, D, tunings{k}, Mn);
            fm = wv_peakfreq(y, 1 / fs, f / 1.1, f * 1.1);
            CE(k, i) = wv_cents(fm, f);
            T60(k, i) = drv_t60(y, fm, fs);
        end
    end
    out.notes = notes;
    out.f = wv_note(notes);
    out.tunings = tunings;
    out.cents = wv_rows(CE);
    out.t60 = wv_rows(T60);

    % a scale, played with each tuning and exactly
    sc = p.scale(:)';
    step = round(p.step * fs);
    ring = round(p.ring * fs);
    Mt = step * (numel(sc) - 1) + ring;
    Y = zeros(Mt, nt + 1);
    for i = 1:numel(sc)
        f = wv_note(sc(i));
        x0 = drv_seed(f, fs, p.beta);
        idx = (i - 1) * step + (1:ring);
        for k = 1:nt
            Y(idx, k) = Y(idx, k) + ks(x0, fs / f, tunings{k}, ring);
        end
        % the reference: the same pluck, modally, at exactly the target pitch,
        % with each partial given the decay the integer loop gives it
        t = (0:ring - 1)' / fs;
        Y(idx, nt + 1) = Y(idx, nt + 1) + drv_modal(t, f, fs, p.beta);
    end
    s = max(abs(Y(:)));
    out.clipInt = wv_audio(Y(:, 1), fs, s);
    out.clipLin = wv_audio(Y(:, 2), fs, s);
    out.clipAp = wv_audio(Y(:, 3), fs, s);
    out.clipExact = wv_audio(Y(:, 4), fs, s);
    out.scale = sc;
end

wv_write(out);

function y = drv_lf(u, uold, M)
% leapfrog at lam = 1 on a fixed string, recording the first node off the end
n1 = numel(u);
y = zeros(M, 1);
unew = zeros(n1, 1);
for n = 1:M
    for j = 2:n1 - 1
        unew(j) = u(j + 1) + u(j - 1) - uold(j);
    end
    uold = u;
    u = unew;
    y(n) = u(2);
end
end

function x0 = drv_seed(f, fs, beta)
% The loop's initial content, the same for every tuning: the band-limited
% bridge force of a string plucked at beta, for a little over one period,
% because the loop reads up to two samples further back than its length.
P = ceil(fs / f) + 2;
t = (0:P - 1)' / fs;
x0 = zeros(P, 1);
for n = 1:floor(fs / 2 / f)
    x0 = x0 + 2 * sin(n * pi * beta) / (n * pi * beta * (1 - beta)) * cos(2 * pi * n * f * t);
end
% A loop with a two-point average in it passes a constant unchanged for
% ever, so a seed with a nonzero mean leaves a DC offset that never decays.
% It is inaudible, but it would sit under the fundamental's decay and spoil
% the measurement of it, so the seed is made to have none.
x0 = x0 - mean(x0);
end

function y = drv_modal(t, f, fs, beta)
% the same pluck by modal synthesis, each partial decaying as the two-point
% average makes it decay in a loop of exactly the right length
y = zeros(size(t));
for n = 1:floor(fs / 2 / f)
    g = abs(cos(pi * n * f / fs));
    s = -log(max(g, 1e-300)) * f;       % amplitude lost per second
    y = y + 2 * sin(n * pi * beta) / (n * pi * beta * (1 - beta)) * exp(-s * t) .* cos(2 * pi * n * f * t);
end
end

function T = drv_t60(y, f, fs)
% decay time of the partial at f: its amplitude in successive tapered frames,
% a straight line fitted to the logarithm, extrapolated to -60 dB
nw = round(8 * fs / f);
nw = max(nw, 256);
hop = nw;
nfr = floor((numel(y) - nw) / hop);
w = wv_window(nw);
tt = (0:nw - 1)' / fs;
e = exp(-2i * pi * f * tt);
a = zeros(nfr, 1);
tc = zeros(nfr, 1);
for k = 1:nfr
    seg = y((k - 1) * hop + (1:nw));
    a(k) = abs(sum(w .* seg .* e));
    tc(k) = ((k - 1) * hop + nw / 2) / fs;
end
% only the frames within 100 dB of the loudest, where the partial is still
% far above rounding error and the window's leakage from anything else
use = find(a > 1e-5 * max(a));
use = use(use > 1);
if numel(use) < 3
    use = 2:min(4, nfr);
end
c = polyfit(tc(use), log(a(use)), 1);
T = -log(1000) / c(1);
if c(1) >= 0
    T = Inf;
end
end