% traveling.m
%
% This Matlab script creates a simple waveguide string model that can
% be used to view traveling-wave components as they move along the
% string.
%
% by Gary P. Scavone
% Music 150: Musical Acoustics
% CCRMA, Stanford University, 1998-2001.

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

fs = 44100;   % sampling rate
N = 300;      % number of iterations to display
T = 1 / fs;   % sampling period

fprintf('\nSelect waveform type:\n\n');
fprintf('  1. Pluck\n');
fprintf('  2. Sinusoid\n');
fprintf('  3. Impulse\n\n');
waveform = input('Enter option number [1]: ');
if isempty(waveform),
  waveform = 1;
end
if waveform < 1 | waveform > 3,
  waveform = 1;
end

if waveform == 1,
  pp = input('\nSet pluck point (0-1.0) [0.5]:\n\n');
  if isempty(pp),
    pp = 0.5;
  end
  if pp < 0.0 | pp > 1.0,
    pp = 0.5;  % pluck point between bridge (0) and nut (1.0)
  end
end

if waveform == 2,
  harmonic = input('\nSet harmonic number [1]:\n\n');
  if isempty(harmonic),
    harmonic = 1;
  end
  if harmonic < 1,
    harmonic = 1;
  end
end

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

strLength = 0.33;                     % real violin string length (meters)
c = [196, 294, 440, 660]*2*strLength; % velocities for 4 violin strings
strLength = 0.3;                      % synthesized string length (meters)
D = round(strLength*fs/c(1));         % length of string in samples
if ~mod(D,2),                         % make D odd, better visual
  D = D+1;
end

% ******* Simple String Dispersion Filter ******* %

b_disp = -0.97*[0.5 0.5];
a_disp = 1;
z_disp = 0;

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

% Offset string samples by 1/2 sample to avoid a double zero end condition.
% In this case, the end of the string is actually 1/2 sample to the right
% and left of the end sample values seen.  However, the displayed
% plots include the end points by incorporating an extra (repeated)
% sample on both the upper and lower delay lines.

if waveform == 1,    % triangular pluck
  ppsamp = round(pp*D);
  dlines = [[0.5:ppsamp-1+0.5]/(ppsamp-1), [(D-ppsamp-0.5):-1:0.5]/(D-ppsamp)];
  dlines = 0.45*[fliplr(dlines); dlines];

elseif waveform == 2, % sinusoidal
  dlines = sin(harmonic*pi*[0.5:D-0.5]/D);
  dlines = 0.45*[fliplr(dlines); dlines];

elseif waveform == 3, % pulse
  plen = round(0.2*D);
  if ~mod(plen,2),
    plen = plen+1;
  end
  dlines = [0.5*sin(pi*[0:plen-1]/(plen-1)), zeros(1,D-plen)];
  dlines = [fliplr(dlines); zeros(1,D)];
end

ptr = 1;

% ******* Plot initial string shape and traveling-wave components ******* %

x = [0:D]/(D);
upline = fliplr([dlines(1,ptr:D), dlines(1,1:ptr-1), -dlines(2,ptr)]);
dwnline = [dlines(2,ptr:D), dlines(2,1:ptr-1), -dlines(1,ptr)];
plot(x,upline+dwnline)
hold on
plot(x,upline,'r')
plot(x,dwnline,'g')
hold off
axis([0 1 -1 1])
grid
ylabel('Vertical Displacement')
xlabel('Normalized String Length')

disp('Paused ... hit any key');
pause

% *************************************** %
%                                         %
%             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 i = 1:N,

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

  % Uncomment next line for commuted losses.
  %[dlines(1,ptr),z_disp] = filter(b_disp, a_disp, temp, z_disp);

  % Uncomment next line for lossless wave propagation.
  dlines(1,ptr) = -temp;
  
  % ****** Increment Pointers & Check Limits ****** %

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

  % ****** Plot string shape and traveling-wave components  ****** %
  if mod(i,1) == 0, 
    upline = fliplr([dlines(1,ptr:D), dlines(1,1:ptr-1), -dlines(2,ptr)]);
    dwnline = [dlines(2,ptr:D), dlines(2,1:ptr-1), -dlines(1,ptr)];
    plot(x,upline+dwnline)
    hold on
    plot(x,upline,'r')
    plot(x,dwnline,'g')
    hold off
    axis([0 1 -1 1])
    grid
    ylabel('Vertical Displacement')
    xlabel('Normalized String Length')
    
    % Uncomment to trigger each tick with key hit
    pause
    % Uncomment to have a uniform pause between each tick
    %pause(0.05)
  end

end
