|
wolffd@0
|
1 function fcoefs=MakeERBFilters(fs,numChannels,lowFreq)
|
|
wolffd@0
|
2 % function [fcoefs]=MakeERBFilters(fs,numChannels,lowFreq)
|
|
wolffd@0
|
3 % This function computes the filter coefficients for a bank of
|
|
wolffd@0
|
4 % Gammatone filters. These filters were defined by Patterson and
|
|
wolffd@0
|
5 % Holdworth for simulating the cochlea.
|
|
wolffd@0
|
6 %
|
|
wolffd@0
|
7 % The result is returned as an array of filter coefficients. Each row
|
|
wolffd@0
|
8 % of the filter arrays contains the coefficients for four second order
|
|
wolffd@0
|
9 % filters. The transfer function for these four filters share the same
|
|
wolffd@0
|
10 % denominator (poles) but have different numerators (zeros). All of these
|
|
wolffd@0
|
11 % coefficients are assembled into one vector that the ERBFilterBank
|
|
wolffd@0
|
12 % can take apart to implement the filter.
|
|
wolffd@0
|
13 %
|
|
wolffd@0
|
14 % The filter bank contains "numChannels" channels that extend from
|
|
wolffd@0
|
15 % half the sampling rate (fs) to "lowFreq". Alternatively, if the numChannels
|
|
wolffd@0
|
16 % input argument is a vector, then the values of this vector are taken to
|
|
wolffd@0
|
17 % be the center frequency of each desired filter. (The lowFreq argument is
|
|
wolffd@0
|
18 % ignored in this case.)
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 % Note this implementation fixes a problem in the original code by
|
|
wolffd@0
|
21 % computing four separate second order filters. This avoids a big
|
|
wolffd@0
|
22 % problem with round off errors in cases of very small cfs (100Hz) and
|
|
wolffd@0
|
23 % large sample rates (44kHz). The problem is caused by roundoff error
|
|
wolffd@0
|
24 % when a number of poles are combined, all very close to the unit
|
|
wolffd@0
|
25 % circle. Small errors in the eigth order coefficient, are multiplied
|
|
wolffd@0
|
26 % when the eigth root is taken to give the pole location. These small
|
|
wolffd@0
|
27 % errors lead to poles outside the unit circle and instability. Thanks
|
|
wolffd@0
|
28 % to Julius Smith for leading me to the proper explanation.
|
|
wolffd@0
|
29
|
|
wolffd@0
|
30 % Execute the following code to evaluate the frequency
|
|
wolffd@0
|
31 % response of a 10 channel filterbank.
|
|
wolffd@0
|
32 % fcoefs = MakeERBFilters(16000,10,100);
|
|
wolffd@0
|
33 % y = ERBFilterBank([1 zeros(1,511)], fcoefs);
|
|
wolffd@0
|
34 % resp = 20*log10(abs(fft(y')));
|
|
wolffd@0
|
35 % freqScale = (0:511)/512*16000;
|
|
wolffd@0
|
36 % semilogx(freqScale(1:255),resp(1:255,:));
|
|
wolffd@0
|
37 % axis([100 16000 -60 0])
|
|
wolffd@0
|
38 % xlabel('Frequency (Hz)'); ylabel('Filter Response (dB)');
|
|
wolffd@0
|
39
|
|
wolffd@0
|
40 % Rewritten by Malcolm Slaney@Interval. June 11, 1998.
|
|
wolffd@0
|
41 % (c) 1998 Interval Research Corporation
|
|
wolffd@0
|
42
|
|
wolffd@0
|
43 T = 1/fs;
|
|
wolffd@0
|
44 if length(numChannels) == 1
|
|
wolffd@0
|
45 cf = ERBSpace(lowFreq, fs/2, numChannels);
|
|
wolffd@0
|
46 else
|
|
wolffd@0
|
47 cf = numChannels(1:end);
|
|
wolffd@0
|
48 if size(cf,2) > size(cf,1)
|
|
wolffd@0
|
49 cf = cf';
|
|
wolffd@0
|
50 end
|
|
wolffd@0
|
51 end
|
|
wolffd@0
|
52
|
|
wolffd@0
|
53 % Change the followFreqing three parameters if you wish to use a different
|
|
wolffd@0
|
54 % ERB scale. Must change in ERBSpace too.
|
|
wolffd@0
|
55 EarQ = 9.26449; % Glasberg and Moore Parameters
|
|
wolffd@0
|
56 minBW = 24.7;
|
|
wolffd@0
|
57 order = 1;
|
|
wolffd@0
|
58
|
|
wolffd@0
|
59 ERB = ((cf/EarQ).^order + minBW^order).^(1/order);
|
|
wolffd@0
|
60 B=1.019*2*pi*ERB;
|
|
wolffd@0
|
61
|
|
wolffd@0
|
62 A0 = T;
|
|
wolffd@0
|
63 A2 = 0;
|
|
wolffd@0
|
64 B0 = 1;
|
|
wolffd@0
|
65 B1 = -2*cos(2*cf*pi*T)./exp(B*T);
|
|
wolffd@0
|
66 B2 = exp(-2*B*T);
|
|
wolffd@0
|
67
|
|
wolffd@0
|
68 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
|
69 exp(B*T))/2;
|
|
wolffd@0
|
70 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
|
71 exp(B*T))/2;
|
|
wolffd@0
|
72 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
|
73 exp(B*T))/2;
|
|
wolffd@0
|
74 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
|
75 exp(B*T))/2;
|
|
wolffd@0
|
76
|
|
wolffd@0
|
77 gain = abs((-2*exp(4*i*cf*pi*T)*T + ...
|
|
wolffd@0
|
78 2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
|
|
wolffd@0
|
79 (cos(2*cf*pi*T) - sqrt(3 - 2^(3/2))* ...
|
|
wolffd@0
|
80 sin(2*cf*pi*T))) .* ...
|
|
wolffd@0
|
81 (-2*exp(4*i*cf*pi*T)*T + ...
|
|
wolffd@0
|
82 2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
|
|
wolffd@0
|
83 (cos(2*cf*pi*T) + sqrt(3 - 2^(3/2)) * ...
|
|
wolffd@0
|
84 sin(2*cf*pi*T))).* ...
|
|
wolffd@0
|
85 (-2*exp(4*i*cf*pi*T)*T + ...
|
|
wolffd@0
|
86 2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
|
|
wolffd@0
|
87 (cos(2*cf*pi*T) - ...
|
|
wolffd@0
|
88 sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) .* ...
|
|
wolffd@0
|
89 (-2*exp(4*i*cf*pi*T)*T + 2*exp(-(B*T) + 2*i*cf*pi*T).*T.* ...
|
|
wolffd@0
|
90 (cos(2*cf*pi*T) + sqrt(3 + 2^(3/2))*sin(2*cf*pi*T))) ./ ...
|
|
wolffd@0
|
91 (-2 ./ exp(2*B*T) - 2*exp(4*i*cf*pi*T) + ...
|
|
wolffd@0
|
92 2*(1 + exp(4*i*cf*pi*T))./exp(B*T)).^4);
|
|
wolffd@0
|
93
|
|
wolffd@0
|
94 allfilts = ones(length(cf),1);
|
|
wolffd@0
|
95 fcoefs = [A0*allfilts A11 A12 A13 A14 A2*allfilts B0*allfilts B1 B2 gain];
|
|
wolffd@0
|
96
|
|
wolffd@0
|
97 if (0) % Test Code
|
|
wolffd@0
|
98 A0 = fcoefs(:,1);
|
|
wolffd@0
|
99 A11 = fcoefs(:,2);
|
|
wolffd@0
|
100 A12 = fcoefs(:,3);
|
|
wolffd@0
|
101 A13 = fcoefs(:,4);
|
|
wolffd@0
|
102 A14 = fcoefs(:,5);
|
|
wolffd@0
|
103 A2 = fcoefs(:,6);
|
|
wolffd@0
|
104 B0 = fcoefs(:,7);
|
|
wolffd@0
|
105 B1 = fcoefs(:,8);
|
|
wolffd@0
|
106 B2 = fcoefs(:,9);
|
|
wolffd@0
|
107 gain= fcoefs(:,10);
|
|
wolffd@0
|
108 chan=1;
|
|
wolffd@0
|
109 x = [1 zeros(1, 511)];
|
|
wolffd@0
|
110 y1=filter([A0(chan)/gain(chan) A11(chan)/gain(chan) ...
|
|
wolffd@0
|
111 A2(chan)/gain(chan)],[B0(chan) B1(chan) B2(chan)], x);
|
|
wolffd@0
|
112 y2=filter([A0(chan) A12(chan) A2(chan)], ...
|
|
wolffd@0
|
113 [B0(chan) B1(chan) B2(chan)], y1);
|
|
wolffd@0
|
114 y3=filter([A0(chan) A13(chan) A2(chan)], ...
|
|
wolffd@0
|
115 [B0(chan) B1(chan) B2(chan)], y2);
|
|
wolffd@0
|
116 y4=filter([A0(chan) A14(chan) A2(chan)], ...
|
|
wolffd@0
|
117 [B0(chan) B1(chan) B2(chan)], y3);
|
|
wolffd@0
|
118 semilogx((0:(length(x)-1))*(fs/length(x)),20*log10(abs(fft(y4))));
|
|
wolffd@0
|
119 end
|
