% The Karplus-Strong plucked string: a delay loop with a filter in it.
%
% Chapter 2 showed that at lam = 1 the leapfrog scheme is exact, and that the
% signal at the first node off the bridge repeats itself every 2N steps.  A
% sequence that repeats itself needs no grid: it needs a buffer of one period
% and a rule for what comes next,
%
%   y(n) = y(n - D)
%
% which is the lossless digital waveguide, at a cost of one sample per sample
% instead of N.  Karplus and Strong's observation (1983) was that putting a
% two-point average inside the loop,
%
%   y(n) = ( y(n - P) + y(n - P - 1) ) / 2,
%
% makes each partial lose a little on every trip round, the high ones more
% than the low ones, which is what a real string does.  The average also
% delays the loop by half a sample, so the loop length is P + 1/2.
%
% The pitch is fs / (loop length), and a loop made of whole samples can only
% be so long.  The `tuning` argument chooses how the missing fraction eta of a
% sample is supplied:
%
%   'lossless'  y(n) = y(n - D), D a whole number: the waveguide itself
%   'integer'   round the loop to the nearest whole length and accept the error
%   'linear'    interpolate linearly between two taps, a delay of eta that is
%               also a low-pass filter
%   'allpass'   a first-order allpass, coefficient (1 - eta)/(1 + eta), whose
%               delay is eta at low frequency and whose gain is exactly one
%
% x0 is the loop's initial content, played out first; D is the desired loop
% length in samples, fs / f0; M is the number of samples to return.

function y = ks(x0, D, tuning, M)
y = zeros(M, 1);
n0 = numel(x0);
y(1:n0) = x0(:);
if strcmp(tuning, 'lossless')
    P = round(D);
    for n = n0 + 1:M
        y(n) = y(n - P);
    end
elseif strcmp(tuning, 'integer')
    P = round(D - 0.5);
    for n = n0 + 1:M
        y(n) = 0.5 * (y(n - P) + y(n - P - 1));
    end
elseif strcmp(tuning, 'linear')
    P = floor(D - 0.5);
    eta = D - 0.5 - P;
    for n = n0 + 1:M
        y(n) = 0.5 * ((1 - eta) * y(n - P) + y(n - P - 1) + eta * y(n - P - 2));
    end
else
    % keep eta in [0.1, 1.1): an allpass with eta near zero has its pole near
    % the unit circle and takes a long time to settle
    P = floor(D - 0.6);
    eta = D - 0.5 - P;
    C = (1 - eta) / (1 + eta);
    wprev = 0;
    for n = n0 + 1:M
        w = 0.5 * (y(n - P) + y(n - P - 1));
        y(n) = C * w + wprev - C * y(n - 1);
        wprev = w;
    end
end
end