% Driver for Chapter 1: an exact string, and the sound of sampling it.
%
% The editable script supplies
%
%     y = pluck(t, f0, beta, where, xp, fmax, T60)
%
% the modal sum for a plucked string, restricted to partials below fmax.
%
% The driver's own ground truth is d'Alembert's solution, which shares no
% algebra with the modal sum: the string's shape is the average of two copies
% of the initial triangle travelling in opposite directions, so the bridge
% force is a square-cornered pulse train and the pickup displacement is
% piecewise linear in time.  Both are evaluated exactly, in closed form, at
% any time, and neither is band-limited.
%
% mode 'pluck'  one note: the first three periods drawn finely, the harmonic
%               amplitudes, and the sound with a listening loss applied.
% mode 'alias'  the same note sampled two ways: d'Alembert's signal read off
%               at the sampling instants (every partial, most of them above
%               the Nyquist frequency and folded back below it), and the
%               modal sum stopped at the Nyquist frequency.  Their difference
%               is exactly the folded part, and its energy is also known in
%               closed form, which the driver prints alongside the measurement.
p = wv_params();
fs = p.fs;
f0 = p.f0;
beta = p.beta;
out = struct();
out.fs = fs;
out.f0 = f0;
out.nyq = floor(fs / 2 / f0);

if strcmp(p.mode, 'pluck')
    % three periods, finely, for the waveform panel
    tf = linspace(0, 3 / f0, 900)';
    ts = (0:floor(3 * fs / f0))' / fs;
    out.tf = tf' * 1000;
    out.exact = drv_dalembert(tf, f0, beta, p.where, p.xp)';
    out.ts = ts' * 1000;
    out.synth = pluck(ts, f0, beta, p.where, p.xp, fs / 2, Inf)';

    % harmonic amplitudes, from the exact coefficients
    nh = min(out.nyq, 60);
    n = 1:nh;
    c = drv_coef(n, beta, p.where, p.xp);
    out.n = n;
    out.cdb = drv_db(c, max(abs(c)));
    out.missing = sum(abs(sin(n * pi * beta)) < 1e-9);

    % the sound: band-limited, with the listening loss
    t = (0:round(p.T * fs) - 1)' / fs;
    y = pluck(t, f0, beta, p.where, p.xp, fs / 2, p.T60);
    [fr, db] = wv_spectrum(y(1:min(end, round(0.5 * fs))), fs, p.fmax, 700);
    out.fr = fr;
    out.db = db;
    out.clip = wv_audio(y, fs);

    % how far the band-limited sum is from d'Alembert at the sample instants
    yb = pluck(t(1:round(0.25 * fs)), f0, beta, p.where, p.xp, fs / 2, Inf);
    ye = drv_dalembert(t(1:round(0.25 * fs)), f0, beta, p.where, p.xp);
    out.reldiff = sqrt(sum((yb - ye) .^ 2) / sum(yb .^ 2));
else
    t = (0:round(p.T * fs) - 1)' / fs;
    env = exp(-6.908 * t / p.T60);          % the same envelope on both clips
    ye = drv_dalembert(t, f0, beta, p.where, p.xp);
    yb = pluck(t, f0, beta, p.where, p.xp, fs / 2, Inf);
    d = ye - yb;
    out.measured = 10 * log10(sum(d .^ 2) / sum(yb .^ 2));
    out.predicted = drv_predicted(f0, fs, beta, p.where, p.xp);
    nw = round(0.5 * fs);
    [fr, dbe] = wv_spectrum(ye(1:nw), fs, fs / 2, 900);
    [~, dbb] = wv_spectrum(yb(1:nw), fs, fs / 2, 900);
    out.fr = fr;
    out.dbnaive = dbe;
    out.dbband = dbb;
    % where the first few folded partials land
    nmax = 3 * out.nyq;
    nn = out.nyq + 1:nmax;
    fa = mod(nn * f0, fs);
    fa(fa > fs / 2) = fs - fa(fa > fs / 2);
    cc = drv_coef(nn, beta, p.where, p.xp);
    keep = abs(cc) > 1e-9 * max(abs(drv_coef(1:out.nyq, beta, p.where, p.xp)));
    out.falias = fa(keep);
    out.harm = (1:out.nyq) * f0;
    sc = max(abs(ye));
    out.clipNaive = wv_audio(ye .* env, fs, sc);
    out.clipBand = wv_audio(yb .* env, fs, sc);
    out.clipDiff = wv_audio(d .* env, fs, sc);
    out.ratio = fs / f0;
end

wv_write(out);

function y = drv_dalembert(t, f0, beta, where, xp)
% d'Alembert: u(x, t) = (phi(x - ct) + phi(x + ct)) / 2, with phi the odd,
% period-2 extension of the initial triangle and c = 2 f0 on a string of
% length 1.
s = 2 * f0 * t;
if strcmp(where, 'bridge')
    % u_x(0, t) = phi'(ct), and phi' is even: 1/beta on (0, beta),
    % -1/(1 - beta) on (beta, 2 - beta), 1/beta again up to 2.  A sample that
    % lands exactly on a jump takes the midpoint, as a Fourier series would.
    r = mod(s, 2);
    hi = 1 / beta;
    lo = -1 / (1 - beta);
    y = lo * ones(size(r));
    y(r < beta | r > 2 - beta) = hi;
    y(abs(r - beta) < 1e-12 | abs(r - (2 - beta)) < 1e-12) = 0.5 * (hi + lo);
else
    y = 0.5 * (drv_phi(xp - s, beta) + drv_phi(xp + s, beta));
end
end

function v = drv_phi(x, beta)
% the initial triangle, extended oddly about 0 and periodically with period 2
r = mod(x, 2);
sg = ones(size(r));
sg(r > 1) = -1;
r(r > 1) = 2 - r(r > 1);
v = min(r / beta, (1 - r) / (1 - beta));
v = sg .* v;
end

function c = drv_coef(n, beta, where, xp)
A = 2 * sin(n * pi * beta) ./ (n .^ 2 * pi ^ 2 * beta * (1 - beta));
if strcmp(where, 'bridge')
    c = n * pi .* A;
else
    c = A .* sin(n * pi * xp);
end
end

function d = drv_db(c, ref)
d = 20 * log10(max(abs(c) / ref, 1e-8));
end

function r = drv_predicted(f0, fs, beta, where, xp)
% Energy above the Nyquist frequency over energy below it, from the exact
% coefficients.  For the bridge force the total is known in closed form by
% Parseval, mean(F^2) = 1 / (beta (1 - beta)), so the tail is the total less
% the head; for the pickup the coefficients fall like 1/n^2 and the tail is
% summed directly to n = 2e5, beyond which it is below rounding.
K = floor(fs / 2 / f0);
head = drv_coef(1:K, beta, where, xp);
eh = sum(head .^ 2) / 2;
if strcmp(where, 'bridge')
    et = 1 / (beta * (1 - beta)) - eh;
else
    tail = drv_coef(K + 1:200000, beta, where, xp);
    et = sum(tail .^ 2) / 2;
end
r = 10 * log10(et / eh);
end