| der(X,DT,W,DEG,PAD,DIM)
|
function [DXDT,DXDTstd] = der(X,DT,W,DEG,PAD,DIM)
%DER Estimates central first derivative from local slope ignoring NaNs.
%
% SYNTAX:
% DXDT = der(X);
% DXDT = der(X,DT);
% DXDT = der(X,DT,W);
% DXDT = der(X,DT,W,DEG);
% DXDT = der(X,DT,W,DEG,PAD);
% DXDT = der(X,DT,W,DEG,PAD,DIM);
% DXDT = der(X,T,...); % Instead of DT
% [DXDT,DXDTstd] = der(...);
%
% INPUT:
% X - Data array with/without NaNs. Must have length N in the
% specified DIM. For example N-by-M.
% DT - * An scalar with the sampling interval (argument is then
% T=(0:N-1)*DT).
% * A vector T of length N with the X argument: X(T).
% DEFAULT: 1 (then T = 0:N-1)
% W - * An scalar with the window width to be used at estimating the
% linear regression slope in units of DT. W/DT is forced to be
% odd. If T was given, W is the number of elements in the
% window.
% * A vector with the M values for each column on X (if DIM=1).
% DEFAULT: 3*DT (or 3 if T was given)
% DEG - Polynomial degree to be used on the linear regression.
% DEFAULT: 1 (fits a line on the small window)
% PAD - Padding edge option. One of (see NOTE below):
% 'circular' - Treats series as periodic.
% 'replicate' - Repeats border elements.
% 'symmetric' - Pads with mirror reflection of X.
% 'mean' - Pads with mean series value.
% 'zero' - Pads with zeros.
% 'nan' - Pads with NaNs (do not pads!)
% An scalar - Padding value.
% A INTERP1 method - Extrapolates by using 'linear', 'cubic'
% 'spline', 'nearest' or 'pchip' method.
% A function handle - Extrapolation function with same syntax
% as x0=INTERP1(x,t,t0).
% DEFAULT: NaN (do not pads at all)
% DIM - Specifies the location dimension of the time series in the X
% array.
% DEFAULT: first non-singleton (generally in columns, DIM=1)
%
% OUTPUT:
% DXDT - First central derivative of X with respect to T at T. Its
% has the same size as X.
% DXDTstd - Estimation of the error on the DXDT estimation. Same size
% as DXDT.
%
% DESCRIPTION:
% This program estimate the first central derivative of an array by
% fitting locally a polynomial (linear regression) and then estimate
% its slope at T.
%
% The locality is achieved by using only the surrounding elements of
% each T element no further than half of the specified rectangular
% window width.
%
% The result is an smoothed derivative where the degree of smoothing
% increases when the window width is increase. The degree of the
% polynomial may be specified too, greater than 3 is not recommended,
% nut is cannot be lower than the number of elements in the window. The
% latter must be less than the length of the array.
%
% The error estimation of the each estimated derivative is obtained by
% assuming X is a random variable. See POLYFIT for details.
%
% NOTE:
% * Optional inputs use its DEFAULT value when not given or [].
% * Optional outputs may or not be called.
% * Property names may be shortened, as long as they are unambiguous.
% Similarly for string values. Case is ignored.
% * Padding options is important for the derivatives at edges values.
% If your data has a constant derivative on the edges 'replicate' or
% 'symmetric' are recommended. If your data is periodic then use
% 'circular'. The best way is to try different options and keep the
% "best" for your data.
%
% EXAMPLE:
% % Compares DEG 1 and 3 for window width (in number of elements):
% w = 11;
% % Data.
% pnoise = 30;
% n = 100;
% t = linspace(0,4*pi,n)';
% x = sin(t) + (pnoise/100)*(rand(n,1)*2-1); x(end)=[]; t(end)=[];
% dt = t(2)-t(1);
% % Linear derivative.
% [y1,y1std] = der(x,dt,dt*w,1,'re');
% % Cubic derivative.
% [y3,y3std] = der(x,dt,dt*w,3,'re');
% % Draws.
% figure(1),clf
% subplot(311)
% plot(t,sin(t),t,x,'.:'), legend('sin(t)','x')
% title 'Data'
% subplot(312)
% plot(t,cos(t),t,y1,'r.:',t,y3,'g--.',...
% t,[y1+y1std y1-y1std],'r',t,[y3+y3std y3-y3std],'g')
% legend('dsin(t)/dt','der(x), DEG=1','der(x), DEG=3')
% title(sprintf(['Derivatives from %.0f surrounding elements, '...
% 'plus error intervals'],w))
% subplot(313)
% stem(t,y1-cos(t),'r')
% hold on
% stem(t,y3-cos(t),'g')
% hold off
% legend('DEG=1','DEG=3')
% title 'Derivative error with respect to cos(t)'
%
% SEE ALSO:
% DIFF, GRADIENT
% and
% MOVINGSLOPE by John D'Errico and NANDERIVATIVE by Carlos Vargas
% at http://www.mathworks.com/matlabcentral/fileexchange
%
%
% ---
% MFILE: der.m
% VERSION: 1.0 (Nov 12, 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. (Nov 12, 2009)
% DISCLAIMER:
% der.m is provided "as is" without warranty of any kind, under the
% revised BSD license.
% Copyright (c) 2009 Carlos Adrian Vargas Aguilera
% INPUTS CHECK-IN
% -------------------------------------------------------------------------
% Defaults.
DXDT = [];
DXDTstd = [];
dDT = 1;
dW = 3; % *DT
dDEG = 1;
dPAD = NaN;
dDIM = 1; % or 2 if X is row vector.
% Checks arguments.
if (nargin<1)
error('CVARGAS:der:notEnoughInputs',...
'At least 1 input is required')
end
if (nargin>6)
error('CVARGAS:der:tooManyInputs',...
'At most 6 inputs are allowed')
end
if (nargout>2)
error('CVARGAS:der:tooManyOutputs',...
'At most 2 outputs are allowed')
end
% Checks X.
if isempty(X), return, end
Xsize = size(X);
% Checks PAD.
if (nargin<5) || isempty(PAD)
PAD = dPAD;
elseif ischar(PAD)
l = length(PAD);
if strncmpi(PAD,'circular',max(l,2))
PAD = 'circular';
elseif strncmpi(PAD,'replicate',l)
PAD = 'replicate';
elseif strncmpi(PAD,'symmetric',l)
PAD = 'symmetric';
elseif strncmpi(PAD,'mean',l)
PAD = 'mean';
elseif strncmpi(PAD,'zero',l)
PAD = 0;
elseif strncmpi(PAD,'nan',max(l,2))
PAD = NaN;
elseif strncmpi(PAD,'linear',l)
PAD = @(x,y,xo) interp1(x,y,xo,'linear');
elseif strncmpi(PAD,'cubic',max(l,2))
PAD = @(x,y,xo) interp1(x,y,xo,'cubic');
elseif strncmpi(PAD,'spline',l)
PAD = @(x,y,xo) interp1(x,y,xo,'spline');
elseif strncmpi(PAD,'nearest',max(l,2))
PAD = @(x,y,xo) interp1(x,y,xo,'nearest');
elseif strncmpi(PAD,'pchip',l)
PAD = @(x,y,xo) interp1(x,y,xo,'pchip');
elseif (exist(PAD,'file')==2)
PAD = str2func(PAD);
else
warning('CVARGAS:der:incorrectPadString',...
'Unrecognize %s padding. Default used.',PAD)
end
elseif (isnumeric(PAD) || islogical(PAD)) && ...
(length(PAD)==1)
% Continue
elseif (length(PAD)==1) && isa(PAD,'function_handle')
% Continue
else
warning('CVARGAS:der:incorrectPadType',...
'Incorrect padding type. Default used.')
PAD = dPAD;
end
% Checks DEG.
if (nargin<4) || isempty(DEG)
DEG = dDEG;
end
% Checks DIM.
if (nargin<6) || isempty(DIM)
DIM = find(Xsize~=1,1); % DIM = min(find(S~=1));
if isempty(DIM), DIM = dDIM; end
end
% Series length and dimension number.
m = Xsize(DIM);
Ndim = length(Xsize);
% Reshape series in columns.
n = prod(Xsize)/m; % Number of series:
X = reshape(permute(X,[DIM 1:DIM-1 DIM+1:Ndim]),m,n);
% Checks DT and T.
if (nargin<2) || isempty(DT)
DT = dDT;
T = repmat((0:m-1)'*DT,1,n);
else
nDT = numel(DT);
if nDT==1
T = repmat((0:m-1)'*DT,1,n);
elseif nDT==m
T = repmat(DT(:),1,n);
DT = [];
elseif nDT==n
DT = DT(:).';
T = repmat((0:m-1)',1,n).*repmat(DT,m,1);
else
T = reshape(permute(DT,[DIM 1:DIM-1 DIM+1:Ndim]),m,n);
DT = [];
end
end
% Checks W.
if (nargin<3) || isempty(W)
W = repmat(dW,1,n);
else
W = W(:).';
if ~isempty(DT)
W = W./DT;
if any(abs(W-round(W))>eps*1e4)
warning('CVARGAS:der:incorrectWMultiple',...
'W must be a multiple of sampling time DT. ROUND used.')
end
W = round(W);
end
ind = (W<3);
if any(ind)
warning('CVARGAS:der:incorrectWValue',...
'W/DT too low. Used default.')
W(ind) = dW;
end
ind = ~rem(W,2);
if any(ind)
warning('CVARGAS:der:incorrectWOdd',...
'W/DT forced to be odd.')
W(ind) = W(ind)+ind;
end
if length(W)==1
W = repmat(W,1,n);
end
end
% Checks DEG and W.
if DEG>=W
warning('CVARGAS:der:incorrectDegW',...
'DEG must be lower than W. NaNs output.')
DXDT = NaN(Xsize);
if nargout>1, DXDTstd = DXDT; end
return
end
% Checks W and X.
if W>=m
warning('CVARGAS:der:tooLargeW',...
'W must be less than the array length. Continue anyway.')
end
% -------------------------------------------------------------------------
% MAIN
% -------------------------------------------------------------------------
% Window semi width.
W2 = floor(W/2);
% Pads W2 elements at each edges.
T = [interp1((1:2)',T(1:2,:) ,repmat((-W2+1:0)' ,1,n),[],'extrap'); T;
interp1((1:2)',T(m-1:m,:),repmat((0:W2-1)'+3,1,n),[],'extrap')];
if ischar(PAD)
switch PAD(1)
case 'r'
% 'replicates'
X = [repmat(X(1,:),W2,1); X; repmat(X(m,:),W2,1)];
case 'c'
% 'circular'
X = [X(mod((m-W2+1:m)-1,m)+1,:); X; X(mod((1:W2)-1,m)+1,:)];
case 's'
% 'symmetric'
temp = [1:m m-1:1];
X = [X(temp(mod((W2:-1: 1)-1,2*m)+1),:); X; ...
X(temp(mod((m :-1:m-W2+1)-1,2*m)+1),:)];
clear temp
case 'm'
% 'mean'
X = [repmat(nanmean(X,1),W2,1); X; repmat(nanmean(X,1),W2,1)];
end
elseif isa(PAD,'function_handle')
try
X = [PAD(T(W2+(1:m),:),X,T( 1:W2 ,:)); X; ...
PAD(T(W2+(1:m),:),X,T(W2+m+(1:W2),:))];
catch
warning('CVARGAS:der:errorDuringPadFunctionEvaluation',...
'An error ocurred during "%s" evaluation. Returned NaNs.',func2str(PAD))
DXDT = NaN(Xsize);
if nargout>1, DXDTstd = DXDT; end
return
end
else
% An scalar
X = [repmat(PAD,W2,n); X; repmat(PAD,W2,n)];
end
% Gets derivatives..
DXDT = NaN(m,n);
inan = isnan(X);
if nargout==1
% Do not estimate slope error.
for col = 1:n
for row = 1:m
r = (-W2(col):W2(col))' + row + W2;
bad = inan(r,col);
if sum(~bad)>DEG
slope = polyfit(T(r(~bad),col)-T(row+W2,col),X(r(~bad),col),DEG);
DXDT(row,col) = slope(DEG);
end
end
end
else
% Estimate slope error.
DXDTstd = DXDT;
for col = 1:n
for row = 1:m
r = (-W2(col):W2(col))' + row + W2;
bad = inan(r,col);
if sum(~bad)>DEG
[slope,E] = polyfit(T(r(~bad),col)-T(row+W2,col),X(r(~bad),col),DEG);
DXDT(row,col) = slope(DEG);
E = sqrt(diag(abs((inv(E.R)*inv(E.R')))*E.normr^2/E.df));
DXDTstd(row,col) = E(DEG);
end
end
end
end
% OUTPUTS CHECK-OUT
% -------------------------------------------------------------------------
% Returns to original shape.
DXDT = ipermute(reshape(DXDT,...
[size(DXDT,1) Xsize([DIM+1:Ndim 1:DIM-1])]),[DIM 1:DIM-1 DIM+1:Ndim]);
if nargout>1
DXDTstd = ipermute(reshape(DXDTstd,...
[size(DXDTstd,1) Xsize([DIM+1:Ndim 1:DIM-1])]),[DIM 1:DIM-1 DIM+1:Ndim]);
end
% [EOF] der.m |
|
