DFT(離散フーリエ変換)係数を直接計算するC++プログラムの実装方法

離散フーリエ変換(DFT)とは、関数を等間隔でサンプリングして得られた有限個のデータ列を、複素正弦波の有限線形結合における係数のリストへ変換する手法です。係数は周波数順に並べられ、元のサンプル値と同じ値を持つように構成されます。これにより、サンプリングされた関数を元の領域(多くの場合、時間や直線上の位置)から周波数領域へと変換することができます。

本記事では、一次関数 ax + by = c のサンプル列に対してDFT係数を実装で直接計算するC++プログラムを紹介します。

アルゴリズム

処理の流れは以下の擬似コードの通りです。まず入力された係数からサンプル列を生成し、次に各周波数kに対応する余弦・正弦の値を求め、それらを畳み込むことでDFT係数の実部と虚部を計算します。

開始
   線形方程式の係数3つ(a, b, c)と最大値k用の変数を宣言する
   変数を読み込む
   実部(real)と虚部(img)の2つのメンバ変数を持つクラスを定義する
   コンストラクタを作成し、realとimgを0.0で初期化する
   変数Mを宣言し、適切な整数値で初期化する(サンプル数)
   配列 function[M] を作成する
   i = 0 から M-1 まで繰り返す
      function[i] = (((a * (double)i) + (b * (double)i)) - c)
   配列 sine[M] を宣言する
   配列 cosine[M] を宣言する
   i = 0 から M-1 まで繰り返す
      cosine[i] = cos((2 * i * k * PI) / M)
      sine[i] = sin((2 * i * k * PI) / M)
   i = 0 から M-1 まで繰り返す
      dft_value.real += function[i] * cosine[i]
      dft_value.img += function[i] * sine[i]
   計算結果を出力する
終了

サンプルコード(C++)

以下が実際のC++による実装例です。DFT係数の実部と虚部を保持するクラス DFT_Coeff を定義し、コンストラクタで両者を0に初期化している点がポイントです。

#include<iostream>
#include<math.h>
using namespace std;
#define PI 3.14159265
class DFT_Coeff {
    public:
    double real, img;
    DFT_Coeff() {
        real = 0.0;
        img = 0.0;
    }
};
int main(int argc, char **argv) {
    int M = 10;
    cout << "1次関数の係数を入力してください:\n";
    cout << "ax + by = c\n";
    double a, b, c;
    cin >> a >> b >> c;
    double function[M];
    for (int i = 0; i < M; i++) {
        function[i] = (((a * (double) i) + (b * (double) i)) - c);
    }
    cout << "最大K値を入力してください: ";
    int k;
    cin >> k;
    double cosine[M];
    double sine[M];
    for (int i = 0; i < M; i++) {
        cosine[i] = cos((2 * i * k * PI) / M);
        sine[i] = sin((2 * i * k * PI) / M);
    }
    DFT_Coeff dft_value;
    cout << "DFT係数は: ";
    for (int i = 0; i < M; i++) {
        dft_value.real += function[i] * cosine[i];
        dft_value.img += function[i] * sine[i];
    }
    cout << "(" << dft_value.real << ") - " << "(" << dft_value.img << " i)";
}

出力結果

係数 a=4、b=6、c=7、最大K値 k=4 を入力した場合の実行結果は以下のようになります。

1次関数の係数を入力してください:
ax + by = c
4 6 7
最大K値を入力してください:
4
DFT係数は: (-50) - (-16.246 i)

プログラムの解説

このプログラムの動作を簡単に整理すると、以下の3段階で構成されています。

1. サンプル列の生成

入力された係数 a、b、c をもとに、i = 0〜M-1 の範囲で f(i) = ai + bi − c を計算し、長さMのサンプル列を作成します。

2. 回転因子(三角関数)の計算

指定された周波数kに対して、cos(2πik/M) と sin(2πik/M) をそれぞれ配列に格納します。これはDFTの核となる複素指数関数 e^-j2πik/M の実部と虚部に相当します。

3. 係数の総和計算

サンプル列と三角関数の値を要素ごとに掛け合わせて総和を取ることで、DFT係数の実部(real)と虚部(img)が求まります。最終的な結果は「実部 − 虚部i」の形式で出力されます。

なお、本コードでは可変長配列(VLA)を使用しているため、厳密にはC++標準規格外となります。移植性を高めたい場合は、std::vector<double> の利用をおすすめします。