X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fnan-2.5.5%2Fdetrend.m;fp=octave_packages%2Fnan-2.5.5%2Fdetrend.m;h=a49ac725e8707a03b7ec0570364ae324d4f86147;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/nan-2.5.5/detrend.m b/octave_packages/nan-2.5.5/detrend.m
new file mode 100644
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;
+
+
+