Python NumPy SciPy サンプルコード: フーリエ変換処理その 1

かくあき

Python の fft 関数でのデータ処理法について、何回かに分けてまとめていきます。

Python の fft 関数

時系列データのフーリエ変換処理は、データの周波数領域での特徴抽出のために様々な分野で利用されています。機械工学の分野では、加速度計で構造物の加速度データを取得し、テータを周波数解析したりすることが多いと思います。

fft 関数でのデータ処理をやろうとした場合、時系列データと周波数データとの関係を理解しておかないと適切なデータ処理ができません。以下のような簡単なプログラムで fft 関数の使い方を説明していきます。

時系列のサンプルデータとして、データ数 512 点、サンプリング間隔 dt=0.01[sec]、周波数 f=20[Hz]の sin 波を作成し、それを fft 関数で離散フーリエ変換しています。

import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 時系列のサンプルデータ作成
N = 512                         # データ数
dt = 0.01                       # サンプリング間隔
f = 20                          # 周波数
t = np.linspace(1, N, N)*dt-dt
y = np.sin(2*np.pi*f*t)

# 離散フーリエ変換
yf = fft(y)

# プロット
# 時系列データ
plt.figure(1)
plt.plot(t, y)
plt.xlim(0, 1)
plt.xlabel("time")
plt.ylabel("y")

# 離散フーリエ変換の結果
fig = plt.figure(2)
ax = fig.add_subplot(111, projection='3d')
ax.plot3D(np.linspace(1, N, N), np.real(yf), np.imag(yf))
ax.set_xlabel("data number")
ax.set_ylabel("real part")
ax.set_zlabel("imaginary part")

# 大きさ、位相
plt.figure(3)
plt.subplot(211)
plt.plot(np.linspace(1, N, N), np.abs(yf))
plt.axis('tight')
plt.ylabel("amplitude")
plt.subplot(212)
plt.plot(np.linspace(1, N, N), np.degrees(np.angle(yf)))
plt.axis('tight')
plt.xlabel("data number")
plt.ylabel("phase[deg]")

# plt.show()

図 1 は時系列の波形、図 2 は fft 処理されたデータが複素数であることを示すために実部、虚部を各軸に表示させた 3 次元プロットです。 Python で 3 次元プロットしてみると、数式や 2 次元表現だけではイメージしにくかった複素数も理解しやすくなると思います。図 3 は fft 関数で処理されたデータの大きさと位相を表示しています。 NumPy の abs, angle 関数を使って、複素数の大きさと位相を求めることができます。図 3 を見るとデータに対称性があることがわかりますが、次回はこの特徴について説明していきます。

図1: 時系列データ

図2: fft 処理されたデータ

図3: fft 処理されたデータの大きさと位相

Python の fft 関数

コメント