Implement detrend1.sci in Scilab

1 files changed, 71 insertions, 22 deletions

diff --git a/macros/detrend1.sci b/macros/detrend1.sci
index 52ac26f..29b0201 100644
--- a/macros/detrend1.sci
+++ b/macros/detrend1.sci
@@ -1,23 +1,72 @@
-function y = detrend1(x, varargin)
-//This function removes the best fit of a polynomial of order P from the data X
-//Calling Sequence
-//detrend1(X,P)
-//Parameters 
-//X: Input vecor or matrix.
-//P: The order of polnomial
-//Description
-//If X is a vector, 'detrend1(X, P)' removes the best fit of apolynomial of order P from the data X.If X is a matrix, 'detrend1(X, P)' does the same for each column in X.
-//
-//The second argument P is optional.  If it is not specified, a value of 1 is assumed.  This corresponds to removing a linear trend.
-//The order of the polynomial can also be given as a string, in which case P must be either "constant" (corresponds to 'P=0') or "linear"(corresponds to 'P=1')
-	rhs= argn(2);
-	if(rhs<1 | rhs> 2)
-		error("Wrong number of input arguments");
-	end
-	select(rhs)
-	case 1 then
-		y= callOctave("detrend", x);
-	case 2 then
-		y= callOctave("detrend", x , varargin(1));
-	end
+function y = detrend1 (x, p)
+//Remove the best fit of a polynomial of order p from the data x.
+//Calling Sequence:
+//detrend1(x,p)
+//Parameters: 
+//x: Input vecor or matrix
+//p: The order of polnomial
+//Description:
+//If X is a vector, 'detrend1(X, P)' removes the best fit of apolynomial of order P from the data X.
+//If X is a matrix, 'detrend1(X, P)' does the same for each column in X.
+//The second argument p is optional.  If it is not specified, a value of 1 is assumed.  This corresponds to removing a linear trend.
+//The order of the polynomial can also be given as a string, in which case p must be either "constant" (corresponds to 'P=0') or "linear" (corresponds to 'P=1')
+//Example:
+//detrend1([1, 6, 9])
+//ans = [ -0.3333, 0.6667, -0.3333]
+
+  funcprot(0);
+  rhs = argn(2);
+  if (rhs < 1 | rhs > 2)
+    error("detrend1: wrong number of input arguments");
+  end
+  
+  if (rhs == 1)
+    p = 1;
+  end
+
+  if (~or(type(x)==[1 5 8]) | ndims (x) > 2)
+    error ("detrend1: X must be a numeric vector or matrix");
+  end
+
+  if (type(p) == 10 & ~strcmp(p, "constant", 'i'))
+    p = 0;
+  elseif (type(p) == 10 & ~strcmp(p, "linear", 'i'))
+    p = 1;
+  elseif (~isscalar(p) | p < 0 | p ~= fix (p))
+    error ("detrend1: P must be constant, linear, or a positive integer");
+  end
+
+  [m, n] = size (x);
+  if (m == 1)
+    x = x.';
+  end
+
+  r = size(x, 'r');
+  b = ((1 : r).' * ones (1, p + 1)) .^ (ones (r, 1) * (0 : p));
+  y = x - b * (b \ x);
+
+  if (m == 1)
+    y = y.';
+  end
+
 endfunction
+//
+//input validation:
+//assert_checkerror("detrend1()", "detrend1: wrong number of input arguments");
+//a = "string";
+//assert_checkerror("detrend1(a)", "detrend1: X must be a numeric vector or matrix");
+//assert_checkerror("detrend1(%T)", "detrend1: X must be a numeric vector or matrix");
+//assert_checkerror("detrend1(1, -1)", "detrend1: P must be constant, linear, or a positive integer");
+//assert_checkerror("detrend1(1, 1.25)", "detrend1: P must be constant, linear, or a positive integer");
+
+//tests:
+//N = 32;
+//x = (0:1:N-1)/N + 2;
+//y = detrend1(x);
+//assert_checktrue(abs (y(:)) < 20*%eps);
+//assert_checkequal(detrend1([2, 5, 8]), detrend1([2. 5, 8], "linear"));
+//assert_checkequal(detrend1([2, 5, 8], 0), detrend1([2. 5, 8], "constant"));
+//assert_checkalmostequal(detrend1([1; 6; 9], "constant"), [-4.33333; 0.66666; 3.66666], 5*10^-5);
+//assert_checkalmostequal(detrend1([5, 12, 14; 8, 16, 14; 5, 10, 12]), [-1, -1.6666, -0.3333; 2, 3.3333, 0.6666; -1, -1.6666, -0.3333], 5*10^-4);
+//assert_checkalmostequal(detrend1([-5-5*%i; 2+%i; -4+3*%i], "linear"), [-2.1667-0.6667*%i; 4.3333+1.3333*%i; -2.1667-0.6667*%i], 5*10^-4);
+//assert_checkalmostequal(detrend1([5*%i, 1+2*%i,; -8, -1-6*%i], 0), [4+2.5*%i, 1+4*%i; -4-2.5*%i, -1-4*%i], 5*10^-4);