Mercurial > hg > camir-aes2014
view 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 source
% [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