ToolGammatoneFb

computes a gammatone filterbank
> see function MakeERBFilters.m from Slaneys Auditory Toolbox
>
> @param afAudioData: time domain sample data, dimension channels X samples
> @param f_s: sample rate of audio data
> @param iNumBands: number of filter bands
> @param f_low: start frequency
>
> @retval X filtered signal (dimension bands X samples)
> @retval f_c filterbank mid frequencies
 ======================================================================

0001 %computes a gammatone filterbank
0002 %> see function MakeERBFilters.m from Slaneys Auditory Toolbox
0003 %>
0004 %> @param afAudioData: time domain sample data, dimension channels X samples
0005 %> @param f_s: sample rate of audio data
0006 %> @param iNumBands: number of filter bands
0007 %> @param f_low: start frequency
0008 %>
0009 %> @retval X filtered signal (dimension bands X samples)
0010 %> @retval f_c filterbank mid frequencies
0011 % ======================================================================
0012 function [X, f_c] = ToolGammatoneFb(afAudioData, f_s, iNumBands, f_low)
0013     
0014     if (nargin < 4)
0015         f_low = 100;
0016     end
0017     if (nargin < 3)
0018         iNumBands = 20;
0019     end
0020         
0021     %initialization
0022     fEarQ   = 9.26449;                
0023     fBW     = 24.7;
0024     iOrder  = 1;
0025     T       = 1/f_s;
0026     
0027     % compute the mid frequencies
0028     f_c = GetMidFrequencies_I (f_low, f_s/2, iNumBands, fEarQ, fBW);
0029     
0030     % compute the coefficients
0031     [afCoeffB, afCoeffA] = GetCoeffs_I (  f_c, ...
0032         1.019*2*pi*(((f_c/fEarQ).^iOrder + fBW^iOrder).^(1/iOrder)), ...
0033         T);
0034  
0035     % allocate output memory
0036     X = zeros(iNumBands, length(afAudioData));
0037  
0038     % pre-processing: down-mixing
0039     if (size(afAudioData,2)> 1)
0040         afAudioData = mean(afAudioData,2);
0041     end
0042  
0043     % do the (cascaded) filter process
0044     for (k = 1:iNumBands)
0045         X(k,:) = afAudioData;
0046         for (j = 1:4)
0047             X(k,:) = filter(afCoeffB(j,:,k), afCoeffA(j,:,k), X(k,:));
0048         end
0049     end
0050 end
0051 
0052 %> see function ERBSpace.m from Slaneys Auditory Toolbox
0053 function [f_c] = GetMidFrequencies_I (f_low, f_hi, iNumBands, fEarQ, fBW)
0054     f_c = -(fEarQ*fBW) + exp((1:iNumBands)'*(-log(f_hi + fEarQ*fBW) + ...
0055             log(f_low + fEarQ*fBW))/iNumBands) * (f_hi + fEarQ*fBW);
0056 end
0057 
0058 
0059 %> see function MakeERBFilters.m from Slaneys Auditory Toolbox
0060 function [afCoeffB, afCoeffA] = GetCoeffs_I (f_c, B, T)
0061 
0062     fCos    = cos(2*f_c*pi*T);
0063     fSin    = sin(2*f_c*pi*T);
0064     fExp    = exp(B*T);
0065     fSqrtA  = 2*sqrt(3+2^1.5);
0066     fSqrtS  = 2*sqrt(3-2^1.5);
0067 
0068     A0      = T;
0069     A2      = 0;
0070     B0      = 1;
0071     B1      = -2*fCos./fExp;
0072     B2      = exp(-2*B*T);
0073 
0074     A11     = -(2*T*fCos./fExp + fSqrtA*T*fSin./ fExp)/2;
0075     A12     = -(2*T*fCos./fExp - fSqrtA*T*fSin./ fExp)/2;
0076     A13     = -(2*T*fCos./fExp + fSqrtS*T*fSin./ fExp)/2;
0077     A14     = -(2*T*fCos./fExp - fSqrtS*T*fSin./ fExp)/2;
0078 
0079     fSqrtA  = sqrt(3 + 2^(3/2));
0080     fSqrtS  = sqrt(3 - 2^(3/2));
0081     fArg    = i*f_c*pi*T;
0082     
0083     afGain    = ...
0084     abs((-2*exp(4*fArg)*T+2*exp(-(B*T)+2*fArg).*T.*(fCos-fSqrtS*fSin)).* ...
0085         (-2*exp(4*fArg)*T+2*exp(-(B*T)+2*fArg).*T.*(fCos+fSqrtS*fSin)).* ...
0086         (-2*exp(4*fArg)*T+2*exp(-(B*T)+2*fArg).*T.*(fCos-fSqrtA*fSin)).* ...
0087         (-2*exp(4*fArg)*T+2*exp(-(B*T)+2*fArg).*T.*(fCos+fSqrtA*fSin))./ ...
0088         (-2 ./ exp(2*B*T) - 2*exp(4*fArg) + 2*(1 + exp(4*fArg))./fExp).^4);
0089 
0090     % this is Slaneys compact format - now resort into 3D Matrices
0091     %fcoefs = [A0*ones(length(f_c),1) A11 A12 A13 A14 A2*ones(length(f_c),1) B0*ones(length(f_c),1) B1 B2 afGain];
0092     
0093     afCoeffB = zeros(4,3,length(B));    
0094     afCoeffA = zeros(4,3,length(B));
0095     
0096     for (k = 1:length(B))
0097         afCoeffB(1,:,k)     = [A0 A11(k) A2]/afGain(k);
0098         afCoeffA(1,:,k)     = [B0 B1(k) B2(k)];
0099         
0100         afCoeffB(2,:,k)     = [A0 A12(k) A2];
0101         afCoeffA(2,:,k)     = [B0 B1(k) B2(k)];
0102         
0103         afCoeffB(3,:,k)     = [A0 A13(k) A2];
0104         afCoeffA(3,:,k)     = [B0 B1(k) B2(k)];
0105         
0106         afCoeffB(4,:,k)     = [A0 A14(k) A2];
0107         afCoeffA(4,:,k)     = [B0 B1(k) B2(k)];
0108         
0109     end
0110 end