octave_packages/nan-2.5.5/detrend.m

   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)
  17 % [...]=detrend(X,'constant')
  18 %       removes the mean
  19 %
  20 % [...]=detrend(X,p)
  21 %       removes polynomial of order p (default p=1)
  22 %
  23 % [...]=detrend(X,1) - default
  24 % [...]=detrend(X,'linear')
  25 %       removes linear trend
  26 %
  27 % [X,T]=detrend(...)
  28 %
  29 % X is the detrended data
  30 % T is the removed trend
  31 %
  32 % see also: SUMSKIPNAN, ZSCORE
  33
  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.
  38 %
  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.
  43 %
  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/>.
  46
  47
  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/
  53
  54
  55 if (nargin == 1)
  56         p = 1;
  57         X = t;
  58         t = [];
  59 elseif (nargin == 2)
  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),
  71                 p = X;
  72                 X = t;
  73                 t = [];
  74         else
  75                 p = 1;
  76         end;
  77 elseif (nargin == 3)
  78         if ischar(X),
  79                 warning('input arguments are not supported');
  80         end;
  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