Mercurial > hg > camir-aes2014
diff 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 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/MIRtoolbox1.3.2/AuditoryToolbox/DesignLyonFilters.m Tue Feb 10 15:05:51 2015 +0000 @@ -0,0 +1,102 @@ +% [filters, freqs] = DesignLyonFilters(fs,EarQ,StepFactor) +% Design the cascade of second order filters and the front filters +% (outer/middle and compensator) needed for Lyon's Passive Short Wave +% (Second Order Sections) cochlear model. The variables used here come +% from Apple ATG Technical Report #13 titled "Lyon's Cochlear Model". +% +% Most of the parameters are hardwired into this m-function. The user +% settable parameters are the digital sampling rate, the basic Q of the +% each stage (usually 8 or 4), and the spacing factor between channels +% (usually .25 and .125, respectively.) +% +% The result is returned as rows of second order filters; three coefficients +% for the numerator and two for the denomiator. The coefficients are +% [A0 A1 A2 B1 B2]..........................(Malcolm 1 Feb. 1993) +% + +% (c) 1998 Interval Research Corporation +function [filters, CenterFreqs, gains] = DesignLyonFilters(fs,EarQ,StepFactor) + +if nargin < 2 + EarQ = 8; +end + +if nargin < 3 + StepFactor=EarQ/32; +end + +Eb = 1000.0; +EarZeroOffset = 1.5; +EarSharpness = 5.0; +EarPremphCorner = 300; + +% Find top frequency, allowing space for first cascade filter. +topf = fs/2.0; +topf = topf - (sqrt(topf^2+Eb^2)/EarQ*StepFactor*EarZeroOffset)+ ... + sqrt(topf^2+Eb^2)/EarQ*StepFactor; + +% Find place where CascadePoleQ < .5 +lowf = Eb/sqrt(4*EarQ^2-1); +NumberOfChannels = floor((EarQ*(-log(lowf + sqrt(lowf^2 + Eb^2)) + ... + log(topf + sqrt(Eb^2 + topf^2))))/StepFactor); + +% Now make an array of CenterFreqs..... This expression was derived by +% Mathematica by integrating 1/EarBandwidth(cf) and solving for f as a +% function of channel number. +cn = 1:NumberOfChannels; +CenterFreqs = (-((exp((cn*StepFactor)/EarQ)*Eb^2)/ ... + (topf + sqrt(Eb^2 + topf^2))) + ... + (topf + sqrt(Eb^2 + topf^2))./exp((cn*StepFactor)/EarQ))/2; + +% OK, now we can figure out the parameters of each stage filter. +EarBandwidth = sqrt(CenterFreqs.^2+Eb^2)/EarQ; +CascadeZeroCF = CenterFreqs + EarBandwidth * StepFactor * EarZeroOffset; +CascadeZeroQ = EarSharpness*CascadeZeroCF./EarBandwidth; +CascadePoleCF = CenterFreqs; +CascadePoleQ = CenterFreqs./EarBandwidth; + +% Now lets find some filters.... first the zeros then the poles +zerofilts = SecondOrderFilter(CascadeZeroCF, CascadeZeroQ, fs); +polefilts = SecondOrderFilter(CascadePoleCF, CascadePoleQ, fs); +filters = [zerofilts polefilts(:,2:3)]; + +% Now we can set the DC gain of each stage. +dcgain(2:NumberOfChannels)=CenterFreqs(1:NumberOfChannels-1)./ ... + CenterFreqs(2:NumberOfChannels); +dcgain(1) = dcgain(2); +for i=1:NumberOfChannels + filters(i,:) = SetGain(filters(i,:), dcgain(i), 0, fs); +end + +% Finally, let's design the front filters. +front(1,:) = SetGain([0 1 -exp(-2*pi*EarPremphCorner/fs) 0 0], 1, fs/4, fs); +topPoles = SecondOrderFilter(topf,CascadePoleQ(1), fs); +front(2,:) = SetGain([1 0 -1 topPoles(2:3)], 1, fs/4, fs); + +% Now, put them all together. +filters = [front; filters]; + +if nargout > 2 + % Compute the gains needed so that everything is flat + % when we do the filter inversion. + [channels, width] = size(filters); + len = length(CenterFreqs); + clear cfs2; + cfs2(1) = fs/2; + cfs2(1:2:2*len) = CenterFreqs; + cfs2(2:2:2*len-1) = (CenterFreqs(1:len-1) + CenterFreqs(2:len))/2; +% cfs2 = CenterFreqs; + + resp = zeros(length(cfs2), channels); + for i=1:channels + resp(:,i) = FreqResp(filters(i,:), cfs2, fs)'; + end + % Each column vector of resp contains one channel's response + % across all frequencies (in dB). + cumresp = cumsum(resp')'; + % cumresp now contains the total gain at output of i'th channel at + % frequency f + a=10.^(cumresp(:,3:channels)/20); + a = a.* a; % Have to go through each filter twice. + gains = [0; 0; a \ ones(length(cfs2),1)]; +end