--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/ltpda_dft/polyregz.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,149 @@
+(*
+
+/home/ghh/math/polyreg/polyregz.m  17.01.2008 Gerhard Heinzel AEI
+
+This is version 1.1
+
+history:
+ version 1.0 17.1.08 initial release
+ version 1.1 18.1.08 fixed two typos in comment
+
+Mathematica code to produce C code for efficient polynomial detrending
+
+The polynomials are defined in terms of z = 2i/(n-1)-1, i=0...n
+such that the range 0...n-1 is mapped to -1...+1 
+
+Thereby the high dynamic range of i^n is reduced to -1..1,
+and furthermore every second entry in the matrix of
+the normal equations and its inverse vanishes due to symmetry.
+
+The regression is done via explicit solution of the normal equations.
+
+These are of dimension (p+1) for degree p:
+
+( s_0 s_1 s_2 .. s_p    )  ( a_0 )  = ( b_0 ) 
+( s_1 s_2 s_3 .. s_p+1  )  ( a_1 )    ( b_1 )
+( ...                   )  ( ... )    ( ... )
+( s_p s_p+1 .. s_2p     )  ( a_p )    ( b_p )
+
+with 
+s_k = Sum(i=0...n-1) z^k,
+b_k = Sum(i=0...n-1) x_i z^k,
+a_k = The unknown coefficient of z^k in the best-fit regression.
+
+The Matrix with s depends only on p and n, not on the data.
+
+The Mathematica code below finds exact explicit expressions for 
+the Matrix and its Inverse and writes the corresponding C code.
+
+If any modifications of the C code are desired apart from "indent"ing,
+this file below should be altered but not the C code.
+
+*)
+
+<<clear.m
+
+(*
+<<polyregz.m
+*)
+
+z[i_,n_] := 2 i /(n-1) -1;
+
+s[k_,n_] := Sum[z[i,n]^k,{i,0,n-1}];
+
+m[n_,p_] := Table[s[i+j,n],{i,0,p},{j,0,p}];
+
+sol[p_] := sol[p] = Simplify[Inverse[m[n, p]]];
+
+cln[x_] := CoefficientList[Numerator[x],n];
+cld[x_] := CoefficientList[Denominator[x],n];
+
+max[x_] := Max[cln[x],cld[x]];
+
+pn[x_] := "0" /; x== 0;
+pn[x_] := ToString[CForm[x]] <> "L";
+
+ncln[x_] := Map[pn, N[cln[x] / max[x], 22]];
+ncld[x_] := Map[pn, N[cld[x] / max[x], 22]];
+
+polyprint[p_, x_] := p[[1]] /; Length[p]==1;
+polyprint[p_, x_] := "(" <> polyprint[Drop[p,1], x] <> " * " <> x <> " + " <> p[[1]] <> ")" /; Length[p]>1 ;
+
+cform[x_] := "0" /; FreeQ[x,n];
+cform[x_] := polyprint[ncln[x], "n"] <> " / " <> polyprint[ncld[x], "n"];
+
+wlist[sbeg_, list_, sep_, send_] := 
+  Write[fp, sbeg <> Apply[StringJoin, Riffle[list, sep]] <> send];
+
+out[p_] := Module[ {fn, s, i, j, s1, list},
+  fn="polyreg"<>ToString[p] <> ".c";
+  fp=OpenWrite[fn,FormatType->OutputForm, PageWidth->1024];
+  Write[fp,"void"];
+  Write[fp,"polyreg",p," (double *x, int nn, double *y, double *a)"];
+  Write[fp,"{"];
+  Write[fp,""];
+  Write[fp,"/*"];
+  Write[fp,"polynominal detrending of a time series, order ",p];
+  Write[fp,"machine-generated file, do not edit!"];
+  Write[fp,"made by polyregz.m (Mathematica )"];
+  Write[fp,"Gerhard Heinzel AEI 17.01.2008"];
+  Write[fp,""];
+  Write[fp,"x[]: input, read-only: time series to be detrended"];
+  Write[fp,"nn: input, read-only: length of x[]"];
+  Write[fp,"y[]: output: time series with trend subtracted"];
+  Write[fp,"a[]: fitting coefficients in for z^0, z^1, z^2,... with"];
+  Write[fp,"     z = 2*i/(nn-1)-1 ; i=0,...,nn-1"];
+  Write[fp,"*/"];
+  Write[fp,""];
+ Write[fp,"  long double n = nn, n1 = 2.L / (n - 1), z;"];
+  For[i=0, i<=p, i++, 
+    Write[fp,"  long double a",i,", temp",i,", sum",i," = 0;"];
+  ];  
+  Write[fp,"  int i;"];
+  Write[fp,"  for (i = 0; i < nn; i++)"];
+  Write[fp,"    {"];
+  Write[fp,"      z = n1 * i - 1.L;"];
+  list={"x[i]"};
+  For[i=0, i<=p, i++,
+     s1="      sum" <> ToString[i] <> " += ";
+     wlist[s1, list, " * ", ";"];
+     list=Append[list,"z"];
+  ];
+  Write[fp,"/* the above code is efficiently optimized by GCC4 */"];
+  Write[fp,""];
+  Write[fp,"    }"];
+
+  s=sol[p];
+  For[i=0, i<=p, i++,
+    list={};
+    For[j=0, j<=p, j++,
+      If[!FreeQ[s[[i+1,j+1]],n],
+	Write[fp, "temp", j, " = ", cform[s[[i+1,j+1]]], ";"];
+	list=Append[list,"temp" <> ToString[j] <> " * sum" <> ToString[j]];
+      ]  
+    ];  
+    s1 = "a" <> ToString[i] <> " = ";
+    wlist[s1, list, " + ", ";"];
+    ];
+  Write[fp,"  for (i = 0; i < nn; i++)"];
+  Write[fp,"    {"];
+  Write[fp,"      z = n1 * i - 1.L;"];
+  list={};
+  For[i=0, i<=p, i++,
+     list=Append[list, "a" <> ToString[i]];
+  ];
+  
+  Write[fp,"      y[i] = x[i] - ", polyprint[list, "z"], ";"];
+  Write[fp,"    }"];
+  For[i=0, i<=p, i++, 
+    Write[fp,"  a[",i,"] = a",i,";"];
+  ];
+  Write[fp,"}"];
+  Close[fp];
+];
+
+For[i=0, i<=10, i++, 
+  out[i];
+  Print["degree ",i," done."];
+]  
+