+
−% 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
+
−%
+
−%
+
−% <a href="matlab:utils.helper.displayMethodInfo('ao', 'firwhiten')">Parameters Description</a>
+
−%
+
−% VERSION:     $Id: firwhiten.m,v 1.29 2011/11/11 15:21:19 luigi Exp $
+
−%
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−
+
−function varargout = firwhiten(varargin)
+
−
+
−  callerIsMethod = utils.helper.callerIsMethod;
+
−
+
−  % Check if this is a call for parameters
+
−  if utils.helper.isinfocall(varargin{:})
+
−    varargout{1} = getInfo(varargin{3});
+
−    return
+
−  end
+
−
+
−  import utils.const.*
+
−  utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
+
−
+
−  % Collect input variable names
+
−  in_names = cell(size(varargin));
+
−  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+
−
+
−  % Collect all AOs and plists
+
−  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+
−  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+
−
+
−  % Decide on a deep copy or a modify
+
−  bs = copy(as, nargout);
+
−  inhists = copy([as.hist],1);
+
−
+
−  % combine plists
+
−  pl = parse(pl, getDefaultPlist());
+
−
+
−  % Extract necessary parameters
+
−  iNfft   = find(pl, 'Nfft');
+
−  bw      = find(pl, 'bw');
+
−  hc      = find(pl, 'hc');
+
−  swin    = find(pl, 'win');
+
−  order   = find(pl, 'order');
+
−  fwin    = find(pl, 'FIRwin');
+
−  Ntaps   = find(pl, 'Ntaps');
+
−
+
−  % Loop over input AOs
+
−  filts    = [];
+
−  nfs      = [];
+
−
+
−  for j=1:numel(bs)
+
−    if ~isa(bs(j).data, 'tsdata')
+
−      warning('!!! %s expects ao/tsdata objects. Skipping AO %s', mfilename, ao_invars{j});
+
−      bs(j) = [];
+
−    else
+
−      % get Nfft
+
−      if iNfft < 0 || isempty(iNfft)
+
−        Nfft = length(bs(j).data.y);
+
−      else
+
−        Nfft = iNfft;
+
−      end
+
−      utils.helper.msg(msg.PROC1, 'building spectrum');
+
−      % Make spectrum
+
−      axx = psd(bs(j), plist('Nfft', Nfft, 'Win', swin, 'Order', order, 'Scale', 'ASD'));
+
−      % make noise floor estimate
+
−      utils.helper.msg(msg.PROC1, 'estimating noise-floor');
+
−      nxx = smoother(axx, plist('width', bw, 'hc', hc));
+
−      % collect noise-floor estimates for output
+
−      nfs = [nfs nxx];
+
−      % invert and make weights
+
−      w = 1./nxx;
+
−      % Make mfir object
+
−      utils.helper.msg(msg.PROC1, 'building filter');
+
−      ff = mfir(w, plist('Win', fwin, 'N', Ntaps));
+
−      % collect filters for output
+
−      filts = [filts ff];
+
−      % Filter data
+
−      utils.helper.msg(msg.PROC1, 'filter data');
+
−      filter(bs(j), ff);
+
−      % Set name
+
−      bs(j).name = sprintf('firwhiten(%s)', ao_invars{j});
+
−      % add history
+
−      if ~callerIsMethod
+
−        bs(j).addHistory(getInfo('None'), pl, ao_invars(j), inhists(j));
+
−      end
+
−      % clear errors
+
−      bs(j).clearErrors;
+
−      
+
−    end
+
−  end
+
−  
+
−  % Any errors are meaningless after this process, so clear them on both
+
−  % axes.
+
−  bs.clearErrors(plist('axis', 'xy'));
+
−
+
−  % Set outputs
+
−  if nargout > 0
+
−    varargout{1} = bs;
+
−  end
+
−  if nargout > 1
+
−    varargout{2} = filts;
+
−  end
+
−  if nargout > 2
+
−    varargout{3} = nfs;
+
−  end
+
−end
+
−
+
−%--------------------------------------------------------------------------
+
−% Get Info Object
+
−%--------------------------------------------------------------------------
+
−function ii = getInfo(varargin)
+
−  if nargin == 1 && strcmpi(varargin{1}, 'None')
+
−    sets = {};
+
−    pl   = [];
+
−  else
+
−    sets = {'Default'};
+
−    pl   = getDefaultPlist;
+
−  end
+
−  % Build info object
+
−  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: firwhiten.m,v 1.29 2011/11/11 15:21:19 luigi Exp $', sets, pl);
+
−end
+
−
+
−%--------------------------------------------------------------------------
+
−% Get Default Plist
+
−%--------------------------------------------------------------------------
+
−function plout = getDefaultPlist()
+
−  persistent pl;  
+
−  if exist('pl', 'var')==0 || isempty(pl)
+
−    pl = buildplist();
+
−  end
+
−  plout = pl;  
+
−end
+
−
+
−function pl = buildplist()
+
−  
+
−  pl = plist();
+
−  
+
−  % Nfft
+
−  p = param({'Nfft', ['The number of points in the FFT used to estimate<br>'...
+
−                      'the power spectrum. If unspecified, this is calculated as Ndata/4.']}, paramValue.DOUBLE_VALUE(-1));
+
−  pl.append(p);
+
−  
+
−  % BW
+
−  p = param({'bw', ['The bandwidth of the running median filter used to<br>'...
+
−                    'estimate the noise-floor.']}, {1, {20}, paramValue.OPTIONAL});
+
−  pl.append(p);                
+
−                  
+
−  % HC
+
−  p = param({'hc', 'The cutoff used to reject outliers (0-1).'}, {1, {0.8}, paramValue.OPTIONAL});
+
−  pl.append(p);
+
−  
+
−  % Win
+
−  p = param({'Win', 'Spectral window used in spectral estimation.'}, paramValue.WINDOW);
+
−  pl.append(p);
+
−  
+
−  % Order
+
−  p = param({'Order',['The order of segment detrending:<ul>', ...
+
−    '<li>-1 - no detrending</li>', ...
+
−    '<li>0 - subtract mean</li>', ...
+
−    '<li>1 - subtract linear fit</li>', ...
+
−    '<li>N - subtract fit of polynomial, order N</li></ul>']}, paramValue.DETREND_ORDER);
+
−  pl.append(p);
+
−  
+
−  % FIR win
+
−  p = param({'FIRwin', 'The window to use in the filter design.'}, paramValue.WINDOW);
+
−  pl.append(p);
+
−  
+
−  % Ntaps
+
−  p = param({'Ntaps', 'The length of the FIR filter to build.'}, {1, {256}, paramValue.OPTIONAL});
+
−  pl.append(p);
+
−  
+
−end
+
−
+
−% 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]
+
−