% karplus.m
%
% This Matlab script implements the Karplus-Strong algorithm.
%
% by Gary P. Scavone
% Mcgill University, 2006.

clear all
close all

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

fs = 44100;   % sampling rate
N = 80000;    % length of vector to compute
D = 200;      % delay line (or wavetable) length

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

b = -0.99*[0.5 0.5];
z = 0;

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

y = zeros(1,N);                  % initialize output vector
dline = 2 * rand(1, D) - 1.0;
ptr = 1;


figure(1)
subplot(3,1,1)
plot(dline)
set(gca,'fontsize',18);

title('Original delayline');


subplot(3,1,2)
plot(dline)
set(gca,'fontsize',18);

title('Filter delayline step n');
loopsound(dline,fs,fs/D);

subplot(3,1,3)
plot(y);
title('Waveform Step n');
set(gca,'fontsize',18);

figure(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,

  
  y(n) = dline(ptr);
  [dline(ptr), z] = filter(b, 1, y(n), z);
  
  % Increment Pointers & Check Limits

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

 
  
  
  if  mod(n,2000) == 0
      
      subplot(3,1,2)
      plot(dline)
      str = sprintf('Filter delayline step %d',n);
      title(str);

      subplot(3,1,3)
      plot(y);
      str = sprintf('Waveform Step %d',n);
      title(str);
      figure(1);

  end

end

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

figure(2)
% Clear Figure Window and Plot
clf
plot(y);
title('Final Step');
set(gca,'fontsize',18);

playsound(y',fs);