tmtf.cpp
Go to the documentation of this file.
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 


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:13