code by jitendra

1 files changed, 24 insertions, 14 deletions

diff --git a/macros/sgolay.sci b/macros/sgolay.sci
index 1156d40..09062a7 100644
--- a/macros/sgolay.sci
+++ b/macros/sgolay.sci
@@ -10,7 +10,6 @@ function F = sgolay (p, n, m, ts)
 //m: positive integer less than 2^31 or logical
 //ts: real or complex value
 //Description
-//This is an Octave function.
 //This function computes the filter coefficients for all Savitzsky-Golay smoothing filters of order p for length n (odd). 
 //m can be used in order to get directly the mth derivative; ts is a scaling factor.
 //Examples
@@ -20,21 +19,32 @@ function F = sgolay (p, n, m, ts)
 //   0.33333   0.33333   0.33333
 //  -0.16667   0.33333   0.83333
 
-funcprot(0);
-rhs = argn(2)
-
-if(rhs<2 | rhs>4)
+if(argn(2)<2 | argn(2)>4) then
 error("Wrong number of input arguments.")
+elseif ((n-fix(n./2).*2)~=1)  then
+error ("sgolay needs an odd filter length n");
+elseif (p>=n) then
+     error ("sgolay needs filter length n larger than polynomial order p");
+ end
+
+if (argn(2)==2) then
+    m=0; ts=1;
+end
+if (argn(2)==3) then
+    ts=1;
 end
 
-	select(rhs)
-	case 2 then
-	F = callOctave("sgolay",p,n)
-	case 3 then
-	F = callOctave("sgolay",p,n,m)
-	case 4 then
-	F = callOctave("sgolay",p,n,m,ts)
-	end
-endfunction
+if length(m) > 1, error("weight vector unimplemented"); end
+
+ F = zeros (n, n);
+ k = floor (n/2);
+    for row = 1:k+1
+      C = ( [(1:n)-row]'*ones(1,p+1) ) .^ ( ones(n,1)*[0:p] );   
+      A = pinv(C);   
+      F(row,:) = A(1+m,:);
+    end
 
+  F(k+2:n,:) = (-1)^m*F(k:-1:1,n:-1:1);
+  F =  F * ( prod(1:m) / (ts^m) );
 
+endfunction