wolffd@0: function fcoefs=MakeERBFilters(fs,numChannels,lowFreq)
wolffd@0: % function [fcoefs]=MakeERBFilters(fs,numChannels,lowFreq)
wolffd@0: % This function computes the filter coefficients for a bank of 
wolffd@0: % Gammatone filters.  These filters were defined by Patterson and 
wolffd@0: % Holdworth for simulating the cochlea.  
wolffd@0: % 
wolffd@0: % The result is returned as an array of filter coefficients.  Each row 
wolffd@0: % of the filter arrays contains the coefficients for four second order 
wolffd@0: % filters.  The transfer function for these four filters share the same
wolffd@0: % denominator (poles) but have different numerators (zeros).  All of these
wolffd@0: % coefficients are assembled into one vector that the ERBFilterBank 
wolffd@0: % can take apart to implement the filter.
wolffd@0: %
wolffd@0: % The filter bank contains "numChannels" channels that extend from
wolffd@0: % half the sampling rate (fs) to "lowFreq".  Alternatively, if the numChannels
wolffd@0: % input argument is a vector, then the values of this vector are taken to
wolffd@0: % be the center frequency of each desired filter.  (The lowFreq argument is
wolffd@0: % ignored in this case.)
wolffd@0: 
wolffd@0: % Note this implementation fixes a problem in the original code by
wolffd@0: % computing four separate second order filters.  This avoids a big
wolffd@0: % problem with round off errors in cases of very small cfs (100Hz) and
wolffd@0: % large sample rates (44kHz).  The problem is caused by roundoff error
wolffd@0: % when a number of poles are combined, all very close to the unit
wolffd@0: % circle.  Small errors in the eigth order coefficient, are multiplied
wolffd@0: % when the eigth root is taken to give the pole location.  These small
wolffd@0: % errors lead to poles outside the unit circle and instability.  Thanks
wolffd@0: % to Julius Smith for leading me to the proper explanation.
wolffd@0: 
wolffd@0: % Execute the following code to evaluate the frequency
wolffd@0: % response of a 10 channel filterbank.
wolffd@0: %	fcoefs = MakeERBFilters(16000,10,100);
wolffd@0: %	y = ERBFilterBank([1 zeros(1,511)], fcoefs);
wolffd@0: %	resp = 20*log10(abs(fft(y')));
wolffd@0: %	freqScale = (0:511)/512*16000;
wolffd@0: %	semilogx(freqScale(1:255),resp(1:255,:));
wolffd@0: %	axis([100 16000 -60 0])
wolffd@0: %	xlabel('Frequency (Hz)'); ylabel('Filter Response (dB)');
wolffd@0: 
wolffd@0: % Rewritten by Malcolm Slaney@Interval.  June 11, 1998.
wolffd@0: % (c) 1998 Interval Research Corporation  
wolffd@0: 
wolffd@0: T = 1/fs;
wolffd@0: if length(numChannels) == 1
wolffd@0: 	cf = ERBSpace(lowFreq, fs/2, numChannels);
wolffd@0: else
wolffd@0: 	cf = numChannels(1:end);
wolffd@0: 	if size(cf,2) > size(cf,1)
wolffd@0: 		cf = cf';
wolffd@0: 	end
wolffd@0: end
wolffd@0: 
wolffd@0: % Change the followFreqing three parameters if you wish to use a different
wolffd@0: % ERB scale.  Must change in ERBSpace too.
wolffd@0: EarQ = 9.26449;				%  Glasberg and Moore Parameters
wolffd@0: minBW = 24.7;
wolffd@0: order = 1;
wolffd@0: 
wolffd@0: ERB = ((cf/EarQ).^order + minBW^order).^(1/order);
wolffd@0: B=1.019*2*pi*ERB;
wolffd@0: 
wolffd@0: A0 = T;
wolffd@0: A2 = 0;
wolffd@0: B0 = 1;
wolffd@0: B1 = -2*cos(2*cf*pi*T)./exp(B*T);
wolffd@0: B2 = exp(-2*B*T);
wolffd@0: 
wolffd@0: A11 = -(2*T*cos(2*cf*pi*T)./exp(B*T) + 2*sqrt(3+2^1.5)*T*sin(2*cf*pi*T)./ ...
wolffd@0: 		exp(B*T))/2;
wolffd@0: A12 = -(2*T*cos(2*cf*pi*T)./exp(B*T) - 2*sqrt(3+2^1.5)*T*sin(2*cf*pi*T)./ ...
wolffd@0: 		exp(B*T))/2;
wolffd@0: A13 = -(2*T*cos(2*cf*pi*T)./exp(B*T) + 2*sqrt(3-2^1.5)*T*sin(2*cf*pi*T)./ ...
wolffd@0: 		exp(B*T))/2;
wolffd@0: A14 = -(2*T*cos(2*cf*pi*T)./exp(B*T) - 2*sqrt(3-2^1.5)*T*sin(2*cf*pi*T)./ ...
wolffd@0: 		exp(B*T))/2;
wolffd@0: 
wolffd@0: gain = abs((-2*exp(4*i*cf*pi*T)*T + ...
wolffd@0:                  2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
wolffd@0:                          (cos(2*cf*pi*T) - sqrt(3 - 2^(3/2))* ...
wolffd@0:                           sin(2*cf*pi*T))) .* ...
wolffd@0:            (-2*exp(4*i*cf*pi*T)*T + ...
wolffd@0:              2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
wolffd@0:               (cos(2*cf*pi*T) + sqrt(3 - 2^(3/2)) * ...
wolffd@0:                sin(2*cf*pi*T))).* ...
wolffd@0:            (-2*exp(4*i*cf*pi*T)*T + ...
wolffd@0:              2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
wolffd@0:               (cos(2*cf*pi*T) - ...
wolffd@0:                sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) .* ...
wolffd@0:            (-2*exp(4*i*cf*pi*T)*T + 2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
wolffd@0:            (cos(2*cf*pi*T) + sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) ./ ...
wolffd@0:           (-2 ./ exp(2*B*T) - 2*exp(4*i*cf*pi*T) +  ...
wolffd@0:            2*(1 + exp(4*i*cf*pi*T))./exp(B*T)).^4);
wolffd@0: 	
wolffd@0: allfilts = ones(length(cf),1);
wolffd@0: fcoefs = [A0*allfilts A11 A12 A13 A14 A2*allfilts B0*allfilts B1 B2 gain];
wolffd@0: 
wolffd@0: if 0						% Test Code
wolffd@0:     figure
wolffd@0: 	A0  = fcoefs(:,1);
wolffd@0: 	A11 = fcoefs(:,2);
wolffd@0: 	A12 = fcoefs(:,3);
wolffd@0: 	A13 = fcoefs(:,4);
wolffd@0: 	A14 = fcoefs(:,5);
wolffd@0: 	A2  = fcoefs(:,6);
wolffd@0: 	B0  = fcoefs(:,7);
wolffd@0: 	B1  = fcoefs(:,8);
wolffd@0: 	B2  = fcoefs(:,9);
wolffd@0: 	gain= fcoefs(:,10);	
wolffd@0: 	chan=1;
wolffd@0: 	x = [1 zeros(1, 511)];
wolffd@0: 	y1=filter([A0(chan)/gain(chan) A11(chan)/gain(chan) ...
wolffd@0: 		A2(chan)/gain(chan)],[B0(chan) B1(chan) B2(chan)], x);
wolffd@0: 	y2=filter([A0(chan) A12(chan) A2(chan)], ...
wolffd@0: 			[B0(chan) B1(chan) B2(chan)], y1);
wolffd@0: 	y3=filter([A0(chan) A13(chan) A2(chan)], ...
wolffd@0: 			[B0(chan) B1(chan) B2(chan)], y2);
wolffd@0: 	y4=filter([A0(chan) A14(chan) A2(chan)], ...
wolffd@0: 			[B0(chan) B1(chan) B2(chan)], y3);
wolffd@0: 	semilogx((0:(length(x)-1))*(fs/length(x)),20*log10(abs(fft(y4))));
wolffd@0: end