% 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 = 31 ; %; %length of delay line
tune = 0.8 ; % [-1/11, 2/3] maps to [0.2 to 1.2] sample delay

% sample freq and time to simulate
Fs = 8000;
tmax = 1;
t = linspace(0, tmax, Fs*tmax);
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

out = out + 0.0*randn(size(out))  ;
soundsc(out,Fs)
figure(1); clf
subplot(3,1,1)
plot(out)

subplot(3,1,2)
pwelch(out,[4096],[],[],Fs);
set(gca, 'xlim', [0 2])

subplot(3,1,3)

% do the crosscorrelation with each sine wave
% middle C 261.62 Hz
for i=260:0.2:263
    si = sin(2*pi*i*t);
    ci = cos(2*pi*i*t);
    % running totals to lowpass
    outXsi = cumsum(out.*si);
    outXci = cumsum(out.*ci);
    % magnitude to plot
    ssq = outXsi.^2 + outXci.^2;
    plot(t, ssq);
    hold on
    % print the total (average) of the ssq
    fprintf('%6.1f %8.0f \n',i, ssq(end) );
end



