Home > voicebox > gammabank.m

gammabank

PURPOSE ^

GAMMABANK gammatone filter bank [b,a,fx,bx,gd]=(n,fs,w,fc,bw,ph,k)

SYNOPSIS ^

function [b,a,fx,bx,gd,ph]=gammabank(n,fs,w,fc,bw,ph,k)

DESCRIPTION ^

GAMMABANK gammatone filter bank [b,a,fx,bx,gd]=(n,fs,w,fc,bw,ph,k)

 Usage:
          (1) [b,a,fx,bx,gd,ph]=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
        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 teh 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 response

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [b,a,fx,bx,gd,ph]=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,ph]=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 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 teh 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()" by "exp(1i*)".
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 38093812, 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: gammabank.m 10028 2017-08-11 16:39:37Z 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=frq2bark(fk);
0156         case 'l'
0157             fk=log10(fk);
0158         case 'e'
0159             fk=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=bark2frq(fx);
0182         case 'm'
0183             fx=mel2frq(fx);
0184         case 'l'
0185             fx=10.^(fx);
0186         case 'e'
0187             fx=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=frq2bark(fx);
0201             case 'm'
0202                 fx=frq2mel(fx);
0203             case 'l'
0204                 fx=log10(fx);
0205             case 'e'
0206                 fx=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=bark2frq(fx);
0231             case 'm'
0232                 fx=mel2frq(fx);
0233             case 'l'
0234                 fx=10.^(fx);
0235             case 'e'
0236                 fx=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]=frq2bark(fx);
0249     case 'm'
0250         [dum,bwf]=frq2mel(fx);
0251     case 'l'
0252         bwf=fx*log(10);
0253     case 'e'
0254         [dum,bwf]=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 (' 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 (' 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

Generated on Tue 19-Sep-2017 12:07:31 by m2html © 2003