V_GAMMABANK gammatone filter bank [b,a,fx,bx,gd]=(n,fs,w,fc,bw,ph,k)
Usage:
(1) [b,a,fx,bx,gd,ph]=v_gammabank(0.35,fs,'',[100 6000]);
y = v_filterbank(b,a,s,gd);
Will create an erb-spaced v_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 v_filterbank or the filter spacing in bark/erb/k-mel.
Set n=0 if fc lists centre frequencies explicitly.
fs sample rate in Hz
w any sensible combination of the following:
(1) Filters are uniformly spaced in
'e' = erb scale (aka erb-rate scale) [default]
'b' = bark scale
'm' = mel scale
'l' = log10 scale
'h' = Hz
(2) Centre frequencies, fc, are in units of
'f' = Hz [default]
'F' = the same units as the filter spacing (see (1) above)
(3) Bandwidths, bw, are specified in units of
'w' = Hz [default]
'W' = the same units as the filter spacing (see (1) above)
'E' = erb scale (aka erb-rate scale)
'B' = bark scale
'M' = mel scale
'L' = log10 scale
'H' = Hz
(4) Number of filters
'n' = n input gives number of filters [default if n>=1]
'N' = n input gives filter spacing [default if n<1]
(5) Filter alignment
'k' = force a filter at 1kHz
(6) Output filter
'c' = Use complex gammatone filters whose impulse responses uses a complex exponential
instead of a cosine; the filter order is k instead of 2k. Take the real part of
the filter output to obtain the same signal as using the real-valued filter.
The imaginary part is approximately the Hilbert transform of the real part and
so the magnitude gives the envelope.
'q' = biquad filter; b(k,6,n) has the sos coefficients in the same form as sosfilt.m.
Use y=sosfilt(b,x) to apply the filter.
(6) Future possible options
['a' = use all-pole gammtone funtion: see [1]]
['s' = use Slaney gammatone approximation: see [2]]
['z' = use one-zero gammatone function: see [1]]
(7) Graph plotting
'g' = plot filter responses [default if no output arguments present]
'G' = plot frequency responses on a log axis
fc frequency range or, if n=0, list of centre frequencies [default: [100 6000] ]
bw bandwidth for all filters or a vector of bandwidths: one per filter [default = 1.019 erb ]
ph phase of all gammatone impulse repsonses or a vector of phases: one per filter
[default is chosen to give zero gain at DC and a positive real part for the gain at low frequencies]
k gammatone filter order; the filters have order 2k [default: 4]
Outputs:
b/a filter coefficients: one filter per row. The gain of each
filter is scaled to have unit magnitude at the centre frequency.
Dimensions are b(n,k+max(1,k-1)) and a(n,2*k+1) normally or else
b(n,max(1,k-1)) and a(n,k+1) if the 'c' option is specified.
If the q option is given then b(k,6,n) contains the biquad coefficients for filter n.
fx,bx centre frequencies and bandwidths in Hz
gd group delay at the centre frequencies (in samples)
ph phase of each gammatone impulse response0001 function [b,a,fx,bx,gd,ph]=v_gammabank(n,fs,w,fc,bw,ph,k) 0002 %V_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,ph]=v_gammabank(0.35,fs,'',[100 6000]); 0006 % y = v_filterbank(b,a,s,gd); 0007 % 0008 % Will create an erb-spaced v_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 v_filterbank or the filter spacing in bark/erb/k-mel. 0014 % Set n=0 if fc lists centre frequencies explicitly. 0015 % fs sample rate in Hz 0016 % w any sensible combination of the following: 0017 % (1) Filters are uniformly spaced in 0018 % 'e' = erb scale (aka erb-rate scale) [default] 0019 % 'b' = bark scale 0020 % 'm' = mel scale 0021 % 'l' = log10 scale 0022 % 'h' = Hz 0023 % 0024 % (2) Centre frequencies, fc, are in units of 0025 % 'f' = Hz [default] 0026 % 'F' = the same units as the filter spacing (see (1) above) 0027 % 0028 % (3) Bandwidths, bw, are specified in units of 0029 % 'w' = Hz [default] 0030 % 'W' = the same units as the filter spacing (see (1) above) 0031 % 'E' = erb scale (aka erb-rate scale) 0032 % 'B' = bark scale 0033 % 'M' = mel scale 0034 % 'L' = log10 scale 0035 % 'H' = Hz 0036 % 0037 % (4) Number of filters 0038 % 'n' = n input gives number of filters [default if n>=1] 0039 % 'N' = n input gives filter spacing [default if n<1] 0040 % 0041 % (5) Filter alignment 0042 % 'k' = force a filter at 1kHz 0043 % 0044 % (6) Output filter 0045 % 'c' = Use complex gammatone filters whose impulse responses uses a complex exponential 0046 % instead of a cosine; the filter order is k instead of 2k. Take the real part of 0047 % the filter output to obtain the same signal as using the real-valued filter. 0048 % The imaginary part is approximately the Hilbert transform of the real part and 0049 % so the magnitude gives the envelope. 0050 % 'q' = biquad filter; b(k,6,n) has the sos coefficients in the same form as sosfilt.m. 0051 % Use y=sosfilt(b,x) to apply the filter. 0052 % 0053 % (6) Future possible options 0054 % ['a' = use all-pole gammtone funtion: see [1]] 0055 % ['s' = use Slaney gammatone approximation: see [2]] 0056 % ['z' = use one-zero gammatone function: see [1]] 0057 % 0058 % (7) Graph plotting 0059 % 'g' = plot filter responses [default if no output arguments present] 0060 % 'G' = plot frequency responses on a log axis 0061 % fc frequency range or, if n=0, list of centre frequencies [default: [100 6000] ] 0062 % bw bandwidth for all filters or a vector of bandwidths: one per filter [default = 1.019 erb ] 0063 % ph phase of all gammatone impulse repsonses or a vector of phases: one per filter 0064 % [default is chosen to give zero gain at DC and a positive real part for the gain at low frequencies] 0065 % k gammatone filter order; the filters have order 2k [default: 4] 0066 % 0067 % Outputs: 0068 % b/a filter coefficients: one filter per row. The gain of each 0069 % filter is scaled to have unit magnitude at the centre frequency. 0070 % Dimensions are b(n,k+max(1,k-1)) and a(n,2*k+1) normally or else 0071 % b(n,max(1,k-1)) and a(n,k+1) if the 'c' option is specified. 0072 % If the q option is given then b(k,6,n) contains the biquad coefficients for filter n. 0073 % fx,bx centre frequencies and bandwidths in Hz 0074 % gd group delay at the centre frequencies (in samples) 0075 % ph phase of each gammatone impulse response 0076 % 0077 0078 % For k>1, the impulse response of filter i is proportional to: 0079 % h[n]=(((n+1)/fs).^(k-1))*cos(2*pi*fx(i)*(n+1)/fs+ph(i))*exp(-2*pi*bx(i)*(n+1)/fs) 0080 % where n>=0. For k=1 replace "(n+1)" by "n". If the 'c' option is used, replace "cos(x)" by "exp(1i*x)". 0081 % Note that the DC gain is only equal to zero for two particular value of ph(i) that differ by pi. 0082 % Each filter is normalized to have unity gain at the centre frequency, fx(i). 0083 % 0084 % Derivation: 0085 % The cos() term in h[n] can be decomposed as a sum of complex exponentials resulting in 0086 % h[n] = a*g[n+1;p] + a'*g[n+1;p]' where p=2*pi*(-bx(i) + j*fx(i))/fs, a=fs^(1-k)*exp(j ph(i)) 0087 % and g[n;p]=n^(k-1)*p^n. 0088 % The z-transform, G(z)=B(z)/A(z) where A(z)=(1 - p/z)^k. The numerator coefficents, b[r], can 0089 % be obtained by convolving a[r] and g[r] where a[r]=nchoosek(k,r)*(-p)^n. 0090 % 0091 % References 0092 % [1] R. F. Lyon, A. G. Katsiamis, and E. M. Drakakis. 0093 % History and future of auditory filter models. 0094 % In Proc Intl Symp Circuits and Systems, pages 3809�3812, 2010. 0095 % doi: 10.1109/ISCAS.2010.5537724. 0096 % [2] M. Slaney. 0097 % An efficient implementation of the patterson-holdsworth auditory filter bank. 0098 % Technical report, Apple Computer, Perception Group, Tech. Rep, 1993. 0099 0100 % Copyright (C) Mike Brookes 2009-2017 0101 % Version: $Id: v_gammabank.m 10865 2018-09-21 17:22:45Z dmb $ 0102 % 0103 % VOICEBOX is a MATLAB toolbox for speech processing. 0104 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0105 % 0106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0107 % This program is free software; you can redistribute it and/or modify 0108 % it under the terms of the GNU General Public License as published by 0109 % the Free Software Foundation; either version 2 of the License, or 0110 % (at your option) any later version. 0111 % 0112 % This program is distributed in the hope that it will be useful, 0113 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0114 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0115 % GNU General Public License for more details. 0116 % 0117 % You can obtain a copy of the GNU General Public License from 0118 % http://www.gnu.org/copyleft/gpl.html or by writing to 0119 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0120 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0121 0122 if nargin<7 0123 k=[]; 0124 if nargin<6 0125 ph=[]; 0126 if nargin<5 0127 bw=[]; 0128 if nargin<4 0129 fc=[]; 0130 if nargin<3 0131 w=''; 0132 end 0133 end 0134 end 0135 end 0136 end 0137 fx=fc(:); 0138 bx=bw(:); 0139 if ~numel(k) 0140 k=4; 0141 end 0142 if ~numel(fx) 0143 fx=[100; 6000]; % default 0144 end 0145 wr='e'; % default frequency warping is erb 0146 for i=1:length(w) 0147 if any(w(i)=='bmlef'); 0148 wr=w(i); 0149 end 0150 end 0151 if any(w=='k') 0152 fk=1000; 0153 switch wr % convert 1kHz to spacing units 0154 case 'b' 0155 fk=v_frq2bark(fk); 0156 case 'l' 0157 fk=log10(fk); 0158 case 'e' 0159 fk=v_frq2erb(fk); 0160 end 0161 else 0162 fk=0; 0163 end 0164 if any(w=='W') 0165 wb=wr; 0166 else 0167 wb='h'; % default bandwidth units are Hz 0168 end 0169 for i=1:length(w) 0170 if any(w(i)=='BMLEF'); 0171 wb=w(i)+'a'-'A'; % convert to lower case 0172 end 0173 end 0174 if ~numel(bx) 0175 bx=1.019; 0176 wb='e'; 0177 end 0178 if any(w=='F') % convert centre frequencies to Hz 0179 switch wr 0180 case 'b' 0181 fx=v_bark2frq(fx); 0182 case 'm' 0183 fx=v_mel2frq(fx); 0184 case 'l' 0185 fx=10.^(fx); 0186 case 'e' 0187 fx=v_erb2frq(fx); 0188 end 0189 end 0190 0191 % now sort out the centre frequencies 0192 0193 if n>0 % n>0: filter end points specified 0194 bx=bx(1); % only use the first bx element 0195 if n==1 % only one filter requested 0196 fx=fx(1); % just use the first frequency 0197 else 0198 switch wr % convert end frequencies to spacing units 0199 case 'b' 0200 fx=v_frq2bark(fx); 0201 case 'm' 0202 fx=v_frq2mel(fx); 0203 case 'l' 0204 fx=log10(fx); 0205 case 'e' 0206 fx=v_frq2erb(fx); 0207 end 0208 if n<1 || any(w=='N') % n = filter spacing 0209 if fk % force filter to 1 kHz 0210 f0=fk-n*floor((fk-fx(1))/n); 0211 else % centre filters in range 0212 f0=(fx(2)+fx(1)-n*floor((fx(2)-fx(1))/n))/2; 0213 end 0214 fx=(f0:n:fx(2))'; 0215 0216 else % n = number of filters specified 0217 % Multiple filters - evenly spaced 0218 fx=linspace(fx(1),fx(2),n)'; % centre frequencies in spacing units 0219 if fk % force a filter at 1kHz 0220 ik=1+ceil((fk-fx(1))*(n-1)/(fx(n)-fx(1))); % index of centre freq immediately above 1 kHz 0221 if ik>n || ik>1 && ((fk-fx(1))*(fx(n)-fx(ik-1))>(fx(n)-fk)*(fx(ik)-fx(1))) 0222 fx=fx(1)+(fx-fx(1))*(fk-fx(1))/(fx(ik)-fx(1)); 0223 else 0224 fx=fx(n)+(fx-fx(n))*(fx(n)-fk)/(fx(n)-fx(ik-1)); 0225 end 0226 end 0227 end 0228 switch wr % convert back to Hz 0229 case 'b' 0230 fx=v_bark2frq(fx); 0231 case 'm' 0232 fx=v_mel2frq(fx); 0233 case 'l' 0234 fx=10.^(fx); 0235 case 'e' 0236 fx=v_erb2frq(fx); 0237 end 0238 end 0239 0240 end 0241 % now sort out the bandwidths 0242 nf=numel(fx); 0243 if numel(bx)==1 0244 bx=bx(ones(nf,1)); % replicate if necessary 0245 end 0246 switch wb % convert bandwidth to Hz 0247 case 'b' 0248 [dum,bwf]=v_frq2bark(fx); 0249 case 'm' 0250 [dum,bwf]=v_frq2mel(fx); 0251 case 'l' 0252 bwf=fx*log(10); 0253 case 'e' 0254 [dum,bwf]=v_frq2erb(fx); 0255 case 'h' 0256 bwf=ones(nf,1); 0257 end 0258 bx=bx.*bwf; 0259 zph=0; % flag to indicate missing phase vector 0260 if ~numel(ph) 0261 ph=zeros(nf,1); % default phase is zero 0262 zph=1; % set missing phase flag 0263 elseif numel(ph)==1 0264 ph=ph(ones(nf,1)); % replicate a scalar value 0265 else 0266 ph=ph(:); % force ph to be a column vector 0267 end 0268 % 0269 % t=(0:ceil(10*fs/(2*pi*bnd)))/fs; % five time constants 0270 % gt=t.^(n-1).*cos(2*pi*cfr*t+phi).*exp(-2*pi*bnd*t); 0271 % gt=gt/sqrt(mean(gt.^2)); % normalize 0272 % figure(1); 0273 % plot(t,gt); 0274 % title('Desired Impulse response'); 0275 % xlabel(['Time (' v_xticksi 's)']); 0276 % 0277 zp=exp((1i*fx-bx)*2*pi/fs); % pole position in top half of z-plane 0278 a=round([1 cumprod((-k:-1)./(1:k))]); % create alternating sign binomial coefficients [1,k+1] 0279 wwp=repmat(zp,1,k+1).^repmat(0:k,nf,1); % powers of pole position [nf,k+1] 0280 denc=repmat(a,nf,1).*wwp; % complex denominator coefficients [nf,k+1] 0281 if k>1 0282 b=conv(a,(1:k-1).^(k-1)); % convolve with integers to the power k-1 [1,2*k-1] 0283 b=exp(1i*ph)*b(1:k-1); % correct for starting phase shift [nf,k-1] 0284 numc=b.*wwp(:,2:k); % complex numerator coefficients [nf,k-1] 0285 else % for the special case of k=1 0286 numc=exp(1i*ph); % complex numerator coefficients [nf,1] 0287 end 0288 kk=size(numc,2); % = max(1,k-1) 0289 % if phase is unspecified choose phase to give zero gain at DC and positive real(gain) for small w 0290 if zph 0291 sd=sum(denc,2); 0292 sn=sum(numc,2); 0293 snn=numc*(0:kk-1)'; 0294 sg=sign(real(conj(sn).*snn)); 0295 ph=angle(1i*sg.*conj(sn./sd)); 0296 numc=numc.*repmat(exp(1i*ph),1,kk); 0297 end 0298 b=zeros(nf,k+kk); % space for numerators 0299 a=zeros(nf,2*k+1); % space for denominators 0300 gd=zeros(nf,1); % space for group delay 0301 ww=exp(2i*fx*pi/fs); % exp(j*centre-freq) 0302 for i=1:nf 0303 b(i,:)=real(conv(numc(i,:),conj(denc(i,:)))); 0304 a(i,:)=real(conv(denc(i,:),conj(denc(i,:)))); % denominator has k repeated poles at ww and ww' 0305 u=polyval(b(i,:),ww(i)); % numerator gain @ centre freqs 0306 v=polyval(a(i,:),ww(i)); % denominator gain @ centre freqs 0307 ud=polyval(b(i,:).*(0:k+kk-1),ww(i)); 0308 vd=polyval(a(i,:).*(0:2*k),ww(i)); 0309 vu=abs(v/u); % scale factor to give unity gain 0310 b(i,:)=b(i,:)*vu; % force unity gain @ centre freqs 0311 numc(i,:)=numc(i,:)*vu; % scale the complex numerator coefficients also 0312 gd(i)=real((v*ud-u*vd)/(u*v)); % group delay at centre freq in samples 0313 end 0314 0315 % now plot graph 0316 0317 if ~nargout || any(w=='g') || any(w=='G') 0318 ng=200; %number of points to plot 0319 if any(w=='G') 0320 fax=logspace(log10(fx(1)/4),log10(fs/2),ng); 0321 else 0322 fax=linspace(0,fs/2,ng); 0323 end 0324 ww=exp(2i*pi*fax/fs); 0325 gg=zeros(nf,ng); 0326 for i=1:nf 0327 gg(i,:)=10*log10(abs(polyval(b(i,:),ww)./polyval(a(i,:),ww))); 0328 end 0329 if any(w=='G') 0330 semilogx(fax,gg','-b'); 0331 set(gca,'xlim',[fax(1) fax(end)]); 0332 else 0333 plot(fax,gg','-b'); 0334 end 0335 xlabel(['Frequency (' v_xticksi 'Hz)']); 0336 set(gca,'ylim',[-50 1]); 0337 title(sprintf('Order-%d Gammatone Filterbank (N=%d, Opt=%s)',k,nf,w)); 0338 end 0339 0340 % sort out output format 0341 0342 if any(w=='c') % if complex filter output 0343 b=numc; % output complex filter of order k 0344 a=denc; 0345 elseif any(w=='q') % if biquad output 0346 bb=zeros(k,6,nf); % space for sos 0347 for i=1:nf 0348 bb(:,:,i)=tf2sos(b(i,:),a(i,:)); % calculate the biquads 0349 bb(:,5,i)=-2*real(zp(i)); % force denominators to be exact 0350 bb(:,6,i)=zp(i)*conj(zp(i)); 0351 end 0352 a=[]; 0353 b=bb; 0354 end