THE AUDITORY MODELING TOOLBOX
Go to function
MIDDLEEARFILTER - Middle ear filter
Program code:
function [b,a] = middleearfilter(fs,varargin)
%MIDDLEEARFILTER Middle ear filter
% Usage: [b,a]=middleearfilter(fs,varargin);
% [b,a]=middleearfilter(fs);
% [b,a]=middleearfilter;
%
% MIDDLEEARFILTER(fs) computes the filter coefficients of a FIR (or IIR)
% filter approximating the effect of the middle ear.
%
% The following parameter and flags can be specified additionally:
%
% 'order',order Sets the filter order of the computed FIR filter.
% Default value is 512.
%
% 'minimum' Calculates a minimum phase filter. This is the default.
%
% 'zero' returns a filter with zero phase. Since Matlab shifts the
% symmetric impulse response due to no negative indices.
% This results in a linear phase and hence a delay in the
% signal chain.
%
% 'lopezpoveda2001' Use data from Lopez-Poveda and Meddis (2001). These
% data are in turn derived from Goode et al. (1994).
% This is the default.
%
% 'jepsen2008' Use the data originally used in the Jepsen et al. (2008).
%
% 'verhulst2012' IIR filter approximating the middle ear transfer function
% based on Puria2003 (M1 filter) as used by Verhulst 2012.
%
% 'verhulst2015' IIR filter approximating the middle ear transfer function
% based on Puria2003 (M1 filter) as used by Verhulst 2015.
%
% 'verhulst2018' IIR filter approximating the middle ear transfer function
% based on Puria2003 (M1 filter) as used by Verhulst 2018.
%
% 'zilany2009' Second-order cascade IIR filters approximating the
% middle ear transfer function described by Ibrahim
% (2012, Appendix).
%
% 'zilany2009cat' Second-order cascade IIR filters approximating the
% middle ear transfer function described by Zilany et al.
% (2006, their Eq. 1--3)
%
% MIDDLEEARFILTER without any input arguments returns a table describing
% the frequency response of the middle ear filter. First column of the
% table contain frequencies and the second column contains the amplitude
% (stapes peak velocity in m/s at 0dB SPL) of the frequency like in figure
% 2b) of Lopez-Poveda and Meddis (2001).
%
% MIDDLEEARFILTER is meant to be used in conjunction with the LOPEZPOVEDA2001
% function, as the output is scaled to make lopezpoveda2001 work. If you are not
% using the lopezpoveda2001, you probably do not want to call this function. The
% following code displays the magnitude response of the filter:
%
% fs=16000;
% x=erbspace(0,fs/2,100);
% b=middleearfilter(fs);
% H=freqz(b,1,x,fs);
% semiaudplot(x,10*log10(abs(H).^2));
% xlabel('Frequency (Hz)');
% ylabel('Magnitude (dB)');
%
% See also: data_lopezpoveda2001, lopezpoveda2001
%
% References:
% R. Ibrahim. The role of temporal fine structure cues in speech
% perception. Ph.d., McMaster University, 2012.
%
% R. Goode, M. Killion, K. Nakamura, and S. Nishihara. New knowledge
% about the function of the human middle ear: development of an improved
% analog model. The American journal of otology, 15(2):145--154, 1994.
%
% E. Lopez-Poveda and R. Meddis. A human nonlinear cochlear filterbank.
% J. Acoust. Soc. Am., 110:3107--3118, 2001.
%
% M. Zilany and I. Bruce. Modeling auditory-nerve responses for high
% sound pressure levels in the normal and impaired auditory periphery. J.
% Acoust. Soc. Am., 120:1446--1466, 2006.
%
%
% Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/common/middleearfilter.php
% Copyright (C) 2009-2020 Piotr Majdak and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 1.0.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% AUTHOR: Peter L. Søndergaard, Katharina Egger, Alejandro Osses
%% ------ Check input options --------------------------------------------
a = 1; % if a == 1 then it is a FIR filter (not the case for the Verhulst et
% et al middle-ear filters, based on Puria2003)
% Define input flags
definput.flags.filtertype = {'lopezpoveda2001','jepsen2008', ...
'verhulst2012','verhulst2015','verhulst2018', ...
'zilany2009','zilany2009cat'};
definput.flags.phase = {'minimum','zero'};
definput.keyvals.order = 512;
% Parse input options
[flags,kv] = ltfatarghelper({},definput,varargin);
if flags.do_lopezpoveda2001
data = data_lopezpoveda2001('fig2b', 'no_plot');
if nargin==0
b = data;
else
if fs<=20000
% In this case, we need to cut the table because the sampling
% frequency is too low to accomodate the full range.
indx=find(data(:,1)<fs/2);
data = data(1:indx(end),:);
else
% otherwise the table will be extrapolated towards fs/2
% data point added every 1000Hz
lgth = size(data,1);
for ii = 1:floor((fs/2-data(end,1))/1000)
data(lgth+ii,1) = data(lgth+ii-1,1) + 1000;
% 1.1 corresponds to the decay of the last amplitude values = approx. ratio
% between amplitudes of frequency values seperated by 1000Hz
data(lgth+ii,2) = data(lgth+ii-1,2) / 1.1;
end
end;
% for the function fir2 the last data point has to be at fs/2
lgth = size(data,1);
if data(lgth,1) ~= fs/2
data(lgth+1,1) = fs/2;
data(lgth+1,2) = data(lgth,2) / (1+(fs/2-data(lgth,1))*0.1/1000);
end
% Extract the frequencies and amplitudes, and put them in the format
% that fir2 likes.
freq=[0;...
data(:,1).*(2/fs);...
];
ampl=[0;...
data(:,2);...
];
b = fir2(kv.order,freq,ampl);
b = b / 20e-6; % scaling for SPL in dB re 20uPa
if flags.do_minimum
X = fft(b);
Xmin = abs(X) .* exp(-1j*imag(hilbert(log(abs(X)))));
b = real(ifft(Xmin));
end
end
end;
if flags.do_jepsen2008
stapes_data = [...
50, 48046.39731;...
100, 24023.19865;...
200, 12011.59933;...
400, 6005.799663;...
600, 3720.406871;...
800, 2866.404385;...
1000, 3363.247811;...
1200, 4379.228921;...
1400, 4804.639731;...
1600, 5732.808769;...
1800, 6228.236688;...
2000, 7206.959596;...
2200, 9172.494031;...
2400, 9554.681282;...
2600, 10779.64042;...
2800, 12011.59933;...
3000, 14013.53255;...
3500, 16015.46577;...
4000, 18017.39899;...
4500, 23852.82136;...
5000, 21020.29882;...
5500, 22931.23508;...
6000, 28027.06509;...
6500, 28745.70779;...
7000, 32098.9;...
7500, 34504.4;...
8000, 36909.9;...
8500, 39315.4;...
9000, 41720.9;...
9500, 44126.4;...
10000,46531.9;...
];
% We need to find inverse because original data is stapes impedance and we
% need stapes velocity.
stapes_data (:,2) = 1./stapes_data(:,2);
if nargin==0
b = stapes_data;
else
if fs<=20000
% In this case, we need to cut the table because the sampling
% frequency is too low to accomodate the full range.
indx=find(stapes_data(:,1)<fs/2);
stapes_data=stapes_data(1:indx(end),:);
end;
% Extract the frequencies and amplitudes, and put them in the format
% that fir2 likes.
freq=[0;...
stapes_data(:,1).*(2/fs);...
1];
ampl=[0;...
stapes_data(:,2);...
0];
b = fir2(kv.order,freq,ampl);
% See the figure text for figure 1, Lopez (2001).
b = b/max(abs(fft(b)))*1e-8*10^(104/20);
end;
end;
if flags.do_verhulst2012 || flags.do_verhulst2015 || flags.do_verhulst2018
%%% PuriaM1 filter (Puria2003, Fig. 2A and 3A)
if flags.do_verhulst2012
% Source: AMT implementation and originally in fortran code
fc1 = 100; % as in cochlear_model.py, ../model2012-fortran/SourceCode/PuriaM1.f90
fc2 = 3000; % as in cochlear_model.py, ../model2012-fortran/SourceCode/PuriaM1.f90
puria_gain = 2*scaletodbspl(1,18); % puria_gain=10**(18./20.)*2. % The factor of 2 is explained
% as: 'VoltageDivisionGain=2d0 !2 resistor matching network(Zme+Zoch)'
end
if flags.do_verhulst2015
% Source: ../model2015/cochlear_model.py
fc1 = 600; % Hz
fc2 = 3000; % Hz
puria_gain = 2*scaletodbspl(1,18); % same gain as in verhulst2012
end
if flags.do_verhulst2018
fc1 = 600; % as in cochlear_model2018.py
fc2 = 4000; % as in cochlear_model2018.py
puria_gain = scaletodbspl(1,18); % puria_gain=10**(18./20.)
end
[b,a] = butter(1,[fc1 fc2]/(fs/2),'bandpass');
b = b*puria_gain;
end
if flags.do_zilany2009
% The following humanised middle ear filter is a digital implementation
% described by Rasha Ibrahim's thesis from 2012. She based this implementation
% on the work described by Pascal et al. (JASA 1998) */
%
% There are three resulting (2nd-order) filters. For an fs=100 kHz, the
% transfer functions to be obtained are (see Ibrahim, 2012, her Eq. A1, A2, and A3):
%
% 0.9979 - 1.9408 z^-1 + 0.9429 z^-2
% H1(z-1) = ------------------------------------
% 1.0000 - 1.9395 z^-1 + 0.9420 z^-2
%
% 0.9984 - 1.9226 z^-1 + 0.9415 z^-2
% H2(z-1) = ------------------------------------
% 1.0000 - 1.9244 z^-1 + 0.9379 z^-2
%
% 0.0286 + 0.0302 z^-1 + 0.0016 z^-2
% H3(z-1) = 0.5 x ------------------------------------
% 1.0000 - 1.6748 z^-1 + 0.7847 z^-2
% Note that the 0.5 is not defined in Rasha's thesis, but in the code
% that was included in the AMT toolbox version of zilany2009's model.
fp = 1e3; % prewarping frequency 1 kHz
C = 2*pi*fp/tan(pi*fp/fs);
m11=1/(C^2+5.9761e3*C+2.5255e7);
m12=(-2*(C^2)+2*2.5255e7);
m13=(C^2-5.9761e3*C+2.5255e7);
m14=(C^2+5.6665e3*C);
m15=-2*(C^2);
m16=(C^2-5.6665e3*C);
m21=1/(C^2+6.4255e3*C+1.3975e8);
m22=(-2*(C^2)+2*1.3975e8);
m23=(C^2-6.4255e3*C+1.3975e8);
m24=(C^2+5.8934e3*C+1.7926e8);
m25=-2*(C^2)+2*1.7926e8;
m26=C^2-5.8934e3*C+1.7926e8;
m31=1/(C^2+2.4891e4*C+1.2700e9);
m32=(-2*(C^2)+2*1.27e9);
m33=(C^2-2.4891e4*C+1.27e9);
m34=(3.1137e3*C+6.9768e8);
m35=2*6.9768e8;
m36=(-3.1137e3*C+6.9768e8);
megainmax=2;
a(1,1:3) = [1 m11*m12 m11*m13];
b(1,1:3) = m11*[m14 m15 m16];
a(2,1:3) = [1 m21*m22 m21*m23];
b(2,1:3) = m21*[m24 m25 m26];
a(3,1:3) = [1 m31*m32 m31*m33];
b(3,1:3) = m31*[m34 m35 m36]/megainmax;
end
if flags.do_zilany2009cat
% The following middle ear filter of the cat is a digital implementation
% described by Rasha Ibrahim's thesis from 2012. She based this implementation
% on the work described by Pascal et al. (JASA 1998) */
%
% There are three resulting (2nd-order) filters. The obtained 2nd order
% transfer functions at an fs=500 kHz can be found in (Zilany et al, 2006,
% their Eq. 1--3):
fp = 1e3; % prewarping frequency 1 kHz
C = 2*pi*fp/tan(pi*fp/fs);
m11 = C/(C + 693.48);
m12 = (693.48-C)/C;
m13 = 0.0;
m14 = 1.0;
m15 = -1.0;
m16 = 0.0;
m21 = 1/(C^2 + 11053*C + 1.163e8);
m22 = -2*(C^2) + 2.326e8;
m23 = C^2 - 11053*C + 1.163e8;
m24 = C^2 + 1356.3*C + 7.4417e8;
m25 = -2*(C^2) + 14.8834e8;
m26 = C^2 - 1356.3*C + 7.4417e8;
m31 = 1/(C^2 + 4620*C + 909059944);
m32 = -2*(C^2) + 2*909059944;
m33 = C^2 - 4620*C + 909059944;
m34 = 5.7585e5*C + 7.1665e7;
m35 = 14.333e7;
m36 = 7.1665e7 - 5.7585e5*C;
megainmax=41.1405;
a(1,1:3) = [1 m11*m12 m11*m13];
b(1,1:3) = m11*[m14 m15 m16];
a(2,1:3) = [1 m21*m22 m21*m23];
b(2,1:3) = m21*[m24 m25 m26];
a(3,1:3) = [1 m31*m32 m31*m33];
b(3,1:3) = m31*[m34 m35 m36]/megainmax;
end
if nargout == 0
if flags.do_zilany2009 || flags.do_zilany2009cat
N = round(fs/2);
insig = [1; zeros(N-1,1)];
Nr_cascaded = size(b,1);
for k = 1:Nr_cascaded
insig = filter(b(k,:),a(k,:),insig);
end
[h,w] = freqz(insig,1,N);
f = (w/pi)*fs/2;
figure;
semilogx(f,20*log10(abs(h)));
xlim([20 fs/2]);
ylim([-50 10])
grid on
xlabel('Frequency (Hz)');
xlabel('IIR middleearfilter');
end
end
% if flags.do_plot
% % Manually calculate the frequency response
% fmid = abs(fftreal(b));
% % Half the filter length.
% n2=length(fmid);
% % x-values for plotting.
% xplot=linspace(0,fs/2,n2);
% loglog(xplot/1000,fmid);
% xlabel('Frequency (kHz)');
% ylabel('FIR middleearfilter');
% end
Build with Bootstrap