File content as of revision 0:c740f515e0b7:

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