% To tune the string:
% change the length, p, with bigger being lower freq
% AND fine tune by changing the 'tune' parameter
% which sets the phase shift of an allpass filter.
% see:
% https://ccrma.stanford.edu/~jos/pasp/First_Order_Allpass_Interpolation.html

%   C2   65.4 p=122 tune=0.12
%   C3   130.8 p=61 tune=0.12
%   C4   262 Hz (middle C)   p=31 tune=0.8
%   D4   294
%	E4   330
%	F4   349
%	G4   392
%	A4   440  p=18 tune=0.15
%	B4   494
%	C5   523  p=15 tune=0.12


p = 122 ; %; %length of delay line
tune = 0.12 ; % [-1/11, 2/3] maps to [0.2 to 1.2] sample delay

% sample freq and time to simulate
Fs = 8000;
tmax = 3;

out = zeros(1,Fs*tmax);
ptrout = 2; %pointer to sr
lastptrout = 1 ;
ptrin = 1; %pointer to sr
last_tune_out = 0 ;
last_tune_in = 0 ;

% decay factor 0<factor<=0.5
factor = 0.5; 

% pluck the string
sr = rand(1,p)-0.5; %initial shift register contents
[b,a] = butter(1,0.105); %modify pluck -- lower freq is softer
sr = filter(b,a,sr);

for i=1:Fs*tmax
    out(i) = sr(ptrin) ; 
    % get low pass for string dynamics
    lowpass_out = (sr(ptrin)+sr(ptrout))*factor; 
    % all-pass to tune string and feed back to shift reg
   % sr(ptrin) = lowpass_out;
    sr(ptrin) = tune*(lowpass_out - last_tune_out) + last_tune_in ;
    % all-pass state vars
    last_tune_out = sr(ptrin);
    last_tune_in = lowpass_out;
    
    %update and wrap pointers
    if (ptrin==p) ptrin=1;
    else ptrin=ptrin+1;
    end
    %
    if (ptrout==p) ptrout=1; 
    else ptrout=ptrout+1;
    end  
end

soundsc(out,Fs)

figure(2)
plot(out)

figure(1); clf
pwelch(out,[],[],[],Fs)
set(gca, 'xlim', [0 .2])
