% 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]