function GFB = may2011_gammatoneinit(fs,lowFreq,upFreq,nFilter,bUseEar,bAlign,bInfo)
%may2011_gammatoneINIT Initialize gammatone filterbank structure
%
% Usage:
% GFB = may2011_gammatoneinit(FS,FLOW,FUP,NFILTER,BEAR,BALIGN,BINFO)
%
% Input parameters:
% FS : sampling frequency in Hz
% FLOW : center frequency of lowest auditory filter
% (default, FLOW = 80)
% FUP : center frequency of highest auditory filter
% (default, FUP = 8e3)
% NFILTER : number of auditory filters which will be spaced linear on the
% ERB scale. If NFILTER is not a scalar but a vector, the first
% value is assumed to represent the number of auditory channels
% and the following values represents the indices of the
% filters which should be processed. This can be useful if a
% large number of channels is required as MATLAB might run out
% of memory for a longer exerpt if all filters are computed in
% one step. (default, NFILTER = round(freq2erb(FS/2))
% BEAR : Adjust gain coefficients of the auditory channels to
% incorporate middle ear effects. Note that this feature can
% only be used if the center frequencies of the auditory
% channels are above 20 hertz and below 12500 hertz.
% (default, BEAR = true)
% BALIGN : time-aligned gammatone output (non-causal output)
% (default, BALIGN = false)
% BINFO : info flag printing gammatone parameters on the screen
% (default, BINFO = false)
%
% Output parameters:
% GFB : gammatone parameter structure which can be passed as second
% input argument to the function gammatone
%
% Examples:
% nSamples = 500;
% % Initialize gammatone parameter structure
% GFB = gammatoneInit(20e3);
% % Filter impulse with gammatone filtering
% bm = gammatone([1; zeros(nSamples-1,1)],GFB);
% % Plot result
% waveplot(1:nSamples,GFB.cf,bm);
%
% See also: gammatone
%
% Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/modelstages/may2011_gammatoneinit.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/>.
% Developed with Matlab 7.4.0.287 (R2007a). Please send bug reports to:
%
% Author : Tobias May, 2008
% TUe Eindhoven and Philips Research
% t.may@tue.nl tobias.may@philips.com
%
% History :
% v.1.0 2008/01/21
% v.1.1 2008/02/04 added phase alignment
% v.1.2 2008/02/08 fixed bug related to center frequency calculation
% ***********************************************************************
% Check for proper input arguments
if nargin < 1 || nargin > 7
help(mfilename);
error('Wrong number of input arguments!');
end
% Set default values ...
if nargin < 2 || isempty(lowFreq); lowFreq = 80; end
if nargin < 3 || isempty(upFreq); upFreq = min(5e3,fs/2); end
if nargin < 4 || isempty(nFilter); nFilter = round(freq2erb(fs/2)); end
if nargin < 5 || isempty(bUseEar); bUseEar = false; end
if nargin < 6 || isempty(bAlign); bAlign = false; end
if nargin < 7 || isempty(bInfo); bInfo = false; end
% First value corresponds to the total number of auditory filters
nTotalFilter = nFilter(1);
% Figure out how many auditory filters should be processed
if length(nFilter) == 1
filter2Process = 1:nTotalFilter;
else
filter2Process = nFilter(2:end);
end
% Transform frequencies to erb domain
lowerERB = freq2erb(lowFreq);
upperERB = freq2erb(upFreq);
% Calculate center frequencies being linear spaced on the ERB scale
cf = erb2freq(linspace(lowerERB,upperERB,nTotalFilter));
% Bandwidth correction factor (related to 4rth order gamma function)
bandwidthCorrection = 1.019;
% Bandwidth
bw = erb(cf) * bandwidthCorrection;
% Compute gammatone envelope delay in samples (for 4th order)
delay = (3 * fs)./(bw * 2 * pi);
% Phase compensation factor
pc = -cf*3./bw;
% Store gammatone-related parameter
GFB = struct('object','4th order gammatone filterbank', ...
'fcnHandle','comp_may2011_gammatone','fs',fs, ...
'nFilter',nTotalFilter, ...
'filter2Process',[nTotalFilter filter2Process], ...
'lowerFreq',lowFreq,'upperFreq',upFreq,'cf',cf, ...
'bw',bw,'delay',delay,'phaseCorrection',pc, ...
'bOuterMiddleEar',bUseEar,'bPhaseAlign',bAlign,'bInfo',bInfo);
function out = isGFB(in)
%isGFB Check if input is a gammatone filterbank structure.
% This is a small helper function in order to check if gammatone
% filterbank is initialized properly.
%
%USAGE
% OUT = isGFB(IN)
%
%INPUT ARGUMENTS
% IN : input
%
%OUTPUT ARGUMENTS
% OUT : true/false depending on whether IN is a gammatone structure
%
%EXAMPLE
% % Create gammatone structure
% p = gammatoneInit(44.1e3);
% % Check
% isGFB(p)
%ans =
% 1
%
% See also gammatone.
% Developed with Matlab 7.4.0.287 (R2007a). Please send bug reports to:
%
% Author : Tobias May, 2008
% TUe Eindhoven and Philips Research
% t.may@tue.nl tobias.may@philips.com
%
% History :
% v.0.1 2008/05/11
% ***********************************************************************
% Check for proper input arguments
if nargin ~= 1
help(mfilename);
error('Wrong number of input arguments!');
end
% Initialize output
out = false;
% Check if IN is a structure
if isstruct(in)
% Required structure fields
reqFields = {'fcnHandle' 'fs' 'lowerFreq' 'upperFreq' ...
'filter2Process' 'bOuterMiddleEar' 'bPhaseAlign' 'bInfo'};
% Check if all required fields are present
if all(isfield(in,reqFields))
% Set flag to true
out = true;
end
end
function freqHz = erb2freq(erbf)
%erb2freq Convert ERB rate to frequency in Hz.
%
% Transformation is done according to Moore and Glasberg.
%
%USAGE
% FREQHZ = erb2freq(ERBF);
%
%INPUT ARGUMENTS
% ERBF : ERB-warped frequencies which should be transformned to
% frequency in Hz
%
%OUTPUT ARGUMENTS
% FREQHZ : frequency vector in Hz
%
% See also erb and freq2erb.
% Developed with Matlab 7.4.0.287 (R2007a). Please send bug reports to:
%
% Author : Tobias May, 2008
% TUe Eindhoven and Philips Research
% t.may@tue.nl tobias.may@philips.com
%
% History :
% v.1.0 2008/04/02
% ***********************************************************************
% Check for proper input arguments
if nargin ~= 1
help(mfilename);
error('Wrong number of input arguments!');
end
% Convert erb to frequency
freqHz = (10.^(erbf/21.4)-1)/4.37e-3;
function erbf = freq2erb(freqHz)
%freq2erb Convert frequency in Hz to ERB rate.
%
% Transformation is done according to Moore and Glasberg.
%
%USAGE
% ERBF = freq2erb(FREQHZ);
%
%INPUT ARGUMENTS
% FREQHZ : frequencies in Hz for which the erb width should be computed
%
%OUTPUT ARGUMENTS
% ERBF : auditory filter width
%
% See also erb and erb2freq.
% Developed with Matlab 7.4.0.287 (R2007a). Please send bug reports to:
%
% Author : Tobias May, 2008
% TUe Eindhoven and Philips Research
% t.may@tue.nl tobias.may@philips.com
%
% History :
% v.1.0 2008/05/08
% ***********************************************************************
% Check for proper input arguments
if nargin ~= 1
help(mfilename);
error('Wrong number of input arguments!');
end
% Convert frequency to erb
erbf = 21.4*log10(freqHz*0.00437 + 1.0);
function erbWidth = erb(freqHz)
%erb Compute Equivalent Rectangular Bandwidth (ERB).
%
% The ERB is calculated according to Moore and Glasberg at each given
% frequency.
%
%USAGE
% ERBWIDTH = ERB(FREQHZ);
%
%INPUT ARGUMENTS
% FREQHZ : frequencies in hertz for which the erb width should be
% computed
%
%OUTPUT ARGUMENTS
% ERBWIDTH : auditory filter width
%
% See also freq2erb and erb2freq.
% Developed with Matlab 7.4.0.287 (R2007a). Please send bug reports to:
%
% Author : Tobias May, 2008
% TUe Eindhoven and Philips Research
% t.may@tue.nl tobias.may@philips.com
%
% History :
% v.1.0 2008/05/08
% ***********************************************************************
% Check for proper input arguments
if nargin ~= 1
help(mfilename);
error('Wrong number of input arguments!');
end
% Compute ERB width
erbWidth = 24.7*(4.37e-3*freqHz+1);
% ***********************************************************************
% 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/>.
% ***********************************************************************