File content as of revision 4:0b4975fffc90:

//------------------------------------------------------------------------------
//  データが実数の場合の FFT class
//      FFT     複素共役になる部分は計算しない
//      IFFT    複素共役の部分を除いた値から計算する
//
//  2020/12/12, Copyright (c) 2020 MIKAMI, Naoki
//------------------------------------------------------------------------------

#include "fftReal.hpp"

namespace Mikami
{
    // コンストラクタ
    FftReal::FftReal(int16_t n)
            : N_FFT_(n), N_INV_(1.0f/n), wTable_(n/2), bTable_(n), u_(n)
    {
        // __clz(): リーディング 0 の数を取得する命令に対応
        uint32_t shifted = n << (__clz(n) + 1);
        MBED_ASSERT(shifted == 0);  // n が 2 のべき乗であることをチェック

        // 回転子の表を作成
        Complex arg = Complex(0, -6.283185f/N_FFT_);
        for (int k=0; k<N_FFT_/2; k++)
            wTable_[k] = exp(arg*(float)k);

        // ビット逆順のための表を作成
        uint16_t nShift = __clz(n) + 1;
        for (int k=0; k<n; k++)
            // __rbit(k): ビットの並びを逆にする命令に対応
            bTable_[k] = __rbit(k) >> nShift;
    }

    // FFT の実行
    //      有効な部分 y[0] ～ y[N_FFT_/2]
    //      それ以外の部分は複素共役になるので，計算をしていない
    void FftReal::Execute(const float x[], Complex y[])
    {
        for (int n=0; n<N_FFT_; n++) u_[n] = x[n];

        // 最終ステージを除いた部分
        ExcludeLastStage();

        // 最終ステージ
        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];

        // ビット逆順の並べ替え
        for (int k=1; k<N_FFT_/2; k++) y[k] = u_[bTable_[k]];
    }

    // IFFT の実行
    //      このクラスで計算された FFT の結果の IFFT を計算する
    //      実行結果    x[0] ～ x[N_FFT_-1]
    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]);      // 後半は複素共役になっているものとする

        // 最終ステージを除いた部分
        ExcludeLastStage();

        // 最終ステージとビット逆順の並べ替え処理
        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);
        }
    }

    // 最終ステージを除いた処理
    void FftReal::ExcludeLastStage()
    {
        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++)
                {
                    // バタフライ演算
                    Complex uTmp = u_[k+nHalf];
                    u_[k+nHalf] = (u_[k] - uTmp)*wTable_[kx];
                    u_[k] = u_[k] + uTmp;
                    kx = kx + stg;
                }
            }
            nHalf = nHalf/2;
        }
    }
}