|
2414
|
1 function [X,T]=detrend(t,X,p) |
|
|
2 % DETREND removes the trend from data, NaN's are considered as missing values |
|
|
3 % |
|
|
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 |
|
|
9 % |
|
|
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) |
|
|
14 % p must be a scalar |
|
|
15 % |
|
|
16 % [...]=detrend(X,0) |
|
2881
|
17 % [...]=detrend(X,'constant') |
|
2414
|
18 % removes the mean |
|
|
19 % |
|
|
20 % [...]=detrend(X,p) |
|
|
21 % removes polynomial of order p (default p=1) |
|
|
22 % |
|
|
23 % [...]=detrend(X,1) - default |
|
2881
|
24 % [...]=detrend(X,'linear') |
|
2414
|
25 % removes linear trend |
|
|
26 % |
|
|
27 % [X,T]=detrend(...) |
|
|
28 % |
|
|
29 % X is the detrended data |
|
|
30 % T is the removed trend |
|
|
31 % |
|
2881
|
32 % see also: SUMSKIPNAN, ZSCORE |
|
2414
|
33 |
|
4907
|
34 % Copyright (C) 1995, 1996 Kurt Hornik <Kurt.Hornik@ci.tuwien.ac.at> |
|
|
35 % $Id$ |
|
12580
|
36 % Copyright (C) 2001,2007 by Alois Schloegl <alois.schloegl@gmail.com> |
|
4907
|
37 % This function is part of the TSA-toolbox |
|
11507
|
38 % http://pub.ist.ac.at/~schloegl/matlab/tsa/ |
|
12580
|
39 % Copyright (C) 1997, 1998, 2008 by Alois Schloegl <alois.schloegl@gmail.com> |
|
4907
|
40 % |
|
|
41 % This program is free software: you can redistribute it and/or modify |
|
2881
|
42 % it under the terms of the GNU General Public License as published by |
|
4907
|
43 % the Free Software Foundation, either version 3 of the License, or |
|
2881
|
44 % (at your option) any later version. |
|
2414
|
45 % |
|
2881
|
46 % This program is distributed in the hope that it will be useful, |
|
|
47 % but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
48 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
49 % GNU General Public License for more details. |
|
2414
|
50 % |
|
2881
|
51 % You should have received a copy of the GNU General Public License |
|
4907
|
52 % along with this program. If not, see <http://www.gnu.org/licenses/>. |
|
2414
|
53 |
|
|
54 |
|
|
55 if (nargin == 1) |
|
|
56 p = 1; |
|
|
57 X = t; |
|
|
58 t = []; |
|
|
59 elseif (nargin == 2) |
|
2881
|
60 if strcmpi(X,'constant'), |
|
|
61 p = 0; |
|
|
62 X = t; |
|
|
63 t = []; |
|
|
64 elseif strcmpi(X,'linear'), |
|
|
65 p = 1; |
|
|
66 X = t; |
|
|
67 t = []; |
|
|
68 elseif ischar(X) |
|
|
69 error('unknown 2nd input argument'); |
|
|
70 elseif all(size(X)==1), |
|
2414
|
71 p = X; |
|
|
72 X = t; |
|
|
73 t = []; |
|
|
74 else |
|
|
75 p = 1; |
|
|
76 end; |
|
|
77 elseif (nargin == 3) |
|
2881
|
78 if ischar(X), |
|
|
79 warning('input arguments are not supported'); |
|
|
80 end; |
|
2414
|
81 |
|
|
82 elseif (nargin > 3) |
|
|
83 fprintf (1,'usage: detrend (x [, p])\n'); |
|
|
84 end; |
|
|
85 |
|
|
86 % check data, must be in culomn order |
|
|
87 [m, n] = size (X); |
|
|
88 if (m == 1) |
|
|
89 X = X'; |
|
|
90 r=n; |
|
|
91 else |
|
|
92 r=m; |
|
|
93 end |
|
|
94 % check time scale |
|
|
95 if isempty(t), |
|
|
96 t = (1:r).'; % make time scale |
|
|
97 elseif ~all(size(t)==size(X)) |
|
|
98 t = t(:); |
|
|
99 end; |
|
|
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'); |
|
|
103 end; |
|
|
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'); |
|
|
107 end |
|
|
108 |
|
|
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 |
|
|
112 |
|
|
113 |
|
|
114 if (size(t,2)>1), % for multiple time scales |
|
|
115 for k=1:size(X,2), |
|
|
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)); |
|
|
119 end; |
|
|
120 |
|
|
121 else % if only one time scale is used |
|
|
122 b = (t * ones (1, p + 1)) .^ (ones (length(t),1) * (0 : p)); |
|
|
123 for k=1:size(X,2), |
|
|
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 |
|
|
128 end; |
|
|
129 end; |
|
|
130 X = X-T; % 3nd alternative |
|
|
131 |
|
|
132 if (m == 1) |
|
|
133 X = X'; |
|
|
134 T = T'; |
|
|
135 end |
|
|
136 else % needs less memory |
|
|
137 if (size(t,2)>1), % for multiple time scales |
|
|
138 for k = 1:size(X,2), |
|
|
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)); |
|
|
142 end; |
|
|
143 else % if only one time scale is used |
|
|
144 b = (t * ones (1, p + 1)) .^ (ones (length(t),1) * (0 : p)); |
|
|
145 for k = 1:size(X,2), |
|
|
146 idx = find(~isnan(X(:,k))); |
|
|
147 X(idx,k) = X(idx,k) - b(idx,:) * (b(idx,:) \ X(idx,k)); |
|
|
148 end; |
|
|
149 end; |
|
|
150 |
|
|
151 if (m == 1) |
|
|
152 X = X'; |
|
|
153 end |
|
|
154 end; |
|
|
155 |
|
|
156 |
|
|
157 |
