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 |