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 0139 switch(params.gainType) 0140 case 'constant' 0141 gainFactor = params.gain; 0142 case 'weighted' 0143 gainFactor = gainVector(iReson); 0144 case 'prenorm_weighted' 0145 %normalizing before multiplying by beta gain 0146 gainFactor = (1-R^2)/2.*gainVector(iReson); 0147 case 'normalized_response' 0148 switch(params.method) 0149 %normalization method depends on the type of reson filter used 0150 case 'reson_r' 0151 gainFactor = 1-R; 0152 case 'reson_z' 0153 gainFactor = (1-R^2)/2; 0154 case 'reson' 0155 gainFactor = (1 - R^2)*sin(theta); 0156 end 0157 end 0158 0159 % approximating actual peak magnitude response from theta 0160 % (Steiglitz, 1994) 0161 0162 switch(params.method) 0163 case 'reson_z' 0164 B_unscaled = [1 0 -1]; 0165 psi = acos((2*R)./(1+R^2).*cos(theta)); 0166 case 'reson' 0167 B_unscaled = [1]; 0168 psi = acos((1+R^2)./(2*R)*cos(theta)); 0169 case 'reson_r' 0170 psi = theta; %we have no peak response correction equation for reson_r 0171 B_unscaled = [1 0 -R]; 0172 end 0173 0174 realFreq = psi./(2*pi).*params.Fs; 0175 resonatorFreqs(iReson) = realFreq; 0176 0177 B = B_unscaled.*gainFactor; 0178 A = [1 -2*R*cos(theta) R^2]; 0179 0180 wvfOut(ibank,iReson,:) = filter(B,A,wvfIn(ibank,:)); 0181 0182 %the filter poles, zeros, and Q-factors can be returned for testing/tuning 0183 BList{iReson} = B; 0184 AList{iReson} = A; 0185 QList(iReson) = resonatorFreqs(iReson)./bw; 0186 RList(iReson) = R; 0187 end 0188 end 0189 return 0190 0191 0192 0193 function [wvfOut,resonatorFreqs,BList,AList,QList] = comb(wvfIn,params) 0194 disp('Filtering signal(s) through comb filterbank(s)...'); 0195 0196 freqLimits = params.freqLimits; 0197 numFilters = params.numFilters; 0198 resonatorFreqs = logspace2(log2(freqLimits(1)),log2(freqLimits(2)),numFilters); 0199 0200 tSamps = params.halfEnergyTimeSecs*params.Fs; 0201 0202 resonatorPeriodSecs = 1./resonatorFreqs; 0203 resonatorPeriodSamps = floor(resonatorPeriodSecs*params.Fs); 0204 0205 %calculation for alpha (gain) 0206 alpha = 0.5.^(resonatorPeriodSamps./tSamps); 0207 0208 if(strcmp(params.gainType,'weighted')) 0209 gainFunc = str2func(params.gainFunction.method); 0210 params.gainFunction.resonatorFreqs = resonatorFreqs; 0211 gainVector = gainFunc(params.gainFunction); 0212 end 0213 0214 for ibank = 1:size(wvfIn,1) 0215 disp(['channel ' num2str(ibank) ' ']); 0216 for iComb = 1:numFilters 0217 A = [1 zeros(1,resonatorPeriodSamps(iComb)-1) -1*alpha(iComb)]; 0218 B = 1-alpha(iComb); 0219 0220 if(strcmp(params.gainType,'weighted')) 0221 B = B.*gainVector(iComb); 0222 end 0223 0224 wvfOut(ibank,iComb,:) = filter(B,A,wvfIn(ibank,:)); 0225 0226 %the filter poles and zeros can be returned for testing/tuning 0227 BList{iComb} = B; 0228 AList{iComb} = A; 0229 QList(iComb) = NaN; 0230 end 0231 end 0232 0233 return 0234 0235 0236 function gainVector = partialHanning(params) 0237 %weighting function for resonator gains 0238 gainVector = zeros(size(params.resonatorFreqs)); 0239 gainVector(find(params.resonatorFreqs >= params.maxFreq)) = params.maxGain; 0240 gainCurveLength = sum(params.resonatorFreqs < params.maxFreq); 0241 hanningCurve = hanning(gainCurveLength*2); 0242 gainVector(find(params.resonatorFreqs < params.maxFreq)) = params.minGain + (params.maxGain-params.minGain).*hanningCurve(1:gainCurveLength); 0243 0244 return 0245 0246 function gainVector = betaDistribution(params) 0247 0248 freqRange = params.resonatorFreqs - min(params.resonatorFreqs); 0249 xvals = freqRange./max(freqRange); 0250 betaDist = betapdf(xvals,params.betaDistribution.alpha,params.betaDistribution.beta); 0251 0252 if(isfield(params.betaDistribution,'flattenAt') && ~isempty(params.betaDistribution.flattenAt)) 0253 [smallestDiff,closestFreqIdx] = min(abs(params.resonatorFreqs - params.betaDistribution.flattenAt)); 0254 scaleFactor = (params.maxGain - params.minGain)./betaDist(closestFreqIdx); 0255 gainVector = betaDist.*scaleFactor + params.minGain; 0256 gainVector(find(gainVector > params.maxGain)) = params.maxGain; 0257 else 0258 gainVector = betaDist./max(betaDist).*(params.maxGain - params.minGain) + params.minGain; 0259 end 0260 0261 return