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