% Modal synthesis of a plucked string.
%
% An ideal string of length 1, fixed at both ends and released from rest in
% the shape of a triangle with its apex at x = beta, moves as
%
%   u(x, t) = sum_n  A_n sin(n pi x) cos(2 pi n f0 t),
%   A_n = 2 sin(n pi beta) / (n^2 pi^2 beta (1 - beta)),
%
% which is exact: it is the Fourier sine series of the initial shape, and each
% term is a standing wave that solves the wave equation on its own.  What the
% listener hears is not u itself but something it drives, and this script
% offers two choices:
%
%   'bridge'   the force on the bridge, proportional to the slope u_x(0, t),
%              whose partials fall off like 1/n.  This is what an acoustic
%              instrument's body is shaken by, and it sounds bright.
%   'pickup'   the displacement at x = xp, whose partials fall off like 1/n^2
%              and are notched wherever sin(n pi xp) vanishes.  This is what a
%              magnetic pickup reports, and it sounds darker.
%
% Only partials below fmax are summed.  With fmax at the Nyquist frequency the
% result is the band-limited signal, which is the right thing to play; see the
% driver for what happens otherwise.
%
% T60 sets a loss for listening: each partial decays by 60 dB in
% T60 / (1 + (f/2000)^2) seconds, so that high partials die first, as they do
% on a real string.  With T60 = Inf the synthesis is lossless and is compared
% exactly against d'Alembert's solution.

function y = pluck(t, f0, beta, where, xp, fmax, T60)
t = t(:);
y = zeros(size(t));
nmax = floor(fmax / f0);
for n = 1:nmax
    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
    f = n * f0;
    if isinf(T60)
        y = y + c * cos(2 * pi * f * t);
    else
        s = 6.908 * (1 + (f / 2000) ^ 2) / T60;
        y = y + c * exp(-s * t) .* cos(2 * pi * f * t);
    end
end
end