comb_filter

 processes a signal through a comb filter

 INPUTS (passed in as tag,value pairs)
   inputSig:       input signal to process with the filter
   resonantFreq:   lowest resonant frequency of the filter
                   (fundamental frequency of the comb filter)
   Fs:             sampling rate of the signal
   halfEnergyTime: Only configured for feedback comb filter (combType='normal').
                   Tunes the filter so that it has a  half energy
                   decay time specified in seconds (see Scheirer, 1998)
   combType:       either 'normal' or 'inverse'. 'inverse'
                   implements a feedfoward comb filter with a magnitude response
                   that is the inverse of the feedback ('normal') version. Default
                   is 'normal'.
   gain:           constant (between 0 and 1) that is multiplied by
                   numerator of transfer function. Default = 1.
   poleRadius:     If 'halfEnergyTime' is set, this argument is
                   ignored. Radius of the poles for the feedback comb filter (combType='normal').
                   This affects the overall gain of the comb filter.
   zeroRadius:     Radius of the zeros for the feedforward comb filter (combType='inverse')
                   This affects the overall gain of the comb filter.

 OUTPUTS
 If inputSig is passed in, outputs are [outputSig,B,A], where B and A are the filter coefficients of
 the numerator and denominator. If inputSig is not passed in, then
 the output simply contains the filter coefficients [B,A].


 1/26/08 Stefan Tomic. Took code out of BTB into a more general
                       purpose function

0001 function varargout = comb_filter(varargin)
0002 %
0003 % processes a signal through a comb filter
0004 %
0005 % INPUTS (passed in as tag,value pairs)
0006 %   inputSig:       input signal to process with the filter
0007 %   resonantFreq:   lowest resonant frequency of the filter
0008 %                   (fundamental frequency of the comb filter)
0009 %   Fs:             sampling rate of the signal
0010 %   halfEnergyTime: Only configured for feedback comb filter (combType='normal').
0011 %                   Tunes the filter so that it has a  half energy
0012 %                   decay time specified in seconds (see Scheirer, 1998)
0013 %   combType:       either 'normal' or 'inverse'. 'inverse'
0014 %                   implements a feedfoward comb filter with a magnitude response
0015 %                   that is the inverse of the feedback ('normal') version. Default
0016 %                   is 'normal'.
0017 %   gain:           constant (between 0 and 1) that is multiplied by
0018 %                   numerator of transfer function. Default = 1.
0019 %   poleRadius:     If 'halfEnergyTime' is set, this argument is
0020 %                   ignored. Radius of the poles for the feedback comb filter (combType='normal').
0021 %                   This affects the overall gain of the comb filter.
0022 %   zeroRadius:     Radius of the zeros for the feedforward comb filter (combType='inverse')
0023 %                   This affects the overall gain of the comb filter.
0024 %
0025 % OUTPUTS
0026 % If inputSig is passed in, outputs are [outputSig,B,A], where B and A are the filter coefficients of
0027 % the numerator and denominator. If inputSig is not passed in, then
0028 % the output simply contains the filter coefficients [B,A].
0029 %
0030 %
0031 % 1/26/08 Stefan Tomic. Took code out of BTB into a more general
0032 %                       purpose function
0033  
0034 for iarg = 1:2:nargin
0035   switch(varargin{iarg})
0036    
0037    case 'inputSig'
0038     inSig = varargin{iarg+1};
0039   
0040    case 'resonantFreq'
0041     Fr = varargin{iarg+1};
0042     
0043    case 'gain'
0044     gain = varargin{iarg+1};
0045     
0046    case 'halfEnergyTime'
0047     tHalf = varargin{iarg+1};
0048     
0049    case 'combType'
0050     combType = varargin{iarg+1};
0051     
0052    case 'Fs'
0053     Fs = varargin{iarg+1};
0054     
0055    case 'poleRadius'
0056     poleRadius = varargin{iarg+1};
0057     
0058    case 'zeroRadius'
0059     zeroRadius = varargin{iarg+1};
0060     
0061   end
0062 end
0063 
0064 try
0065   combType;
0066 catch
0067   combType = 'normal';
0068 end
0069 
0070 try
0071   gain;
0072 catch
0073   gain = 1;
0074 end
0075      
0076 switch(combType)
0077   
0078   case 'normal'
0079 
0080    if(exist('tHalf','var') && exist('poleRadius','var'))
0081      warning(['Both half energy decay time and pole radius were specified.' ...
0082           ' Ignoring the pole radius value.']);
0083    end
0084 
0085    Tr = 1./Fr;
0086    trSamps = Tr*Fs;
0087 
0088    if(exist('tHalf','var'))
0089      tHalfSamps = tHalf*Fs;
0090      alpha = 0.5.^(trSamps./tHalfSamps);
0091      A = [1 zeros(1,trSamps-1) -1*alpha];
0092      B = 1-alpha;
0093    elseif(exist('poleRadius','var'))
0094      A = [1 zeros(1,trSamps-1) -1*poleRadius^trSamps];
0095      B = gain;
0096    else
0097      error('comb_filter needs either a half energy time or pole radius');
0098    end
0099 
0100  case 'inverse'
0101   
0102   %for the inverse comb filter, the first zero is actually at twice
0103   %the angle of the lowest resonance, so the delay is half of the
0104   %period corresponding to the first resonance.
0105   Tr = 1./Fr;
0106   trSamps = Tr*Fs./2;
0107   
0108   B =  gain.*[1 zeros(1,trSamps-1) -1*zeroRadius^trSamps];
0109   A = 1;
0110   
0111 end
0112     
0113 varargout = {};
0114 if(exist('inSig','var'))
0115   varargout{end+1} = filter(B,A,inSig);
0116 end
0117 
0118 varargout{end+1} = B;
0119 varargout{end+1} = A;
0120     
0121 return