Filters each channel of a multi channel signal through a resonator bank
hfunc = resonatorBanks(select)
returns a function handle to one of the two subfunctions for
implementing comb filters or reson filters.
h = resonatorBanks('reson') returns a handle to the reson subfunction.
h = resonatorBands('comb') returns a handle to the comb subfunction.
SUBFUNCTION RESON:
[wvfOut,resonatorFreqs,BList,AList,QList] = reson(wvfIn,params)
INPUTS
wvfIn: waveform used as input to the filters. Expects a matrix
where the rows are individual band signals and time is the x-axis.
params.freqLimits: Two element vector. Frequency limits of reson filters
params.resonFreqSpacing: 'log2' uses log base 2 spacing.
params.numFilters determines the
number of filters
'linear' uses linear spacing. The
number of filters are determined by params.spacingHz.
'filterQ' spaces filters apart by half
a bandwidth. Bandwidth is determined by
the filters' Q factor.
params.numFilters: Specifies the number of filters to
use. This parameter is only used for 'log2' spacing.
params.spacingHz: Used for spacing the filters when using
linear spacing.
params.numPeriodDecay: Filter Q-factor (also the number of
periods the filters decay by approx. -27 dB, see BTB paper)
params.gainType: 'constant'
'weighted'
'normalized_response'
params.gain: gain value if using constant gain.
params.gainFunction.method: determines shape of rise of gain function
across resonators. Currently
params.gainFunction.minGain: minimum gain value.
params.gainFunction.maxGain: maximum gain value.
params.gainFunction.maxFreq: frequency at which gain rise
reaches the maximum gain value.
params.Fs: sampling rate of input signal
(also sampling rate of output signal)
params.bw_measurement: 'constant' - uses a constant
bandwidth (in Hz) across all filters
'halfPowerBandwidth' - determines filter bandwidth so that filter magnitude
meet at their half power points.
'periodBasedDecay' - A constant Q factor across all the filters is
used to determine their bandwidths.
params.spacingFactor: when params.resonFreqSpacing is set
to 'filterQ', this is multiplied by the bandwidth of the
filters to determine the spacing of filters.
params.bandwidth: when using a constant bandwidth in Hz (params.bw_measurement = 'constant'),
this is the bandwidth used in Hz.
Copyright (c) 2006 The Regents of the University of California
All Rights Reserved
Author(s):
Stefan Tomic 12/060001 function hfunc = resonatorBanks(select) 0002 % Filters each channel of a multi channel signal through a resonator bank 0003 % 0004 % hfunc = resonatorBanks(select) 0005 % 0006 % returns a function handle to one of the two subfunctions for 0007 % implementing comb filters or reson filters. 0008 % h = resonatorBanks('reson') returns a handle to the reson subfunction. 0009 % h = resonatorBands('comb') returns a handle to the comb subfunction. 0010 % 0011 % SUBFUNCTION RESON: 0012 % [wvfOut,resonatorFreqs,BList,AList,QList] = reson(wvfIn,params) 0013 % 0014 % INPUTS 0015 % wvfIn: waveform used as input to the filters. Expects a matrix 0016 % where the rows are individual band signals and time is the x-axis. 0017 % params.freqLimits: Two element vector. Frequency limits of reson filters 0018 % params.resonFreqSpacing: 'log2' uses log base 2 spacing. 0019 % params.numFilters determines the 0020 % number of filters 0021 % 'linear' uses linear spacing. The 0022 % number of filters are determined by params.spacingHz. 0023 % 'filterQ' spaces filters apart by half 0024 % a bandwidth. Bandwidth is determined by 0025 % the filters' Q factor. 0026 % params.numFilters: Specifies the number of filters to 0027 % use. This parameter is only used for 'log2' spacing. 0028 % params.spacingHz: Used for spacing the filters when using 0029 % linear spacing. 0030 % params.numPeriodDecay: Filter Q-factor (also the number of 0031 % periods the filters decay by approx. -27 dB, see BTB paper) 0032 % params.gainType: 'constant' 0033 % 'weighted' 0034 % 'normalized_response' 0035 % params.gain: gain value if using constant gain. 0036 % params.gainFunction.method: determines shape of rise of gain function 0037 % across resonators. Currently 0038 % params.gainFunction.minGain: minimum gain value. 0039 % params.gainFunction.maxGain: maximum gain value. 0040 % params.gainFunction.maxFreq: frequency at which gain rise 0041 % reaches the maximum gain value. 0042 % params.Fs: sampling rate of input signal 0043 % (also sampling rate of output signal) 0044 % params.bw_measurement: 'constant' - uses a constant 0045 % bandwidth (in Hz) across all filters 0046 % 'halfPowerBandwidth' - determines filter bandwidth so that filter magnitude 0047 % meet at their half power points. 0048 % 'periodBasedDecay' - A constant Q factor across all the filters is 0049 % used to determine their bandwidths. 0050 % params.spacingFactor: when params.resonFreqSpacing is set 0051 % to 'filterQ', this is multiplied by the bandwidth of the 0052 % filters to determine the spacing of filters. 0053 % params.bandwidth: when using a constant bandwidth in Hz (params.bw_measurement = 'constant'), 0054 % this is the bandwidth used in Hz. 0055 % 0056 % Copyright (c) 2006 The Regents of the University of California 0057 % All Rights Reserved 0058 % 0059 % Author(s): 0060 % Stefan Tomic 12/06 0061 0062 switch select 0063 case {'reson','reson_r','reson_z'} 0064 hfunc=@reson; 0065 case 'comb' 0066 hfunc=@comb; 0067 case 'partialHanning' 0068 hfunc=@partialHanning; 0069 case 'betaDistribution' 0070 hfunc=@betaDistribution; 0071 otherwise 0072 hfunc=[]; 0073 error('Unknown Input'); 0074 end 0075 return 0076 0077 function [wvfOut,resonatorFreqs,BList,AList,QList,RList] = reson(wvfIn,params) 0078 disp('Filtering signal(s) through reson filterbank(s)...'); 0079 0080 freqLimits = params.freqLimits; 0081 0082 switch(params.resonFreqSpacing) 0083 case 'log2' 0084 numFilters = params.numFilters; 0085 resonatorFreqs = logspace2(log2(freqLimits(1)),log2(freqLimits(2)),numFilters); 0086 0087 case 'linear' 0088 spacingHz = params.spacingHz; 0089 resonatorFreqs = [freqLimits(1):spacingHz:freqLimits(2)]; 0090 numFilters = length(resonatorFreqs); 0091 0092 case 'filterQ' 0093 Q = params.numPeriodDecay; 0094 freqIdx = 1; 0095 resonatorFreqs(freqIdx) = freqLimits(1); 0096 while(resonatorFreqs(freqIdx) <= freqLimits(2)) 0097 freqIdx = freqIdx + 1; 0098 lastBW = resonatorFreqs(freqIdx-1)./Q; 0099 resonatorFreqs(freqIdx) = resonatorFreqs(freqIdx-1) + params.spacingFactor*lastBW; 0100 end 0101 numFilters = length(resonatorFreqs); 0102 0103 end 0104 0105 if(ismember(params.gainType,{'weighted','prenorm_weighted'})) 0106 gainFunc = str2func(params.gainFunction.method); 0107 params.gainFunction.resonatorFreqs = resonatorFreqs; 0108 gainVector = gainFunc(params.gainFunction); 0109 end 0110 0111 if(~isfield(params,'bw_measurement')) 0112 params.bw_measurement = 'constant'; 0113 end 0114 0115 for ibank = 1:size(wvfIn,1) 0116 disp(['channel ' num2str(ibank) ' ']); 0117 for iReson = 1:numFilters 0118 0119 switch(params.bw_measurement) 0120 case 'constant' 0121 bw = params.bandwidth; 0122 0123 case 'halfPowerBandwidth' 0124 if(iReson == numFilters) 0125 percentBw = (resonatorFreqs(iReson) - resonatorFreqs(iReson-1))./resonatorFreqs(iReson-1); 0126 bw = percentBw*resonatorFreqs(iReson); 0127 else 0128 bw = (resonatorFreqs(iReson+1) - resonatorFreqs(iReson)); 0129 end 0130 0131 case 'periodBasedDecay' 0132 Q = params.numPeriodDecay; 0133 bw = resonatorFreqs(iReson)./Q; 0134 end 0135 0136 R = 1 - bw./params.Fs*pi; 0137 %theta = resonatorFreqs(iReson)./params.Fs*(2*pi); 0138 psi = resonatorFreqs(iReson)./params.Fs*(2*pi); 0139 0140 switch(params.gainType) 0141 case 'constant' 0142 gainFactor = params.gain; 0143 case 'weighted' 0144 gainFactor = gainVector(iReson); 0145 case 'prenorm_weighted' 0146 %normalizing before multiplying by beta gain 0147 gainFactor = (1-R^2)/2.*gainVector(iReson); 0148 case 'normalized_response' 0149 switch(params.method) 0150 %normalization method depends on the type of reson filter used 0151 case 'reson_r' 0152 gainFactor = 1-R; 0153 case 'reson_z' 0154 gainFactor = (1-R^2)/2; 0155 case 'reson' 0156 gainFactor = (1 - R^2)*sin(theta); 0157 end 0158 end 0159 0160 % approximating actual peak magnitude response from theta 0161 % (Steiglitz, 1994) 0162 0163 switch(params.method) 0164 case 'reson_z' 0165 B_unscaled = [1 0 -1]; 0166 theta = acos((1+R^2)./(2*R).*cos(psi)); 0167 case 'reson' 0168 B_unscaled = [1]; 0169 theta = acos((2*R)./(1+R^2).*cos(psi)); 0170 case 'reson_r' %we have no peak response correction equation for reson_r 0171 theta = psi; 0172 B_unscaled = [1 0 -R]; 0173 end 0174 0175 realFreq = psi./(2*pi).*params.Fs; 0176 resonatorFreqs(iReson) = realFreq; 0177 0178 B = B_unscaled.*gainFactor; 0179 A = [1 -2*R*cos(theta) R^2]; 0180 0181 wvfOut(ibank,iReson,:) = filter(B,A,wvfIn(ibank,:)); 0182 0183 %the filter poles, zeros, and Q-factors can be returned for testing/tuning 0184 BList{iReson} = B; 0185 AList{iReson} = A; 0186 QList(iReson) = resonatorFreqs(iReson)./bw; 0187 RList(iReson) = R; 0188 end 0189 end 0190 return 0191 0192 0193 0194 function [wvfOut,resonatorFreqs,BList,AList,QList] = comb(wvfIn,params) 0195 disp('Filtering signal(s) through comb filterbank(s)...'); 0196 0197 freqLimits = params.freqLimits; 0198 numFilters = params.numFilters; 0199 resonatorFreqs = logspace2(log2(freqLimits(1)),log2(freqLimits(2)),numFilters); 0200 0201 tSamps = params.halfEnergyTimeSecs*params.Fs; 0202 0203 resonatorPeriodSecs = 1./resonatorFreqs; 0204 resonatorPeriodSamps = floor(resonatorPeriodSecs*params.Fs); 0205 0206 %calculation for alpha (gain) 0207 alpha = 0.5.^(resonatorPeriodSamps./tSamps); 0208 0209 if(strcmp(params.gainType,'weighted')) 0210 gainFunc = str2func(params.gainFunction.method); 0211 params.gainFunction.resonatorFreqs = resonatorFreqs; 0212 gainVector = gainFunc(params.gainFunction); 0213 end 0214 0215 for ibank = 1:size(wvfIn,1) 0216 disp(['channel ' num2str(ibank) ' ']); 0217 for iComb = 1:numFilters 0218 A = [1 zeros(1,resonatorPeriodSamps(iComb)-1) -1*alpha(iComb)]; 0219 B = 1-alpha(iComb); 0220 0221 if(strcmp(params.gainType,'weighted')) 0222 B = B.*gainVector(iComb); 0223 end 0224 0225 wvfOut(ibank,iComb,:) = filter(B,A,wvfIn(ibank,:)); 0226 0227 %the filter poles and zeros can be returned for testing/tuning 0228 BList{iComb} = B; 0229 AList{iComb} = A; 0230 QList(iComb) = NaN; 0231 end 0232 end 0233 0234 return 0235 0236 0237 function gainVector = partialHanning(params) 0238 %weighting function for resonator gains 0239 gainVector = zeros(size(params.resonatorFreqs)); 0240 gainVector(find(params.resonatorFreqs >= params.maxFreq)) = params.maxGain; 0241 gainCurveLength = sum(params.resonatorFreqs < params.maxFreq); 0242 hanningCurve = hanning(gainCurveLength*2); 0243 gainVector(find(params.resonatorFreqs < params.maxFreq)) = params.minGain + (params.maxGain-params.minGain).*hanningCurve(1:gainCurveLength); 0244 0245 return 0246 0247 function gainVector = betaDistribution(params) 0248 0249 freqRange = params.resonatorFreqs - min(params.resonatorFreqs); 0250 xvals = freqRange./max(freqRange); 0251 betaDist = betapdf(xvals,params.betaDistribution.alpha,params.betaDistribution.beta); 0252 0253 if(isfield(params.betaDistribution,'flattenAt') && ~isempty(params.betaDistribution.flattenAt)) 0254 [smallestDiff,closestFreqIdx] = min(abs(params.resonatorFreqs - params.betaDistribution.flattenAt)); 0255 scaleFactor = (params.maxGain - params.minGain)./betaDist(closestFreqIdx); 0256 gainVector = betaDist.*scaleFactor + params.minGain; 0257 gainVector(find(gainVector > params.maxGain)) = params.maxGain; 0258 else 0259 gainVector = betaDist./max(betaDist).*(params.maxGain - params.minGain) + params.minGain; 0260 end 0261 0262 return