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