--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/delayEstimate.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,166 @@
+% DELAYESTIMATE estimates the delay between two AOs
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: DELAYESTIMATE returns an estimate of the delay between the two
+%              input analysis objects. Different weights in frequency
+%              domain can be used.
+%
+% CALL:        bs = delayEstimate(a1,a2,pl)
+%
+% INPUTS:      a1    - input analysis objects
+%              a2    - delayed analysis objects
+%              pl    - input parameter list
+%
+% OUTPUTS:     bs    - analysis object with the delay (as cdata)
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'delayEstimate')">Parameters Description</a>
+%
+% VERSION:    $Id: delayEstimate.m,v 1.6 2011/04/08 08:56:15 hewitson Exp $
+%
+% HISTORY:    15-10-09 M Nofrarias
+%             Creation
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+function varargout = delayEstimate(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
+  [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);
+  
+  % Combine plists
+  pl = parse(pl, getDefaultPlist);
+  
+  % get frequency weight
+  weight = find(pl, 'weight');
+  
+  N = len(bs(1));
+  T = bs(1).x(end);
+  fs = bs(1).fs;
+  % zeropad to avoid circular convolution in ifft
+  bs(1).zeropad;
+  bs(2).zeropad;
+  % compute spectral density, scale applied after to avoid round-off errors
+  scale = fs/N;
+  Gx11 = conj(fft(bs(1).y)).*fft(bs(1).y);
+  Gx11 = Gx11*scale;
+  Gx12 = conj(fft(bs(1).y)).*fft(bs(2).y);
+  Gx12 = Gx12*scale;
+  Gx22 = conj(fft(bs(2).y)).*fft(bs(2).y);
+  Gx22 = Gx22*scale;
+  % frequencies for two-sided spectrum
+  f = linspace(-fs,fs,2*N);
+  % select weight
+  if strcmp(weight,'roth')
+    weight = Gx11;
+  elseif strcmp(weight,'scot')
+    weight = sqrt(Gx11.* Gx22);
+  elseif strcmp(weight,'scc')
+    weight = 1;
+  elseif strcmp(weight,'phat')
+    weight = abs(Gx12);
+  elseif strcmp(weight,'eckart')
+    weight = 1;
+  elseif strcmp(weight,'ML')
+    weight = 1;
+  else
+    error('###  Unknown value for ''weight'' parameter '  )
+  end
+  % compute unscaled correlation function
+  Ru = ifft(Gx12./weight);
+  % lag= 0:\deltat*(N-1)
+  % n= 1:N/2-1
+  r = linspace(0,T-1/fs,length(Ru)/2-1);
+  % scaling to correct zeropad bias
+  R = (N./(N-r))'.*Ru(1:length(Ru)/2-1);
+  % get maximum
+  [m,ind] = max(R);
+  del = r(ind);
+  % plot if required
+  plt = find(pl, 'plot');
+  if plt
+    Rxy = ao(xydata(r,R));
+    Rxy.setName(sprintf('Correlation(%s,%s)', ao_invars{1},ao_invars{2}))
+    iplot(Rxy)
+  end
+  
+  % create new output cdata
+  cd = cdata(del);
+  % add unitss
+  cd.yunits = 's';
+  % update AO
+  c = ao(cd);
+  % add error
+  %   bs(jj).data.dy = dev;
+  % Add name
+  c.name = sprintf('delayEst(%s,%s)', ao_invars{1},ao_invars{2});
+  % Add history
+  c.addHistory(getInfo('None'), pl, [ao_invars(:)], [bs(:).hist]);
+  
+  % set output
+  varargout{1} = c;
+  
+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: delayEstimate.m,v 1.6 2011/04/08 08:56:15 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();
+  % method
+  p = param({'weight',['scaling of output. Choose from:<ul>', ...
+    '<li>scc  - </li>', ...
+    '<li>roth  - </li>', ...
+    '<li>scot  - </li>', ...
+    '<li>phat  - </li>', ...
+    '<li>eckart  - </li>', ...
+    '<li>ML - </li></ul>']}, {1, {'scc','roth', 'scot','phat','eckart','ML'}, paramValue.SINGLE});
+  pl.append(p);
+  % Plot
+  p = param({'Plot', 'Plot the correlation function'}, ...
+    {2, {'true', 'false'}, paramValue.SINGLE});
+  pl.append(p);
+end
+% END
+