--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/fpsder.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,173 @@
+% FPSDER performs the numeric time derivative
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: FPSDER (Five Points Stencil Derivative) performs the numeric
+% time derivative using the method of five points stencil.
+% The function can perform first, second, third and fourth derivetive
+% of a series of input data.
+% Given a discrete series of data points, the five-point-stencil method for
+% the derivative approximation, at a given time t0, is calculated by means
+% of finite differences between the element at t0 with its four neighbors.
+% The n-order derivative at a certain time can be approximated by a five
+% point difference equation:
+%
+% d^{n}y[k]
+% --------- = (1/T^n) * {a*y[k-2] + b*y[k-1] + c*y[k] + d*y[k+1] + e*y[k+2]}
+%  dt^{n}
+%
+% It can be demonstrated [1,2] that the five coefficients [a, b, c, d, e] can
+% be written in terms of only one of them. In fpsder the independent
+% coefficient is fixed to be the first and is called m. It can be input as
+% a parameter when the function is called.
+%
+%
+% CALL:                 Deriv = fpsder(data, params)
+%
+% INPUTS:
+%
+%         - a is a vector containing the data to be differentiated.
+%         - params is a struct with the input parameters:
+%
+%         - 'ORDER' set the derivative order. Its allowed options are:
+%             - 'ZERO' perform data smoothing using the couefficients
+%             vector d0 = [m -4*m 1+6*m -4*m m].
+%             - 'FIRST' perform the first derivative using the
+%             couefficients vector d1 = [m -(0.5+2*m) 0 (0.5+2*m) m]./T.
+%             Recomended values of m are in the interval [-0.1, 0.1].
+%             - 'SECOND' perform the second derivative using the
+%             coefficients vector d2 = [m 1-4*m 6*m-2 1-4*m m]./(T^2).
+%             Recomended values of m are in the interval [-0.11, 0.3].
+%             - 'THIRD' perform the third derivative using the
+%             coefficients vector d3 = []./(T^3)
+%             - 'FOURTH' perform the third derivative using the
+%             coefficients vector d4 = []./(T^4)
+%
+%         - 'COEFF' set m coefficient values.
+%           In case of data smoothing: m = -3/35 correspond to the
+%           parabolic fit approximation.
+%           In case of first order derivative: m = -1/5 correspond to the
+%           parabolic fit approximation, m = 1/12 correspond to the
+%           Taylor series approximation.
+%           In case of second order derivative: m = 2/7 corresponds to
+%           the parabolic fit approximation, m = -1/12 corresponds to the
+%           Taylor series approximation and m = 1/4 gives the notch
+%           feature at the Nyquist frequency
+%
+%         - 'FS' set the data sampling frequency in Hz
+%
+%         NOTE1: T is the sampling period
+%         NOTE2: The default option for 'ORDER' is 'SECOND'
+%         NOTE3: The default option for 'COEFF' is 2/7
+%         NOTE4: The default option for 'FS' is 10
+%
+% OUTPUTS:
+%           - D is a vector containing the resulting data after the
+%           differentiation procedure
+%
+% 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>
+%
+% EXAMPLES:
+%           - Performing the second order derivative of a series of data, m
+%           coefficient is fixed to 2/7 and data sampling frequency is
+%           fixed to 10 Hz.
+%           params = struct('ORDER', 'SECOND', 'COEFF', 2/7, 'FS', 10);
+%           Deriv = fpsder(data, params);
+%
+% VERSION: '$Id: fpsder.m,v 1.8 2010/04/09 09:59:59 mauro Exp $';
+%
+% HISTORY:     18-03-2008 L Ferraioli
+%                 Creation
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+function   Deriv = fpsder(a, params)
+
+  % Getting input parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Collect inputs
+
+  % Default input struct
+  defaultparams = struct('ORDER','SECOND',...
+    'COEFF',2/7,...
+    'FS',10);
+
+  names = {'ORDER','COEFF','FS'};
+
+  % collecting input and default params
+  if ~isempty(params)
+    for jj=1:length(names)
+      if isfield(params, names(jj)) && ~isempty(params.(names{1,jj}))
+       defaultparams.(names{1,jj}) = params.(names{1,jj});
+      end
+    end
+  end
+
+  % values for input variables
+  order = defaultparams.ORDER;
+  m = defaultparams.COEFF;
+  fs = defaultparams.FS;
+  
+  % willing to work with columns
+  if size(a,2)>1
+    a = a.';
+  end
+
+  % Assigning coefficients values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  % Assigning coefficients values based on the input options
+  switch upper(order)
+    case 'ZERO'
+      Coeffs = [m -4*m 1+6*m -4*m m];
+    case 'FIRST'
+      Coeffs = [m -(0.5+2*m) 0 (0.5+2*m) -m];
+    case 'SECOND'
+      Coeffs = [m 1-4*m 6*m-2 1-4*m m];
+    case 'THIRD'
+      Coeffs = [0 0 0 0 0];
+      disp('Not yet implemented, sorry!');
+    case 'FOURTH'
+      Coeffs = [0 0 0 0 0];
+      disp('Not yet implemented, sorry!');
+    otherwise
+      error('### Unknown order %s', order);
+  end
+
+  % Sampling period
+  T = 1/fs;
+
+  % Building vectors for calculation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  % Building the 'extended' vector for calculation
+  % a_temp = [a(1);a(1);a(1);a(1);a;a(end);a(end);a(end);a(end)];
+  a_temp = [2*a(1)-a((4+1):-1:2);a;2*a(end)-a((end-1):-1:end-4)];
+
+  % Switching between the input options differentiate
+  switch upper(order)
+    case 'ZERO'
+      Deriv = (Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
+      Deriv = Deriv(3:end-2);
+    case 'FIRST'
+      Deriv = (1/T).*(Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
+      Deriv = Deriv(3:end-2);
+    case 'SECOND'
+      Deriv = (1/T^2).*(Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
+      Deriv = Deriv(3:end-2);
+    case 'THIRD'
+      Deriv = (1/T^3).*(Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
+      Deriv = Deriv(3:end-2);
+    case 'FOURTH'
+      Deriv = (1/T^4).*(Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
+      Deriv = Deriv(3:end-2);
+    otherwise
+      error('### Unknown order %s', order);
+  end
+
+end