[cf,fcoefs]=MakeERBFilters(fs,numChannels,lowFreq,highFreq)

function [cf,fcoefs]=MakeERBFilters(fs,numChannels,lowFreq,highFreq)
% function [fcoefs]=MakeERBFilters(fs,numChannels,lowFreq)
% This function computes the filter coefficients for a bank of
% Gammatone filters.  These filters were defined by Patterson and
% Holdworth for simulating the cochlea.
%
% The result is returned as an array of filter coefficients.  Each row
% of the filter arrays contains the coefficients for four second order
% filters.  The transfer function for these four filters share the same
% denominator (poles) but have different numerators (zeros).  All of these
% coefficients are assembled into one vector that the ERBFilterBank
% can take apart to implement the filter.
%
% The filter bank contains "numChannels" channels that extend from
% half the sampling rate (fs) to "lowFreq".  Alternatively, if the numChannels
% input argument is a vector, then the values of this vector are taken to
% be the center frequency of each desired filter.  (The lowFreq argument is
% ignored in this case.)

% Note this implementation fixes a problem in the original code by
% computing four separate second order filters.  This avoids a big
% problem with round off errors in cases of very small cfs (100Hz) and
% large sample rates (44kHz).  The problem is caused by roundoff error
% when a number of poles are combined, all very close to the unit
% circle.  Small errors in the eigth order coefficient, are multiplied
% when the eigth root is taken to give the pole location.  These small
% errors lead to poles outside the unit circle and instability.  Thanks
% to Julius Smith for leading me to the proper explanation.

% Execute the following code to evaluate the frequency
% response of a 10 channel filterbank.
% 	fcoefs = MakeERBFilters(16000,10,100);
% 	y = ERBFilterBank([1 zeros(1,511)], fcoefs);
% 	resp = 20*log10(abs(fft(y')));
% 	freqScale = (0:511)/512*16000;
%     figure();
% 	semilogx(freqScale(1:255),resp(1:255,:));
% 	axis([100 16000 -60 0])
% 	xlabel('Frequency (Hz)'); ylabel('Filter Response (dB)');

% Rewritten by Malcolm Slaney@Interval.  June 11, 1998.
% (c) 1998 Interval Research Corporation

T = 1/fs;
% Change the followFreqing three parameters if you wish to use a different
% ERB scale.  Must change in ERBSpace too.
EarQ = 9.26449;				%  Glasberg and Moore Parameters
minBW = 24.7;
order = 1;
if length(numChannels) == 1
    % 	cf = ERBSpace(lowFreq, fs/2, numChannels);
    cf = -(EarQ*minBW) + exp((1:numChannels)'*(-log(highFreq + EarQ*minBW) + ...
        log(lowFreq + EarQ*minBW))/numChannels) * (highFreq + EarQ*minBW);
else
    cf = numChannels(1:end);
    if size(cf,2) > size(cf,1)
        cf = cf';
    end
end


ERB = ((cf/EarQ).^order + minBW^order).^(1/order);
B=1.019*2*pi*ERB;

A0 = T;
A2 = 0;
B0 = 1;
B1 = -2*cos(2*cf*pi*T)./exp(B*T);
B2 = exp(-2*B*T);

A11 = -(2*T*cos(2*cf*pi*T)./exp(B*T) + 2*sqrt(3+2^1.5)*T*sin(2*cf*pi*T)./ ...
    exp(B*T))/2;
A12 = -(2*T*cos(2*cf*pi*T)./exp(B*T) - 2*sqrt(3+2^1.5)*T*sin(2*cf*pi*T)./ ...
    exp(B*T))/2;
A13 = -(2*T*cos(2*cf*pi*T)./exp(B*T) + 2*sqrt(3-2^1.5)*T*sin(2*cf*pi*T)./ ...
    exp(B*T))/2;
A14 = -(2*T*cos(2*cf*pi*T)./exp(B*T) - 2*sqrt(3-2^1.5)*T*sin(2*cf*pi*T)./ ...
    exp(B*T))/2;

gain = abs((-2*exp(4*i*cf*pi*T)*T + ...
    2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
    (cos(2*cf*pi*T) - sqrt(3 - 2^(3/2))* ...
    sin(2*cf*pi*T))) .* ...
    (-2*exp(4*i*cf*pi*T)*T + ...
    2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
    (cos(2*cf*pi*T) + sqrt(3 - 2^(3/2)) * ...
    sin(2*cf*pi*T))).* ...
    (-2*exp(4*i*cf*pi*T)*T + ...
    2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
    (cos(2*cf*pi*T) - ...
    sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) .* ...
    (-2*exp(4*i*cf*pi*T)*T + 2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
    (cos(2*cf*pi*T) + sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) ./ ...
    (-2 ./ exp(2*B*T) - 2*exp(4*i*cf*pi*T) +  ...
    2*(1 + exp(4*i*cf*pi*T))./exp(B*T)).^4);

allfilts = ones(length(cf),1);
fcoefs = [A0*allfilts A11 A12 A13 A14 A2*allfilts B0*allfilts B1 B2 gain];