home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / fft.lha / fft / vfft.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-08-09  |  12.4 KB  |  472 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *        Verify the Fast Fourier Transform Package
  6.  *
  7.  ************************************************************************
  8.  */
  9.  
  10. #include "fft.h"
  11. #include <math.h>
  12. #include <builtin.h>
  13. #include <ostream.h>
  14.  
  15. /*
  16.  *------------------------------------------------------------------------
  17.  *             Timing the program execution
  18.  */
  19.  
  20. #include <time.h>
  21.  
  22. static tms clock_acc;
  23.  
  24. static void start_timing(void)
  25. {
  26.   times(&clock_acc);
  27. }
  28.  
  29. static void print_timing(const char * header)
  30. {
  31.   register int old_tick = clock_acc.tms_utime;
  32.   register double timing = (times(&clock_acc),clock_acc.tms_utime
  33.                 - old_tick )/60.;     // In secs  
  34.   printf("\nIt took %.2f sec to perform %s\n",timing,header);
  35. }
  36.  
  37. /*
  38.  *-----------------------------------------------------------------------
  39.  */
  40.  
  41.                 // Simplified printing a vector
  42. static void print_seq(const char * header, const Vector& v)        
  43. {
  44.   register int i;
  45.   cout << "\n" << header << "\t";
  46.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  47.     cout << form("%7.4f ",v(i));
  48.   cout << "\n";
  49. }
  50.  
  51.  
  52. static void verify_vector_identity(const Vector& v1, const Vector& v2)
  53. {
  54.   register imax = 0;
  55.   register double max_dev = 0;
  56.   register int i;
  57.   are_compatible(v1,v2);
  58.   for(i=v1.q_lwb(); i<=v1.q_upb(); i++)
  59.   {
  60.     register double dev = abs(v1(i)-v2(i));
  61.     if( dev >= max_dev )
  62.       imax = i, max_dev = dev;
  63.   }
  64.  
  65.   if( max_dev == 0 )
  66.     return;
  67.   if( max_dev < 1e-5 )
  68.     message("Two #%d elements of the vectors with values %g and %g\n"
  69.         "differ the most, though the deviation %g is small\n",
  70.         imax,v1(imax),v2(imax),max_dev);
  71.   else
  72.     _error("A significant difference between the vectors encountered\n"
  73.        "at #%d element, with values %g and %g",
  74.        imax,v1(imax),v2(imax));
  75. }
  76.  
  77. /*
  78.  *-----------------------------------------------------------------------
  79.  *          Check FFT of the Arithmetical Progression sequence
  80.  *                x[j] = j
  81.  * Analytical transform is
  82.  *    SUM{ j*W^(kj) } = N/(W^k - 1), k > 0,
  83.  *                  N*(N-1)/2,   k = 0
  84.  */
  85.  
  86. void test_ap_series(const int N)
  87. {
  88.   cout << "\n\nVerify the computed FFT of the AP series x[j]=j\n";
  89.   cout << "j = 0.." << N-1 << "\n";
  90.  
  91.   Vector xre(0,N-1);
  92.   Vector xim(xre);
  93.  
  94.   register int j;
  95.   for(j=0; j<N; j++)
  96.     xre(j) = j;
  97.   xim = 0;
  98.  
  99.  
  100.   Vector xfe_re(xre), xfe_im(xre), xfe_abs(xre); // Exact transform
  101.   xfe_re(0) = xfe_abs(0) = N*(N-1)/2;
  102.   xfe_im(0) = 0;
  103.   register int k;
  104.   for(k=1; k<N; k++)
  105.   {
  106.     Complex arg(0,-2*PI/N * k);
  107.     Complex t = N / ( exp(arg) - 1 );
  108.     xfe_re(k)  = t.real();
  109.     xfe_im(k)  = t.imag();
  110.     xfe_abs(k) = abs(t);
  111.   }
  112.  
  113.   FFT fft(N);
  114.   Vector xf_re(xre), xf_im(xre), xf_abs(xre);
  115.  
  116.   cout << "\nPerforming Complex FFT of AP series (IM part being set to 0)\n";
  117.   start_timing();
  118.   fft.input(xre,xim);
  119.   fft.real(xf_re);
  120.   fft.imag(xf_im);
  121.   print_timing("Complex Fourier transform");
  122.   fft.abs(xf_abs);
  123.  
  124.   cout << "Verifying the Re part of the transform ...\n";
  125.   verify_vector_identity(xfe_re,xf_re);
  126.   cout << "Verifying the Im part of the transform ...\n";
  127.   verify_vector_identity(xfe_im,xf_im);
  128.   cout << "Verifying the power spectrum ...\n";
  129.   verify_vector_identity(xfe_abs,xf_abs);
  130.  
  131.   Vector xfr_re(xre), xfr_im(xre), xfr_abs(xre);
  132.   cout << "\nPerforming FFT of a REAL AP sequence\n";
  133.   start_timing();
  134.   fft.input(xre);
  135.   fft.real(xfr_re);
  136.   fft.imag(xfr_im);
  137.   print_timing("\"Real\" Fourier transform");
  138.   fft.abs(xfr_abs);
  139.  
  140.   cout << "Check out that \"Real\" and Complex FFT give identical results";
  141.   verify_vector_identity(xfr_re,xf_re);
  142.   verify_vector_identity(xfr_im,xf_im);
  143.   verify_vector_identity(xfr_abs,xf_abs);
  144.  
  145.   cout << "\nDone\n";
  146. }
  147.  
  148. /*
  149.  *-----------------------------------------------------------------------
  150.  *          Check out the orthogonality of the basis functions of FFT
  151.  *            
  152.  * x[j] = W^(-l*j)
  153.  * SUM{ x[j] * W^(kj) } = 0,   k <> l
  154.  *                  N,   k = l
  155.  */
  156.  
  157. void test_orth(const int N, const int l)
  158. {
  159.   cout << "\n\nVerify the computed FFT for x[j] = W^(-l*j)\n";
  160.   cout << "j = 0.." << N-1 << ", l=" << l << "\n";
  161.  
  162.   Vector xre(0,N-1);
  163.   Vector xim(xre);
  164.  
  165.   register int j;
  166.   for(j=0; j<N; j++)
  167.   {
  168.     Complex arg(0, 2*PI/N * l * j);
  169.     Complex t = exp(arg);
  170.     xre(j) = t.real(), xim(j) = t.imag();
  171.   }
  172.  
  173.   Vector xfe_re(xre), xfe_im(xre);        // Exact transform
  174.   xfe_re(l) = N; xfe_im = 0;
  175.  
  176.   FFT fft(N);
  177.   Vector xf_re(xre), xf_im(xre), xf_abs(xre);
  178.  
  179.   cout << "\nPerforming Complex FFT\n";
  180.   start_timing();
  181.   fft.input(xre,xim);
  182.   fft.real(xf_re);
  183.   fft.imag(xf_im);
  184.   print_timing("Complex Fourier transform");
  185.   fft.abs(xf_abs);
  186.  
  187.   cout << "Verifying the Re part of the transform ...\n";
  188.   verify_vector_identity(xfe_re,xf_re);
  189.   cout << "Verifying the Im part of the transform ...\n";
  190.   verify_vector_identity(xfe_im,xf_im);
  191.   cout << "Verifying the power spectrum ...\n";
  192.   verify_vector_identity(xfe_re,xf_abs);
  193.  
  194.   cout << "\nDone\n";
  195. }
  196.  
  197.  
  198. /*
  199.  *-----------------------------------------------------------------------
  200.  *      Check FFT of the truncated Arithmetical Progression sequence
  201.  *                 x[j] = j, j=0..N/2
  202.  * Analytical transform is
  203.  *    SUM{ j*W^(kj) } = N/2 (W^k - 1), k > 0 and even
  204.  *              2*W^k/(W^k - 1)^2 - N/2 * 1/(W^k - 1), k being odd
  205.  *                  N/2 * (N/2-1)/2,   k = 0
  206.  */
  207.  
  208. void test_ap_series_pad0(const int N)
  209. {
  210.   cout << "\n\nVerify the computed FFT of the truncated AP sequence x[j]=j\n";
  211.   cout << "j = 0.." << N/2-1 << ", with N=" << N << "\n";
  212.  
  213.   Vector xre(0,N/2-1);
  214.   Vector xim(xre);
  215.  
  216.   register int j;
  217.   for(j=0; j<N/2; j++)
  218.     xre(j) = j;
  219.   xim = 0;
  220.  
  221.  
  222.   Vector xfe_re(0,N-1), xfe_im(0,N-1), xfe_abs(0,N-1); // Exact transform
  223.   xfe_re(0) = xfe_abs(0) = N/2 * (N/2 - 1)/2;
  224.   xfe_im(0) = 0;
  225.   register int k;
  226.   for(k=1; k<N; k++)
  227.   {
  228.     Complex wk = exp(Complex(0,-2*PI/N * k));
  229.     Complex t = 1 / ( wk - 1 );
  230.     if( k & 1 )
  231.       t = 2*wk*t*t - N/2 * t;            // for k odd
  232.     else
  233.       t *= N/2;                    // for k even
  234.     xfe_re(k)  = t.real();
  235.     xfe_im(k)  = t.imag();
  236.     xfe_abs(k) = abs(t);
  237.   }
  238.  
  239.   FFT fft(N);
  240.   Vector xf_re(0,N-1), xf_im(0,N-1), xf_abs(0,N-1);
  241.  
  242.   cout << "\nPerforming Complex FFT (with IM part being set to 0)\n";
  243.   start_timing();
  244.   fft.input_pad0(xre,xim);
  245.   fft.real(xf_re);
  246.   fft.imag(xf_im);
  247.   print_timing("Complex Fourier transform");
  248.   fft.abs(xf_abs);
  249.  
  250.   if( N <= 16 )
  251.   {
  252.     print_seq("Source Vector         ",xre);
  253.     print_seq("Computed cos transform",xf_re);
  254.     print_seq("Computed sin transform",xf_im);
  255.   }
  256.  
  257.   cout << "Verifying the Re part of the transform ...\n";
  258.   verify_vector_identity(xfe_re,xf_re);
  259.   cout << "Verifying the Im part of the transform ...\n";
  260.   verify_vector_identity(xfe_im,xf_im);
  261.   cout << "Verifying the power spectrum ...\n";
  262.   verify_vector_identity(xfe_abs,xf_abs);
  263.  
  264.   Vector xfr_re(xf_re), xfr_im(xf_re), xfr_abs(xf_re);
  265.   cout << "\nPerforming FFT of a REAL AP sequence\n";
  266.   start_timing();
  267.   fft.input_pad0(xre);
  268.   fft.real(xfr_re);
  269.   fft.imag(xfr_im);
  270.   print_timing("\"Real\" Fourier transform");
  271.   fft.abs(xfr_abs);
  272.  
  273.   cout << "Check out that \"Real\" and Complex FFT give identical results\n";
  274.   verify_vector_identity(xfr_re,xf_re);
  275.   verify_vector_identity(xfr_im,xf_im);
  276.   verify_vector_identity(xfr_abs,xf_abs);
  277.  
  278.   Vector xfh_re(xre), xfh_im(xre), xfh_abs(xre);
  279.   cout << "Check out the functions returning the half of the transform\n";
  280.   fft.real_half(xfh_re);
  281.   fft.imag_half(xfh_im);
  282.   fft.abs_half(xfh_abs);
  283.   Vector xfhr_re(xre), xfhr_im(xre), xfhr_abs(xre);
  284.   for(j=0; j<N/2; j++)
  285.     xfhr_re(j) = xfh_re(j), xfhr_im(j) = xfh_im(j), xfhr_abs(j) = xfh_abs(j);
  286.   verify_vector_identity(xfhr_re,xfh_re);
  287.   verify_vector_identity(xfhr_im,xfh_im);
  288.   verify_vector_identity(xfhr_abs,xfh_abs);
  289.  
  290.   cout << "\nDone\n";
  291. }
  292.  
  293.  
  294. /*
  295.  *-----------------------------------------------------------------------
  296.  *             Verify the sin/cos Fourier transform
  297.  *
  298.  *    r*exp( -r/a )    <=== sin-transform ===> 2a^3 k / (1 + (ak)^2)^2
  299.  *    r*exp( -r/a )    <=== cos-transform ===> a^2 (1 - (ak)^2)/(1 + (ak)^2)^2
  300.  */
  301.  
  302. static void test_lorentzian(void)
  303. {
  304.   const double R = 20;            // Cutoff distance
  305.   const int n = 512;            // No. of grids
  306.   const double a = 4;            // Constant a, see above
  307.   const double dr = R/n,         // Grid meshes
  308.            dk = PI/R;    
  309.  
  310.   cout << "\n\nVerify the sin/cos transform for the following example\n";
  311.   cout << "\tr*exp( -r/a )\t<=== sin-transform ===>\t2a^3 k/(1 + (ak)^2)^2\n";
  312.   cout << "\tr*exp( -r/a )\t<=== cos-transform ===>\t"
  313.           "a^2 (1-(ak)^2)/(1 + (ak)^2)^2\n";
  314.  
  315.   cout << "\nParameter a is " << form("%.2f",a);
  316.   cout << "\nNo. of grids   " << n;
  317.   cout << "\nGrid mesh in the r-space dr = " << form("%.3f",dr);
  318.   cout << "\nGrid mesh in the k-space dk = " << form("%.3f",dk);
  319.  
  320.   FFT fft(2*n,dr);
  321.  
  322.   cout << "\nCheck out the inquires to FFT package about N, dr, dk, cutoffs\n";
  323.   assert( fft.q_N() == 2*n );
  324.   assert( fft.q_dr() == dr );
  325.   assert( fft.q_dk() == dk );
  326.   assert( fft.q_r_cutoff() == R );
  327.   assert( fft.q_k_cutoff() == dk*n );
  328.  
  329.  
  330.   Vector xr(0,n-1);            // Tabulate the source function
  331.   register int j;
  332.   for(j=0; j<n; j++)
  333.   {
  334.     double r = j*dr;
  335.     xr(j) = r * exp( -r/a );
  336.   }
  337.  
  338.   Vector xs(xr), xc(xr);
  339.   start_timing();
  340.   fft.sin_transform(xs,xr);
  341.   print_timing("Sine transform");
  342.   start_timing();
  343.   fft.cos_transform(xc,xr);
  344.   print_timing("Cosine transform");
  345.  
  346.   Vector xs_ex(xs), xc_ex(xc);        // Compute exact transforms
  347.   for(j=0; j<n; j++)
  348.   {
  349.     double k = j * dk;
  350.     Complex t = 1/Complex(1/a,-k);
  351.     t = t*t;
  352.     xc_ex(j)  = t.real();
  353.     xs_ex(j)  = t.imag();
  354.   }
  355.  
  356.   compare(xs,xs_ex,"Computed and Exact sin-transform");
  357.   compare(xc,xc_ex,"Computed and Exact cos-transform");
  358.  
  359.   Vector xt(xc); xt = xc; xt -= xc(xc.q_upb());
  360.   compare(xt,xc_ex,"Computed cos-transform with DC component removed, and "
  361.       "exact result");
  362.  
  363.   Vector xs_inv(xr), xc_inv(xr);    // Compute Inverse transforms
  364.   start_timing();
  365.   fft.sin_inv_transform(xs_inv,xs);
  366.   print_timing("Inverse sine transform");
  367.   start_timing();
  368.   fft.cos_inv_transform(xc_inv,xc);
  369.   print_timing("Inverse cosine transform");
  370.  
  371.   compare(xs_inv,xr,"Computed inverse sin-transform vs the original function");
  372.   compare(xc_inv,xr,"Computed inverse cos-transform vs the original function");
  373.  
  374.   xt = xc_inv; xt -= xc_inv(xc_inv.q_upb());
  375.   compare(xt,xr,"Computed inverse cos-transform with DC component removed,\n"
  376.       "and the original function");
  377.  
  378.   cout << "\nDone\n";
  379. }
  380.  
  381. /*
  382.  *-----------------------------------------------------------------------
  383.  *             Verify the cos Fourier transform
  384.  *
  385.  *    exp( -r^2/4a )    <=== cos-transform ===> sqrt(a*pi) exp(-a*k^2)
  386.  */
  387.  
  388. static void test_gaussian(void)
  389. {
  390.   const double R = 20;            // Cutoff distance
  391.   const int n = 512;            // No. of grids
  392.   const double a = 4;            // Constant a, see above
  393.   const double dr = R/n,         // Grid meshes
  394.            dk = PI/R;    
  395.  
  396.   cout << "\n\nVerify the sin/cos transform for the following example\n";
  397.   cout << "\texp( -r^2/4a )\t<=== cos-transform ===> sqrt(a*pi) exp(-a*k^2)\n";
  398.  
  399.   cout << "\nParameter a is " << form("%.2f",a);
  400.   cout << "\nNo. of grids   " << n;
  401.   cout << "\nGrid mesh in the r-space dr = " << form("%.3f",dr);
  402.   cout << "\nGrid mesh in the k-space dk = " << form("%.3f",dk);
  403.  
  404.   FFT fft(2*n,dr);
  405.  
  406.   cout << "\nCheck out the inquires to FFT package about N, dr, dk, cutoffs\n";
  407.   assert( fft.q_N() == 2*n );
  408.   assert( fft.q_dr() == dr );
  409.   assert( fft.q_dk() == dk );
  410.   assert( fft.q_r_cutoff() == R );
  411.   assert( fft.q_k_cutoff() == dk*n );
  412.  
  413.  
  414.   Vector xr(0,n-1);            // Tabulate the source function
  415.   register int j;
  416.   for(j=0; j<n; j++)
  417.   {
  418.     double r = j*dr;
  419.     xr(j) = exp(-r*r/4/a );
  420.   }
  421.  
  422.   Vector xc(xr);
  423.   start_timing();
  424.   fft.cos_transform(xc,xr);
  425.   print_timing("Cosine transform");
  426.  
  427.   Vector xc_ex(xr);            // Compute exact transforms
  428.   for(j=0; j<n; j++)
  429.   {
  430.     double k = j * dk;
  431.     xc_ex(j)  = sqrt(PI*a) * exp( -a*k*k );
  432.   }
  433.  
  434.   compare(xc,xc_ex,"Computed and Exact cos-transform");
  435.   xc -= xc(xc.q_upb());
  436.   compare(xc,xc_ex,"Computed with DC removed, and Exact cos-transform");
  437.  
  438.   Vector xc_inv(xr);            // Compute Inverse transforms
  439.   start_timing();
  440.   fft.cos_inv_transform(xc_inv,xc);
  441.   print_timing("Inverse cosine transform");
  442.  
  443.   compare(xc_inv,xr,"Computed inverse cos-transform vs the original function");
  444.   xc_inv -= xc_inv(xc_inv.q_upb());
  445.   compare(xc_inv,xr,"Computed inverse with DC removed vs the original");
  446.  
  447.   cout << "\nDone\n";
  448. }
  449.  
  450. /*
  451.  *------------------------------------------------------------------------
  452.  *                Root module
  453.  */
  454.  
  455. main()
  456. {
  457.   cout << "\n\n" << _Equals << 
  458.        "\n\n\t\tVerify Fast Fourier Transform Package\n\n";
  459.  
  460.   test_ap_series(8);
  461.   test_ap_series(1024);
  462.  
  463.   test_orth(1024,1);
  464.  
  465.   test_ap_series_pad0(16);
  466.   test_ap_series_pad0(1024);
  467.  
  468.   test_lorentzian();
  469.   test_gaussian();
  470. }
  471.  
  472.