
function Karplus(p1,tune)

% 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


%  C4   p=170
%  D4   p=151
%  E4   p=134
%  F4   p=127
%  G4   p=113
%  A4   p=101
%  B4   p=90
%  C5   p=85
%  D5   p=76
%  E5   p=67


%p1 = 170 ; %; %length of delay lin
offset = 5;
p2 = p1 + offset;
p3 = p1 - offset;
%tune = 0.775 ; % [-1/11, 2/3] maps to [0.2 to 1.2] sample delay

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

out = zeros(1,Fs*tmax);
ptrout1 = 2; %pointer to sr
lastptrout = 1 ;
ptrin1 = 1; %pointer to sr
ptrout2 = 2;
ptrin2 = 1; %pointer to sr
ptrout3 = 2;
ptrin3 = 1; %pointer to sr


last_tune_out1 = 0 ;
last_tune_in1 = 0 ;
last_tune_out2 = 0 ;
last_tune_in2 = 0 ;
last_tune_out3 = 0 ;
last_tune_in3 = 0 ;

% decay factor 0<factor<=0.5
factor1 = 0.495; 
factor2 = 0.465; 
factor3 = 0.445; 

% pluck the string
sr1 = zeros(1,p1);
sr2 = zeros(1,p2);
sr3 = zeros(1,p3);
pluckLen = 20;

sr1(1:pluckLen) = rand(1,pluckLen)-2.75; %initial shift register contents
[b,a] = butter(1,[0.01 0.02]); %modify pluck -- lower freq is softer
sr1 = filter(b,a,sr1);

sr2(1:pluckLen) = rand(1,pluckLen)-0.50; %initial shift register contents
sr2 = filter(b,a,sr2);

sr3(1:pluckLen) = rand(1,pluckLen)-0.25; %initial shift register contents
sr3 = filter(b,a,sr3);


for i=1:Fs*tmax
    out(i) = 0.8.*sr1(ptrin1) + 0.1.*sr2(ptrin2) + 0.1.*sr3(ptrin3); 
    % get low pass for string dynamics
    lowpass_out1 = (sr1(ptrin1)+sr1(ptrout1))*factor1; 
    lowpass_out2 = (sr2(ptrin2)+sr2(ptrout2))*factor2; 
    lowpass_out3 = (sr3(ptrin3)+sr3(ptrout3))*factor3; 
    
    % all-pass to tune string and feed back to shift reg
   % sr(ptrin) = lowpass_out;
    sr1(ptrin1) = tune*(lowpass_out1 - last_tune_out1) + last_tune_in1 ;
    sr2(ptrin2) = tune*(lowpass_out2 - last_tune_out2) + last_tune_in2 ;
    sr3(ptrin3) = tune*(lowpass_out3 - last_tune_out3) + last_tune_in3 ;
    
    % all-pass state vars
    last_tune_out1 = sr1(ptrin1);
    last_tune_in1 = lowpass_out1;
    
    last_tune_out2 = sr2(ptrin2);
    last_tune_in2 = lowpass_out2;
    
    last_tune_out3 = sr3(ptrin3);
    last_tune_in3 = lowpass_out3;
    
    %update and wrap pointers
    if (ptrin1==p1) ptrin1=1;
    else ptrin1=ptrin1+1;
    end
    if (ptrin2==p2) ptrin2=1;
    else ptrin2=ptrin2+1;
    end
    if (ptrin3==p3) ptrin3=1;
    else ptrin3=ptrin3+1;
    end
    %
    if (ptrout1==p1) ptrout1=1; 
    else ptrout1=ptrout1+1;
    end
    if (ptrout2==p2) ptrout2=1;
    else ptrout2=ptrout2+1;
    end
    if (ptrout3==p3) ptrout3=1;
    else ptrout3=ptrout3+1;
    end  
%      if (i == Fs/100) 
%          sr1 = rand(1,p1)-0.5; %initial shift register contents
%          sr1 = filter(b,a,sr1);
%      end
end
echoPt = 290;
out = out + .85.*[out(echoPt:end) zeros(1,echoPt-1)];
L = length(out);
A = linspace(0.2, 1, 0.1*(L));
D1 = linspace(1, 0.7, 0.1*(L));
D2 = linspace(0.7, 0.5,0.25*(L)); 
S = linspace(0.5, 0.45, 0.25*(L));
R = linspace(0.45, 0, 0.3*(L));
ADSR = [A D1 D2 S R]; 
out = ADSR .* out;

soundsc(out,Fs)

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

%figure(3)
Y = fft(fftOut);
Y = abs(Y);
%plot(Y);