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