PythonのSciPyを使って離散フーリエ変換（DFT）を実行する方法

離散フーリエ変換（DFT：Discrete Fourier Transform）は、空間データを周波数データへ変換するための数学的手法です。

一方、高速フーリエ変換（FFT：Fast Fourier Transform）は、空間データの離散フーリエ変換を効率的に計算するために設計されたアルゴリズムです。DFTを直接計算すると計算量が大きくなりますが、FFTを利用することで大幅に処理を高速化できます。

ここでいう空間データとは、通常多次元配列の形式で表されるデータを指します。また、周波数データとは、特定の時間内に含まれる信号や波長の数に関する情報を持つデータのことです。

この記事では、Pythonの「SciPy」ライブラリを使用して、このDFTをどのように実現できるかを見ていきましょう。

グラフの作成にはmatplotlibライブラリを使用し、データの生成にはNumpyライブラリを使用します。

サンプルコード

from matplotlib import pyplot as plt
import numpy as np

my_freq = 6
freq_samp = 70
time_val = np.linspace(0, 3, 3 * freq_samp, endpoint=False)
amp_val = np.sin(my_freq * 3 * np.pi * time_val)

figure, axis = plt.subplots()
axis.plot(time_val, amp_val)
axis.set_xlabel('Time (in seconds)')
axis.set_ylabel('Amplitude of signal')
plt.show()

from scipy import fftpack
A = fftpack.fft(amp_val)
frequency = fftpack.fftfreq(len(amp_val)) * freq_samp

figure, axis = plt.subplots()
axis.stem(frequency, np.abs(A))
axis.set_xlabel('Frequency in Hz')
axis.set_ylabel('Frequency Spectrum Magnitude')
axis.set_xlim(-freq_samp / 2, freq_samp / 2)
axis.set_ylim(-7, 125)
plt.show()

出力結果

コードの解説

• まず、必要なパッケージ（matplotlib、numpy、scipy）をインポートします。
• Numpyライブラリのlinspace関数とsin関数を使って、正弦波のサンプルデータを生成します。ここでは周波数6Hzの信号を、サンプリング周波数70Hzで3秒間分作成しています。
• 生成したデータをmatplotlibライブラリを使って正弦波としてコンソール上にプロットし、時間領域での波形を可視化します。
• 次に、SciPyの「fftpack」パッケージを使用して、生成したデータの高速フーリエ変換（FFT）を計算します。fftfreq関数により、各周波数成分に対応する周波数の値も取得しています。
• 最後に、得られた周波数スペクトルの絶対値をステムプロットとしてグラフに描画し、信号に含まれる周波数成分を確認します。

このように、SciPyのfftpackモジュールを利用すれば、わずか数行のコードで離散フーリエ変換を実行し、時間領域の信号を周波数領域に変換して分析することができます。