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;
    }        
}