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 function0001 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