Diff: fftReal.cpp

Revision:

0:c740f515e0b7

diff -r 000000000000 -r c740f515e0b7 fftReal.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/fftReal.cpp	Tue Jul 15 10:14:37 2014 +0000
@@ -0,0 +1,76 @@
+//------------------------------------------------------------------------------
+// FFT class for real data
+// Copyright (c) 2014 MIKAMI, Naoki,  2014/06/29
+//------------------------------------------------------------------------------
+
+#include "fftReal.hpp"
+
+namespace Mikami
+{
+    // Constructor
+    FftReal::FftReal(uint16_t n) : N_FFT_(n)
+    {
+        wTable_ = new Complex[n/2];
+        bTable_ = new uint16_t[n/2];
+        u_ = new Complex[n];
+
+        // calculation of twiddle factor
+        Complex arg = Complex(0, -6.283185f/N_FFT_);
+        for (int k=0; k<N_FFT_/2; k++)
+            wTable_[k] = exp(arg*(float)k);
+
+        // for bit reversal
+        int nHalf = N_FFT_/2;
+        bTable_[0] = 0;
+        for (int m=1; m<N_FFT_/2; m*=2)
+        {
+            for (int k=0; k<m; k++)
+                bTable_[k+m] = bTable_[k] + nHalf;
+            nHalf = nHalf/2;
+        }
+    }
+
+    // Destructor
+    FftReal::~FftReal()
+    {
+        delete[] wTable_;
+        delete[] bTable_;
+        delete[] u_;    
+    }
+
+    // Execute FFT
+    void FftReal::Execute(const float x[], Complex y[])
+    {
+        for (int n=0; n<N_FFT_; n++) u_[n] = x[n];
+
+    // except for last stage
+        uint16_t nHalf = N_FFT_/2;
+        for (int stg=1; stg<N_FFT_/2; stg*=2)
+        {
+            uint16_t nHalf2 = nHalf*2;
+            for (int kp=0; kp<N_FFT_; kp+=nHalf2)
+            {
+                uint16_t kx = 0;
+                for (int k=kp; k<kp+nHalf; k++)
+                {
+                    // Butterfly operation
+                    Complex uTmp = u_[k+nHalf];
+                    u_[k+nHalf] = (u_[k] - uTmp)*wTable_[kx];
+                    u_[k] = u_[k] + uTmp;
+                    kx = kx + stg;
+                }
+            }
+            nHalf = nHalf/2;
+        }
+
+        // Last stage
+        y[0] = u_[0] + u_[1];
+        y[N_FFT_/2] = u_[0] - u_[1];
+        for (int k=2; k<N_FFT_; k+=2)
+            u_[k] = u_[k] + u_[k+1];
+
+        // Reorder according to bit reversal
+        for (int k=1; k<N_FFT_/2; k++)
+            y[k] = u_[bTable_[k]];
+    }
+}