% SPSD implements the smoothed (binned) PSD algorithm for analysis objects.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: SPSD implements the smoothed PSD algorithm for analysis objects.
%
% CALL:        bs = spsd(a1,a2,a3,...,pl)
%              bs = spsd(as,pl)
%              bs = as.spsd(pl)
%
% INPUTS:      aN   - input analysis objects
%              as   - input analysis objects array
%              pl   - input parameter list
%
% OUTPUTS:     bs   - array of analysis objects, one for each input
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'spsd')">Parameters Description</a>
%
% VERSION:     $Id: spsd.m,v 1.21 2011/07/11 10:43:35 adrien Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = spsd(varargin)

  import utils.const.*

  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end

  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, rest] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  [pl, pl_invars, rest] = utils.helper.collect_objects(rest(:), 'plist', in_names);

  % Decide on a deep copy or a modify
  bs = copy(as, nargout);

  % Combine plists
  pl = combine(pl, plist(rest(:)),  getDefaultPlist); 

  inhists = [];

  %% Go through each input AO
  for jj = 1 : numel(bs)
    % gather the input history objects
    inhists = [inhists bs(jj).hist]; %#ok<AGROW>

    % check this is a time-series object
    if ~isa(bs(jj).data, 'tsdata')
      warning('!!! spsd requires tsdata (time-series) inputs. Skipping AO %s', ao_invars{jj}); %#ok<WNTAG>
    else

      % Check the time range.
      time_range = find(pl, 'times');
      if ~isempty(time_range)
        bs(jj) = split(bs(jj), plist('method', 'times', 'times', time_range));
      end
      % Check the length of the object
      if bs(jj).len <= 0
        error('### The object is empty! Please revise your settings ...');
      end

      % pl = utils.helper.process_spectral_options(pl, 'log');
      pl = pl.combine(getDefaultPlist());

      % getting data
      y = bs(jj).y;

      % Window function
      Win = find(pl, 'Win');
      nfft = length(y);
      Win = ao( combine(plist('win', Win , 'length', nfft), pl) );

      % detrend
      order = find(pl,'order');
      if ~(order < 0)
        y = ltpda_polyreg(y,  order).';
      else 
        y = reshape(y, 1, nfft);
      end

      % computing PSD
      window = Win.data.y;
      window = window/norm(window)*sqrt(nfft);
      yASD = real(fft(y.*window, nfft)).^2 + imag(fft(y.*window, nfft)).^2;
      pow = [yASD(1) yASD(2:floor(nfft/2))*2];
      pow = pow /  ( bs(jj).data.fs * nfft);
      Freqs = linspace(0, bs(jj).data.fs/2, nfft/2);

      % smoothing PSD
      if ~isempty(find(pl,'frequencies'))
        error('the option "frequencies" is deprecated, frequencies are "removed" by default')
      end
      [Freqs, pow, nFreqs, nDofs] = ltpda_spsd(Freqs, pow, find(pl,'linCoef'), find(pl,'logCoef') );
      % create new output fsdata
      scale = find(pl, 'Scale');
      switch lower(scale)
        case 'asd'
          fsd = fsdata(Freqs, sqrt(pow), bs(jj).data.fs);
          fsd.setYunits(bs(jj).data.yunits / unit('Hz^0.5'));
          %           stdDev = 0.5 * sqrt( pow ./ nDofs ); % linear approximation of the sqrt of a distribution
          % approximation knowing the STD of the PSD
          % STD assuming amplitude samples are independent, Chi^1_2 distibuted
          % (with both variables of powe expectancy pow/2), and of different
          % magnitude
          stdDev = 2 * sqrt(pow./nDofs) .* ( nDofs - 2*exp( 2*(gammaln((nDofs+1)/2)-gammaln(nDofs/2)) ) ); % std of the chi_2N^1 
        case 'psd'
          fsd = fsdata(Freqs, pow, bs(jj).data.fs);
          fsd.setYunits(bs(jj).data.yunits.^2/unit('Hz'));
          % STD assuming power samples are independent, Chi^2_2 distibuted
          % (with both variables of expectancy pow/2), and of different
          % magnitude
          stdDev = sqrt(2) * (pow./nDofs) .* sqrt(2*nDofs);  % std of the chi_2N^2 
        otherwise
          error(['### Unknown scaling:' scale]);
      end

      fsd.setXunits('Hz');
      fsd.setDx(nFreqs*Freqs(2)/2);
      fsd.setEnbw(1);% WARNING HERE!!!
      fsd.setT0(bs(jj).data.t0);
      % make output analysis object
      bs(jj).data = fsd;
      % set name
      bs(jj).name = ['SPSD(', ao_invars{jj}, ') ' upper(scale)];
      % Add standard deviation
      bs(jj).data.dy = stdDev;
      % Add history
      bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), inhists(jj));

    end % End tsdata if/else
  end % End AO loop

  %% 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   = [];
  else
    sets = {'Default'};
    pl   = getDefaultPlist;
  end
  % Build info object
  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: spsd.m,v 1.21 2011/07/11 10:43:35 adrien Exp $', sets, pl);
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------
function pl = getDefaultPlist()

  % Plist for Welch-based, log-scale spaced spectral estimators.
  pl = plist;

  % Win
  p = param({'Win',['the window to be applied to the data to remove the ', ...
    'discontinuities at edges of segments. [default: taken from user prefs] <br>', ...
    'Only the design parameters of the window object are used. Enter either: <ul>', ...
    '<li> a specwin window object OR</li>', ...
    '<li> a string value containing the window name</li></ul>', ...
    'e.g., <tt>plist(''Win'', ''Kaiser'', ''psll'', 200)</tt>']}, paramValue.WINDOW);
  pl.append(p);

  % Psll
  p = param({'Psll',['the peak sidelobe level for Kaiser windows.<br>', ...
    'Note: it is ignored for all other windows']}, paramValue.DOUBLE_VALUE(200));
  pl.append(p);

  % Psll
  p = param({'levelOrder','the contracting order for levelledHanning window'}, paramValue.DOUBLE_VALUE(2));
  pl.append(p);

  % Order
  p = param({'Order',['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);
  p.val.setValIndex(2);
  pl.append(p);

  % Times
  p = param({'Times','time range. If not empty, sets the restricted interval to analyze'}, paramValue.DOUBLE_VALUE([]));
  pl.append(p);

  % Scale
  p = param({'Scale',['scaling of output. Choose from:<ul>', ...
    '<li>PSD - Power Spectral Density</li>', ...
    '<li>ASD - Amplitude (linear) Spectral Density</li>'...
    ]}, {1, {'PSD', 'ASD', 'PS', 'AS'}, paramValue.SINGLE});
  pl.append(p);

  p = param( {'lincoef', 'Linear scale smoothing coefficent (freq. bins)'}, 1);
  pl.append(p);

  p = param( {'logcoef', ['Logarithmic scale smoothing coefficent<br>', 'Best compromise for both axes is 2/3']}, 2/3);
  pl.append(p);
end