FFT(Fast Fourier transform)をC++で実装する

はじめに

今回は、波形を解析する際に重要となるフーリエ変換について、FFTを具体的に実装することで学んでいく事を目的としました。任意個数のデータをFFTできるように本当はしたかったのですが、難しかったので諦めました。

環境

Visual Studio 2019 communityを用いてコーディングを行いました。

FFTについて

FFTとは、すごく大まかに言うと、時間tによって変化する関数f(t)を周期ωによる関数F(ω)に変換するフーリエ変換について、t、ω両方を離散的に考えて変換を行うDFTを改良し、計算のオーダーをO(N²)からO(NlogN)にしたものです。
もし、FFTを勉強する上で必要になるフーリエ変換という計算の詳細が分からないけど興味がある人がいればこのページでフーリエ変換について勉強する事をお勧めします。

FFTの性質

今回のFFTを実装する上で、以下の3つの性質が重要になります。

1. バタフライ演算
   。DFTの計算を並び替えて、同じ計算や値を纏めて用いる事が出来るようにする演算で、この手法によって計算のオーダーが減るようになります。

2. ビット逆順
   。バタフライ演算を行う上で重要になるインデックスの参照について、ビットを逆順にして上から並べるだけでバタフライの配置を再現できる性質があります。(ビット逆順の例：1000 -> 0001 、1010 -> 0101)

3. 計算に必要な変数格納領域はN個
   。バタフライ演算を行う上で、計算する値を格納するための領域が当然必要になります。しかし、その際にバタフライ演算の途中の値を保管する領域が必要になるように思えます。ですが、バタフライ演算を行う際に、一度使って加工されたデータは、その演算の中で、もう使われることがないという性質があります。つまり、最初のデータが格納されているであろうN個の保存領域があれば、その領域を使いまわしても計算可能であるという性質が考えられます。

実際に書いたコード

まず、初めにFFTを行う関数をラップするクラスを宣言しました。(個人的に使える機能を増やしてるので、最低限以外の機能があります)

クラスの要素

template<class T> class Fourier {
public:
    //T型のデータを受け取ってそれを使う（g(n)を受け取る）
    Fourier(const T* data, int size);
    Fourier(int size);
    Fourier();
    ~Fourier();

    //フーリエ変換するデータの更新
    void update(const T* data, int size);

    //FFT関数
    void FFT(const T* data, int size);
    //DFT関数
    void DFT(const T* data, int size);

    //位置指定をして特定のデータを受け取る
    double CkValue(int) const;
    double CkRange(int) const;

    //ポインタでごっそり受け取る
    const double* const CkValue() const {
        return mCkValue;
    }
    const double* const CkRange() const {
        return mCkRange;
    }
    const int difference() const {
        return mDifference;
    }

    void extractionValue(double point) {
        for (int i = 0; i < mSize; ++i) {
            mCkValue[i] *= (mCkValue[i] > point);
        }
    }

    //「=」の処理
    void operator=(Fourier& f);

    int size() {
        return mSize;
    }

private:
    //データに関する変数
    int mSize;
    //2のべき乗と受け取ったデータの差分
    int mDifference;

    //mCkの絶対値成分と角度成分が入る
    double* mCkValue;
    double* mCkRange;

    //ローカルで宣言するとFFTを何回も呼び出す場合遅くなるので、バッファを計算時間短縮のためにここで宣言しておく
    double* CkSin;      //虚数成分(絶対値成分も含まれる)
    double* CkCos;      //実数成分(絶対値成分も含まれる)
    double* Wsin;       //e^jθの虚数成分であるsinθの値を入れる
    double* Wcos;       //e^jθの実数成分であるcosθの値を入れる
};

次にFFTを行う関数の説明...と行きたいところですが、その前に、バタフライ演算を行う際に必要になる特殊なインデックス順である逆順ビットを再現するための、ビットを逆順にする関数を以下に示します。

受け取ったビットを逆順にする関数

///iは普通のインデックスで、mはビット数
int reverseIndex(int index, int m) {
    int t = 0;
    //ビット反転処理
    for (int i = m - 1; i >= 0; --i) {
        t |= (0x00000001 & (index >> i)) << ((m - 1) - i);  //一番左のビットから論理和して右に詰めていく
    }
    return t;
}

そしてFFTを行う関数をメンバとして実装しました。(waveファイルのデータを受ける事を想定して、shortでもcharでも受け取れるようにするためにtemplateを用いてます)

FFTの中身

//FFTの実体（2のべき乗以外の数は0で埋めるので注意）
template<class T> void Fourier<T>::FFT(const T* data, int size) {
    //Mはデータ数上限の時のビット数。-1によって、例えば128個の要素があった時に0から127の数字だと出来る
    unsigned int M = log2(mSize - 1) + 1; int r = 0;
    //dataの中身をクラスの中に取り入れる(インデックスをビット逆順で格納)
    for (int i = 0; i < mSize; ++i) {
        CkSin[i] = 0;
        if ((r = reverseIndex(i, M)) > size) {
            CkCos[i] = 0;
        }
        else {
            CkCos[i] = data[r];
        }
        //初期値はビット逆順な値
    }

    //stepは2をステップ数分乗算したもの
    int step = 1;
    //MはNを表現するのに必要なビット数。このループ表現がO(logN * N)になる所以である
    for (int i = 0; i < M; ++i) {
        //このタイミングでstepを倍にすることで係数の数を一致させつつ、最初が2、次は4、8、16、...となる表現を実現してる
        step *= 2;
        //予め必要な係数を計算しておく。係数の数はstepの1/2倍分ある
        for (int j = 0; j < step / 2; ++j) {
            //計算 W(kn N) = exp(-j2π/N * kn)を参考にW(j step)を実装
            Wsin[j] = sin(2 * PI * (double)j / step);
            Wcos[j] = cos(2 * PI * (double)j / step);
        }
        //ここでfor文を2重にしているので、o(N^2)の処理をしている様に見えるが、
        //実際は段階を縦に区切る処理が入ってるだけなので、O(N)分の計算量である
        for (int k = 0; k < mSize; k += step) {
            for (int j = step / 2; j < step; ++j) {
                //バタフライの上側を示すインデックス値
                int upBuff = k + j - (step / 2);
                //バタフライの下側を示すインデックス値
                int downBuff = k + j;
                //掛ける係数を示すインデックス
                int Wbuff = j - (step / 2);

                //＊バタフライ的な計算の開始地点 (cosは実成分、sinは虚数成分を示している)
                //4回の掛け算と2回の加算
                double WXcosBuf = CkCos[downBuff] * Wcos[Wbuff] - CkSin[downBuff] * Wsin[Wbuff];
                double WXsinBuf = CkCos[downBuff] * Wsin[Wbuff] + CkSin[downBuff] * Wcos[Wbuff];

                //2つの要素を用いて2つの要素を更新する
                CkSin[downBuff] = CkSin[upBuff] - WXsinBuf;
                CkCos[downBuff] = CkCos[upBuff] - WXcosBuf;
                CkSin[upBuff] = CkSin[upBuff] + WXsinBuf;
                CkCos[upBuff] = CkCos[upBuff] + WXcosBuf;
            }
        }
    }

    for (int i = 0; i < mSize; ++i) {
        mCkValue[i] = sqrt(CkSin[i] * CkSin[i] + CkCos[i] * CkCos[i]);  //√(Cksin^2 + Ckcos^2)
        mCkRange[i] = atan2(CkSin[i], CkCos[i]);    //tan^-1(y/x)
    }
}

大体の内容はコメントに書かれていますが、一応大まかに説明をしたいと思います。

1. バタフライ演算を行う下準備
   今回のデータが何ビットで表現可能かを獲得->受け取った離散波形データを、計算に用いる実数領域に格納（この際、データはi番目の時、iの逆順ビット番目のデータがi番目のCkCosに入る）という流れで処理をおこなっています。

2. バタフライ演算
   まず初めに、バタフライ演算が2のべき乗ごとのインデックス参照を行うことが考えられたので、stepという、今のステップの際に2の何乗のデータを要するかを格納する変数を用意しました。
   そしてそれを用いてバタフライの演算を行っています。バタフライの計算については、それぞれのステップの段階に応じて、計算する組み合わせを考慮しつつ、積と和を行うといったものですが、結構ややこしい参照順になっているので、もしバタフライ演算の全ての動作を追うのなら、バタフライを表した図を参考にしながらプログラムを追っていけば理解が出来ると思います。（図の参考やメモなしでやるとかなりきついと思います）この演算を考える時に重要になるのは、インデックスの参照の仕方と、虚数成分を計算に組み込むので、一度の掛け算で4回の掛け算と2回の足し算が行われる事に注意することです。

3. 結果の格納
   最後に、バタフライ演算終了後のCkSin、CkCosの各領域はフーリエ変換が終わった後のデータ（周波数がnの時にどのようなパルスが出るかのデータ）が、実成分と虚数成分で分かれているので、それを絶対値成分と角度成分として格納します。

以上の主な3つの処理によってFFTが構成されました。

おまけ

データ数が2のべき乗数でない時には、アホみたいに計算がかかりますが、DFTが動く様になっているので、おまけとして一応DFTのプログラムも載せておきます。

DFT

template<class T>
void Fourier<T>::DFT(const T* data, int size) {
    mSize = size;
    for (int i = 0; i < mSize; ++i) {
        CkSin[i] = 0;
        CkCos[i] = 0;
        for (int j = 0; j < mSize; ++j) {
            CkSin[i] += data[j] * -(sin((2 * PI / mSize) * i * j));
            CkCos[i] += data[j] * cos((2 * PI / mSize) * i * j);
        }
        CkSin[i] /= mSize;
        CkCos[i] /= mSize;
        mCkValue[i] = sqrt(CkSin[i] * CkSin[i] + CkCos[i] * CkCos[i]);  //√(Cksin^2 + Ckcos^2)
        mCkRange[i] = atan2(CkSin[i], CkCos[i]);    //tan^-1(y/x)
    }
}

FFTと比べるとすごく単純ですよね！でもやっぱりFFTと比べるとヤバイほど重いので、なるべく使わない様に気を付けましょう。

おわりに

ここまで、FFTの実装について長々と書いていましたが、FFTを行いたい場合は用意されているライブラリを用いるべきです。理由としては、論文などで自作のFFTを使うと、ちゃんとFFTが出来ているかの保証が出来ないからです。
また、今回の内容では、任意個数でのFFTが出来ないです。もしそちらに興味があれば、この記事に書いてあったので、そこで確認すると良いと思います。