% WHITEN1D whitens the input time-series.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: WHITEN1D whitens the input time-series. The filter is built
%              by fitting to the model provided. If no model is provided, a
%              fit is made to a spectral-density estimate of the
%              time-series (made using psd+bin_data or lpsd).
%              Note: The function assumes that the input model corresponds
%              to the one-sided psd of the data to be whitened.
%
% ALGORITHM:
%            1) If no model provided, make psd+bin_data or lpsd
%               of time-series and take it as a model
%               for the data power spectral density
%            2) Fit a set of partial fraction z-domain filters using
%               utils.math.psd2wf. The fit is automatically stopped when
%               the accuracy tolerance is reached.
%            3) Convert to array of MIIR filters
%            4) Assemble into a parallel filterbank object
%            5) Filter time-series in parallel
%
%
% CALL:         b = whiten1D(a, pl)
%               [b1,b2,...,bn] = whiten1D(a1,a2,...,an, pl);
%
% INPUT:
%               - as are time-series analysis objects or a vector of
%               analysis objects
%               - pl is a plist with the input parameters
%
% OUTPUT:
%               - bs "whitened" time-series AOs. The whitening filters used
%               are stored in the objects procinfo field under the
%               parameter 'Filter'.
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'whiten1D')">Parameters Description</a>
%
% VERSION:     $Id: whiten1D.m,v 1.43 2011/04/08 08:56:12 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = whiten1D(varargin)
  
  % 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 = [as.hist];
  
  % combine plists
  if isempty(pl)
    model = 'psd';
  else
    model = find(pl, 'model');
    if isempty(model)
      model = 'psd';
    end
  end
  
  if ischar(model)
    pl = parse(pl, getDefaultPlist(model));
  else
    pl = parse(pl, getDefaultPlist('Default'));
  end
  pl.getSetRandState();
  
  scale = find(pl, 'scaleOut');
  flim = find(pl, 'flim');
  
  
  % Loop over input AOs
  for jj = 1:numel(as)
    if ~isa(as(jj).data, 'tsdata')
      utils.helper.msg(msg.IMPORTANT, '%s expects ao/tsdata objects. Skipping AO %s', mfilename, ao_invars{jj});
    else
      
      %-------------- Whiten this AO
      
      % 1) Build whitening filterbank
      switch class(model)
        case 'char'
          % Model is to be evaluated from data
          in = as(jj);
          pl.pset('model', model);
        case 'ao'
          % Model was provided as fsdata
          in = model;
          pl.pset('model', []);          
      end
      wf = buildWhitener1D(in, pl);
            
      % 1.5) Scale the date if demanded
      if (scale)
        spsd = lpsd(as(jj));
        freqs = spsd.x;
        if isempty(flim)
          error('Please specify a flim field, to know the analysis band.');
        elseif (flim(2) < flim(1))
          error('flim should go from the smaller frequency to the bigger frequency. Please reverse them!')
        else
          index = find((freqs > flim(1)) & (freqs < flim(2)));
        end
        
        v1 = spsd.y(index(1):index(end-1));
        v2 = spsd.y(index(2):index(end));
        m = (v1 + v2) /2;
        p = sum(m.* diff(freqs(index(1):index(end)))); 
      end
      
      % 2) Filter data and scale it if necessary
      bs(jj).filter(wf);
      if (scale)
        bs(jj) = bs(jj) * sqrt(p);
      end

      
      % 3) Output data
      % name for this object
      bs(jj).name = sprintf('whiten1D(%s)', ao_invars{jj});
      % Collect the filters into procinfo
      bs(jj).procinfo = combine(plist('Filter', wf.filters), as(jj).procinfo);
      if(scale)
        bs(jj).procinfo = combine(plist('ScaleFactor', p, 'Filter', wf.filters), as(jj).procinfo);
      end
      % Make sure that the output yunits are empty
      if ~eq(bs(jj).yunits, unit(''))
        utils.helper.msg(msg.PROC1, 'Resetting output yunits to empty');
        bs(jj).setYunits(unit(''));
      end
      % add history
      bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), inhists(jj));
      % clear errors
      bs(jj).clearErrors;
      
      
    end
  end
  
  
  
  % Set output
  if nargout == numel(bs)
    % List of outputs
    for ii = 1:numel(bs)
      varargout{ii} = bs(ii);
    end
  else
    % Single output
    varargout{1} = bs;
  end
end

%--------------------------------------------------------------------------
% Get Info Object
%--------------------------------------------------------------------------
function ii = getInfo(varargin)
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pl   = [];
  elseif nargin == 1 && ~isempty(varargin{1}) && ischar(varargin{1})
    sets{1} = varargin{1};
    pl = getDefaultPlist(sets{1});
  else
    sets = SETS();
    % get plists
    pl(size(sets)) = plist;
    for kk = 1:numel(sets)
      pl(kk) =  getDefaultPlist(sets{kk});
    end
  end
  % Build info object
  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: whiten1D.m,v 1.43 2011/04/08 08:56:12 hewitson Exp $', sets, pl);
end


%--------------------------------------------------------------------------
% Defintion of Sets
%--------------------------------------------------------------------------

function out = SETS()
  out = ao.getInfo('buildWhitener1D').sets;
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------
function plout = getDefaultPlist(set)
  persistent pl;  
  persistent lastset;
  if ~exist('pl', 'var') || isempty(pl) || ~strcmp(lastset, set)
    pl = buildplist(set);
    lastset = set;
  end
  plout = pl;  
end

function pl = buildplist(set)
  
  pl = plist();
  
  % Append sets of parameters according to the chosen spectral estimator
  if ~utils.helper.ismember(lower(SETS), lower(set))
    error('### Unknown set [%s]', set);
  else
    pl = ao.getInfo('buildWhitener1D', lower(set)).plists;
  end
  
   switch lower(set)
    case 'default'
      % Model
      p = param({'model', ['A frequency-series AO describing the model<br>'...
        'response to build the filter from. <br>' ...
        'As an alternative, the user '...
        'can choose a model estimation technique:<br>'...
        '<li>PSD - using <tt>psd</tt> + <tt>bin_data</tt></li>'...
        '<li>LPSD - using <tt>lpsd</tt></li>']}, paramValue.EMPTY_DOUBLE);
      pl = combine(plist(p), pl);
    otherwise
   end
  
   p = param({'scaleOut', ['Scale your output by the inband power']},paramValue.FALSE_TRUE);
   pl = combine(plist(p), pl);

   p = param({'flim', ['Band to calculate the scaling power']},[1e-3 30e-3]);
   pl = combine(plist(p), pl);

   
end