Home > utils > signal > reson_filter.m

reson_filter

PURPOSE ^

SYNOPSIS ^

function varargout = reson_filter(varargin)

DESCRIPTION ^

 processes a signal through a reson filter

 INPUTS (passed in as tag,value pairs)
   inputSig:      input signal to process with the filter
   resonantFreq:  resonant frequency of the filter
   bandwidth:     half-power bandwidth of the filter in Hz.
   Fs:            sampling rate of the signal
   gain:          constant (between 0 and 1) that is multiplied by
                  numerator of transfer function. Default = 1.
   resonType:     either 'reson' or 'reson_z'. reson_z adds two zeros
                  to the filter in order to improve the low frequency rolloff
                  (for low frequency filters), or high frequency rolloff (for
                  high frequency filters). Default is 'reson_z' if omitted.
   normalize:     logical 1 or 0 for true or false. Default = 1.
   Qfactor:       Ratio of resonantFreq/bandwidth. Either Qfactor or
                  bandwidth may be passed in as arguments. If both are passed in,
                  make sure they are compatible (i.e. Q=Fr/bw)


 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/23/08 Stefan Tomic. Took code out of BTB into a more general
                       purpose function

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function varargout = reson_filter(varargin)
0002 %
0003 % processes a signal through a reson filter
0004 %
0005 % INPUTS (passed in as tag,value pairs)
0006 %   inputSig:      input signal to process with the filter
0007 %   resonantFreq:  resonant frequency of the filter
0008 %   bandwidth:     half-power bandwidth of the filter in Hz.
0009 %   Fs:            sampling rate of the signal
0010 %   gain:          constant (between 0 and 1) that is multiplied by
0011 %                  numerator of transfer function. Default = 1.
0012 %   resonType:     either 'reson' or 'reson_z'. reson_z adds two zeros
0013 %                  to the filter in order to improve the low frequency rolloff
0014 %                  (for low frequency filters), or high frequency rolloff (for
0015 %                  high frequency filters). Default is 'reson_z' if omitted.
0016 %   normalize:     logical 1 or 0 for true or false. Default = 1.
0017 %   Qfactor:       Ratio of resonantFreq/bandwidth. Either Qfactor or
0018 %                  bandwidth may be passed in as arguments. If both are passed in,
0019 %                  make sure they are compatible (i.e. Q=Fr/bw)
0020 %
0021 %
0022 % OUTPUTS
0023 % If inputSig is passed in, outputs are [outputSig,B,A], where B and A are the filter coefficients of
0024 % the numerator and denominator. If inputSig is not passed in, then
0025 % the output simply contains the filter coefficients [B,A].
0026 %
0027 %
0028 % 1/23/08 Stefan Tomic. Took code out of BTB into a more general
0029 %                       purpose function
0030 %
0031  
0032 for iarg = 1:2:nargin
0033   switch(varargin{iarg})
0034    
0035    case 'inputSig'
0036     inSig = varargin{iarg+1};
0037   
0038    case 'resonantFreq'
0039     Fr = varargin{iarg+1};
0040     
0041    case 'bandwidth'
0042     bw = varargin{iarg+1};
0043     
0044    case 'gain'
0045     gain = varargin{iarg+1};
0046     
0047    case 'resonType'
0048     resonType = varargin{iarg+1};
0049     
0050    case 'normalize'
0051     normalize = varargin{iarg+1};
0052     
0053    case 'Fs'
0054     Fs = varargin{iarg+1};
0055     
0056    case 'Qfactor'
0057     Qfactor = varargin{iarg+1};
0058     
0059   end
0060 end
0061 
0062 try
0063   normalize;
0064 catch
0065   normalize = 1;
0066 end
0067 
0068 try
0069   gain;
0070 catch
0071   gain = 1;
0072 end
0073 
0074 if(resonantFreq > Fs/2)
0075   error('Resonant freqency is above Nyquist frequency');
0076 end
0077 
0078 if(gain < 0 || gain > 1)
0079   error('Gain is not between 0 and 1.');
0080 end
0081 
0082 %if Qfactor was passed instead of bw, calculate bw from Qfactor.
0083 %if both Qfactor and bw were passed in, make sure that they are compatible.
0084 %if only bw was passed in, then no need to calculate Qfactor, since
0085 %we only need bw from here on.
0086 if(~exist('bw','var') && exist('Qfactor','var'))
0087   bw = Fr/Qfactor;
0088 elseif(exist('bw','var') && exist('Qfactor','var'))
0089   if(bw ~= Fr/Qfactor)
0090     error('Incompatible bandwidth and Qfactor were passed in');
0091   end
0092 end
0093   
0094 R = 1 - bw./Fs*pi;  
0095 
0096 %psi is the desired resonant frequency. Note that with reson
0097 %filters that this differs slightly from the location of the pole
0098 %angle. So we need to calculate the pole angle theta from the
0099 %desired resonant frequency angle, psi.
0100 psi = Fr/Fs*(2*pi);    
0101 
0102 if(normalize)
0103   %we will set a value for normFactor here, which will normalize
0104   %the filter's magnitude response. The gain (passed in as a param)
0105   %can also be imposed on the normalized response, so that we can
0106   %scale the magnitude response effectively. Since non-normalized
0107   %reson filters' peak magnitude response is dependent on the
0108   %resonant frequency, normalization gives us better control over
0109   %the magnitude response.
0110   
0111   switch(resonType)
0112     %normalization method depends on the type of reson filter used
0113    case 'reson_r'
0114     normFactor = 1-R; 
0115    case 'reson_z'
0116     normFactor = (1-R^2)/2;
0117    case 'reson'
0118     normFactor = (1-R^2)*sin(theta);
0119   end
0120 else
0121   normFactor = 1;
0122 end
0123 
0124 switch(resonType)
0125  case 'reson_z'
0126   B_unscaled = [1 0 -1];
0127   % approximating actual pole angle theta from desired peak response
0128   % (Steiglitz, 1994)
0129 
0130   theta = acos( (1+R^2)./(2*R).*cos(psi) );
0131  case 'reson'
0132   B_unscaled = [1];
0133   theta = acos( (2*R)./(1+R^2).*cos(psi) );
0134  case 'reson_r'
0135   psi = theta; %we have no peak response correction equation for reson_r
0136   B_unscaled = [1 0 -R];
0137 end
0138 
0139 
0140 B = B_unscaled.*gain.*normFactor;
0141 A = [1 -2*R*cos(theta) R^2];
0142 
0143 if(exist('inSig','var'))
0144   varargout{1} = filter(B,A,inSig);   
0145   varargout{2} = B;
0146   varargout{3} = A;
0147 else
0148   varargout{1} = B;
0149   varargout{2} = A;
0150 end

Generated on Wed 20-Sep-2023 04:00:50 by m2html © 2003