tmtf.cpp
Go to the documentation of this file.
1 
6 
7 
8 //#define WANT_STREAM
9 #define WANT_MATH
10 
11 #include "include.h"
12 
13 #include "newmatap.h"
14 
15 //#include "newmatio.h"
16 
17 #include "tmt.h"
18 
19 #ifdef use_namespace
20 using namespace NEWMAT;
21 #endif
22 
23 
24 
25 static void SlowFT(const ColumnVector& a, const ColumnVector&b,
27 {
28  int n = a.Nrows();
29  x.ReSize(n); y.ReSize(n);
30  Real f = 6.2831853071795864769/n;
31  for (int j=1; j<=n; j++)
32  {
33  Real sumx = 0.0; Real sumy = 0.0;
34  for (int k=1; k<=n; k++)
35  {
36  Real theta = - (j-1) * (k-1) * f;
37  Real c = cos(theta); Real s = sin(theta);
38  sumx += c * a(k) - s * b(k); sumy += s * a(k) + c * b(k);
39  }
40  x(j) = sumx; y(j) = sumy;
41  }
42 }
43 
44 static void SlowDTT_II(const ColumnVector& a, ColumnVector& c, ColumnVector& s)
45 {
46  int n = a.Nrows(); c.ReSize(n); s.ReSize(n);
47  Real f = 6.2831853071795864769 / (4*n);
48  int k;
49 
50  for (k=1; k<=n; k++)
51  {
52  Real sum = 0.0;
53  const int k1 = k-1; // otherwise Visual C++ 5 fails
54  for (int j=1; j<=n; j++) sum += cos(k1 * (2*j-1) * f) * a(j);
55  c(k) = sum;
56  }
57 
58  for (k=1; k<=n; k++)
59  {
60  Real sum = 0.0;
61  for (int j=1; j<=n; j++) sum += sin(k * (2*j-1) * f) * a(j);
62  s(k) = sum;
63  }
64 }
65 
66 static void SlowDTT(const ColumnVector& a, ColumnVector& c, ColumnVector& s)
67 {
68  int n1 = a.Nrows(); int n = n1 - 1;
69  c.ReSize(n1); s.ReSize(n1);
70  Real f = 6.2831853071795864769 / (2*n);
71  int k;
72 
73  int sign = 1;
74  for (k=1; k<=n1; k++)
75  {
76  Real sum = 0.0;
77  for (int j=2; j<=n; j++) sum += cos((j-1) * (k-1) * f) * a(j);
78  c(k) = sum + (a(1) + sign * a(n1)) / 2.0;
79  sign = -sign;
80  }
81 
82  for (k=2; k<=n; k++)
83  {
84  Real sum = 0.0;
85  for (int j=2; j<=n; j++) sum += sin((j-1) * (k-1) * f) * a(j);
86  s(k) = sum;
87  }
88  s(1) = s(n1) = 0;
89 }
90 
91 void SlowFT2(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y)
92 {
93  Tracer trace("SlowFT2");
94  int m = U.Nrows(); int n = U.Ncols();
95  if (m != V.Nrows() || n != V.Ncols() || m == 0 || n == 0)
96  Throw(ProgramException("Matrix dimensions unequal or zero", U, V));
97  X.ReSize(U); Y.ReSize(V);
98  const Real pi2 = atan(1.0) * 8.0;
99  for (int i = 0; i < m; ++i) for (int j = 0; j < n; ++j)
100  {
101  Real sumr = 0.0, sumi = 0.0;
102  for (int k = 0; k < m; ++k) for (int l = 0; l < n; ++l)
103  {
104  Real a = -pi2 * ( (Real)k * (Real)i / (Real)m
105  + (Real)j * (Real)l / (Real)n );
106  Real cs = cos(a); Real sn = sin(a);
107  sumr += cs * U(k+1,l+1) - sn * V(k+1,l+1);
108  sumi += cs * V(k+1,l+1) + sn * U(k+1,l+1);
109  }
110  X(i+1,j+1) = sumr; Y(i+1,j+1) = sumi;
111  }
112 }
113 
114 
115 static void test(int n)
116 {
117  Tracer et("Test FFT");
118  MultWithCarry mwc;
119 
120  ColumnVector A(n), B(n), X, Y;
121  for (int i=0; i<n; i++)
122  { A.element(i) = mwc.Next(); B.element(i) = mwc.Next(); }
123  FFT(A, B, X, Y); FFTI(X, Y, X, Y);
124  X = X - A; Y = Y - B;
125  Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
126 }
127 
128 static void test1(int n)
129 {
130  Tracer et("Test RealFFT");
131  MultWithCarry mwc;
132 
133  ColumnVector A(n), B(n), X, Y;
134  for (int i=0; i<n; i++) A.element(i) = mwc.Next();
135  B = 0.0;
136  FFT(A, B, X, Y);
137  B.CleanUp(); // release some space
138  int n2 = n/2+1;
139  ColumnVector X1,Y1,X2,Y2;
140  RealFFT(A, X1, Y1);
141  X2 = X1 - X.Rows(1,n2); Y2 = Y1 - Y.Rows(1,n2);
142  Clean(X2,0.000000001); Clean(Y2,0.000000001); Print(X2); Print(Y2);
143  X2.CleanUp(); Y2.CleanUp(); // release some more space
144  RealFFTI(X1,Y1,B);
145  B = A - B;
146  Clean(B,0.000000001); Print(B);
147 }
148 
149 static void test2(int n)
150 {
151  Tracer et("cf FFT and slow FT");
152  MultWithCarry mwc;
153 
154  ColumnVector A(n), B(n), X, Y, X1, Y1;
155  for (int i=0; i<n; i++)
156  { A.element(i) = mwc.Next(); B.element(i) = mwc.Next(); }
157  FFT(A, B, X, Y);
158  SlowFT(A, B, X1, Y1);
159  X = X - X1; Y = Y - Y1;
160  Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
161 }
162 
163 static void test3(int n)
164 {
165  Tracer et("cf slow and fast DCT_II and DST_II");
166  MultWithCarry mwc;
167 
168  ColumnVector A(n), X, Y, X1, Y1;
169  for (int i=0; i<n; i++) A.element(i) = mwc.Next();
170  DCT_II(A, X); DST_II(A, Y);
171  SlowDTT_II(A, X1, Y1);
172  X -= X1; Y -= Y1;
173  Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
174 }
175 
176 static void test4(int n)
177 {
178  Tracer et("Test DCT_II");
179  MultWithCarry mwc;
180 
181  ColumnVector A1(n);
182  for (int i=0; i<n; i++) A1.element(i) = mwc.Next();
183  // do DCT II by ordinary FFT
184  ColumnVector P(2*n), Q(2*n);
185  P = 0.0; Q = 0.0; P.Rows(1,n) = A1;
186  FFT(P, Q, P, Q);
187  ColumnVector B1(n);
188  for (int k=0; k<n; k++)
189  {
190  Real arg = k * 6.2831853071795864769 / (4 * n);
191  B1(k+1) = P(k+1) * cos(arg) + Q(k+1) * sin(arg);
192  }
193  // use DCT_II routine
194  ColumnVector B2;
195  DCT_II(A1,B2);
196  // test inverse
197  ColumnVector A2;
198  DCT_II_inverse(B2,A2);
199  A1 -= A2; B1 -= B2;
200  Clean(A1,0.000000001); Clean(B1,0.000000001); Print(A1); Print(B1);
201 }
202 
203 static void test5(int n)
204 {
205  Tracer et("Test DST_II");
206  MultWithCarry mwc;
207 
208  ColumnVector A1(n);
209  for (int i=0; i<n; i++) A1.element(i) = mwc.Next();
210  // do DST II by ordinary FFT
211  ColumnVector P(2*n), Q(2*n);
212  P = 0.0; Q = 0.0; P.Rows(1,n) = A1;
213  FFT(P, Q, P, Q);
214  ColumnVector B1(n);
215  for (int k=1; k<=n; k++)
216  {
217  Real arg = k * 6.2831853071795864769 / (4 * n);
218  B1(k) = P(k+1) * sin(arg) - Q(k+1) * cos(arg);
219  }
220  // use DST_II routine
221  ColumnVector B2;
222  DST_II(A1,B2);
223  // test inverse
224  ColumnVector A2;
225  DST_II_inverse(B2,A2);
226  A1 -= A2; B1 -= B2;
227  Clean(A1,0.000000001); Clean(B1,0.000000001); Print(A1); Print(B1);
228 }
229 
230 static void test6(int n)
231 {
232  Tracer et("Test DST");
233  MultWithCarry mwc;
234 
235  ColumnVector A1(n+1);
236  A1(1) = A1(n+1) = 0;
237  for (int i=1; i<n; i++) A1.element(i) = mwc.Next();
238 
239  // do DST by ordinary FFT
240  ColumnVector P(2*n), Q(2*n); P = 0.0; Q = 0.0; P.Rows(1,n+1) = A1;
241  FFT(P, Q, P, Q);
242  ColumnVector B1 = -Q.Rows(1,n+1);
243  // use DST routine
244  ColumnVector B2;
245  DST(A1,B2);
246  // test inverse
247  ColumnVector A2;
248  DST_inverse(B2,A2);
249  A1 -= A2; B1 -= B2;
250  Clean(A1,0.000000001); Clean(B1,0.000000001); Print(A1); Print(B1);
251 }
252 
253 
254 
255 static void test7(int n)
256 {
257  Tracer et("Test DCT");
258  MultWithCarry mwc;
259 
260  ColumnVector A1(n+1);
261  for (int i=0; i<=n; i++) A1.element(i) = mwc.Next();
262 
263  // do DCT by ordinary FFT
264  ColumnVector P(2*n), Q(2*n); P = 0.0; Q = 0.0; P.Rows(1,n+1) = A1;
265  P(1) /= 2.0; P(n+1) /= 2.0;
266  FFT(P, Q, P, Q);
267  ColumnVector B1 = P.Rows(1,n+1);
268  // use DCT routine
269  ColumnVector B2;
270  DCT(A1,B2);
271  // test inverse
272  ColumnVector A2;
273  DCT_inverse(B2,A2);
274  A1 -= A2; B1 -= B2;
275  Clean(A1,0.000000001); Clean(B1,0.000000001); Print(A1); Print(B1);
276 }
277 
278 static void test8(int n)
279 {
280  Tracer et("cf slow and fast DCT and DST");
281  MultWithCarry mwc;
282 
283  ColumnVector A(n+1), X, Y, X1, Y1;
284  for (int i=0; i<=n; i++) A.element(i) = mwc.Next();
285 
286  DCT(A, X); DST(A, Y);
287  SlowDTT(A, X1, Y1);
288  X -= X1; Y -= Y1;
289  Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
290 }
291 
292 static void test9(int m, int n)
293 {
294  Tracer et("cf FFT2 and slow FT2");
295  MultWithCarry mwc;
296 
297  Matrix A(m,n), B(m,n), X, Y, X1, Y1;
298  for (int i=0; i<m; i++) for (int j=0; j<n; j++)
299  { A.element(i,j) = mwc.Next(); B.element(i,j) = mwc.Next(); }
300  FFT2(A, B, X, Y);
301  SlowFT2(A, B, X1, Y1);
302  X = X - X1; Y = Y - Y1;
303  Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
304  FFT2I(X1, Y1, X1, Y1);
305  X1 -= A; Y1 -= B;
306  Clean(X1,0.000000001); Clean(Y1,0.000000001); Print(X1); Print(Y1);
307 }
308 
309 
310 
311 
312 void trymatf()
313 {
314  Tracer et("Fifteenth test of Matrix package");
316 
317  int i;
318  ColumnVector A(12), B(12);
319  for (i = 1; i <=12; i++)
320  {
321  Real i1 = i - 1;
322  A(i) = .7
323  + .2 * cos(6.2831853071795864769 * 4.0 * i1 / 12)
324  + .3 * sin(6.2831853071795864769 * 3.0 * i1 / 12);
325  B(i) = .9
326  + .5 * sin(6.2831853071795864769 * 2.0 * i1 / 12)
327  + .4 * cos(6.2831853071795864769 * 1.0 * i1 / 12);
328  }
329  FFT(A, B, A, B);
330  A(1) -= 8.4; A(3) -= 3.0; A(5) -= 1.2; A(9) -= 1.2; A(11) += 3.0;
331  B(1) -= 10.8; B(2) -= 2.4; B(4) += 1.8; B(10) -= 1.8; B(12) -= 2.4;
332  Clean(A,0.000000001); Clean(B,0.000000001); Print(A); Print(B);
333 
334 
335  // test FFT
336  test(2048); test(2000); test(27*81); test(2310); test(49*49);
337  test(13*13*13); test(43*47);
338  test(16*16*3); test(16*16*5); test(16*16*7);
339  test(8*8*5);
340 
341  // test realFFT
342  test1(2); test1(98); test1(22); test1(100);
343  test1(2048); test1(2000); test1(35*35*2);
344 
345  // compare FFT and slowFFT
346  test2(1); test2(13); test2(12); test2(9); test2(16); test2(30); test2(42);
347  test2(24); test2(8); test2(40); test2(48); test2(4); test2(3); test2(5);
348  test2(32); test2(2);
349 
350  // compare DCT_II, DST_II and slow versions
351  test3(2); test3(26); test3(32); test3(18);
352 
353  // test DCT_II and DST_II
354  test4(2); test5(2);
355  test4(202); test5(202);
356  test4(1000); test5(1000);
357 
358  // test DST and DCT
359  test6(2); test7(2);
360  test6(274); test7(274);
361  test6(1000); test7(1000);
362 
363  // compare DCT, DST and slow versions
364  test8(2); test8(26); test8(32); test8(18);
365 
366  // compare FFT2 with slow version
367  test9(1,16); test9(16,1);
368  test9(4,3); test9(4,4); test9(4,5); test9(5,3);
369 
370 
371  // now do the same thing all over again forcing use of old FFT
373 
374  for (i = 1; i <=12; i++)
375  {
376  Real i1 = i - 1;
377  A(i) = .7
378  + .2 * cos(6.2831853071795864769 * 4.0 * i1 / 12)
379  + .3 * sin(6.2831853071795864769 * 3.0 * i1 / 12);
380  B(i) = .9
381  + .5 * sin(6.2831853071795864769 * 2.0 * i1 / 12)
382  + .4 * cos(6.2831853071795864769 * 1.0 * i1 / 12);
383  }
384  FFT(A, B, A, B);
385  A(1) -= 8.4; A(3) -= 3.0; A(5) -= 1.2; A(9) -= 1.2; A(11) += 3.0;
386  B(1) -= 10.8; B(2) -= 2.4; B(4) += 1.8; B(10) -= 1.8; B(12) -= 2.4;
387  Clean(A,0.000000001); Clean(B,0.000000001); Print(A); Print(B);
388 
389 
390  // test FFT
391  test(2048); test(2000); test(27*81); test(2310); test(49*49);
392  test(13*13*13); test(43*47);
393  test(16*16*3); test(16*16*5); test(16*16*7);
394  test(8*8*5);
395 
396  // test realFFT
397  test1(2); test1(98); test1(22); test1(100);
398  test1(2048); test1(2000); test1(35*35*2);
399 
400  // compare FFT and slowFFT
401  test2(1); test2(13); test2(12); test2(9); test2(16); test2(30); test2(42);
402  test2(24); test2(8); test2(40); test2(48); test2(4); test2(3); test2(5);
403  test2(32); test2(2);
404 
405  // compare DCT_II, DST_II and slow versions
406  test3(2); test3(26); test3(32); test3(18);
407 
408  // test DCT_II and DST_II
409  test4(2); test5(2);
410  test4(202); test5(202);
411  test4(1000); test5(1000);
412 
413  // test DST and DCT
414  test6(2); test7(2);
415  test6(274); test7(274);
416  test6(1000); test7(1000);
417 
418  // compare DCT, DST and slow versions
419  test8(2); test8(26); test8(32); test8(18);
420 
421  // compare FFT2 with slow version
422  // don't redo these
423 
425 
426 
427 
428 }
429 
void trymatf()
Definition: tmtf.cpp:312
void DST_inverse(const ColumnVector &V, ColumnVector &U)
Definition: fft.cpp:406
void DST_II(const ColumnVector &U, ColumnVector &V)
Definition: fft.cpp:308
static void test5(int n)
Definition: tmtf.cpp:203
void RealFFT(const ColumnVector &U, ColumnVector &X, ColumnVector &Y)
Definition: fft.cpp:130
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
void ReSize(int m)
Definition: newmat.h:1034
void SlowFT2(const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
Definition: tmtf.cpp:91
static void SlowFT(const ColumnVector &a, const ColumnVector &b, ColumnVector &x, ColumnVector &y)
Definition: tmtf.cpp:25
static void test4(int n)
Definition: tmtf.cpp:176
static void test(int n)
Definition: tmtf.cpp:115
virtual void ReSize(int m, int n)
Definition: newmat.h:662
static void test7(int n)
Definition: tmtf.cpp:255
double Real
Definition: include.h:307
Real & element(int, int)
Definition: newmat6.cpp:682
int Nrows() const
Definition: newmat.h:494
void DCT_II(const ColumnVector &U, ColumnVector &V)
Definition: fft.cpp:253
void FFT(const ColumnVector &U, const ColumnVector &V, ColumnVector &X, ColumnVector &Y)
Definition: fft.cpp:201
Real Next()
Definition: tmt.cpp:187
void FFT2(const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
Definition: fft.cpp:444
void DCT_II_inverse(const ColumnVector &V, ColumnVector &U)
Definition: fft.cpp:281
void FFT2I(const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
Definition: fft.cpp:465
void Clean(Matrix &A, Real c)
Definition: tmt.cpp:128
void FFTI(const ColumnVector &U, const ColumnVector &V, ColumnVector &X, ColumnVector &Y)
Definition: fft.cpp:120
void DST(const ColumnVector &U, ColumnVector &V)
Definition: fft.cpp:434
static void test2(int n)
Definition: tmtf.cpp:149
static void test8(int n)
Definition: tmtf.cpp:278
#define Throw(E)
Definition: myexcept.h:191
static void PrintTrace()
Definition: myexcept.cpp:109
The usual rectangular matrix.
Definition: newmat.h:625
static void test9(int m, int n)
Definition: tmtf.cpp:292
Real & element(int)
Definition: newmat6.cpp:778
FloatVector FloatVector * a
int Ncols() const
Definition: newmat.h:495
Real trace(const BaseMatrix &B)
Definition: newmat.h:2099
void CleanUp()
Definition: newmat.h:396
void RealFFTI(const ColumnVector &A, const ColumnVector &B, ColumnVector &U)
Definition: fft.cpp:166
static void SlowDTT_II(const ColumnVector &a, ColumnVector &c, ColumnVector &s)
Definition: tmtf.cpp:44
static bool OnlyOldFFT
Definition: newmatap.h:145
Real sum(const BaseMatrix &B)
Definition: newmat.h:2105
static void test3(int n)
Definition: tmtf.cpp:163
void DCT(const ColumnVector &U, ColumnVector &V)
Definition: fft.cpp:397
void DCT_inverse(const ColumnVector &V, ColumnVector &U)
Definition: fft.cpp:363
void Print(const Matrix &X)
Definition: tmt.cpp:42
Column vector.
Definition: newmat.h:1008
static void SlowDTT(const ColumnVector &a, ColumnVector &c, ColumnVector &s)
Definition: tmtf.cpp:66
static void test1(int n)
Definition: tmtf.cpp:128
void DST_II_inverse(const ColumnVector &V, ColumnVector &U)
Definition: fft.cpp:336
static void test6(int n)
Definition: tmtf.cpp:230
GetSubMatrix Rows(int f, int l) const
Definition: newmat.h:2151


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:17