GAMMABANK gammatone filter bank [b,a,fx,bx,gd]=(n,fs,w,fc,bw,ph,k)
Usage:
(1) [b,a,fx,bx,gd]=gammabank(0.35,fs,'',[100 6000]);
y = filterbank(b,a,s,gd);
Will create an erb-spaced filterbank between 100 Hz and 6kHz
with a filter spacing of 0.35 erb and a default bandwidth
of 1.019 erb. Omitting the "y =" from the second line will plot
a spectrogram.
Inputs:
n number of filters in filterbank or the filter spacing in
bark/erb/k-mel. Set n=0 if fc lists centre frequencies explicitly.
fs sample rate in Hz
fc centre frequencies [default = [100 6000] ]
bw bandwidths [default = 1.019*erb(fc) ]
ph phase of gammatone impulse repsonse [default = 0]
k filter ordere [default = 4]
w any sensible combination of the following:
'e' = erb scale for filter spacing, frequencies ('F' option)
and bandwidths ('W' option but can be overridden by 'EMBLH')
'm','b','l','h' = mel, bark, log10 and Hz scale
'E' = erb scale for bandwidths ('W' option)
'M','B','L','H' = mel, bark, log10 and Hz scale
'n' = n input gives number of filters [default if n>=1]
'N' = n input gives filter spacing [default if n<1]
'f' = fc is in Hz [default]
'F' = fc is in mel/erb-rate/bark/log10
'w' = bw inputs are in Hz [default]
'W' = bw inputs are in multiples of df/dx where x=bark/erb/mel etc
'k' = force a filter at 1kHz
['d' = choose ph() so that all filters have zero DC gain]
['a' = use all-pole gammtone funtion: see [1]]
['s' = use Slaney gammatone approximation: see [2]]
['z' = use one-zero gammatone function: see [1]]
'g' = plot filter responses [default if no output arguments present]
'G' = plot frequency responses on a log axis
Outputs:
b/a filter coefficients: one filter per row
fx,bx centre frequencies and bandwidths in Hz
gd group delay at the centre frequencies (in samples)0001 function [b,a,fx,bx,gd]=gammabank(n,fs,w,fc,bw,ph,k) 0002 %GAMMABANK gammatone filter bank [b,a,fx,bx,gd]=(n,fs,w,fc,bw,ph,k) 0003 % 0004 % Usage: 0005 % (1) [b,a,fx,bx,gd]=gammabank(0.35,fs,'',[100 6000]); 0006 % y = filterbank(b,a,s,gd); 0007 % 0008 % Will create an erb-spaced filterbank between 100 Hz and 6kHz 0009 % with a filter spacing of 0.35 erb and a default bandwidth 0010 % of 1.019 erb. Omitting the "y =" from the second line will plot 0011 % a spectrogram. 0012 % Inputs: 0013 % n number of filters in filterbank or the filter spacing in 0014 % bark/erb/k-mel. Set n=0 if fc lists centre frequencies explicitly. 0015 % fs sample rate in Hz 0016 % fc centre frequencies [default = [100 6000] ] 0017 % bw bandwidths [default = 1.019*erb(fc) ] 0018 % ph phase of gammatone impulse repsonse [default = 0] 0019 % k filter ordere [default = 4] 0020 % w any sensible combination of the following: 0021 % 'e' = erb scale for filter spacing, frequencies ('F' option) 0022 % and bandwidths ('W' option but can be overridden by 'EMBLH') 0023 % 'm','b','l','h' = mel, bark, log10 and Hz scale 0024 % 'E' = erb scale for bandwidths ('W' option) 0025 % 'M','B','L','H' = mel, bark, log10 and Hz scale 0026 % 0027 % 'n' = n input gives number of filters [default if n>=1] 0028 % 'N' = n input gives filter spacing [default if n<1] 0029 % 0030 % 'f' = fc is in Hz [default] 0031 % 'F' = fc is in mel/erb-rate/bark/log10 0032 % 'w' = bw inputs are in Hz [default] 0033 % 'W' = bw inputs are in multiples of df/dx where x=bark/erb/mel etc 0034 % 0035 % 'k' = force a filter at 1kHz 0036 % ['d' = choose ph() so that all filters have zero DC gain] 0037 % ['a' = use all-pole gammtone funtion: see [1]] 0038 % ['s' = use Slaney gammatone approximation: see [2]] 0039 % ['z' = use one-zero gammatone function: see [1]] 0040 % 0041 % 'g' = plot filter responses [default if no output arguments present] 0042 % 'G' = plot frequency responses on a log axis 0043 % 0044 % Outputs: 0045 % b/a filter coefficients: one filter per row 0046 % fx,bx centre frequencies and bandwidths in Hz 0047 % gd group delay at the centre frequencies (in samples) 0048 % 0049 0050 % The impulse response of filter i is proportional to: 0051 % h(n)=((n/fs).^(k-1))*cos(2*pi*fx(i)*n/fs+ph(i))*exp(-2*pi*bx(i)*n/fs) 0052 % where n=0,1,2,... 0053 % Note that the DC gain is only equal to zero for one particular value of ph(i) 0054 % The filters are normalized to have unity gain at the centre frequencies 0055 % 0056 % References 0057 % [1] R. F. Lyon, A. G. Katsiamis, and E. M. Drakakis. 0058 % History and future of auditory filter models. 0059 % In Proc Intl Symp Circuits and Systems, pages 3809–3812, 2010. 0060 % doi: 10.1109/ISCAS.2010.5537724. 0061 % [2] M. Slaney. 0062 % An efficient implementation of the patterson-holdsworth auditory filter bank. 0063 % Technical report, Apple Computer, Perception Group, Tech. Rep, 1993. 0064 0065 % Copyright (C) Mike Brookes 2009-2010 0066 % Version: $Id: gammabank.m 713 2011-10-16 14:45:43Z dmb $ 0067 % 0068 % VOICEBOX is a MATLAB toolbox for speech processing. 0069 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0070 % 0071 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0072 % This program is free software; you can redistribute it and/or modify 0073 % it under the terms of the GNU General Public License as published by 0074 % the Free Software Foundation; either version 2 of the License, or 0075 % (at your option) any later version. 0076 % 0077 % This program is distributed in the hope that it will be useful, 0078 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0079 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0080 % GNU General Public License for more details. 0081 % 0082 % You can obtain a copy of the GNU General Public License from 0083 % http://www.gnu.org/copyleft/gpl.html or by writing to 0084 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0085 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0086 0087 if nargin<7 0088 k=[]; 0089 if nargin<6 0090 ph=[]; 0091 if nargin<5 0092 bw=[]; 0093 if nargin<4 0094 fc=[]; 0095 if nargin<3 0096 w=''; 0097 end 0098 end 0099 end 0100 end 0101 end 0102 fx=fc(:); 0103 bx=bw(:); 0104 if ~numel(k) 0105 k=4; 0106 end 0107 if ~numel(fx) 0108 fx=[100; 6000]; % default 0109 end 0110 wr='e'; % default frequency warping is erb 0111 for i=1:length(w) 0112 if any(w(i)=='bmlef'); 0113 wr=w(i); 0114 end 0115 end 0116 if any(w=='k') 0117 fk=1000; 0118 switch wr % convert 1kHz to spacing units 0119 case 'b' 0120 fk=frq2bark(fk); 0121 case 'l' 0122 fk=log10(fk); 0123 case 'e' 0124 fk=frq2erb(fk); 0125 end 0126 else 0127 fk=0; 0128 end 0129 if any(w)=='W' 0130 wb=wr; 0131 else 0132 wb='h'; % default bandwidth units are Hz 0133 end 0134 for i=1:length(w) 0135 if any(w(i)=='BMLEF'); 0136 wb=w(i)+'a'-'A'; % convert to lower case 0137 end 0138 end 0139 if ~numel(bx) 0140 bx=1.019; 0141 wb='e'; 0142 end 0143 if any(w)=='F' % convert centre frequencies to Hz 0144 switch wr 0145 case 'b' 0146 fx=bark2frq(fx); 0147 case 'm' 0148 fx=mel2frq(fx); 0149 case 'l' 0150 fx=10.^(fx); 0151 case 'e' 0152 fx=erb2frq(fx); 0153 end 0154 end 0155 0156 % now sort out the centre frequencies 0157 0158 if n>0 % n>0: filter end points specified 0159 bx=bx(1); % only use the first bx element 0160 if n==1 % only one filter requested 0161 fx=fx(1); % just use the first frequency 0162 else 0163 switch wr % convert end frequencies to spacing units 0164 case 'b' 0165 fx=frq2bark(fx); 0166 case 'm' 0167 fx=frq2mrl(fx); 0168 case 'l' 0169 fx=log10(fx); 0170 case 'e' 0171 fx=frq2erb(fx); 0172 end 0173 if n<1 || any(w=='N') % n = filter spacing 0174 if fk % force filter to 1 kHz 0175 f0=fk-n*floor((fk-fx(1))/n); 0176 else % centre filters in range 0177 f0=(fx(2)+fx(1)-n*floor((fx(2)-fx(1))/n))/2; 0178 end 0179 fx=(f0:n:fx(2))'; 0180 0181 else % n = number of filters specified 0182 % Multiple filters - evenly spaced 0183 fx=linspace(fx(1),fx(2),n)'; % centre frequencies in spacing units 0184 if fk % force a filter at 1kHz 0185 ik=1+ceil((fk-fx(1))*(n-1)/(fx(n)-fx(1))); % index of centre freq immediately above 1 kHz 0186 if ik>n || ik>1 && ((fk-fx(1))*(fx(n)-fx(ik-1))>(fx(n)-fk)*(fx(ik)-fx(1))) 0187 fx=fx(1)+(fx-fx(1))*(fk-fx(1))/(fx(ik)-fx(1)); 0188 else 0189 fx=fx(n)+(fx-fx(n))*(fx(n)-fk)/(fx(n)-fx(ik-1)); 0190 end 0191 end 0192 end 0193 switch wr % convert back to Hz 0194 case 'b' 0195 fx=bark2frq(fx); 0196 case 'm' 0197 fx=mel2frq(fx); 0198 case 'l' 0199 fx=10.^(fx); 0200 case 'e' 0201 fx=erb2frq(fx); 0202 end 0203 end 0204 0205 end 0206 % now sort out the bandwidths 0207 nf=numel(fx); 0208 if numel(bx)==1 0209 bx=bx(ones(nf,1)); % replicate if necessary 0210 end 0211 switch wb % convert bandwidth to Hz 0212 case 'b' 0213 [dum,bwf]=frq2bark(fx); 0214 case 'm' 0215 [dum,bwf]=frq2mel(fx); 0216 case 'l' 0217 bwf=fx*log(10); 0218 case 'e' 0219 [dum,bwf]=frq2erb(fx); 0220 case 'h' 0221 bwf=ones(nf,1); 0222 end 0223 bx=bx.*bwf; 0224 if ~numel(ph) 0225 ph=0; 0226 end 0227 if numel(ph)==1 0228 ph=ph(ones(nf,1)); % replicate if necessary 0229 else 0230 ph=ph(:); 0231 end 0232 % 0233 % t=(0:ceil(10*fs/(2*pi*bnd)))/fs; % five time constants 0234 % gt=t.^(n-1).*cos(2*pi*cfr*t+phi).*exp(-2*pi*bnd*t); 0235 % gt=gt/sqrt(mean(gt.^2)); % normalize 0236 % figure(1); 0237 % plot(t,gt); 0238 % title('Desired Impulse response'); 0239 % xlabel(['Time (' xticksi 's)']); 0240 % 0241 ww=exp((1i*fx-bx)*2*pi/fs); 0242 a=round([1 cumprod((-k:-1)./(1:k))]); % create binomial coefficients 0243 b=conv(a,(0:k-1).^(k-1)); 0244 b=exp(1i*ph)*b(1:k); 0245 wwp=repmat(ww,1,k+1).^repmat(0:k,nf,1); 0246 denc=repmat(a,nf,1).*wwp; 0247 numc=b.*wwp(:,1:k); 0248 b=zeros(nf,2*k); 0249 a=zeros(nf,2*k+1); 0250 gd=zeros(nf,1); 0251 ww=exp(2i*fx*pi/fs); 0252 for i=1:nf 0253 b(i,:)=real(conv(numc(i,:),conj(denc(i,:)))); 0254 a(i,:)=real(conv(denc(i,:),conj(denc(i,:)))); 0255 u=polyval(b(i,:),ww(i)); 0256 v=polyval(a(i,:),ww(i)); 0257 ud=polyval(b(i,:).*(0:2*k-1),ww(i)); 0258 vd=polyval(a(i,:).*(0:2*k),ww(i)); 0259 b(i,:)=b(i,:)*abs(v/u); 0260 gd(i)=real((v*ud-u*vd)/(u*v)); % group delay at centre freq in samples 0261 end 0262 0263 % now plot graph 0264 0265 if ~nargout || any(w=='g') || any(w=='G') 0266 ng=200; %number of points to plot 0267 if any(w=='G') 0268 fax=logspace(log10(fx(1)/4),log10(fs/2),ng); 0269 else 0270 fax=linspace(0,fs/2,ng); 0271 end 0272 ww=exp(2i*pi*fax/fs); 0273 gg=zeros(nf,ng); 0274 for i=1:nf 0275 gg(i,:)=10*log10(abs(polyval(b(i,:),ww)./polyval(a(i,:),ww))); 0276 end 0277 if any(w=='G') 0278 semilogx(fax,gg','-b'); 0279 set(gca,'xlim',[fax(1) fax(end)]); 0280 else 0281 plot(fax,gg','-b'); 0282 end 0283 xlabel(['Frequency (' xticksi 'Hz)']); 0284 set(gca,'ylim',[-50 1]); 0285 title(sprintf('%d Gammatone Filters (Opt=%s)',nf,w)); 0286 end