1% The Karplus-Strong plucked string: a delay loop with a filter in it.
2%
3% Chapter 2 showed that at lam = 1 the leapfrog scheme is exact, and that the
4% signal at the first node off the bridge repeats itself every 2N steps. A
5% sequence that repeats itself needs no grid: it needs a buffer of one period
6% and a rule for what comes next,
7%
8% y(n) = y(n - D)
9%
10% which is the lossless digital waveguide, at a cost of one sample per sample
11% instead of N. Karplus and Strong's observation (1983) was that putting a
12% two-point average inside the loop,
13%
14% y(n) = ( y(n - P) + y(n - P - 1) ) / 2,
15%
16% makes each partial lose a little on every trip round, the high ones more
17% than the low ones, which is what a real string does. The average also
18% delays the loop by half a sample, so the loop length is P + 1/2.
19%
20% The pitch is fs / (loop length), and a loop made of whole samples can only
21% be so long. The `tuning` argument chooses how the missing fraction eta of a
22% sample is supplied:
23%
24% 'lossless' y(n) = y(n - D), D a whole number: the waveguide itself
25% 'integer' round the loop to the nearest whole length and accept the error
26% 'linear' interpolate linearly between two taps, a delay of eta that is
27% also a low-pass filter
28% 'allpass' a first-order allpass, coefficient (1 - eta)/(1 + eta), whose
29% delay is eta at low frequency and whose gain is exactly one
30%
31% x0 is the loop's initial content, played out first; D is the desired loop
32% length in samples, fs / f0; M is the number of samples to return.
34function y = ks(x0, D, tuning, M)
35y = zeros(M, 1);
36n0 = numel(x0);
37y(1:n0) = x0(:);
38if strcmp(tuning, 'lossless')
39 P = round(D);
40 for n = n0 + 1:M
41 y(n) = y(n - P);
42 end
43elseif strcmp(tuning, 'integer')
44 P = round(D - 0.5);
45 for n = n0 + 1:M
46 y(n) = 0.5 * (y(n - P) + y(n - P - 1));
47 end
48elseif strcmp(tuning, 'linear')
49 P = floor(D - 0.5);
50 eta = D - 0.5 - P;
51 for n = n0 + 1:M
52 y(n) = 0.5 * ((1 - eta) * y(n - P) + y(n - P - 1) + eta * y(n - P - 2));
53 end
54else
55 % keep eta in [0.1, 1.1): an allpass with eta near zero has its pole near
56 % the unit circle and takes a long time to settle
57 P = floor(D - 0.6);
58 eta = D - 0.5 - P;
59 C = (1 - eta) / (1 + eta);
60 wprev = 0;
61 for n = n0 + 1:M
62 w = 0.5 * (y(n - P) + y(n - P - 1));
63 y(n) = C * w + wprev - C * y(n - 1);
64 wprev = w;
65 end
66end
67end