Mercurial > hg > ltpda

comparison src/ltpda_dft/polyregz.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 (*
2
3 /home/ghh/math/polyreg/polyregz.m 17.01.2008 Gerhard Heinzel AEI
4
5 This is version 1.1
6
7 history:
8 version 1.0 17.1.08 initial release
9 version 1.1 18.1.08 fixed two typos in comment
10
11 Mathematica code to produce C code for efficient polynomial detrending
12
13 The polynomials are defined in terms of z = 2i/(n-1)-1, i=0...n
14 such that the range 0...n-1 is mapped to -1...+1
15
16 Thereby the high dynamic range of i^n is reduced to -1..1,
17 and furthermore every second entry in the matrix of
18 the normal equations and its inverse vanishes due to symmetry.
19
20 The regression is done via explicit solution of the normal equations.
21
22 These are of dimension (p+1) for degree p:
23
24 ( s_0 s_1 s_2 .. s_p ) ( a_0 ) = ( b_0 )
25 ( s_1 s_2 s_3 .. s_p+1 ) ( a_1 ) ( b_1 )
26 ( ... ) ( ... ) ( ... )
27 ( s_p s_p+1 .. s_2p ) ( a_p ) ( b_p )
28
29 with
30 s_k = Sum(i=0...n-1) z^k,
31 b_k = Sum(i=0...n-1) x_i z^k,
32 a_k = The unknown coefficient of z^k in the best-fit regression.
33
34 The Matrix with s depends only on p and n, not on the data.
35
36 The Mathematica code below finds exact explicit expressions for
37 the Matrix and its Inverse and writes the corresponding C code.
38
39 If any modifications of the C code are desired apart from "indent"ing,
40 this file below should be altered but not the C code.
41
42 *)
43
44 <<clear.m
45
46 (*
47 <<polyregz.m
48 *)
49
50 z[i_,n_] := 2 i /(n-1) -1;
51
52 s[k_,n_] := Sum[z[i,n]^k,{i,0,n-1}];
53
54 m[n_,p_] := Table[s[i+j,n],{i,0,p},{j,0,p}];
55
56 sol[p_] := sol[p] = Simplify[Inverse[m[n, p]]];
57
58 cln[x_] := CoefficientList[Numerator[x],n];
59 cld[x_] := CoefficientList[Denominator[x],n];
60
61 max[x_] := Max[cln[x],cld[x]];
62
63 pn[x_] := "0" /; x== 0;
64 pn[x_] := ToString[CForm[x]] <> "L";
65
66 ncln[x_] := Map[pn, N[cln[x] / max[x], 22]];
67 ncld[x_] := Map[pn, N[cld[x] / max[x], 22]];
68
69 polyprint[p_, x_] := p[[1]] /; Length[p]==1;
70 polyprint[p_, x_] := "(" <> polyprint[Drop[p,1], x] <> " * " <> x <> " + " <> p[[1]] <> ")" /; Length[p]>1 ;
71
72 cform[x_] := "0" /; FreeQ[x,n];
73 cform[x_] := polyprint[ncln[x], "n"] <> " / " <> polyprint[ncld[x], "n"];
74
75 wlist[sbeg_, list_, sep_, send_] :=
76 Write[fp, sbeg <> Apply[StringJoin, Riffle[list, sep]] <> send];
77
78 out[p_] := Module[ {fn, s, i, j, s1, list},
79 fn="polyreg"<>ToString[p] <> ".c";
80 fp=OpenWrite[fn,FormatType->OutputForm, PageWidth->1024];
81 Write[fp,"void"];
82 Write[fp,"polyreg",p," (double *x, int nn, double *y, double *a)"];
83 Write[fp,"{"];
84 Write[fp,""];
85 Write[fp,"/*"];
86 Write[fp,"polynominal detrending of a time series, order ",p];
87 Write[fp,"machine-generated file, do not edit!"];
88 Write[fp,"made by polyregz.m (Mathematica )"];
89 Write[fp,"Gerhard Heinzel AEI 17.01.2008"];
90 Write[fp,""];
91 Write[fp,"x[]: input, read-only: time series to be detrended"];
92 Write[fp,"nn: input, read-only: length of x[]"];
93 Write[fp,"y[]: output: time series with trend subtracted"];
94 Write[fp,"a[]: fitting coefficients in for z^0, z^1, z^2,... with"];
95 Write[fp," z = 2*i/(nn-1)-1 ; i=0,...,nn-1"];
96 Write[fp,"*/"];
97 Write[fp,""];
98 Write[fp," long double n = nn, n1 = 2.L / (n - 1), z;"];
99 For[i=0, i<=p, i++,
100 Write[fp," long double a",i,", temp",i,", sum",i," = 0;"];
101 ];
102 Write[fp," int i;"];
103 Write[fp," for (i = 0; i < nn; i++)"];
104 Write[fp," {"];
105 Write[fp," z = n1 * i - 1.L;"];
106 list={"x[i]"};
107 For[i=0, i<=p, i++,
108 s1=" sum" <> ToString[i] <> " += ";
109 wlist[s1, list, " * ", ";"];
110 list=Append[list,"z"];
111 ];
112 Write[fp,"/* the above code is efficiently optimized by GCC4 */"];
113 Write[fp,""];
114 Write[fp," }"];
115
116 s=sol[p];
117 For[i=0, i<=p, i++,
118 list={};
119 For[j=0, j<=p, j++,
120 If[!FreeQ[s[[i+1,j+1]],n],
121 Write[fp, "temp", j, " = ", cform[s[[i+1,j+1]]], ";"];
122 list=Append[list,"temp" <> ToString[j] <> " * sum" <> ToString[j]];
123 ]
124 ];
125 s1 = "a" <> ToString[i] <> " = ";
126 wlist[s1, list, " + ", ";"];
127 ];
128 Write[fp," for (i = 0; i < nn; i++)"];
129 Write[fp," {"];
130 Write[fp," z = n1 * i - 1.L;"];
131 list={};
132 For[i=0, i<=p, i++,
133 list=Append[list, "a" <> ToString[i]];
134 ];
135
136 Write[fp," y[i] = x[i] - ", polyprint[list, "z"], ";"];
137 Write[fp," }"];
138 For[i=0, i<=p, i++,
139 Write[fp," a[",i,"] = a",i,";"];
140 ];
141 Write[fp,"}"];
142 Close[fp];
143 ];
144
145 For[i=0, i<=10, i++,
146 out[i];
147 Print["degree ",i," done."];
148 ]
149