% Finite difference string
% Linear string, gaussian pluck at T=0
% AND bow model from Bilabo
% dv = (u1(x_mid) - u2(x_mid)) - bow_speed  ;
% bow_fn = bow_amp * sign(dv) * (eps + exp(-(abs(dv))/bow_width)); 
% Bruce Land, Cornell University, Jan 2017
% second order scheme from
% http://arxiv.org/PS_cache/physics/pdf/0009/0009068v2.pdf,
% page 14, eqn 2.18

clear all

%linear dimension of membrane -- bigger is lower pitch
%n=50, rho=0.05 111 Hz is close to A2 (110)
n = 50 ;
u = zeros(n); %time t
u1 = zeros(n); %time t-1
u2 = zeros(n); %time t-2
uHit = zeros(n); %input strike
uDist = zeros(n); %input distirbuted
bow_fn = zeros(n); 

%%%% time %%%%
Fs = 48000 ;
time = 0.0:1/Fs:4.0;
%output sound
aud = zeros(size(time)) ;
% output drive wave
drive = zeros(size(time)) ;
% output string velocity at the bow point
u1_minus_u2 = zeros(size(time));
%index for storing outputs
tindex = 1;

%%%% string parameters %%%%
% 0 < rho < 0.5 -- lower rho => lower pitch
% rho = (vel*dt/dx)^2
rho = 0.05 ; % 0.0485
% estimate the frequency
Fs*sqrt(rho)/(2*(n-1))
% eta = damping*dt/2
% higher damping => shorter sound
eta = 0.00005 ;
one_eta = 1 - eta ;

%%%% pluck parameters %%%%
%sets the amplitude of the stick strike
% a very small pluck helps to start the bow
ampIn = .001; %.3;
%sets the position of the stick strike: 0<pos<n
x_mid = fix(n/4); %/4
%sets width of the gaussian strike input -- see figure 1
alpha = .5 ;
%compute the gaussian strike amplitude
for i=2:n-1
    uHit(i) = ampIn*exp(-alpha*(((i-1)-x_mid)^2));
end
%enforce boundary conditions
uHit(1) = 0;
uHit(n) = 0;

%%%% bow parameters %%%%
% with n=50, rho=0.0485 (note A2)
% speed   width  amp   sound
% 0.00005 0.0005  0.0035 good rise, good sustain, nice modulation <***
% 0.00005 0.0005  0.0033 good rise, good sustain, nice modulation <***
% with n=50, rho=0.0485*4 (note A3)
% speed   width  amp   sound
% 0.00005 0.0005  0.0021 good rise, good sustain, nice modulation <***
% 0.00005 0.0005  0.0015 good rise, good sustain, nice modulation <***
bow_speed =  0.00005 ; %0.0001
bow_width =  0.0005 ; % 0.001 ;  0.005 interesting modulation
eps = 0.02 ;
current_bow_amp = 0;
input_amp = .002 ; % 0.01
bow_rise = 20000;
bow_rise_inc = input_amp/bow_rise ;
bow_fall = 1000;
bow_fall_inc = input_amp/bow_fall ;
bow_sustain = 30000 ;

%%%% string location to use as sound putput %%%%
% pick a location about 0.75 of the length of the string
% any other sample or sum of points is also valid to try
sampleLoc = fix(n*3/4);

for t=time
    % audio output
    aud(tindex) = u1(sampleLoc); ; %sum(sum(u))/(n^2); %
    % drive (pluck or bow location) output
    drive(tindex) = bow_fn(x_mid);
    % speed of string at bow location output
    u1_minus_u2(tindex) = u1(x_mid) - u2(x_mid) ;
    
    % calculate the bow envelope
    if (tindex < (bow_rise + bow_fall + bow_sustain))
        if(tindex <= bow_rise)
            current_bow_amp = current_bow_amp + bow_rise_inc ;
        elseif (tindex <= bow_rise + bow_sustain)
            current_bow_amp = current_bow_amp;
        else
            current_bow_amp = current_bow_amp - bow_fall_inc ;
        end
    else
        current_bow_amp = 0;
    end
    bow_amp = current_bow_amp ;
    
    % now the bow model
    % get the relative speed
    dv =((u1(x_mid) - u2(x_mid)) - bow_speed ) ;
    
    % bow_fn(x_mid) = bow_amp * dv * exp(-(dv*dv)/bow_width); 
    % the actual bow dynamics
    bow_fn(x_mid) = bow_amp * sign(dv) * (eps + exp(-(abs(dv))/bow_width)); 
    % reduce spatial aliasing
    bow_fn(x_mid-1) =  bow_fn(x_mid)/2 ;
    bow_fn(x_mid+1) =  bow_fn(x_mid)/2 ;
    bow_fn(x_mid-2) =  bow_fn(x_mid)/4 ;
    bow_fn(x_mid+2) =  bow_fn(x_mid)/4 ;
%     bow_fn(x_mid-3) = bow_fn(x_mid)/8 ;
%     bow_fn(x_mid+3) = bow_fn(x_mid)/8 ;
%     bow_fn(x_mid-4) = bow_fn(x_mid)/16 ;
%     bow_fn(x_mid+4) = bow_fn(x_mid)/16 ;
%     
    %vectorized update of positions, except for boundaries
    u(2:n-1) = one_eta * (...
        rho * (u1(3:n)+u1(1:n-2)-2*u1(2:n-1)) ...
        + 2 * u1(2:n-1) ...
        - one_eta * (u2(2:n-1)) ...
        )  ;
    
    %enforce boundary conditions
    u(1) = 0 ;
    u(n) = 0 ;
    
    %update history
    u2 = u1;
    % update and add bow energy
    u1 = u + bow_fn ;
    
    %apply drum strikes 
    if t==0 ; %| t==1.0
        u1 = u1 + uHit;
    end
   
    tindex = tindex + 1;
end

figure(1); clf;

subplot(2,2,1)
plot(time, aud);
hold on
plot(time, drive*50);
title('sound out and bow drive (red)')

subplot(2,2,2)
plot(time, aud);
hold on
plot(time, drive*50);
plot(time, u1_minus_u2*20);
title('expanded sound out, bow(r), rel speed (y)')
set(gca, 'ylim', [-1.1 1.1]);
set(gca, 'xlim', [0.5 0.55]);

subplot(2,2,3)
fft_size = 2^14;
spectrogram(aud,fft_size,0,fft_size, Fs,'yaxis');
title('spectrogram')
set(gca, 'ylim', [0 1.0]);

subplot(2,2,4)
Hs=spectrum.welch('Hamming',2^16,50);
psd(Hs,aud,'Fs',Fs)
title('String spectrum')
set(gca,'xlim',[0 0.5])

max(aud)
mean(abs(aud))

% get the spectrum

% set(gca,'ylim',[0 10])
%play the sound


soundsc(aud,Fs)

%end
%============================================================