% modal.m
%
% A Matlab script to demonstrate modal synthesis.
%
% by Gary P. Scavone
% MUMT 307, McGill University

fs = 44100;
Ts = 1 / fs;
N = 100000;    % number of samples to compute

% Mode parameters
nModes = 4;
f = [ 440 720 1200 2000 ];             % modal center frequencies
radii = [ 0.999 0.9999 0.9995 0.9995 ]; % modal radii
amps = [ 40 60 50 70 ];                % modal amplitudes

% Calculate the resonance filter coefficients
as = [ ones(nModes, 1) (-2*radii.*cos( 2*pi*f*Ts))' radii.^2' ];

% Choose from several input signal options
sOption = 1;
if sOption == 1
  % Generate a single impulse signal
  x = [1, zeros(1, N-1)];
elseif sOption == 2
    % Generate an impulse train
    x = zeros(1, N);
    x(1:400:N) = 1;
else
  % Generate a noise signal and perhaps pre-filter it
  x = 2 * rand( 1, N ) - 1;
  [b, a] = butter( 2, [0.25 0.4] ); % 2nd-order passband
  x = filter( b, a, x );
end

% Filter signal with parallel filter bank
y = zeros(size(x));
for n = 1:nModes
  y = y + amps(n) * filter( 1, as(n,:), x );
end

soundsc( y, fs )