00001
00002
00003 #define WANT_MATH
00004
00005 #include "include.h"
00006
00007 #include "newmatap.h"
00008
00009
00010
00011 #include "tmt.h"
00012
00013 #ifdef use_namespace
00014 using namespace NEWMAT;
00015 #endif
00016
00017
00018
00019 static void SlowFT(const ColumnVector& a, const ColumnVector&b,
00020 ColumnVector& x, ColumnVector& y)
00021 {
00022 int n = a.Nrows();
00023 x.ReSize(n); y.ReSize(n);
00024 Real f = 6.2831853071795864769/n;
00025 for (int j=1; j<=n; j++)
00026 {
00027 Real sumx = 0.0; Real sumy = 0.0;
00028 for (int k=1; k<=n; k++)
00029 {
00030 Real theta = - (j-1) * (k-1) * f;
00031 Real c = cos(theta); Real s = sin(theta);
00032 sumx += c * a(k) - s * b(k); sumy += s * a(k) + c * b(k);
00033 }
00034 x(j) = sumx; y(j) = sumy;
00035 }
00036 }
00037
00038 static void SlowDTT_II(const ColumnVector& a, ColumnVector& c, ColumnVector& s)
00039 {
00040 int n = a.Nrows(); c.ReSize(n); s.ReSize(n);
00041 Real f = 6.2831853071795864769 / (4*n);
00042 int k;
00043
00044 for (k=1; k<=n; k++)
00045 {
00046 Real sum = 0.0;
00047 const int k1 = k-1;
00048 for (int j=1; j<=n; j++) sum += cos(k1 * (2*j-1) * f) * a(j);
00049 c(k) = sum;
00050 }
00051
00052 for (k=1; k<=n; k++)
00053 {
00054 Real sum = 0.0;
00055 for (int j=1; j<=n; j++) sum += sin(k * (2*j-1) * f) * a(j);
00056 s(k) = sum;
00057 }
00058 }
00059
00060 static void SlowDTT(const ColumnVector& a, ColumnVector& c, ColumnVector& s)
00061 {
00062 int n1 = a.Nrows(); int n = n1 - 1;
00063 c.ReSize(n1); s.ReSize(n1);
00064 Real f = 6.2831853071795864769 / (2*n);
00065 int k;
00066
00067 int sign = 1;
00068 for (k=1; k<=n1; k++)
00069 {
00070 Real sum = 0.0;
00071 for (int j=2; j<=n; j++) sum += cos((j-1) * (k-1) * f) * a(j);
00072 c(k) = sum + (a(1) + sign * a(n1)) / 2.0;
00073 sign = -sign;
00074 }
00075
00076 for (k=2; k<=n; k++)
00077 {
00078 Real sum = 0.0;
00079 for (int j=2; j<=n; j++) sum += sin((j-1) * (k-1) * f) * a(j);
00080 s(k) = sum;
00081 }
00082 s(1) = s(n1) = 0;
00083 }
00084
00085 static void test(int n)
00086 {
00087 Tracer et("Test FFT");
00088
00089 ColumnVector A(n), B(n), X, Y;
00090 for (int i=0; i<n; i++)
00091 {
00092 Real f = 6.2831853071795864769*i/n;
00093 A.element(i) = fabs(sin(7.0*f) + 0.5 * cos(19.0 * f)) + (Real)i/(Real)n;
00094 B.element(i) = fabs(0.25 * cos(31.0 * f)) + (Real)i/(Real)n;
00095 }
00096 FFT(A, B, X, Y); FFTI(X, Y, X, Y);
00097 X = X - A; Y = Y - B;
00098 Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
00099 }
00100
00101 static void test1(int n)
00102 {
00103 Tracer et("Test RealFFT");
00104
00105 ColumnVector A(n), B(n), X, Y;
00106 for (int i=0; i<n; i++)
00107 {
00108 Real f = 6.2831853071795864769*i/n;
00109 A.element(i) = fabs(sin(7.0*f) + 0.5 * cos(19.0 * f)) + (Real)i/(Real)n;
00110 }
00111 B = 0.0;
00112 FFT(A, B, X, Y);
00113 B.CleanUp();
00114 int n2 = n/2+1;
00115 ColumnVector X1,Y1,X2,Y2;
00116 RealFFT(A, X1, Y1);
00117 X2 = X1 - X.Rows(1,n2); Y2 = Y1 - Y.Rows(1,n2);
00118 Clean(X2,0.000000001); Clean(Y2,0.000000001); Print(X2); Print(Y2);
00119 X2.CleanUp(); Y2.CleanUp();
00120 RealFFTI(X1,Y1,B);
00121 B = A - B;
00122 Clean(B,0.000000001); Print(B);
00123 }
00124
00125 static void test2(int n)
00126 {
00127 Tracer et("cf FFT and slow FT");
00128
00129 ColumnVector A(n), B(n), X, Y, X1, Y1;
00130 for (int i=0; i<n; i++)
00131 {
00132 Real f = 6.2831853071795864769*i/n;
00133 A.element(i) = fabs(sin(7.0*f) - 0.5 * cos(19.0 * f)) + (Real)i/(Real)n;
00134 B.element(i) = fabs(0.25 * cos(31.0 * f)) - (Real)i/(Real)n;
00135 }
00136 FFT(A, B, X, Y);
00137 SlowFT(A, B, X1, Y1);
00138 X = X - X1; Y = Y - Y1;
00139 Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
00140 }
00141
00142 static void test3(int n)
00143 {
00144 Tracer et("cf slow and fast DCT_II and DST_II");
00145
00146 ColumnVector A(n), X, Y, X1, Y1;
00147 for (int i=0; i<n; i++)
00148 {
00149 Real f = 6.2831853071795864769*i/n;
00150 A.element(i) = fabs(sin(7.0*f) - 0.55 * cos(19.0 * f)
00151 + .73 * sin(6.0 * f)) + (Real)i/(Real)n;
00152 }
00153 DCT_II(A, X); DST_II(A, Y);
00154 SlowDTT_II(A, X1, Y1);
00155 X -= X1; Y -= Y1;
00156 Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
00157 }
00158
00159 static void test4(int n)
00160 {
00161 Tracer et("Test DCT_II");
00162
00163 ColumnVector A1(n);
00164 for (int i=0; i<n; i++)
00165 {
00166 Real f = 6.2831853071795864769*i/n;
00167 A1.element(i) =
00168 fabs(sin(7.0*f) + 0.7588 * cos(19.0 * f) + (Real)i/(Real)n);
00169 }
00170
00171 ColumnVector P(2*n), Q(2*n);
00172 P = 0.0; Q = 0.0; P.Rows(1,n) = A1;
00173 FFT(P, Q, P, Q);
00174 ColumnVector B1(n);
00175 for (int k=0; k<n; k++)
00176 {
00177 Real arg = k * 6.2831853071795864769 / (4 * n);
00178 B1(k+1) = P(k+1) * cos(arg) + Q(k+1) * sin(arg);
00179 }
00180
00181 ColumnVector B2;
00182 DCT_II(A1,B2);
00183
00184 ColumnVector A2;
00185 DCT_II_inverse(B2,A2);
00186 A1 -= A2; B1 -= B2;
00187 Clean(A1,0.000000001); Clean(B1,0.000000001); Print(A1); Print(B1);
00188 }
00189
00190 static void test5(int n)
00191 {
00192 Tracer et("Test DST_II");
00193
00194 ColumnVector A1(n);
00195 for (int i=0; i<n; i++)
00196 {
00197 Real f = 6.2831853071795864769*i/n;
00198 A1.element(i) =
00199 fabs(sin(11.0*f) + 0.7588 * cos(19.0 * f) + (Real)i/(Real)n);
00200 }
00201
00202 ColumnVector P(2*n), Q(2*n);
00203 P = 0.0; Q = 0.0; P.Rows(1,n) = A1;
00204 FFT(P, Q, P, Q);
00205 ColumnVector B1(n);
00206 for (int k=1; k<=n; k++)
00207 {
00208 Real arg = k * 6.2831853071795864769 / (4 * n);
00209 B1(k) = P(k+1) * sin(arg) - Q(k+1) * cos(arg);
00210 }
00211
00212 ColumnVector B2;
00213 DST_II(A1,B2);
00214
00215 ColumnVector A2;
00216 DST_II_inverse(B2,A2);
00217 A1 -= A2; B1 -= B2;
00218 Clean(A1,0.000000001); Clean(B1,0.000000001); Print(A1); Print(B1);
00219 }
00220
00221 static void test6(int n)
00222 {
00223 Tracer et("Test DST");
00224
00225 ColumnVector A1(n+1);
00226 A1(1) = A1(n+1) = 0;
00227 for (int i=1; i<n; i++)
00228 {
00229 Real f = 6.2831853071795864769*i/n;
00230 A1.element(i) =
00231 fabs(sin(11.0*f) + 0.7588 * cos(19.0 * f) + (Real)i/(Real)n);
00232 }
00233
00234 ColumnVector P(2*n), Q(2*n); P = 0.0; Q = 0.0; P.Rows(1,n+1) = A1;
00235 FFT(P, Q, P, Q);
00236 ColumnVector B1 = -Q.Rows(1,n+1);
00237
00238 ColumnVector B2;
00239 DST(A1,B2);
00240
00241 ColumnVector A2;
00242 DST_inverse(B2,A2);
00243 A1 -= A2; B1 -= B2;
00244 Clean(A1,0.000000001); Clean(B1,0.000000001); Print(A1); Print(B1);
00245 }
00246
00247
00248
00249 static void test7(int n)
00250 {
00251 Tracer et("Test DCT");
00252
00253 ColumnVector A1(n+1);
00254 for (int i=0; i<=n; i++)
00255 {
00256 Real f = 6.2831853071795864769*i/n;
00257 A1.element(i) =
00258 fabs(sin(17.0*f) + 0.6399 * cos(23.0 * f) + 1.32*(Real)i/(Real)n);
00259 }
00260
00261 ColumnVector P(2*n), Q(2*n); P = 0.0; Q = 0.0; P.Rows(1,n+1) = A1;
00262 P(1) /= 2.0; P(n+1) /= 2.0;
00263 FFT(P, Q, P, Q);
00264 ColumnVector B1 = P.Rows(1,n+1);
00265
00266 ColumnVector B2;
00267 DCT(A1,B2);
00268
00269 ColumnVector A2;
00270 DCT_inverse(B2,A2);
00271 A1 -= A2; B1 -= B2;
00272 Clean(A1,0.000000001); Clean(B1,0.000000001); Print(A1); Print(B1);
00273 }
00274
00275 static void test8(int n)
00276 {
00277 Tracer et("cf slow and fast DCT and DST");
00278
00279 ColumnVector A(n+1), X, Y, X1, Y1;
00280 for (int i=0; i<=n; i++)
00281 {
00282 Real f = 6.2831853071795864769*i/n;
00283 A.element(i) = fabs(sin(7.0*f) - 0.5 * cos(19.0 * f) +
00284 0.3 * (Real)i/(Real)n);
00285 }
00286 DCT(A, X); DST(A, Y);
00287 SlowDTT(A, X1, Y1);
00288 X -= X1; Y -= Y1;
00289 Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
00290 }
00291
00292
00293
00294
00295 void trymatf()
00296 {
00297 Tracer et("Fifteenth test of Matrix package");
00298 Tracer::PrintTrace();
00299
00300 int i;
00301 ColumnVector A(12), B(12);
00302 for (i = 1; i <=12; i++)
00303 {
00304 Real i1 = i - 1;
00305 A(i) = .7
00306 + .2 * cos(6.2831853071795864769 * 4.0 * i1 / 12)
00307 + .3 * sin(6.2831853071795864769 * 3.0 * i1 / 12);
00308 B(i) = .9
00309 + .5 * sin(6.2831853071795864769 * 2.0 * i1 / 12)
00310 + .4 * cos(6.2831853071795864769 * 1.0 * i1 / 12);
00311 }
00312 FFT(A, B, A, B);
00313 A(1) -= 8.4; A(3) -= 3.0; A(5) -= 1.2; A(9) -= 1.2; A(11) += 3.0;
00314 B(1) -= 10.8; B(2) -= 2.4; B(4) += 1.8; B(10) -= 1.8; B(12) -= 2.4;
00315 Clean(A,0.000000001); Clean(B,0.000000001); Print(A); Print(B);
00316
00317
00318
00319 test(2048); test(2000); test(27*81); test(2310); test(49*49);
00320 test(13*13*13); test(43*47);
00321 test(16*16*3); test(16*16*5); test(16*16*7);
00322 test(8*8*5);
00323
00324
00325 test1(2); test1(98); test1(22); test1(100);
00326 test1(2048); test1(2000); test1(35*35*2);
00327
00328
00329 test2(1); test2(13); test2(12); test2(9); test2(16); test2(30); test2(42);
00330 test2(24); test2(8); test2(40); test2(48); test2(4); test2(3); test2(5);
00331 test2(32); test2(2);
00332
00333
00334 test3(2); test3(26); test3(32); test3(18);
00335
00336
00337 test4(2); test5(2);
00338 test4(202); test5(202);
00339 test4(1000); test5(1000);
00340
00341
00342 test6(2); test7(2);
00343 test6(274); test7(274);
00344 test6(1000); test7(1000);
00345
00346
00347 test8(2); test8(26); test8(32); test8(18);
00348
00349
00350
00351 FFT_Controller::OnlyOldFFT = true;
00352
00353 for (i = 1; i <=12; i++)
00354 {
00355 Real i1 = i - 1;
00356 A(i) = .7
00357 + .2 * cos(6.2831853071795864769 * 4.0 * i1 / 12)
00358 + .3 * sin(6.2831853071795864769 * 3.0 * i1 / 12);
00359 B(i) = .9
00360 + .5 * sin(6.2831853071795864769 * 2.0 * i1 / 12)
00361 + .4 * cos(6.2831853071795864769 * 1.0 * i1 / 12);
00362 }
00363 FFT(A, B, A, B);
00364 A(1) -= 8.4; A(3) -= 3.0; A(5) -= 1.2; A(9) -= 1.2; A(11) += 3.0;
00365 B(1) -= 10.8; B(2) -= 2.4; B(4) += 1.8; B(10) -= 1.8; B(12) -= 2.4;
00366 Clean(A,0.000000001); Clean(B,0.000000001); Print(A); Print(B);
00367
00368
00369
00370 test(2048); test(2000); test(27*81); test(2310); test(49*49);
00371 test(13*13*13); test(43*47);
00372 test(16*16*3); test(16*16*5); test(16*16*7);
00373 test(8*8*5);
00374
00375
00376 test1(2); test1(98); test1(22); test1(100);
00377 test1(2048); test1(2000); test1(35*35*2);
00378
00379
00380 test2(1); test2(13); test2(12); test2(9); test2(16); test2(30); test2(42);
00381 test2(24); test2(8); test2(40); test2(48); test2(4); test2(3); test2(5);
00382 test2(32); test2(2);
00383
00384
00385 test3(2); test3(26); test3(32); test3(18);
00386
00387
00388 test4(2); test5(2);
00389 test4(202); test5(202);
00390 test4(1000); test5(1000);
00391
00392
00393 test6(2); test7(2);
00394 test6(274); test7(274);
00395 test6(1000); test7(1000);
00396
00397
00398 test8(2); test8(26); test8(32); test8(18);
00399
00400 FFT_Controller::OnlyOldFFT = false;
00401
00402
00403
00404 }