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