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