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