--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/sineParams.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,241 @@
+% 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
+