fftReal.cpp

00001 //------------------------------------------------------------------------------
00002 // FFT class for real data
00003 // Copyright (c) 2014 MIKAMI, Naoki,  2014/06/29
00004 //------------------------------------------------------------------------------
00005 
00006 #include "fftReal.hpp"
00007 
00008 namespace Mikami
00009 {
00010     // Constructor
00011     FftReal::FftReal(uint16_t n) : N_FFT_(n)
00012     {
00013         wTable_ = new Complex[n/2];
00014         bTable_ = new uint16_t[n/2];
00015         u_ = new Complex[n];
00016 
00017         // calculation of twiddle factor
00018         Complex arg = Complex(0, -6.283185f/N_FFT_);
00019         for (int k=0; k<N_FFT_/2; k++)
00020             wTable_[k] = exp(arg*(float)k);
00021 
00022         // for bit reversal
00023         int nHalf = N_FFT_/2;
00024         bTable_[0] = 0;
00025         for (int m=1; m<N_FFT_/2; m*=2)
00026         {
00027             for (int k=0; k<m; k++)
00028                 bTable_[k+m] = bTable_[k] + nHalf;
00029             nHalf = nHalf/2;
00030         }
00031     }
00032 
00033     // Destructor
00034     FftReal::~FftReal()
00035     {
00036         delete[] wTable_;
00037         delete[] bTable_;
00038         delete[] u_;    
00039     }
00040 
00041     // Execute FFT
00042     void FftReal::Execute(const float x[], Complex y[])
00043     {
00044         for (int n=0; n<N_FFT_; n++) u_[n] = x[n];
00045 
00046     // except for last stage
00047         uint16_t nHalf = N_FFT_/2;
00048         for (int stg=1; stg<N_FFT_/2; stg*=2)
00049         {
00050             uint16_t nHalf2 = nHalf*2;
00051             for (int kp=0; kp<N_FFT_; kp+=nHalf2)
00052             {
00053                 uint16_t kx = 0;
00054                 for (int k=kp; k<kp+nHalf; k++)
00055                 {
00056                     // Butterfly operation
00057                     Complex uTmp = u_[k+nHalf];
00058                     u_[k+nHalf] = (u_[k] - uTmp)*wTable_[kx];
00059                     u_[k] = u_[k] + uTmp;
00060                     kx = kx + stg;
00061                 }
00062             }
00063             nHalf = nHalf/2;
00064         }
00065 
00066         // Last stage
00067         y[0] = u_[0] + u_[1];
00068         y[N_FFT_/2] = u_[0] - u_[1];
00069         for (int k=2; k<N_FFT_; k+=2)
00070             u_[k] = u_[k] + u_[k+1];
00071 
00072         // Reorder according to bit reversal
00073         for (int k=1; k<N_FFT_/2; k++)
00074             y[k] = u_[bTable_[k]];
00075     }
00076 }