% DIFF differentiates the data in AO.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: DIFF differentiates the data in AO. The result is a data
%              series the same length as the input series.
%              In case of method 'diff' computes the difference between two samples, in which
%              case the resulting time object has the length of the input
%              series -1 sample.
% CALL:        bs = diff(a1,a2,a3,...,pl)
%              bs = diff(as,pl)
%              bs = as.diff(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,
%                     containing the differentiated data
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'diff')">Parameters Description</a>
%
% REFERENCES:
% [1] L. Ferraioli, M. Hueller and S. Vitale, Discrete derivative
%     estimation in LISA Pathfinder data reduction,
%     <a
%     href="matlab:web('http://www.iop.org/EJ/abstract/0264-9381/26/9/094013/','-browser')">Class. Quantum Grav. 26 (2009) 094013.</a>
% [2] L. Ferraioli, M. Hueller and S. Vitale, Discrete derivative
%     estimation in LISA Pathfinder data reduction
%     <a href="matlab:web('http://arxiv.org/abs/0903.0324v1','-browser')">http://arxiv.org/abs/0903.0324v1</a>
%
% VERSION:     $Id: diff.m,v 1.36 2011/08/03 19:18:56 adrien Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PARAMETERS:  method - the method to use: [default: '3POINT']
%                       'diff'   - like MATLABs diff
%                                  compute difference between each two samples
%                                  xaxis will be integers from 1 to length
%                                  of resulting object
%                       '2POINT' - 2 point derivative computed as
%                                  [y(i+1)-y(i)]./[x(i+1)-x(i)].
%                       '3POINT' - 3 point derivative. Compute derivative
%                                  at i as [y(i+1)-y(i-1)] / [x(i+1)-x(i-1)].
%                                  For i==1, the output is computed as
%                                  [y(2)-y(1)]/[x(2)-x(1)]. The last sample
%                                  is computed as [y(N)-y(N-1)]/[x(N)-x(N-1)].
%                       '5POINT' - 5 point derivative. Compute derivative dx
%                                  at i as
%                                  [-y(i+2)+8*y(i+1)-8*y(i-1)+y(i-2)] /
%                                       [3*(x(i+2)-x(i-2))].
%                                  For i==1, the output is computed as
%                                  [y(2)-y(1)]/[x(2)-x(1)]. The last sample
%                                  is computed as [y(N)-y(N-1)]/[x(N)-x(N-1)].
%                       'ORDER2' - Compute derivative using a 2nd order
%                                  method.
%                 'ORDER2SMOOTH' - Compute derivative using a 2nd order
%                                  method with a parabolic fit to 5
%                                  consecutive samples.
%                       'filter' - applies an IIR filter built from a
%                                  single pole at the chosen frequency. The
%                                  filter is applied forwards and backwards
%                                  (filtfilt) to achieve the desired f^2
%                                  response. This only works for time-series
%                                  AOs. For this method, you can specify the
%                                  pole frequency with an additional parameter
%                                  'f0' [default: 1/Nsecs]
%                          'FPS' - Calculates five points derivative using
%                                  utils.math.fpsder function. If you call
%                                  with this oprtion you may add also the
%                                  parameters:
%                                   'ORDER' derivative order, supperted
%                                   values are:
%                                     'ZERO', 'FIRST', 'SECOND'
%                                   'COEFF' coefficient used for the
%                                   derivation. Refers to the fpsder help
%                                   for further details.
%
%

function varargout = diff(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);
  
  % combine plists
  pl = parse(pl, getDefaultPlist());
  
  % Extract method
  method = find(pl, 'method');
  
  for jj = 1:numel(bs)
    
    % Diff can't work for cdata objects since we need x data
    if isa(bs(jj).data, 'cdata')
      error('### diff doesn''t work with cdata AOs since we need an x-data vector.');
    end
    
    % Compute derivative with selected method
    switch lower(method)
      case 'diff'
        yunit = bs(jj).yunits;
        y      = bs(jj).y;
        x      = bs(jj).x;
        newX   = x(1:end-1); % cut the last sample from the time series to make x and y same length
        dy     = diff(y);
        bs(jj).data.setY(dy);
        bs(jj).data.setX(newX);
        bs(jj).setYunits(yunit);
      case '2point'
        x      = bs(jj).data.getX;
        dx     = diff(x);
        y      = bs(jj).data.getY;
        dy     = diff(y);
        z      = dy./dx;
        bs(jj).data.setY(z);
        bs(jj).data.setX((x(1:end-1)+x(2:end))/2);
        bs(jj).data.setYunits(bs(jj).data.yunits./bs(jj).data.xunits);
      case '3point'
        x  = bs(jj).data.getX;
        dx = diff(x);
        y  = bs(jj).data.getY;
        z  = zeros(size(y));
        z(2:end-1) = (y(3:end)-y(1:end-2)) ./ (dx(2:end)+dx(1:end-1));
        z(1)       = (y(2)-y(1)) ./ (dx(1));
        z(end)     = 2*z(end-1)-z(end-2);
        bs(jj).data.setY(z);
        bs(jj).data.setYunits(bs(jj).data.yunits./bs(jj).data.xunits);
      case '5point'
        x  = bs(jj).data.getX;
        dx = diff(x);
        y  = bs(jj).data.getY;
        z  = zeros(size(y));
        z(1)       = (y(2)-y(1)) ./ (dx(1));
        z(2)       = (y(3)-y(1))./(dx(2)+dx(1));
        z(3:end-2) = (-y(5:end) + 8.*y(4:end-1) - 8.*y(2:end-3) + y(1:end-4)) ./ (3.*(x(5:end)-x(1:end-4)));
        z(end-1)   = 2*z(end-2)-z(end-3);
        z(end)     = 2*z(end-1)-z(end-2);
        bs(jj).data.setY(z);
        bs(jj).data.setYunits(bs(jj).data.yunits./bs(jj).data.xunits);
      case 'order2'
        x     = bs(jj).data.getX;
        dx    = diff(x);
        y     = bs(jj).data.getY;
        z     = zeros(size(y));
        m     = length(y);
        % y'(x1)
        z(1) = (1/dx(1)+1/dx(2))*(y(2)-y(1))+...
          dx(1)/(dx(1)*dx(2)+dx(2)^2)*(y(1)-y(3));
        % y'(xm)
        z(m) = (1/dx(m-2)+1/dx(m-1))*(y(m)-y(m-1))+...
          dx(m-1)/(dx(m-1)*dx(m-2)+dx(m-2)^2)*(y(m-2)-y(m));
        % y'(xi) (i>1 & i<m)
        dx1 = repmat(dx(1:m-2),1,1);
        dx2 = repmat(dx(2:m-1),1,1);
        y1 = y(1:m-2); y2 = y(2:m-1); y3 = y(3:m);
        z(2:m-1) = 1./(dx1.*dx2.*(dx1+dx2)).*...
          (-dx2.^2.*y1+(dx2.^2-dx1.^2).*y2+dx1.^2.*y3);
        bs(jj).data.setY(z);
        bs(jj).data.setYunits(bs(jj).data.yunits./bs(jj).data.xunits);
      case 'order2smooth'
        x  = bs(jj).data.getX;
        y  = bs(jj).data.getY;
        dx = diff(x);
        m  = length(y);
        if max(abs(diff(dx)))>sqrt(eps(max(abs(dx))))
          error('### The x-step must be constant for method ''ORDER2SMOOTH''')
        elseif m<5
          error('### Length of y must be at least 5 for method ''ORDER2SMOOTH''.')
        end
        h = mean(dx);
        z = zeros(size(y));
        % y'(x1)
        z(1) = sum(y(1:5).*[-54; 13; 40; 27; -26])/70/h;
        % y'(x2)
        z(2) = sum(y(1:5).*[-34; 3; 20; 17; -6])/70/h;
        % y'(x{m-1})
        z(m-1) = sum(y(end-4:end).*[6; -17; -20; -3; 34])/70/h;
        % y'(xm)
        z(m) = sum(y(end-4:end).*[26; -27; -40; -13; 54])/70/h;
        % y'(xi) (i>2 & i<(N-1))
        Dc = [2 1 0 -1 -2];
        tmp = convn(Dc,y)/10/h;
        z(3:m-2) = tmp(5:m);
        bs(jj).data.setY(z);
        bs(jj).data.setYunits(bs(jj).data.yunits./bs(jj).data.xunits);
      case 'filter'
        error('### Comming with release 2.5');
      case 'fps'
        order = find(pl, 'ORDER');
        coeff = find(pl, 'COEFF');
        x  = bs(jj).data.getX;
        dx = x(2)-x(1);
        fs = 1/dx;
        y  = bs(jj).data.getY;
        params = struct('ORDER', order, 'COEFF', coeff, 'FS', fs);
        z = utils.math.fpsder(y, params);
        bs(jj).data.setY(z);
        % setting units
        switch lower(order)
          case 'first'
            bs(jj).data.setYunits(bs(jj).data.yunits./bs(jj).data.xunits);
          case 'second'
            bs(jj).data.setYunits(bs(jj).data.yunits.*bs(jj).data.xunits.^(-2));
        end
      otherwise
        error('### Unknown method for computing the derivative.');
    end
    
    % name for this object
    bs(jj).name = sprintf('diff(%s)', ao_invars{jj});
    % add history
    bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), bs(jj).hist);
  end
  
  % Clear the errors since they don't make sense anymore
  clearErrors(bs);
    
  % 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: diff.m,v 1.36 2011/08/03 19:18:56 adrien 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({'method',['The method to use. Choose between:<ul>', ...
    '' ...
    '<li>''2POINT'' - 2 point derivative computed as [y(i+1)-y(i)]./[x(i+1)-x(i)]', ...
    '</li>' ...
    '<li>''3POINT'' - 3 point derivative. Compute derivative dx at i as <br>', ...
    '<tt>[y(i+1)-y(i-1)] / [x(i+1)-x(i-1)]</tt><br>', ...
    'For <tt>i==1</tt>, the output is computed as <tt>[y(2)-y(1)]/[x(2)-x(1)]</tt>.<br>', ...
    'The last sample is computed as <tt>[y(N)-y(N-1)]/[x(N)-x(N-1)]</tt>', ...
    '</li>' ...
    '<li>''5POINT'' - 5 point derivative. Compute derivative dx at i as <br>', ...
    '<tt>[-y(i+2)+8*y(i+1)-8*y(i-1)+y(i-2)] / [3*(x(i+2)-x(i-2))]</tt><br>', ...
    'For <tt>i==1</tt>, the output is computed as <tt>[y(2)-y(1)]/[x(2)-x(1)]</tt><br>', ...
    'The last sample is computed as <tt>[y(N)-y(N-1)]/[x(N)-x(N-1)]</tt>', ...
    '</li>' ...
    '<li>''ORDER2'' - Compute derivative using a 2nd order method', ...
    '</li>' ...
    '<li>''ORDER2SMOOTH'' - Compute derivative using a 2nd order method<br>', ...
    'with a parabolic fit to 5 consecutive samples', ...
    '</li>' ...
    '<li>''filter'' - applies an IIR filter built from a single pole at the chosen frequency.<br>', ...
    'The filter is applied forwards and backwards (filtfilt) to achieve the desired f^2<br>', ...
    'response. This only works for time-series AOs.<br>', ...
    'For this method, you can specify the pole frequency with an additional parameter ''F0'' (see below):', ...
    '</li>'...
    '<li>''FPS'' - Calculates five points derivative using utils.math.fpsder.<br>', ...
    'When calling with this option you may add also the parameters ''ORDER'' (see below)<br>', ...
    'and ''COEFF'' (see below)' ...
    '</li>' ...
    ]},  {1, {'2POINT', '3POINT', '5POINT', 'ORDER2', 'ORDER2SMOOTH', 'FILTER', 'FPS'}, paramValue.SINGLE});
  pl.append(p);
  
  % F0
  p = param({'f0','The pole frequency for the ''filter'' method.'}, {1, {'1/Nsecs'}, paramValue.OPTIONAL});
  pl.append(p);
  
  % Order
  p = param({'ORDER','Derivative order'}, {1, {'ZERO', 'FIRST', 'SECOND'}, paramValue.SINGLE});
  pl.append(p);
  
  % Coeff
  p = param({'COEFF',['Coefficient used for the derivation. <br>', ...
    'Refer to the <a href="matlab:doc(''utils.math.fpsder'')">fpsder help</a> for further details']}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
end