GAMMATONE - Gammatone filter coefficients

Program code:

function [b,a]=gammatone(fc,fs,n,betamul);

%   AUTHOR : Stephan Ewert, Peter L. Soendergaard

% ------ Checking of input parameters ---------


error(nargchk(2,4,nargin));

if ~isnumeric(fs) || ~isscalar(fs) || fs<=0
  error('%s: fs must be a positive scalar.',upper(mfilename));
end;

if ~isnumeric(fc) || ~isvector(fc) || any(fc<0) || any(fc>fs/2)
  error(['%s: fc must be a vector of positive values that are less than half ' ...
         'the sampling rate.'],upper(mfilename));
end;

if nargin==4
  if ~isnumeric(betamul) || ~isscalar(betamul) || betamul<=0
    error('%s: beta must be a positive scalar.',upper(mfilename));
  end;
end;

if nargin>2
  if ~isnumeric(n) || ~isscalar(n) || n<=0 || fix(n)~=n
    error('%s: n must be a positive, integer scalar.',upper(mfilename));
  end;
end;


% ------ Computation --------------------------

if nargin==2
  % Choose a 4th order filter.
  n=4;
end;

if nargin<4
  % Determine the correct multiplier for beta depending on the filter
  % order.
  switch(n)
   case 2
    betamul =  0.637;
   case 4
    betamul = 1.0183;
   otherwise
    error(['GAMMATONE: Default value for beta can only be computed for 2nd ' ...
           'and 4th order filters.']);
  end;
end;

nchannels = length(fc);

b=zeros(nchannels,1);
a=zeros(nchannels,n+1);

% ourbeta is used in order not to mask the beta function.

ourbeta = betamul*audfiltbw(fc);

for ii = 1:nchannels

  btmp=1-exp(-2*pi*ourbeta(ii)/fs);
  atmp=[1, -exp(-(2*pi*ourbeta(ii) + i*2*pi*fc(ii))/fs)];

  b2=1;
  a2=1;

  for jj=1:n
    b2=conv(btmp,b2);
    a2=conv(atmp,a2);
  end

  % Place the result (a row vector) in the output matrices.
  b(ii,:)=b2;
  a(ii,:)=a2;

end;