1 function [X,T]=detrend(t,X,p)
2 % DETREND removes the trend from data, NaN's are considered as missing values
4 % DETREND is fully compatible to previous Matlab and Octave DETREND with the following features added:
5 % - handles NaN's by assuming that these are missing values
6 % - handles unequally spaced data
7 % - second output parameter gives the trend of the data
8 % - compatible to Matlab and Octave
10 % [...]=detrend([t,] X [,p])
11 % removes trend for unequally spaced data
12 % t represents the time points
13 % X(i) is the value at time t(i)
17 % [...]=detrend(X,'constant')
21 % removes polynomial of order p (default p=1)
23 % [...]=detrend(X,1) - default
24 % [...]=detrend(X,'linear')
25 % removes linear trend
29 % X is the detrended data
30 % T is the removed trend
32 % see also: SUMSKIPNAN, ZSCORE
34 % This program is free software; you can redistribute it and/or modify
35 % it under the terms of the GNU General Public License as published by
36 % the Free Software Foundation; either version 2 of the License, or
37 % (at your option) any later version.
39 % This program is distributed in the hope that it will be useful,
40 % but WITHOUT ANY WARRANTY; without even the implied warranty of
41 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
42 % GNU General Public License for more details.
44 % You should have received a copy of the GNU General Public License
45 % along with this program; If not, see <http://www.gnu.org/licenses/>.
48 % Copyright (C) 1995, 1996 Kurt Hornik <Kurt.Hornik@ci.tuwien.ac.at>
49 % $Id: detrend.m 8223 2011-04-20 09:16:06Z schloegl $
50 % Copyright (C) 2001,2007 by Alois Schloegl <alois.schloegl@gmail.com>
51 % This function is part of the NaN-toolbox
52 % http://pub.ist.ac.at/~schloegl/matlab/NaN/
60 if strcmpi(X,'constant'),
64 elseif strcmpi(X,'linear'),
69 error('unknown 2nd input argument');
70 elseif all(size(X)==1),
79 warning('input arguments are not supported');
83 fprintf (1,'usage: detrend (x [, p])\n');
86 % check data, must be in culomn order
96 t = (1:r).'; % make time scale
97 elseif ~all(size(t)==size(X))
100 % check dimension of t and X
101 if ~all(size(X,1)==size(t,1))
102 fprintf (2,'detrend: size(t,1) must same as size(x,1) \n');
104 % check the order of the polynomial
105 if (~(all(size(p)==1) && (p == round (p)) && (p >= 0)))
106 fprintf (2,'detrend: p must be a nonnegative integer\n');
109 if (nargout>1) , % needs more memory
110 T = zeros(size(X))+nan;
111 %T=repmat(nan,size(X)); % not supported by Octave 2.0.16
114 if (size(t,2)>1), % for multiple time scales
116 idx=find(~isnan(X(:,k)));
117 b = (t(idx,k) * ones (1, p + 1)) .^ (ones (length(idx),1) * (0 : p));
118 T(idx,k) = b * (b \ X(idx,k));
121 else % if only one time scale is used
122 b = (t * ones (1, p + 1)) .^ (ones (length(t),1) * (0 : p));
124 idx=find(~isnan(X(:,k)));
125 T(idx,k) = b(idx,:) * (b(idx,:) \ X(idx,k));
126 %X(idx,k) = X(idx,k) - T(idx,k); % 1st alternative implementation
127 %X(:,k) = X(:,k) - T(:,k); % 2nd alternative
130 X = X-T; % 3nd alternative
136 else % needs less memory
137 if (size(t,2)>1), % for multiple time scales
139 idx = find(~isnan(X(:,k)));
140 b = (t(idx,k) * ones (1, p + 1)) .^ (ones (length(idx),1) * (0 : p));
141 X(idx,k) = X(idx,k) - b * (b \ X(idx,k));
143 else % if only one time scale is used
144 b = (t * ones (1, p + 1)) .^ (ones (length(t),1) * (0 : p));
146 idx = find(~isnan(X(:,k)));
147 X(idx,k) = X(idx,k) - b(idx,:) * (b(idx,:) \ X(idx,k));