C++で複素数2次元配列に対して2D FFTをインプレースで実行する方法

2D FFT（二次元高速フーリエ変換）とは

高速フーリエ変換（FFT）は、離散フーリエ変換（DFT）およびその逆変換を効率的に計算するためのアルゴリズムです。フーリエ解析では、時間領域（または空間領域）の信号を周波数領域へ、あるいはその逆方向へと変換します。FFTは、DFT行列を疎行列（ほとんどの要素がゼロ）の積に因数分解することで、変換処理を大幅に高速化します。

本記事では、複素数を含む2次元配列（画像データなど）を対象に、C++で2D FFTをインプレースで実行するプログラムを紹介します。出力として、各周波数成分の実部、虚部、そして振幅（amp）を求めます。

アルゴリズム

Begin
    配列のサイズ n を宣言する
    2次元配列の要素を入力する
    実部用・虚部用・振幅用の3つの配列を宣言する
    height = n、width = n として初期化する
    出力データ（周波数側）を走査する外側の二重ループを作成する
    入力データ（空間側）を走査する内側の二重ループを作成する
    各成分について実部・虚部・振幅を計算する
End

サンプルコード

以下のC++プログラムでは、入力された n×n の2次元配列に対してDFTの定義式に基づく計算を行い、実部（realOut）、虚部（imgOut）、振幅（ampOut）を出力します。

#include <iostream>
#include <math.h>
using namespace std;
#define PI 3.14159265
int n;
int main(int argc, char **argv) {
    cout << "Enter the size: ";
    cin >> n;
    double Data[n][n];
    cout << "Enter the 2D elements ";
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            cin >> Data[i][j];
    double realOut[n][n];
    double imgOut[n][n];
    double ampOut[n][n];
    int height = n;
    int width = n;
    for (int yWave = 0; yWave < height; yWave++) {
        for (int xWave = 0; xWave < width; xWave++) {
            for (int ySpace = 0; ySpace < height; ySpace++) {
                for (int xSpace = 0; xSpace < width; xSpace++) {
                    realOut[yWave][xWave] += (Data[ySpace][xSpace] * cos(2 *
                        PI * ((1.0 * xWave * xSpace / width) + (1.0 * yWave * ySpace /
                        height)))) / sqrt(width * height);
                    imgOut[yWave][xWave] -= (Data[ySpace][xSpace] * sin(2 * PI
                        * ((1.0 * xWave * xSpace / width) + (1.0 * yWave * ySpace / height)))) /
                        sqrt( width * height);
                    ampOut[yWave][xWave] = sqrt(
                        realOut[yWave][xWave] * realOut[yWave][xWave] +
                        imgOut[yWave][xWave] * imgOut[yWave][xWave]);
                }
                cout << realOut[yWave][xWave] << " + " <<
                imgOut[yWave][xWave] << " i (" << ampOut[yWave][xWave] << ")\n";
            }
        }
    }
}

実行結果

Enter the size: 2
Enter the 2D elements
4 5
6 7
4.5 + 6.60611e-310 i (4.5)
11 + 6.60611e-310 i (11)
-0.5 + -8.97448e-09 i (0.5)
-1 + -2.15388e-08 i (1)
4.5 + 6.60611e-310 i (4.5)
-2 + -2.33337e-08 i (2)
-0.5 + -8.97448e-09 i (0.5)
0 + 5.38469e-09 i (5.38469e-09)

コードのポイント

このプログラムでは、外側の二重ループ（yWave、xWave）が出力側の周波数成分を、内側の二重ループ（ySpace、xSpace）が入力側の空間座標を走査します。各要素ごとに cos 項から実部を加算し、sin 項から虚部を減算することでDFTを構成し、最後に実部と虚部の平方和の平方根として振幅を求めています。

なお、この実装はDFTの定義を直接計算する方式のため、計算量は O(n⁴) となります。実際の大規模データ処理では、行方向・列方向に分割して1D FFTを適用する手法や、FFTWなどの最適化ライブラリを利用することで、O(n² log n) への高速化が可能です。