% SINEPARAMS estimates parameters of sinusoids
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: SINEPARAMS estimates the parameters of the sinusoids in the
%              time series. Number of sinusoids needs to be specified.
%
% CALL:        b = sineParams(a,pl)
%
% INPUTS:      a  - input AOs
%              pl - parameter list (see below)
%
% OUTPUTs:     b  - pest object
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'sineParams')">Parameters Description</a>
%
% VERSION:     $Id: sineParams.m,v 1.9 2011/04/08 08:56:13 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = sineParams(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
  [aos, 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(aos, nargout);
  
  % combine plists
  pl = parse(pl, getDefaultPlist());
  
  % check (current version with only one AO)
  if numel(aos) >1
    error('### Current version of ''sineParams'' works only with one AO')
  end
  
  % get parameters
  Nsine = find(pl, 'N');
  realSine  = find(pl, 'real');
  method = find(pl,'method');
  nonlin = find(pl,'non-linear');
  
  % check number of sinusoids
  if isempty(Nsine)
    error('### You need to set the number of sinusoids ''N'' ')
  end
  
  if realSine
    Nexp = 2*Nsine;
  else
    Nexp = Nsine;
  end
  f =[];A=[];
  %   for i = 1:numel(bs)
  
  switch method
    case 'MUSIC'
      % MUSIC algorithm
      [w,pow,w_mse,pow_mse] = utils.math.rootmusic(bs.y,Nexp,bs.fs);
      
      if realSine
        f = w(1:2:end);
        df = w_mse(1:2:end);
        %         A = pow;
        %         dA = pow_mse;
        A = sqrt(2*pow); % from the expression for power in a sinusoid
        dA = 1/sqrt(2*pow)*pow_mse;
      else
        f = w;
        df = w_mse;
        A = sqrt(2*pow); % from the expression for power in a sinusoid
        dA = A*sqrt(2/pow)*pow_mse;
      end
      
    case 'ESPRIT'
      
      % p = 1; % num. sinusoid
      % N = 50;
      %
      % m = xcov(c.y,N-1);
      %
      % clear cmat
      % for i =1:N
      % cmat(i,:) = m(N-(i-1):end-(i-1));
      % end
      %
      % [U,S,UT] = svd(cmat);
      %
      % D = diag(S);
      %
      % Ls = D(1:2*p);
      % Us = U(:,1:2*p);
      %
      % A1 = Us(1:end-1,:);
      % A2 = Us(2:end,:);
      %
      % M = A1\A2;
      %
      % freq = imag(-log(eig(M))/2/pi*fs)
      % error = real(-log(eig(M))/2/pi*fs)
      
    case 'IFFT'
      
      %   % fft
      %     p = psd(phi(i),plist('win','Hanning','scale','AS','Navs',1));
      %
      %     [m,index] = max(p.y);
      %     ratio between max. and next value
      %     alpha = p.y(index+1)/p.y(index);
      %
      %     xm = (2*alpha-1)/(alpha+1);
      %     frequency
      %     f(i) = (i+xm)/phi(i).nsecs;
      %     amplitude
      %     A(i) = abs(2*pi*xm*(1-xm)/sin(pi*xm)*exp(-pi*1i*xm)*(1*xm)*p.y(i));
      %
      %
      %
      %
  end
  %   end
  
  % create output pest
  mdl = [];
  p = pest();
  for i = 1:Nsine
    Cname = sprintf('C%d',i);
    fname = sprintf('f%d',i);
    Aname = sprintf('A%d',i);
    pname = sprintf('phi%d',i);
    % setY
    p.setYforParameter(Cname,0);
    p.setYforParameter(Aname,A(i));
    p.setYforParameter(fname,f(i));
    p.setYforParameter(pname,0);
    % errors for single sinusoid only, error is the sqrt(mse)
    if Nsine == 1
      p.setDyForParameter(Cname,inf);
      p.setDyForParameter(Aname,sqrt(dA(i)));
      p.setDyForParameter(fname,sqrt(df(i)));
      p.setDyForParameter(pname,inf);
    end
    % smodel for each sinusoid
    m = smodel(sprintf('%s + %s*sin(2*pi*%s*t+%s)',Cname,Aname,fname,pname));
    m.setXunits('s');
    m.setXvar('t');
    m.setXvals(bs.x);
    m.setParams({Cname,Aname,fname,pname},[0 A(i) f(i) 0]);
    % optional non-linear
    if nonlin
      pl = plist('Function',m);
      pnl = xfit(bs,pl);
      p.setY(pnl.y);
      p.setDy(pnl.dy);
      p.setCorr(pnl.corr);
      p.setCov(pnl.cov);
      p.setDof(pnl.dof);
      p.setChi2(pnl.chi2);
      % set values in model as well
      m.setValues(pnl.y);
    end
    mdl = [mdl m];
  end
  % set name
  p.name = sprintf('sineParams(%s)', ao_invars{:});
  % set models
  p.setModels(mdl);
  % Add history
  p.addHistory(getInfo('None'), pl, ao_invars(:), bs(:).hist);
  
  % Set outputs
  if nargout > 0
    varargout{1} = p;
  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: sineParams.m,v 1.9 2011/04/08 08:56:13 hewitson 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();
  
  % number of sinusoids
  p = param({'N', 'Number of sinusoids/complex exp. in the time series'},paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % number of sinusoids
  p = param({'real', 'Set to true if working with real sinusoids'},paramValue.TRUE_FALSE);
  pl.append(p);
  
  % number of sinusoids
  p = param({'method', ['Choose one of the following methods:<ul>', ...
    '<li>MUSIC  - MUltiple SIgnal Classification algorithm </li>']}, {1, {'MUSIC'}, paramValue.SINGLE});
  pl.append(p);
  
  % non-linear
  p = param({'non-linear','Set to true to perform non-linear fit'},...
    paramValue.FALSE_TRUE);
  pl.append(p);
  
end
% END