comparison toolboxes/MIRtoolbox1.3.2/AuditoryToolbox/DesignLyonFilters.m @ 0:e9a9cd732c1e tip

first hg version after svn

+% [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