Home > classes > @ao > firwhiten.m

firwhiten

PURPOSE ^

FIRWHITEN whitens the input time-series by building an FIR whitening filter.

SYNOPSIS ^

function varargout = firwhiten(varargin)

DESCRIPTION ^

 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003