C++で離散フーリエ変換と高速フーリエ変換を実装

サウンドプログラマーを目指す人にとってフーリエ変換は「切っても切れない関係」と言っても
過言ではない知識、技術だと思っています。
しかし、「フーリエ変換とは何ぞ？」な人からするとwikiの公式を見たところでプログラムに起こすことも理解することもできないと思います。
今回はそのような人向けにできるだけ内容が把握しやすいコードにしたのでよかったら拝見してみてください。

環境
・Vidual Studio 2019

離散フーリエ変換

#include <cmath>  //数学関数
#include <vector> //可変長配列
#include <complex>//複素数
std::vector<std::complex<double>> DFT(const std::vector<double>& data)
{
    //虚数単位
    const std::complex<double>Imaginary(0, 1);
    //計算用バッファ
    std::vector<std::complex<double>>buf(data.size(), 0);

    for(size_t i = 0; i < data.size(); ++i)
    {
        //畳み込み
        for(size_t n = 0; n < data.size(); ++n)
        {
            /*
            std::acos(-1.0)は円周率を表しています。
            M_PIの使用がめんどくさくなったので代用しています。
            知っておくと便利なのでこの機会に覚えておくのも良いかもしれません。
            */
            //公式をそのまま書いた場合
            buf[i] += data[n] * std::exp(-Imaginary * 2.0 * std::acos(-1.0) * double(i * n) / double(data.size()));
            /*オイラーの公式を利用した場合
            buf[i] += data[n] * (std::cos(2.0 * std::acos(-1.0) * double(i * n) / double(data.size())) - (Imaginary * std::sin(2.0 * std::acos(-1.0) * double(i * n) / double(data.size()))));
            */
        }
    }
    return buf;
}

上記のコードが離散フーリエ変換になります。
実装する中でいろいろな方のサイトを拝見させていただき、参考にしてみた結果
変換→逆変換で元データとの誤差が起きてしまうことがほとんどでした。
最初はコンパイラの問題なのか、計算誤差の問題なのかと考えたのですが解決せず...
公式を見直すと「虚数単位」が自分の中で引っ掛かりました。
そこで今回のコードにもあるように「虚数単位」を考慮した計算式にしています。
結果、まだ誤差はあるが無視できる範囲にまで抑えることができました。

逆離散フーリエ変換は「虚数単位」の符号が「-」から「+」に代わり、
畳み込みが終わった後に「データ数で除算」したら完成です。
(std::complexのreal()「実数」部分が元のデータになります。)
※コードは割愛させていただきます。

おまけに高速フーリエ変換のコードも記載しておきます。

高速フーリエ変換

#include <cmath>  //数学関数
#include <vector> //可変長配列
#include <complex>//複素数
std::vector<std::complex<double>> FFT(const std::vector<double>& data)
{
    //虚数単位
    const std::complex<double>Imaginary(0, 1);
    //累乗を求める
    unsigned int exponent = int(std::ceil(std::log2(data.size())));

    //2の累乗の数
    unsigned int num = int(std::pow(2.0, exponent));

    //計算用バッファ
    std::vector<std::complex<double>>buf(num, 0);

    //ビット反転
    std::vector<unsigned int>index(buf.size());
    for (size_t i = 0; i < index.size(); ++i)
    {
        unsigned int tmp = i;
        for (unsigned int n = 0; n < exponent; ++n)
        {
             index[i] <<= 1;
             index[i] |= (tmp >> n) & 0x0001;
        }
        //渡されたデータの並びを変える
        if (index[i] < data.size())
        {
            buf[i] = data[index[i]];
        }
        else
        {
            //データ数を超えたら0埋め
            buf[i] = 0.0;
        }
    }

    //バタフライ演算
    for (unsigned int stage = 1; stage <= exponent; ++stage)
    {
        for (unsigned int i = 0; i < std::pow(2.0, exponent - stage); ++i)
        {
            for (unsigned int n = 0; n < std::pow(2.0, stage - 1); ++n)
            {
                std::complex<double> corre1 = std::exp(-Imaginary * 2.0 * std::acos(-1.0) * double(n) / std::pow(2.0, stage));
                std::complex<double> corre2 = std::exp(-Imaginary * 2.0 * std::acos(-1.0) * double(n + std::pow(2.0, stage - 1)) / std::pow(2.0, stage));

                unsigned int No1 = int(i * std::pow(2.0, stage) + n);
                unsigned int No2 = int(No1 + std::pow(2.0, stage - 1));

                std::complex<double> tmp1 = buf[No1];
                std::complex<double> tmp2 = buf[No2];

                buf[No1] = tmp1 + (tmp2 * corre1);
                buf[No2] = tmp1 + (tmp2 * corre2);
            }
        }
    }
    return buf;
}

最初にビット反転で並び替えを行い、バタフライ演算で求めていく流れの「時間間引き」を記載しました。
また、std::ceilによって渡されたデータ数以上で近しい2の累乗を求めることができます。

逆高速フーリエ変換は「虚数単位」の符号が「-」から「+」に代わり、
バタフライ演算が終わった後に「データ数で除算」したら完成です。
ここら辺は離散フーリエ変換と一緒なので分かりやすいですね。
※コードは割愛させていただきます。

よければ多くの方の参考になればうれしいです。