離散フーリエ変換の実装と結果解釈【No.2 離散フーリエ変換連載】 - YU2TA7KA's BLOG ～take one step at a time～

はじめに
実装方針
実装コード
サンプルデータの入力と出力
結果解釈
- 周波数分解能
- DFTのパラメーター
おわりに

はじめに

前回は離散フーリエ変換(DFT)を理解するための参考情報を示し、数式を解説しました。今回は、離散フーリエ変換のpython実装と結果解釈を行います。

実装方針

離散フーリエ変換の数式は以下です。

以下の関数を実装し、入力データを与え、出力グラフを得られるようにします。

サンプリングデータ数からDFT行列(W)を得る関数
DFTを演算する関数
DFT結果をグラフ表示する関数

実装コード

import numpy as np
import matplotlib.pyplot as plt

def generate_dft_matrix(size):
    n = size
    w = np.array([(np.exp(-2j * np.pi * k * np.arange(n) / n)) for k in range(n)])
    w.real[abs(w.real) <= 1.0e-10] = 0
    w.imag[abs(w.imag) <= 1.0e-10] = 0
    print(f"DFT matrix:\n{w}")
    return w

def dft(x):
    n = x.size
    w = generate_dft_matrix(n)
    x_dft = np.dot(w, x) / n
    print(f"DFT result: {x_dft}")
    return x_dft

def plot_dft_result(x, dft_result, period, sampling_number):
    frequency_resolution = 1 / period
    k_array = np.arange(0, sampling_number, 1)
    frequency_array = np.array([k * frequency_resolution if k <= sampling_number / 2 else (k - sampling_number) * frequency_resolution for k in range(sampling_number)])
    plt.stem(frequency_array, dft_result.real, linefmt="blue")
    plt.stem(frequency_array, dft_result.imag, linefmt="green")
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Amplitude")
    plt.grid()
    plt.show()

# input data
x = np.array([3,0,-3,0])
print(f"Input data: {x}")
period = 8 # period[sec] T0
sampling_number = 4 # sampling number

dft_result = dft(x)

plot_dft_result(x, dft_result, period, sampling_number)

サンプルデータの入力と出力

サンプルデータはこちらの例題を利用します。
それぞれの結果は以下になり、参照したPDFの解答と一致しています。

No.	サンプルデータ	データ数	サンプリング時間
1	[3,0,-3,0]	4	8[sec]

2	[1,1,-1,-1]	4	1[sec]

3	[2,0,2,0]	4	2[sec]

結果解釈

2	[1,1,-1,-1]	4	1[sec]

上記を例に結果を解釈します。DFTによって入力データに含まれる周波数成分を求めています。変換結果は実部と虚部で構成されます。この例では、DFT resultの1,3番目要素(0番目始まり)とグラフ画像から、実部虚部ともに±1.0Hzに成分が含まれていることが読み取れます。グラフ画像の-1.0Hzは実部(青線)もあるはずですが、虚部(赤線)で上書きされています。

周波数分解能

グラフ画像の横軸に周波数で数値があるため、ここから周波数成分を読み取れます。この軸の求め方の理解に私は時間がかかりました。結論としては入力データのサンプリング時間とサンプリング数から周波数分解能が求まります。この周波数分解能がDFTで得られる周波数成分の最小単位で、この倍数ごとにどの程度周波数成分を含むかがDFTによって解ります。
例えば、サンプリング時間が1.0secでサンプリング数が4個の場合、周波数分解能は1.0Hzで-1.0, 0, 1.0. 2.0Hzの周波数成分が得られます。

DFTのパラメーター

周波数分解について説明しましたが、DFTで考慮するパラメーターを以下にまとめます。これらのパラメーターを意識することで、入出力データを解釈できると思います。

変数名	記号	概要
サンプリング時間	To	入力データの取得開始から終了までの時間
サンプリング数	N	入力データの取得数
サンプリング時間間隔	⊿T(=To/N)	入力データを取得する時間感覚
周波数分解能	⊿f(=1/To)	DFTで得られる周波数成分の最小基準単位、検出能力
サンプリング周波数	fs(=N⊿f)	単位時間あたりにサンプリング(取得)するデータ数