function dYdX = nanderivative(Y,X,N,K,MET)
% NANDERIVATIVE N-order derivative by interpolation.
%
% SYNTAX:
% dYdX = nanderivative(Y);
% dYdX = nanderivative(Y,X);
% dYdX = nanderivative(Y,X,N);
% dYdX = nanderivative(Y,X,N,K);
% dYdX = nanderivative(Y,X,N,K,MET);
%
% INPUT:
% Y - Vector or matrix data with/without Nans.
% X - Argument of the data where the derivative will be estimated. It
% should have the same dimensions as Y, or at least the same
% number of rows, or can be a single scalar specifying the
% increment.
% Default: 1 (unit increment)
% N - Differential order.
% Default: 1.
% K - Interval reduction factor (see DESCRIPTION below).
% Default: 100.
% MET - Method to be use in the interolation.
% Default: 'spline' ('linear' is also good)
%
% OUTPUT:
% dYdX - Array of the same size as Y with its N-order derivative.
%
% DESCRIPTION:
% This program estimates the N-order derivative of a set of data by
% interpolating the data in two nearby points around each argument, and
% then calculates the slope of the line crossing these points.
%
% Both points are at a distance (XO-XL)/K and (XR-XO)/K from the Xo
% point, where XL is the closer point at the left and XR at the right,
% respectively. This explains the K factor. Usually K=1 and
% MET='linear' are used by any numerical derivation method.
%
% Extrema point are used as one of these two nearby points, so
% extrapolation is avoided.
%
% NaNs elements of Y are just ignored.
%
% NOTES:
% * Optional inputs use its DEFAULT value when not given or [].
% * If Y is not smooth, try a bigger K.
% * The programs sets a warning if floating point rounding are
% occurring.
% * Error increase with the derivative order, and near the edges and
% NaNs points (see the example).
%
% EXAMPLE:
% x = linspace(0,10*pi).';
% fx = cos(x); fxnan = fx; fxnan(50:53) = NaN;
% dfx = -sin(x); % Theoretical derivative
% dfdx = nanderivative(fxnan,x,4); % Numerical derivative
% subplot(511)
% plot(x, fx,x,fxnan,'.'), ylabel('f(x) = cos(x)'), axis tight
% legend('Theory','nanderivative')
% subplot(512)
% plot(x, dfx,x,dfdx(:,1),'.'), ylabel('df/dx'), axis tight
% subplot(513)
% plot(x, -fx,x,dfdx(:,2),'.'), ylabel('d^2f/dx^2'), axis tight
% subplot(514)
% plot(x,-dfx,x,dfdx(:,3),'.'), ylabel('d^3f/dx^3'), axis tight
% subplot(515)
% plot(x, fx,x,dfdx(:,4),'.'), ylabel('d^4f/dx^4'), axis tight
%
% SEE ALSO:
% INTERP1, DIFF and GRADIENT.
%
%
% ---
% MFILE: nanderivative.m
% VERSION: 2.1 (Jul 29, 2009) (<a href="matlab:web('http://www.mathworks.com/matlabcentral/fileexchange/authors/11258')">download</a>)
% MATLAB: 7.7.0.471 (R2008b)
% AUTHOR: Carlos Adrian Vargas Aguilera (MEXICO)
% CONTACT: nubeobscura@hotmail.com
% REVISIONS:
% 1.0 Released. (Aug 22, 2008)
% 2.0 Rewritten code. Mayor changes with inputs and link axes. Fixed
% bugs with floating point error, monotonically
% increasing/decreasing X and if N is not 1 do not outputs all
% other derivatives anymore. (Jun 08, 2009)
% 2.1 Fixed bug with INTERP1. (Jul 29, 2009)
% DISCLAIMER:
% nanderivative.m is provided "as is" without warranty of any kind, under
% the revised BSD license.
% Copyright (c) 2008-2009 Carlos Adrian Vargas Aguilera
% INPUTS CHECK-IN
% -------------------------------------------------------------------------
% Checks number of inputs:
if nargin<1
error('CVARGAS:nanderivative:notEnoughInputs',...
'At least 1 input is required.')
elseif nargin>5
error('CVARGAS:nanderivative:tooManyInputs',...
'At most 5 inputs are allowed.')
elseif nargout>1
error('CVARGAS:nanderivative:tooManyOutputs',...
'At most 1 output is allowed.')
end
% Checks Y:
if isempty(Y)
dYdX = [];
return
end
% Checks X:
if nargin<2 || isempty(X)
X = 1;
end
% Checks N:
if nargin<3 || isempty(N)
N = 1;
end
% Checks K:
if nargin<4 || isempty(K)
K = 100;
else
K = abs(K);
end
% Checks MET:
if nargin<5 || isempty(MET)
MET = 'spline';
end
% Checks dimensions:
m = size(Y);
ms = prod(m);
if length(m)>2 || ms==1
error('CVARGAS:nanderivative:incorrectInputDimension',...
'Y should be a matrix or a vector with at least 2 values.')
end
% Forces column vector:
ndata = m(1);
if m(1)==1
Y = Y.';
ndata = m(2);
end
nseries = ms/ndata;
% Checks X:
[mX,nX] = size(X);
if nX==1
if mX==1
X = (1:ndata).';
mX = ndata;
end
X = repmat(X,1,nseries);
nX = nseries;
singleX = true;
else
if mX==1
X = repmat(X(:),1,nseries);
mX = nX;
nX = nseries;
singleX = true;
else
singleX = false;
end
end
if (mX~=ndata) || (nX~=nseries)
error('CVARGAS:nanderivative:incorrectXSize',...
'Incorrect X size.')
end
% Checks order: Fixed bug MAY 2009
Nn = length(N);
% Initialization:
dYdX = NaN(size(Y));
if ms~=ndata % matrix input
dYdX = repmat(dYdX,[1 1 Nn]);
else
dYdX = repmat(dYdX(:),[1 Nn]);
end
% -------------------------------------------------------------------------
% MAIN
% -------------------------------------------------------------------------
% Checks for NaNs:
nnan = ~isnan(X) & ~isnan(Y);
X(~nnan) = NaN;
Y(~nnan) = NaN;
doInterp = (sum(nnan,1)>1); % Fixed BUG JUL 2009
% Checks if series as equal Nans:
if singleX
Nnan = sum(nnan,2);
if any(Nnan(Nnan~=0)~=nseries)
singleX = false;
end
end
% Finds indixes of no NaNs extrema points:
[ini,ini] = max( nnan ,[],1);
[fin,fin] = max(flipud(nnan),[],1); fin = ndata-fin+1;
if singleX
X = X(:,1);
ini = ini(1);
fin = fin(1);
end
% Generates increment and argument for interpolation:
dX = diff(X,1,1);
% Checks monotonically increasing/decreasing:
if any(dX(:)==0)
error('CVARGAS:nanderivative:repeatedXValue',...
'All elements of X must be distint.')
end
ind = any(dX>0,1) & any(dX<0,1);
if any(ind)
error('CVARGAS:nanderivative:repeatedXValue',...
'X must be monotonically increasing or decreasing.')
end
[ind,ind] = min(abs(dX),[],1);
dX = dX(sub2ind(size(dX),ind,1:size(dX,2)));% Fixed BUG MAY 2009
dX = dX/K;
dX = repmat(dX,ndata,1);
Xi = [reshape(X-dX/2,1,[]); reshape(X+dX/2,1,[])];
Xi = reshape(Xi,2*ndata,[]);
Xi(2*ini-1,:) = X(ini,:);
Xi(2*fin,:) = X(fin,:);
% Checks for any floating point problem:
if any(max(10*eps*abs(X(:)))>=min(abs(dX(:)))) % fixed bug, OCT, 2008
warning('CVARGAS:nanderivative:floatingPointWarning',...
['Could be some floating point rounding problems. '...
'Recommend to use a K < ' int2str(K) '.'])
end
if singleX
X = repmat( X,1,nseries);
dX = repmat(dX,1,nseries);
end
for c = 1:Nn
for nn = 1:N(c) % Fixed bug MAY 2009
Yi = NaN(size(Xi));
if singleX
if doInterp(1)
Yi = interp1(X(nnan(:,1)),Y(nnan(:,1),:),Xi,MET);
end
else
for j = 1:nseries
if doInterp(j)
Yi(:,j) = interp1(X(nnan(:,j),j),Y(nnan(:,j),j),Xi(:,j),MET);
end
end
end
Yi = reshape(Yi,2,[]);
Y = diff(Yi,1,1);
Y = reshape(Y,ndata,[]);
Y = Y./dX;
Y(ini,:) = Y(ini,:)*2;
Y(fin,:) = Y(fin,:)*2;
end
if nseries>1
dYdX(:,:,c) = Y;
else
dYdX(:,c) = Y;
end
end
% OUTPUTS CHECK-OUT
% -------------------------------------------------------------------------
% Returns dimension:
if N==1 && m(1)==1
dYdX = dYdX.';
end
% [EOF] nanderivative.m
