Diff: fftReal.cpp

Revision:

1:ba9ce95ec9a4

Parent:

0:982a9acf3a07

diff -r 982a9acf3a07 -r ba9ce95ec9a4 fftReal.cpp
--- a/fftReal.cpp	Fri Dec 19 12:10:34 2014 +0000
+++ b/fftReal.cpp	Thu Oct 15 15:37:24 2015 +0000
@@ -3,75 +3,103 @@
 //      This class can execute FFT and IFFT 
 // Copyright (c) 2014 MIKAMI, Naoki,  2014/12/19
 //------------------------------------------------------------------------------
+// for LPPC1114FN28(Cortex-M0) modified by Takehisa Oneta(ohneta) 4th Aug.2015
 
 #include "fftReal.hpp"
+#include "mbed.h"
 
-namespace Mikami
+extern Serial pc;
+
+// Constructor
+FftReal::FftReal(int16_t n) : N_FFT_(n), N_INV_(1.0f / n)
 {
-    // Constructor
-    FftReal::FftReal(int16_t n)
-            : N_FFT_(n), N_INV_(1.0f/n)
-    {
-        // __clz(): Count leading zeros
-        uint32_t shifted = n << (__clz(n)+1);
-        if (shifted != 0)
-        {
-            fprintf(stderr, "\r\nNot power of 2, in FftReal class.");
-            fprintf(stderr, "\r\nForce to exit the program.");
-            exit(EXIT_FAILURE); // Terminate program
-        }
+    // __clz(): Count leading zeros
+    uint32_t shifted = n << (__clz(n) + 1);
+    if (shifted != 0) {
+        pc.printf("\r\nNot power of 2, in FftReal class.");
+        pc.printf("\r\nForce to exit the program.");
+        exit(EXIT_FAILURE); // Terminate program
+    }
 
-        wTable_ = new Complex[n/2];
-        bTable_ = new uint16_t[n];
-        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);
+    wTable_ = new Complex[n / 2];
+    if (wTable_ == NULL) {
+        pc.printf("wTable_ is NULL.\n");
+        exit(EXIT_FAILURE);
+    }
+    bTable_ = new uint16_t[n];
+    if (bTable_ == NULL) {
+        pc.printf("bTable_ is NULL.\n");
+        exit(EXIT_FAILURE);
+    }
+    u_ = new Complex[n];
+    if (u_ == NULL) {
+        pc.printf("u_ is NULL.\n");
+        exit(EXIT_FAILURE);
+    }
 
-        // for bit reversal table
-        uint16_t nShift = __clz(n) + 1;
-        for (int k=0; k<n; k++)
-            // __rbit(k): Reverse the bit order in a 32-bit word
-            bTable_[k] = __rbit(k) >> nShift;        
+    // 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);
     }
 
-    // Destructor
-    FftReal::~FftReal()
-    {
-        delete[] wTable_;
-        delete[] bTable_;
-        delete[] u_;    
+    // for bit reversal table
+    uint16_t nShift = __clz(n) + 1;
+    for (int k = 0; k < n; k++) {
+        //bTable_[k] = __rbit(k) >> nShift;   // __rbit(k): Reverse the bit order in a 32-bit word
+        unsigned int x = (unsigned int)k;
+        x = ((x & 0x55555555) << 1) | ((x & 0xAAAAAAAA) >> 1);
+        x = ((x & 0x33333333) << 2) | ((x & 0xCCCCCCCC) >> 2);
+        x = ((x & 0x0F0F0F0F) << 4) | ((x & 0xF0F0F0F0) >> 4);
+        x = ((x & 0x00FF00FF) << 8) | ((x & 0xFF00FF00) >> 8);
+        x = ((x & 0x0000FFFF) << 16) | ((x & 0xFFFF0000) >> 16);
+        bTable_[k] = (x >> nShift);
+    }
+}
+
+// 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];
     }
 
-    // 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
-        ExcludeLastTtage();
+    // except for last stage
+    ExcludeLastTtage();
 
-        // 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 to bit reversal
-        for (int k=1; k<N_FFT_/2; k++)
-            y[k] = u_[bTable_[k]];
+    // 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];
     }
 
-    // Execute IFFT
-    void FftReal::ExecuteIfft(const Complex y[], float x[])
+    // Reorder to bit reversal
+    for (int k = 1; k < N_FFT_ / 2; k++) {
+        y[k] = u_[bTable_[k]];
+    }
+}
+
+#if 0
+// Execute IFFT
+void FftReal::ExecuteIfft(const Complex y[], float x[])
     {
         int half = N_FFT_/2;
 
-        for (int n=0; n<=half; n++) u_[n] = y[n];
-        for (int n=half+1; n<N_FFT_; n++)
+        for (int n=0; n<=half; n++) {
+             u_[n] = y[n];
+        }
+        for (int n=half+1; n<N_FFT_; n++) {
             u_[n] = conj(y[N_FFT_-n]);
+        }
 
         // except for last stage
         ExcludeLastTtage();
@@ -87,29 +115,27 @@
             x[Index(n)]   = N_INV_*(un + un1);
             x[Index(n+1)] = N_INV_*(un - un1);
         }
-    }
+}
+#endif
 
-    // Processing except for last stage
-    void FftReal::ExcludeLastTtage()
-    {
-        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;
-                }
+// Processing except for last stage
+void FftReal::ExcludeLastTtage()
+{
+    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;
-        }        
-    }
+        }
+        nHalf = nHalf / 2;
+    }        
 }
 
+