# v_gammabank

## PURPOSE

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

## SYNOPSIS

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

## DESCRIPTION

    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 response

## CROSS-REFERENCE INFORMATION

This function calls:

- v_bark2frq V_BARK2FRQ Convert the BARK frequency scale to Hertz FRQ=(BARK)
- v_erb2frq V_ERB2FRQ Convert ERB frequency scale to Hertz FRQ=(ERB)
- v_frq2bark V_FRQ2BARK Convert Hertz to BARK frequency scale BARK=(FRQ)
- v_frq2erb V_FRQ2ERB Convert Hertz to ERB frequency scale ERB=(FRQ)
- v_frq2mel V_FRQ2ERB Convert Hertz to Mel frequency scale MEL=(FRQ)
- v_mel2frq V_MEL2FRQ Convert Mel frequency scale to Hertz FRQ=(MEL)
- v_xticksi V_XTIXKSI labels the x-axis of a plot using SI multipliers S=(AH)

This function is called by:

## SOURCE CODE

    0001 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 Lesser General Public License as published by
    0109 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0116 %
    0117 %   You can obtain a copy of the GNU Lesser General Public License from
    0118 %   https://www.gnu.org/licenses/ .
    0119 %
    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