FIRWHITEN whitens the input time-series by building an FIR whitening filter. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION: FIRWHITEN whitens the input time-series by building an FIR whitening filter. The algorithm ultimately uses fir2() to build the whitening filter. ALGORITHM: 1) Make ASD of time-series 2) Perform running median to get noise-floor estimate (ao/smoother) 3) Invert noise-floor estimate 4) Call mfir() on noise-floor estimate to produce whitening filter 5) Filter data CALL: b = firwhiten(a, pl) % returns whitened time-series AOs [b, filts] = firwhiten(a, pl) % returns the mfir filters used [b, filts, nfs] = firwhiten(a, pl) % returns the noise-floor % estimates as fsdata AOs PARAMETERS: 'Ntaps' - the length of the FIR filter to build [default: 256]. 'FIRwin' - the window to use in the filter design. Pass a specwin object of the desired type and of any length. [default: Hanning] parameters passed to ltpda_pwelch() 'Nfft' - The number of points in the FFT used to estimate the power spectrum. [default: Ndata/4] 'Win' - Spectral window used in spectral estimation. [default: Kaiser -150dB] 'Order' - order of segment detrending: -1 - no detrending 0 - subtract mean [default] 1 - subtract linear fit N - subtract fit of polynomial, order N (Segment overlap is taken from the window function.) parameters passed to ltpda_nfest() 'bw' - The bandwidth of the running median filter used to estimate the noise-floor. [default: 20 samples] 'hc' - The cutoff used to reject outliers (0-1) [default: 0.8] M-FILE INFO: Get information about this methods by calling >> ao.getInfo('firwhiten') Get information about a specified set-plist by calling: >> ao.getInfo('firwhiten', 'None') VERSION: $Id: firwhiten.m,v 1.10 2008/09/08 08:30:12 hewitson Exp $ HISTORY: 21-04-08 M Hewitson Creation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 % FIRWHITEN whitens the input time-series by building an FIR whitening filter. 0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0003 % 0004 % DESCRIPTION: FIRWHITEN whitens the input time-series by building an FIR 0005 % whitening filter. The algorithm ultimately uses fir2() to 0006 % build the whitening filter. 0007 % 0008 % ALGORITHM: 0009 % 1) Make ASD of time-series 0010 % 2) Perform running median to get noise-floor estimate (ao/smoother) 0011 % 3) Invert noise-floor estimate 0012 % 4) Call mfir() on noise-floor estimate to produce whitening filter 0013 % 5) Filter data 0014 % 0015 % CALL: b = firwhiten(a, pl) % returns whitened time-series AOs 0016 % [b, filts] = firwhiten(a, pl) % returns the mfir filters used 0017 % [b, filts, nfs] = firwhiten(a, pl) % returns the noise-floor 0018 % % estimates as fsdata AOs 0019 % 0020 % PARAMETERS: 0021 % 0022 % 'Ntaps' - the length of the FIR filter to build [default: 256]. 0023 % 'FIRwin' - the window to use in the filter design. Pass a 0024 % specwin object of the desired type and of any length. 0025 % [default: Hanning] 0026 % 0027 % parameters passed to ltpda_pwelch() 0028 % 0029 % 'Nfft' - The number of points in the FFT used to estimate 0030 % the power spectrum. 0031 % [default: Ndata/4] 0032 % 'Win' - Spectral window used in spectral estimation. 0033 % [default: Kaiser -150dB] 0034 % 'Order' - order of segment detrending: 0035 % -1 - no detrending 0036 % 0 - subtract mean [default] 0037 % 1 - subtract linear fit 0038 % N - subtract fit of polynomial, order N 0039 % 0040 % (Segment overlap is taken from the window function.) 0041 % 0042 % parameters passed to ltpda_nfest() 0043 % 0044 % 'bw' - The bandwidth of the running median filter used to 0045 % estimate the noise-floor. 0046 % [default: 20 samples] 0047 % 'hc' - The cutoff used to reject outliers (0-1) 0048 % [default: 0.8] 0049 % 0050 % M-FILE INFO: Get information about this methods by calling 0051 % >> ao.getInfo('firwhiten') 0052 % 0053 % Get information about a specified set-plist by calling: 0054 % >> ao.getInfo('firwhiten', 'None') 0055 % 0056 % VERSION: $Id: firwhiten.m,v 1.10 2008/09/08 08:30:12 hewitson Exp $ 0057 % 0058 % HISTORY: 21-04-08 M Hewitson 0059 % Creation 0060 % 0061 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0062 0063 function varargout = firwhiten(varargin) 0064 0065 % Check if this is a call for parameters 0066 if utils.helper.isinfocall(varargin{:}) 0067 varargout{1} = getInfo(varargin{3}); 0068 return 0069 end 0070 0071 import utils.const.* 0072 utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename); 0073 0074 % Collect input variable names 0075 in_names = cell(size(varargin)); 0076 for ii = 1:nargin,in_names{ii} = inputname(ii);end 0077 0078 % Collect all AOs and plists 0079 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); 0080 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names); 0081 0082 % Decide on a deep copy or a modify 0083 bs = copy(as, nargout); 0084 inhists = copy([as.hist],1); 0085 0086 % combine plists 0087 pl = combine(pl, getDefaultPlist()); 0088 0089 % Extract necessary parameters 0090 iNfft = find(pl, 'Nfft'); 0091 bw = find(pl, 'bw'); 0092 hc = find(pl, 'hc'); 0093 swin = find(pl, 'win'); 0094 order = find(pl, 'order'); 0095 fwin = find(pl, 'FIRwin'); 0096 Ntaps = find(pl, 'Ntaps'); 0097 0098 % Loop over input AOs 0099 filts = []; 0100 nfs = []; 0101 0102 for j=1:length(bs) 0103 if ~isa(bs(j).data, 'tsdata') 0104 warning('!!! %s expects ao/tsdata objects. Skipping AO %s', mfilename, ao_invars{j}); 0105 bs(j) = []; 0106 else 0107 % get Nfft 0108 if iNfft < 0 0109 Nfft = length(bs(j).data.y); 0110 else 0111 Nfft = iNfft; 0112 end 0113 utils.helper.msg(msg.PROC1, 'building spectrum'); 0114 % Make spectrum 0115 axx = psd(bs(j), plist('Nfft', Nfft, 'Win', swin, 'Order', order, 'Scale', 'ASD')); 0116 % make noise floor estimate 0117 utils.helper.msg(msg.PROC1, 'estimating noise-floor'); 0118 nxx = smoother(axx, plist('width', bw, 'hc', hc)); 0119 % collect noise-floor estimates for output 0120 nfs = [nfs nxx]; 0121 % invert and make weights 0122 w = 1./nxx; 0123 % Make mfir object 0124 utils.helper.msg(msg.PROC1, 'building filter'); 0125 ff = mfir(w, plist('Win', fwin, 'N', Ntaps)); 0126 % collect filters for output 0127 filts = [filts ff]; 0128 % Filter data 0129 utils.helper.msg(msg.PROC1, 'filter data'); 0130 filter(bs(j), ff); 0131 % Set name 0132 bs(j).setName(sprintf('firwhiten(%s)', ao_invars{j}), 'internal'); 0133 % add history 0134 bs(j).addHistory(getInfo, pl, ao_invars(j), inhists(j)); 0135 end 0136 end 0137 0138 % Set outputs 0139 if nargout > 0 0140 varargout{1} = bs; 0141 end 0142 if nargout > 1 0143 varargout{2} = filts; 0144 end 0145 if nargout > 2 0146 varargout{3} = nfs; 0147 end 0148 end 0149 0150 %-------------------------------------------------------------------------- 0151 % Get Info Object 0152 %-------------------------------------------------------------------------- 0153 function ii = getInfo(varargin) 0154 if nargin == 1 && strcmpi(varargin{1}, 'None') 0155 sets = {}; 0156 pl = []; 0157 else 0158 sets = {'Default'}; 0159 pl = getDefaultPlist; 0160 end 0161 % Build info object 0162 ii = minfo(mfilename, 'ao', '', utils.const.categories.sigproc, '$Id: firwhiten.m,v 1.10 2008/09/08 08:30:12 hewitson Exp $', sets, pl); 0163 end 0164 0165 %-------------------------------------------------------------------------- 0166 % Get Default Plist 0167 %-------------------------------------------------------------------------- 0168 function pl_default = getDefaultPlist() 0169 pl_default = plist(... 0170 'Nfft', -1, ... 0171 'bw', 20, ... 0172 'hc', 0.8, ... 0173 'Win', specwin('Kaiser', 10, 150), ... 0174 'order', 0, ... 0175 'FIRwin', specwin('Hanning', 10), ... 0176 'Ntaps', 256); 0177 end 0178 0179