home *** CD-ROM | disk | FTP | other *** search
/ Netrunner 2004 October / NETRUNNER0410.ISO / regular / dctm.lzh / DCTM / source.lzh / source / fft.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  2003-02-08  |  3.4 KB  |  150 lines
  1. /*************************************************************************
  2.     fft.cpp
  3.  
  4.     03/02/08    Xiaohong
  5. *************************************************************************/
  6. #include "fft.h"
  7. #include <math.h>
  8. #include <windows.h>
  9.  
  10. #define    PI 3.14159265358979
  11. #define    BLKSIZE_S   256
  12. #define    LOGBLKSIZE_S 8
  13.  
  14. class fft_values
  15. {
  16. public:
  17.     fft_values();
  18.     double w_real[2][LOGBLKSIZE], w_imag[2][LOGBLKSIZE];
  19. };
  20. fft_values::fft_values()
  21. {
  22.     int L,le,le1;
  23.     ZeroMemory(w_real,sizeof(w_real));
  24.     ZeroMemory(w_imag,sizeof(w_imag));
  25.     for(L=0; L<LOGBLKSIZE; L++){
  26.         le = 1 << (LOGBLKSIZE-L);
  27.         le1 = le >> 1;
  28.         w_real[0][L] = cos(PI/((double)le1));
  29.         w_imag[0][L] = -sin(PI/((double)le1));
  30.     }          
  31.     for(L=0; L<LOGBLKSIZE_S; L++){
  32.         le = 1 << (LOGBLKSIZE_S-L);
  33.         le1 = le >> 1;
  34.         w_real[1][L] = cos(PI/((double)le1));
  35.         w_imag[1][L] = -sin(PI/((double)le1));
  36.     }
  37. }
  38.  
  39.  
  40. bool fft(double x_real[BLKSIZE],
  41.          double x_imag[BLKSIZE],
  42.          double energy[BLKSIZE],
  43.          double phi[BLKSIZE],
  44.          const int N)
  45. {
  46.     int     M,MM1;
  47.     int     NV2, NM1, MP;
  48.     int        i,j,k,L;
  49.     int        ip, le,le1;
  50.     double    t_real, t_imag, u_real, u_imag;
  51.     static fft_values stFftValues;
  52.  
  53.     switch(N) {
  54.     case BLKSIZE:
  55.         M = LOGBLKSIZE;
  56.         MP = 0;
  57.         break;
  58.     case BLKSIZE_S:
  59.         M = LOGBLKSIZE_S;
  60.         MP = 1;
  61.         break;
  62.     default:
  63.         return false;
  64.     }
  65.     MM1 = M-1;
  66.     NV2 = N >> 1;
  67.     NM1 = N - 1;
  68.     for(L=0; L<MM1; L++){
  69.         le = 1 << (M-L);
  70.         le1 = le >> 1;
  71.         u_real = 1;
  72.         u_imag = 0;
  73.         for(j=0; j<le1; j++)
  74.         {
  75.             for(i=j; i<N; i+=le)
  76.             {
  77.                 ip = i + le1;
  78.                 t_real = x_real[i] + x_real[ip];
  79.                 t_imag = x_imag[i] + x_imag[ip];
  80.                 x_real[ip] = x_real[i] - x_real[ip];
  81.                 x_imag[ip] = x_imag[i] - x_imag[ip];
  82.                 x_real[i] = t_real;
  83.                 x_imag[i] = t_imag;
  84.                 t_real = x_real[ip];
  85.                 x_real[ip] = x_real[ip]*u_real - x_imag[ip]*u_imag;
  86.                 x_imag[ip] = x_imag[ip]*u_real + t_real*u_imag;
  87.             }
  88.             t_real = u_real;
  89.             u_real = u_real*stFftValues.w_real[MP][L] - u_imag*stFftValues.w_imag[MP][L];
  90.             u_imag = u_imag*stFftValues.w_real[MP][L] + t_real*stFftValues.w_imag[MP][L];
  91.         }
  92.     }
  93.     /* special case: L = M-1; all Wn = 1 */
  94.     for(i=0; i<N; i+=2)
  95.     {
  96.         ip = i + 1;
  97.         t_real = x_real[i] + x_real[ip];
  98.         t_imag = x_imag[i] + x_imag[ip];
  99.         x_real[ip] = x_real[i] - x_real[ip];
  100.         x_imag[ip] = x_imag[i] - x_imag[ip];
  101.         x_real[i] = t_real;
  102.         x_imag[i] = t_imag;
  103.         energy[i] = x_real[i]*x_real[i] + x_imag[i]*x_imag[i];
  104.         if(energy[i] <= 0.0005)
  105.         {
  106.             phi[i] = 0;
  107.             energy[i] = 0.0005;
  108.         }
  109.         else
  110.             phi[i] = atan2((double) x_imag[i],(double) x_real[i]);
  111.         energy[ip] = x_real[ip]*x_real[ip] + x_imag[ip]*x_imag[ip];
  112.         if(energy[ip] == 0)
  113.             phi[ip] = 0;
  114.         else
  115.             phi[ip] = atan2((double) x_imag[ip],(double) x_real[ip]);
  116.     }
  117.     /* this section reorders the data to the correct ordering */
  118.     j = 0;
  119.     for(i=0; i<NM1; i++)
  120.     {
  121.         if(i<j)
  122.         {
  123.             /* use this section only if you need the FFT in complex number form *
  124.             * (and in the correct ordering)                                    */
  125.             t_real = x_real[j];
  126.             t_imag = x_imag[j];
  127.             x_real[j] = x_real[i];
  128.             x_imag[j] = x_imag[i];
  129.             x_real[i] = t_real;
  130.             x_imag[i] = t_imag;
  131.             /* reorder the energy and phase, phi                                        */
  132.             t_real = energy[j];
  133.             energy[j] = energy[i];
  134.             energy[i] = t_real;
  135.             t_real = phi[j];
  136.             phi[j] = phi[i];
  137.             phi[i] = t_real;
  138.         }
  139.         k=NV2;
  140.         while(k<=j)
  141.         {
  142.             j = j-k;
  143.             k = k >> 1;
  144.         }
  145.         j = j+k;
  146.     }
  147.  
  148.     return true;
  149. }
  150.