|
wolffd@0
|
1 % [filters, freqs] = DesignLyonFilters(fs,EarQ,StepFactor)
|
|
wolffd@0
|
2 % Design the cascade of second order filters and the front filters
|
|
wolffd@0
|
3 % (outer/middle and compensator) needed for Lyon's Passive Short Wave
|
|
wolffd@0
|
4 % (Second Order Sections) cochlear model. The variables used here come
|
|
wolffd@0
|
5 % from Apple ATG Technical Report #13 titled "Lyon's Cochlear Model".
|
|
wolffd@0
|
6 %
|
|
wolffd@0
|
7 % Most of the parameters are hardwired into this m-function. The user
|
|
wolffd@0
|
8 % settable parameters are the digital sampling rate, the basic Q of the
|
|
wolffd@0
|
9 % each stage (usually 8 or 4), and the spacing factor between channels
|
|
wolffd@0
|
10 % (usually .25 and .125, respectively.)
|
|
wolffd@0
|
11 %
|
|
wolffd@0
|
12 % The result is returned as rows of second order filters; three coefficients
|
|
wolffd@0
|
13 % for the numerator and two for the denomiator. The coefficients are
|
|
wolffd@0
|
14 % [A0 A1 A2 B1 B2]..........................(Malcolm 1 Feb. 1993)
|
|
wolffd@0
|
15 %
|
|
wolffd@0
|
16
|
|
wolffd@0
|
17 % (c) 1998 Interval Research Corporation
|
|
wolffd@0
|
18 function [filters, CenterFreqs, gains] = DesignLyonFilters(fs,EarQ,StepFactor)
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 if nargin < 2
|
|
wolffd@0
|
21 EarQ = 8;
|
|
wolffd@0
|
22 end
|
|
wolffd@0
|
23
|
|
wolffd@0
|
24 if nargin < 3
|
|
wolffd@0
|
25 StepFactor=EarQ/32;
|
|
wolffd@0
|
26 end
|
|
wolffd@0
|
27
|
|
wolffd@0
|
28 Eb = 1000.0;
|
|
wolffd@0
|
29 EarZeroOffset = 1.5;
|
|
wolffd@0
|
30 EarSharpness = 5.0;
|
|
wolffd@0
|
31 EarPremphCorner = 300;
|
|
wolffd@0
|
32
|
|
wolffd@0
|
33 % Find top frequency, allowing space for first cascade filter.
|
|
wolffd@0
|
34 topf = fs/2.0;
|
|
wolffd@0
|
35 topf = topf - (sqrt(topf^2+Eb^2)/EarQ*StepFactor*EarZeroOffset)+ ...
|
|
wolffd@0
|
36 sqrt(topf^2+Eb^2)/EarQ*StepFactor;
|
|
wolffd@0
|
37
|
|
wolffd@0
|
38 % Find place where CascadePoleQ < .5
|
|
wolffd@0
|
39 lowf = Eb/sqrt(4*EarQ^2-1);
|
|
wolffd@0
|
40 NumberOfChannels = floor((EarQ*(-log(lowf + sqrt(lowf^2 + Eb^2)) + ...
|
|
wolffd@0
|
41 log(topf + sqrt(Eb^2 + topf^2))))/StepFactor);
|
|
wolffd@0
|
42
|
|
wolffd@0
|
43 % Now make an array of CenterFreqs..... This expression was derived by
|
|
wolffd@0
|
44 % Mathematica by integrating 1/EarBandwidth(cf) and solving for f as a
|
|
wolffd@0
|
45 % function of channel number.
|
|
wolffd@0
|
46 cn = 1:NumberOfChannels;
|
|
wolffd@0
|
47 CenterFreqs = (-((exp((cn*StepFactor)/EarQ)*Eb^2)/ ...
|
|
wolffd@0
|
48 (topf + sqrt(Eb^2 + topf^2))) + ...
|
|
wolffd@0
|
49 (topf + sqrt(Eb^2 + topf^2))./exp((cn*StepFactor)/EarQ))/2;
|
|
wolffd@0
|
50
|
|
wolffd@0
|
51 % OK, now we can figure out the parameters of each stage filter.
|
|
wolffd@0
|
52 EarBandwidth = sqrt(CenterFreqs.^2+Eb^2)/EarQ;
|
|
wolffd@0
|
53 CascadeZeroCF = CenterFreqs + EarBandwidth * StepFactor * EarZeroOffset;
|
|
wolffd@0
|
54 CascadeZeroQ = EarSharpness*CascadeZeroCF./EarBandwidth;
|
|
wolffd@0
|
55 CascadePoleCF = CenterFreqs;
|
|
wolffd@0
|
56 CascadePoleQ = CenterFreqs./EarBandwidth;
|
|
wolffd@0
|
57
|
|
wolffd@0
|
58 % Now lets find some filters.... first the zeros then the poles
|
|
wolffd@0
|
59 zerofilts = SecondOrderFilter(CascadeZeroCF, CascadeZeroQ, fs);
|
|
wolffd@0
|
60 polefilts = SecondOrderFilter(CascadePoleCF, CascadePoleQ, fs);
|
|
wolffd@0
|
61 filters = [zerofilts polefilts(:,2:3)];
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 % Now we can set the DC gain of each stage.
|
|
wolffd@0
|
64 dcgain(2:NumberOfChannels)=CenterFreqs(1:NumberOfChannels-1)./ ...
|
|
wolffd@0
|
65 CenterFreqs(2:NumberOfChannels);
|
|
wolffd@0
|
66 dcgain(1) = dcgain(2);
|
|
wolffd@0
|
67 for i=1:NumberOfChannels
|
|
wolffd@0
|
68 filters(i,:) = SetGain(filters(i,:), dcgain(i), 0, fs);
|
|
wolffd@0
|
69 end
|
|
wolffd@0
|
70
|
|
wolffd@0
|
71 % Finally, let's design the front filters.
|
|
wolffd@0
|
72 front(1,:) = SetGain([0 1 -exp(-2*pi*EarPremphCorner/fs) 0 0], 1, fs/4, fs);
|
|
wolffd@0
|
73 topPoles = SecondOrderFilter(topf,CascadePoleQ(1), fs);
|
|
wolffd@0
|
74 front(2,:) = SetGain([1 0 -1 topPoles(2:3)], 1, fs/4, fs);
|
|
wolffd@0
|
75
|
|
wolffd@0
|
76 % Now, put them all together.
|
|
wolffd@0
|
77 filters = [front; filters];
|
|
wolffd@0
|
78
|
|
wolffd@0
|
79 if nargout > 2
|
|
wolffd@0
|
80 % Compute the gains needed so that everything is flat
|
|
wolffd@0
|
81 % when we do the filter inversion.
|
|
wolffd@0
|
82 [channels, width] = size(filters);
|
|
wolffd@0
|
83 len = length(CenterFreqs);
|
|
wolffd@0
|
84 clear cfs2;
|
|
wolffd@0
|
85 cfs2(1) = fs/2;
|
|
wolffd@0
|
86 cfs2(1:2:2*len) = CenterFreqs;
|
|
wolffd@0
|
87 cfs2(2:2:2*len-1) = (CenterFreqs(1:len-1) + CenterFreqs(2:len))/2;
|
|
wolffd@0
|
88 % cfs2 = CenterFreqs;
|
|
wolffd@0
|
89
|
|
wolffd@0
|
90 resp = zeros(length(cfs2), channels);
|
|
wolffd@0
|
91 for i=1:channels
|
|
wolffd@0
|
92 resp(:,i) = FreqResp(filters(i,:), cfs2, fs)';
|
|
wolffd@0
|
93 end
|
|
wolffd@0
|
94 % Each column vector of resp contains one channel's response
|
|
wolffd@0
|
95 % across all frequencies (in dB).
|
|
wolffd@0
|
96 cumresp = cumsum(resp')';
|
|
wolffd@0
|
97 % cumresp now contains the total gain at output of i'th channel at
|
|
wolffd@0
|
98 % frequency f
|
|
wolffd@0
|
99 a=10.^(cumresp(:,3:channels)/20);
|
|
wolffd@0
|
100 a = a.* a; % Have to go through each filter twice.
|
|
wolffd@0
|
101 gains = [0; 0; a \ ones(length(cfs2),1)];
|
|
wolffd@0
|
102 end
|
