THE AUDITORY MODELING TOOLBOX
Go to function
GAMMATONE - Gammatone filter coefficients
Program code:
function [b,a,delay,z,p,k]=gammatone(fc,fs,varargin)
%GAMMATONE Gammatone filter coefficients
% Usage: [b,a] = gammatone(fc,fs,n,betamul);
% [b,a] = gammatone(fc,fs,n);
% [b,a] = gammatone(fc,fs);
%
% Input parameters:
% fc : center frequency in Hz.
% fs : sampling rate in Hz.
% n : filter order.
% beta : bandwidth of the filter.
%
% Output parameters:
% b : nominator coefficients.
% a : denominator coefficients.
%
% GAMMATONE(fc,fs,n,betamul) computes the filter coefficients of a
% digital gammatone filter with center frequency fc, order n, sampling
% rate fs and bandwith determined by betamul. The bandwidth beta of
% each filter is determined as betamul times audfiltbw of the center
% frequency of corresponding filter.
%
% By default, the returned filter coefficients comes from the all-pole
% approximation described in Lyon (1997). The filters are normalized to
% have a 0 dB attenuation at the center frequency (another way of
% stating this is that their impulse responses will have unit area).
%
% GAMMATONE(fc,fs,n) will do the same but choose a filter bandwidth
% according to Glasberg and Moore (1990).
%
% GAMMATONE(fc,fs) will do as above for a 4th order filter.
%
% If fc is a vector, each entry of fc is considered as one center
% frequency, and the corresponding coefficients are returned as row
% vectors in the output.
%
% The inpulse response of the gammatone filter is given by:
%
% g(t) = a*t^(n-1)*cos(2*pi*fc*t)*exp(-2*pi*beta*t)
%
% GAMMATONE takes the following flags at the end of the line of input
% arguments:
%
% 'allpole' Compute the all-pole approximation of Gammatone
% filters by Lyon. This is the default
%
% 'classic' Compute the classical mixed pole-zero approximation of
% gammatone filters.
%
% 'complex' Generate filter coefficients corresponding to a
% complex valued filterbank modulated by exponential
% functions. This is useful for envelope extration
% purposes.
%
% 'real' Generate real-valued filters.
%
% 'casualphase' This makes the phase of each filter start at zero.
% This is the default.
%
% 'peakphase' This makes the phase of each filter be zero when the
% envelope of the impulse response of the filter peaks.
%
% 'exppeakphase' Experimental version of peakphase. In addition to
% peakphase, the output signal is delayed such that the
% maxima of the corresponding Gammatone impulse
% responses are aligned. This option has been created to
% produce some of the figures from Patterson et al.
% (1987).
%
% '0dBforall' This scales the amplitude of each filter to have an
% impulse response of 0dB. This is default.
%
% '6dBperoctave' This scales the amplitude of each filter to have an
% impulse response of +/-6dB per octave.
%
%
% To create the filter coefficients of a 1-erb spaced filter bank using
% gammatone filters use the following construction:
%
% [b,a] = gammatone(erbspacebw(flow,fhigh),fs,'complex');
%
% To apply the (complex valued) filters to an input signal, use
% FILTERBANKZ:
%
% outsig = 2*real(ufilterbankz(b,a,insig));
%
% References:
% A. Aertsen and P. Johannesma. Spectro-temporal receptive fields of
% auditory neurons in the grassfrog. I. Characterization of tonal and
% natural stimuli. Biol. Cybern, 38:223--234, 1980.
%
% R. Lyon. All pole models of auditory filtering. In E. R. Lewis, G. R.
% Long, R. F. Lyon, P. M. Narins, C. R. Steele, and E. Hecht-Poinar,
% editors, Diversity in auditory mechanics, pages 205--211. World
% Scientific Publishing, Singapore, 1997.
%
% R. Patterson, I. Nimmo-Smith, J. Holdsworth, and P. Rice. An efficient
% auditory filterbank based on the gammatone function. APU report, 2341,
% 1987.
%
%
% Url: http://amtoolbox.org/amt-1.1.0/doc/common/gammatone.php
% Copyright (C) 2009-2021 Piotr Majdak, Clara Hollomey, and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 1.1.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 : Stephan Ewert, Peter L. Soendergaard, modified by Christian
% Klemenschitz
% ------ Checking of input parameters ---------
% TODO: The phases of the filters all start at zero. This means that the
% real value of the impulse response of the filters does peak at the same
% time as the absolute value does. Include option to shift the phases so
% all filters have a distinct peak.
if nargin<2
error('%s: Too few input arguments.',upper(mfilename));
end;
if ~isnumeric(fs) || ~isscalar(fs) || fs<=0
error('%s: fs must be a positive scalar.',upper(mfilename));
end;
if ~isnumeric(fc) || ~isvector(fc) || any(fc<0) || any(fc>fs/2)
error(['%s: fc must be a vector of positive values that are less than half ' ...
'the sampling rate.'],upper(mfilename));
end;
definput.keyvals.n=4;
definput.keyvals.betamul=[];
definput.flags.real={'real','complex'};
definput.flags.phase={'causalphase','peakphase','exppeakphase'};
definput.flags.scale={'0dBforall','6dBperoctave'};
definput.flags.filtertype={'allpole','classic'};
[flags,keyvals,n,betamul] = ltfatarghelper({'n','betamul'},definput,varargin);
if ~isnumeric(n) || ~isscalar(n) || n<=0 || fix(n)~=n
error('%s: n must be a positive, integer scalar.',upper(mfilename));
end;
if isempty(betamul)
% This formula comes from patterson1987efficient, but it is easier to
% find in the Hohmann paper.
betamul = (factorial(n-1))^2/(pi*factorial(2*n-2)*2^(-(2*n-2)));
else
if ~isnumeric(betamul) || ~isscalar(betamul) || betamul<=0
error('%s: beta must be a positive scalar.',upper(mfilename));
end;
end;
% ------ Computation --------------------------
% ourbeta is used in order not to mask the beta function.
ourbeta = betamul*audfiltbw(fc);
nchannels = length(fc);
if flags.do_allpole
if flags.do_real
warning(['FIXME: The real-valued allpole filters are not scaled ' ...
'correctly.']);
b=zeros(nchannels,1);
a=zeros(nchannels,2*n+1);
% This is when the function peaks.
delay = 3./(2*pi*ourbeta);
for ii = 1:nchannels
% convert to radians
theta = 2*pi*fc(ii)/fs;
phi = 2*pi*ourbeta(ii)/fs;
alpha = -exp(-phi)*cos(theta);
b1 = 2*alpha;
b2 = exp(-2*phi);
a0 = abs( (1+b1*cos(theta)-1i*b1*sin(theta)+b2*cos(2*theta)-1i*b2*sin(2*theta)) / (1+alpha*cos(theta)-1i*alpha*sin(theta)) );
% Compute the position of the pole
atilde = exp(-phi - 1i*theta);
% Repeat the pole n times, and expand the polynomial
a2=poly([atilde*ones(1,n),conj(atilde)*ones(1,n)]);
if flags.do_6dBperoctave
b2=a0^n *(fs/fc(ii)/n);
else
% Scale to get 0 dB attenuation, FIXME: Does not work, works only
% for fc=fs/4
b2=a0^n;
end
% Signal peaks at envelope maximum
if flags.do_exppeakphase
insig = [1 , zeros(1,8191)];
outsig = 2*real(ufilterbankz(b2,a2,insig));
envmax = find( abs(outsig) == max(abs(outsig)) );
sigmax = find( outsig == max(outsig) );
% Equation 18 from Hohmanns paper, but 45 deg phasedelayed.
phi_delay = fc(ii)*(-2*pi-pi/4)*(envmax - sigmax)/fs;
% Equation 19 from Hohmanns paper
b2 = b2 * exp(1i *phi_delay);
end
if flags.do_peakphase
b2=b2*exp(2*pi*1i*fc(ii)*delay(ii));
end;
% Place the result (a row vector) in the output matrices.
b(ii,:)=b2;
a(ii,:)=a2;
end;
end;
if flags.do_complex
b=zeros(nchannels,1);
a=zeros(nchannels,n+1);
% This is when the function peaks.
delay = 3./(2*pi*ourbeta);
for ii = 1:nchannels
% convert to radians
theta = 2*pi*fc(ii)/fs;
phi = 2*pi*ourbeta(ii)/fs;
% Compute the position of the pole
% The commented line is the old code from the days of yore
atilde = exp(-2*pi*ourbeta(ii)/fs - 1i*2*pi*fc(ii)/fs);
%atilde = exp(-2*pi*ourbeta(ii)/fs + 1i*2*pi*fc(ii)/fs);
% Repeat the pole n times, and expand the polynomial
a2=poly(atilde*ones(1,n));
btmp=1-exp(-2*pi*ourbeta(ii)/fs);
% Amplitude scaling
if flags.do_6dBperoctave
b2=btmp.^n *(fs/fc(ii)/n);
else
b2=btmp.^n;
end
% Signal peaks at envelope maximum
if flags.do_exppeakphase
insig = [1 , zeros(1,8191)];
outsig = 2*real(ufilterbankz(b2,a2,insig));
envmax = find( abs(outsig) == max(abs(outsig)) );
sigmax = find( outsig == max(outsig) );
% Equation 18 from Hohmanns paper, but 45 deg phasedelayed.
phi_delay = fc(ii)*(-2*pi-pi/4)*(envmax - sigmax)/fs;
% Equation 19 from Hohmanns paper
b2 = b2 * exp(1i *phi_delay);
end
if flags.do_peakphase
b2=b2*exp(2*pi*1i*fc(ii)*delay(ii));
end;
% Place the result (a row vector) in the output matrices.
b(ii,:)=b2;
a(ii,:)=a2;
% Compute the z,p,k representation
z=[];
p=atilde*ones(1,n);
k=1;
end;
end;
else
if flags.do_real
b=zeros(nchannels,n+1);
a=zeros(nchannels,2*n+1);
% This is when the function peaks.
delay = 3./(2*pi*ourbeta);
for ii = 1:nchannels
% convert to radians
theta = 2*pi*fc(ii)/fs;
phi = 2*pi*ourbeta(ii)/fs;
alpha = -exp(-phi)*cos(theta);
b1 = 2*alpha;
b2 = exp(-2*phi);
a0 = abs( (1+b1*cos(theta)-1i*b1*sin(theta)+b2*cos(2*theta)-1i*b2*sin(2*theta)) / (1+alpha*cos(theta)-1i*alpha*sin(theta)) );
% Compute the position of the pole
atilde = exp(-phi-1i*theta);
% Repeat the conjugate pair n times, and expand the polynomial
a2 = poly([atilde*ones(1,n),conj(atilde)*ones(1,n)]);
% Compute the position of the zero, just the real value of the pole
btilde = real(atilde);
% Repeat the zero n times, and expand the polynomial
b2 = poly(btilde*ones(1,n));
% Amplitude scaling
if flags.do_6dBperoctave
b2 = b2*(a0^n) *(fs/fc(ii)/n);
else
% Scale to get 0 dB attenuation
b2=b2*(a0^n);
end
% Signal peaks at envelope maximum
if flags.do_exppeakphase
insig = [1 , zeros(1,8191)];
outsig = 2*real(ufilterbankz(b2,a2,insig));
envmax = find( abs(outsig) == max(abs(outsig)) );
sigmax = find( outsig == max(outsig) );
% Equation 18 from Hohmanns paper
phi_delay = fc(ii)*(-2*pi)*(envmax - sigmax)/fs;
% Equation 19 from Hohmanns paper
b2 = b2 * exp(1i *phi_delay);
end
if flags.do_peakphase
b2=b2*exp(2*pi*1i*fc(ii)*delay(ii));
end;
% Place the result (a row vector) in the output matrices.
b(ii,:)=b2;
a(ii,:)=a2;
end;
% Octave produces small imaginary values
b=real(b);
a=real(a);
end;
if flags.do_complex
warning(['FIXME: The complex-valued mixed pole-zero filters are not scaled ' ...
'correctly.']);
b=zeros(nchannels,n+1);
a=zeros(nchannels,n+1);
% This is when the function peaks.
delay = 3./(2*pi*ourbeta);
for ii = 1:nchannels
% convert to radians
theta = 2*pi*fc(ii)/fs;
phi = 2*pi*ourbeta(ii)/fs;
alpha = -exp(-phi)*cos(theta);
b1 = 2*alpha;
b2 = exp(-2*phi);
a0 = abs( (1+b1*cos(theta)-1i*b1*sin(theta)+b2*cos(2*theta)-1i*b2*sin(2*theta)) / (1+alpha*cos(theta)-1i*alpha*sin(theta)) );
% Compute the position of the pole
% The line commented line below is the old code that returned
% negative frequencies.
%atilde = exp(-2*pi*ourbeta(ii)/fs - 1i*2*pi*fc(ii)/fs);
atilde = exp(-2*pi*ourbeta(ii)/fs + 1i*2*pi*fc(ii)/fs);
% Repeat the pole n times, and expand the polynomial
a2=poly(atilde*ones(1,n));
% Compute the position of the zero, just the real value of the pole
btilde = real(atilde);
% Repeat the zero n times, and expand the polynomial
b2 = poly(btilde*ones(1,n));
% Amplitude scaling
if flags.do_6dBperoctave
b2=b2*(a0^n) *(fs/fc(ii)/n);
else
% Scale to get 0 dB attenuation
b2=b2*(a0^n);
end
% Signal peaks at envelope maximum
if flags.do_exppeakphase
insig = [1 , zeros(1,8191)];
outsig = 2*real(ufilterbankz(b2,a2,insig));
envmax = find( abs(outsig) == max(abs(outsig)) );
sigmax = find( outsig == max(outsig) );
% Equation 18 from Hohmanns paper
phi_delay = fc(ii)*(-2*pi)*(envmax - sigmax)/fs;
% Equation 19 from Hohmanns paper
b2 = b2 * exp(1i *phi_delay);
end
if flags.do_peakphase
b2=b2*exp(2*pi*1i*fc(ii)*delay(ii));
end;
% Place the result (a row vector) in the output matrices.
b(ii,:)=b2;
a(ii,:)=a2;
end;
end;
end;
Build with Bootstrap