% 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.
% 
% better lowpass loop filter
% nonlinearity in loop
% dispersion allpass in loop
% wavetable drive with filtering
%
% 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
% see: http://www.phy.mtu.edu/~suits/notefreqs.html
 clear all

% setting f1 to 1.58*string fundamental makes piano sound
% plucking makes plucked string
% sawtooth makes bowed string 
%=== parmaters ===========================================
% "string" parameters
p = 61 ; %; %length of delay line
tune = 0.12 ; % [-1/11, 2/3] maps to [0.2 to 1.2] sample delay
lp_alpha = 0.5 ;
decay_factor = 0.98 ;

% drive parameters
f1 = 130.8*1.57 ; amp1 = 0.1; fall_time1 = .5 ; attack_time1 = .2  ;
pluck = 0 ;
%=========================================================

% sample freq and time to simulate
Fs = 8000;
tmax = 3;
t = 0:1/Fs:tmax ;

a1 = amp1 * (t.^attack_time1) .* exp(-t/fall_time1) ;
drive_in = a1 .* sawtooth(2*pi*f1*t,.90) ;%cos(2*pi*f1*t) ;
%drive_in = a1 .* cos(2*pi*f1*t) ;%cos(2*pi*f1*t) ;
%drive_in = a1 .* square(2*pi*f1*t, 50) ;
%drive_in = a1 ;
%drive_in = 0.2*a1 .* randn(size(t)) + a1.* cos(2*pi*f1*t)  ;

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 ;
last_lowpass_out = 0;

% pluck the string
if pluck
    sr = (rand(1,p)-0.5);
    sr = sawtooth(2*pi*(1:p)/p, 0.15);
else
    sr = 0*(1:p);%(rand(1,p)-0.5); %initial shift register contents 0*(1:p);%
end
[b,a] = butter(1,0.05); %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)*(1-lp_alpha)+sr(ptrout)*lp_alpha)*decay_factor; 
    %lowpass_out = (last_lowpass_out*(1-lp_alpha)+sr(ptrout)*lp_alpha)*decay_factor; 
    last_lowpass_out = lowpass_out ;
    
    % all-pass to tune string and feed back to shift reg
   % sr(ptrin) = lowpass_out;
    tune_out = tune*(lowpass_out - last_tune_out) + last_tune_in ;
    % all-pass state vars
    last_tune_out = tune_out;
    last_tune_in = lowpass_out;
    
    % cubic nonlinearity
%     thresh = 2/3 ;
%     tune_out = tune_out - tune_out^3/3 ;
%     if tune_out > thresh
%         tune_out =2/3;
%     elseif tune_out < -thresh
%         tune_out = -2/3;
%     end
    
    % update shift reg and wrap pointers
    thresh = .1 ;
    sr(ptrin) = tune_out + drive_in(i) ;
    
    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])
