% pluck.m
%
% This Matlab script creates a simple waveguide plucked
% string model and uses commuted synthesis to implement
% the body filter.  The pluck position of the string is
% not accounted for here.  This code uses concepts covered
% by several patents held by Stanford University.
%
% by Gary P. Scavone
% CCRMA, Stanford University, 1998-2001.
% Mcgill University, 2004.

fprintf('\nSelect string initialization:\n\n');
fprintf('  1. Impulse\n');
fprintf('  2. Noise\n');
fprintf('  3. Body/Air response (test.wav)\n\n');
type = input('Enter option number [1]: ');
if isempty(type),
  type = 1;
end

% ******* Constants and Other Parameters ******* %

fs = 22050;   % sampling rate
N = 20000;    % length of vector to compute

% ******* String dimensions (m) ******* %

f0 = [196, 294, 440, 660];  % violin open-string fundamentals
D = round(0.5*fs/f0(2));    % length of string in samples

% ******* Simple String Attenuation Filter ******* %

b_loss = -0.99*[0.5 0.5];
a_loss = 1;
z_loss = 0;

% ******* Initialize delay lines ******* %

y = zeros(1,N);                  % initialize output vector
dlines = zeros(2,N);

if type == 1,
  x = [1, zeros(1, N-1)];
elseif type == 2,
    x = zeros(1, N);
    dlines = rand(2, N);
else
  x = audioread('test.wav');         % recorded body-to-air response
  x = [x', zeros(1,N-length(x))];  % zero remaining length
end

ptr = 1;
readloc = 0.1;                   % 0 = left-end, 1 = right-end
upreadptr = round((1-readloc)*(D-1))+1;
dwnreadptr = D-upreadptr+1;

% *************************************** %
%                                         %
%             Delay Line(d)               %
%    |-------------------------------|    %
%    ^                                    %
%  pointer                                %
%                                         %
%     >>--- pointer increments --->>      %
%                                         %
% *************************************** %
%
% The pointer initially points to the delay line output.
% We can take the output and calculate a new input value
% which is placed where the output was taken from.  The
% pointer is then incremented and the process repeated.

% ******* Run Loop Start ******* %

for n = 1:N,

  temp = dlines(2,ptr);
  y(n) = dlines(1,upreadptr);
  y(n) = y(n) + dlines(2,dwnreadptr);

  dlines(2,ptr) = -dlines(1,ptr);

  [dlines(1,ptr),z_loss] = filter(b_loss, a_loss, temp, z_loss);
  dlines(1,ptr) = dlines(1,ptr) + x(n);
  
  % ****** Increment Pointers & Check Limits ****** %

  ptr = ptr + 1;
  upreadptr = upreadptr + 1;
  dwnreadptr = dwnreadptr + 1;
  if ptr > D
    ptr = 1;
  end
  if upreadptr > D
    upreadptr = 1;
  end
  if dwnreadptr > D
    dwnreadptr = 1;
  end
end

% Scale soundfile if necessary
if max(abs(y)) > 0.95
  y = y./(max(abs(y))+0.1);
  disp('Scaled waveform');
end

% Clear Figure Window and Plot
clf
plot(y);
sound(y',fs);