% Maraca PhISEM model.
% Uses a periodic "shake" or re-energizing of the system K samples.
% By Gary Scavone, McGill University, based on model by Perry Cook.

fs = 44100;  % sample rate
T = 1 / fs;  % sample period
N = 1 * fs;  % number of samples to compute

nBeans = 64;         % number of beans
shellF = 3200.0;     % gourd resonant frequency
shellR = 0.96;       % gourd resonant radius
systemDecay = 0.999; 
soundDecay = 0.95;
K = 3000;            % number of samples between shakes

gain = log( nBeans ) / (log( 4.0 ) * 40.0 * nBeans);

% Resonant filter coefficients
b = [1 -1];
a = [1 -2 * shellR * cos( 2 * pi * shellF * T ) shellR^2];

sndLevel = 0;
z = [0 0];
kcount = 1;  % shake counter
e = 1;  % system energy
x = 0;  % collision sound state

y = zeros(1, N);
for n = 1:N
  
  e = e * systemDecay;

  % Random collision computation
  if ( rand(1) * 1024 < nBeans )
    % a collision occurred
    sndLevel = sndLevel + (gain * e);
    x = x + 1;
  end

  x = sndLevel * (2 * rand(1) - 1.0);
  %x = sndLevel * x;
  [y(n), z] = filter( b, a, x, z );
  
  sndLevel = sndLevel * soundDecay;
  kcount = kcount + 1;
  if ( kcount > K ), e = 1; kcount = 0; end
  
end

ynorm = 0.95 * y / max(abs(y));
sound( ynorm, fs );