Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / nan-2.5.5 / detrend.m

diff --git a/octave_packages/nan-2.5.5/detrend.m b/octave_packages/nan-2.5.5/detrend.m
new file mode 100644 (file)
index 0000000..a49ac72
--- /dev/null
+++ b/octave_packages/nan-2.5.5/detrend.m
@@ -0,0 +1,157 @@
+function [X,T]=detrend(t,X,p)
+% DETREND removes the trend from data, NaN's are considered as missing values
+% 
+% DETREND is fully compatible to previous Matlab and Octave DETREND with the following features added:
+% - handles NaN's by assuming that these are missing values
+% - handles unequally spaced data
+% - second output parameter gives the trend of the data
+% - compatible to Matlab and Octave 
+%
+% [...]=detrend([t,] X [,p])
+%      removes trend for unequally spaced data
+%      t represents the time points
+%      X(i) is the value at time t(i)
+%      p must be a scalar
+%
+% [...]=detrend(X,0)
+% [...]=detrend(X,'constant')
+%      removes the mean
+%
+% [...]=detrend(X,p)
+%      removes polynomial of order p (default p=1)
+%
+% [...]=detrend(X,1) - default
+% [...]=detrend(X,'linear')
+%      removes linear trend 
+%
+% [X,T]=detrend(...) 
+%
+% X is the detrended data
+% T is the removed trend
+% 
+% see also: SUMSKIPNAN, ZSCORE         
+
+%    This program is free software; you can redistribute it and/or modify
+%    it under the terms of the GNU General Public License as published by
+%    the Free Software Foundation; either version 2 of the License, or
+%    (at your option) any later version.
+%
+%    This program is distributed in the hope that it will be useful,
+%    but WITHOUT ANY WARRANTY; without even the implied warranty of
+%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%    GNU General Public License for more details.
+%
+%    You should have received a copy of the GNU General Public License
+%    along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+
+% Copyright (C) 1995, 1996 Kurt Hornik <Kurt.Hornik@ci.tuwien.ac.at>
+%       $Id: detrend.m 8223 2011-04-20 09:16:06Z schloegl $
+%       Copyright (C) 2001,2007 by Alois Schloegl <alois.schloegl@gmail.com>   
+%       This function is part of the NaN-toolbox
+%       http://pub.ist.ac.at/~schloegl/matlab/NaN/
+
+
+if (nargin == 1)
+       p = 1;
+               X = t;
+       t = [];
+elseif (nargin == 2)
+       if strcmpi(X,'constant'),
+               p = 0; 
+               X = t; 
+               t = [];
+       elseif strcmpi(X,'linear'),
+               p = 1; 
+               X = t; 
+               t = [];
+       elseif ischar(X)
+               error('unknown 2nd input argument');    
+        elseif all(size(X)==1), 
+                p = X;
+                X = t;
+                t = [];
+        else
+                p = 1;
+       end;            
+elseif (nargin == 3)
+        if ischar(X),
+               warning('input arguments are not supported');   
+        end; 
+        
+elseif (nargin > 3)
+       fprintf (1,'usage: detrend (x [, p])\n');
+end;
+
+% check data, must be in culomn order
+[m, n] = size (X);
+if (m == 1)
+        X = X';
+        r=n;
+else
+        r=m;
+end
+% check time scale
+if isempty(t),
+       t = (1:r).'; % make time scale 
+elseif ~all(size(t)==size(X)) 
+        t = t(:);
+end;
+% check dimension of t and X
+if ~all(size(X,1)==size(t,1))
+        fprintf (2,'detrend: size(t,1) must same as size(x,1) \n');
+end;
+% check the order of the polynomial 
+if (~(all(size(p)==1) && (p == round (p)) && (p >= 0)))
+       fprintf (2,'detrend:  p must be a nonnegative integer\n');
+end
+
+if (nargout>1)  , % needs more memory
+        T = zeros(size(X))+nan; 
+        %T=repmat(nan,size(X)); % not supported by Octave 2.0.16
+        
+        
+       if (size(t,2)>1),       % for multiple time scales
+                for k=1:size(X,2),
+                       idx=find(~isnan(X(:,k)));
+                        b = (t(idx,k) * ones (1, p + 1)) .^ (ones (length(idx),1) * (0 : p));
+                       T(idx,k) = b * (b \ X(idx,k));
+                end;
+                       
+        else                   % if only one time scale is used
+               b = (t * ones (1, p + 1)) .^ (ones (length(t),1) * (0 : p));
+                for k=1:size(X,2),
+                       idx=find(~isnan(X(:,k)));
+                       T(idx,k) = b(idx,:) * (b(idx,:) \ X(idx,k));
+                       %X(idx,k) = X(idx,k) - T(idx,k); % 1st alternative implementation
+                        %X(:,k) = X(:,k) - T(:,k); % 2nd alternative 
+                end;
+       end;
+        X = X-T;  % 3nd alternative 
+        
+       if (m == 1)
+               X = X';
+               T = T';
+       end
+else % needs less memory
+       if (size(t,2)>1),       % for multiple time scales
+                for k = 1:size(X,2),
+                       idx = find(~isnan(X(:,k)));
+                        b = (t(idx,k) * ones (1, p + 1)) .^ (ones (length(idx),1) * (0 : p));
+                       X(idx,k) = X(idx,k) -  b * (b \ X(idx,k));
+                end;
+        else                   % if only one time scale is used
+               b = (t * ones (1, p + 1)) .^ (ones (length(t),1) * (0 : p));
+                for k = 1:size(X,2),
+                       idx = find(~isnan(X(:,k)));
+                       X(idx,k) = X(idx,k) - b(idx,:) * (b(idx,:) \ X(idx,k));
+                end;
+       end;
+
+       if (m == 1)
+                X = X';
+       end
+end;
+
+
+