自分用に書き換えてしまったので forkします
Fork of UIT_FFT_Real by
fftReal.cpp
- Committer:
- ohneta
- Date:
- 2015-10-15
- Revision:
- 1:148bf3a521a0
- Parent:
- 0:982a9acf3a07
File content as of revision 1:148bf3a521a0:
//------------------------------------------------------------------------------ // FFT class for real data usind decimation-in-frequency algorithm // 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" #include "SoundWS2812B-FFT.h" #ifdef _DEBUG_UART_ extern Serial pc; #endif // 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) { #ifdef _DEBUG_UART_ pc.printf("\r\nNot power of 2, in FftReal class."); pc.printf("\r\nForce to exit the program."); #endif exit(EXIT_FAILURE); // Terminate program } wTable_ = new Complex[n / 2]; if (wTable_ == NULL) { #ifdef _DEBUG_UART_ pc.printf("wTable_ is NULL.\n"); #endif exit(EXIT_FAILURE); } bTable_ = new uint16_t[n]; if (bTable_ == NULL) { #ifdef _DEBUG_UART_ pc.printf("bTable_ is NULL.\n"); #endif exit(EXIT_FAILURE); } u_ = new Complex[n]; if (u_ == NULL) { #ifdef _DEBUG_UART_ pc.printf("u_ is NULL.\n"); #endif exit(EXIT_FAILURE); } // 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 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 // for Cortex-M0 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]; } // 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]]; } } #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++) { u_[n] = conj(y[N_FFT_-n]); } // except for last stage ExcludeLastTtage(); // Last stage including bit reversal x[0] = N_INV_*(u_[0].real() + u_[1].real()); x[half] = N_INV_*(u_[0].real() - u_[1].real()); for (int n=2; n<N_FFT_; n+=2) { float un = u_[n].real(); float un1 = u_[n+1].real(); 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; } } nHalf = nHalf / 2; } }