$search
00001 00002 00003 00006 00007 00008 //#define WANT_STREAM 00009 #define WANT_MATH 00010 00011 #include "include.h" 00012 00013 #include "newmatap.h" 00014 00015 //#include "newmatio.h" 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; // otherwise Visual C++ 5 fails 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(); // release some space 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(); // release some more space 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 // do DCT II by ordinary FFT 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 // use DCT_II routine 00194 ColumnVector B2; 00195 DCT_II(A1,B2); 00196 // test inverse 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 // do DST II by ordinary FFT 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 // use DST_II routine 00221 ColumnVector B2; 00222 DST_II(A1,B2); 00223 // test inverse 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 // do DST by ordinary FFT 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 // use DST routine 00244 ColumnVector B2; 00245 DST(A1,B2); 00246 // test inverse 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 // do DCT by ordinary FFT 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 // use DCT routine 00269 ColumnVector B2; 00270 DCT(A1,B2); 00271 // test inverse 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 // test FFT 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 // test realFFT 00342 test1(2); test1(98); test1(22); test1(100); 00343 test1(2048); test1(2000); test1(35*35*2); 00344 00345 // compare FFT and slowFFT 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 // compare DCT_II, DST_II and slow versions 00351 test3(2); test3(26); test3(32); test3(18); 00352 00353 // test DCT_II and DST_II 00354 test4(2); test5(2); 00355 test4(202); test5(202); 00356 test4(1000); test5(1000); 00357 00358 // test DST and DCT 00359 test6(2); test7(2); 00360 test6(274); test7(274); 00361 test6(1000); test7(1000); 00362 00363 // compare DCT, DST and slow versions 00364 test8(2); test8(26); test8(32); test8(18); 00365 00366 // compare FFT2 with slow version 00367 test9(1,16); test9(16,1); 00368 test9(4,3); test9(4,4); test9(4,5); test9(5,3); 00369 00370 00371 // now do the same thing all over again forcing use of old FFT 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 // test FFT 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 // test realFFT 00397 test1(2); test1(98); test1(22); test1(100); 00398 test1(2048); test1(2000); test1(35*35*2); 00399 00400 // compare FFT and slowFFT 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 // compare DCT_II, DST_II and slow versions 00406 test3(2); test3(26); test3(32); test3(18); 00407 00408 // test DCT_II and DST_II 00409 test4(2); test5(2); 00410 test4(202); test5(202); 00411 test4(1000); test5(1000); 00412 00413 // test DST and DCT 00414 test6(2); test7(2); 00415 test6(274); test7(274); 00416 test6(1000); test7(1000); 00417 00418 // compare DCT, DST and slow versions 00419 test8(2); test8(26); test8(32); test8(18); 00420 00421 // compare FFT2 with slow version 00422 // don't redo these 00423 00424 FFT_Controller::OnlyOldFFT = false; 00425 00426 00427 00428 } 00429