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

前回 に引き続き、Python の fft 関数でのデータ処理法についてまとめていきます。

FFT 処理したデータとサンプリング定理との関係

前回は時系列のサンプルデータを SciPy の fft 関数で離散フーリエ変換し、それを単純にプロットしてみました。
プロットを見るとデータの対称性に気づきますが、今回はこの点について説明をしていきます。
まあ結論から言うと、これはサンプリング定理に関連しています。

まずは以下のような適当な 8 点のデータを用いて fft 処理したデータの結果を見てみます。

# coding: utf-8
import numpy as np
from scipy.fftpack import fft

y = np.array([[1, 2, 3, 4, 5, 6, 7, 8]])
yf = fft(y)

# 見やすいように結果を転置して表示
print(yf.T)
[[ 36.+0.j        ]
 [ -4.+9.65685425j]
 [ -4.+4.j        ]
 [ -4.+1.65685425j]
 [ -4.+0.j        ]
 [ -4.-1.65685425j]
 [ -4.-4.j        ]
 [ -4.-9.65685425j]]

前回説明し忘れましたが、データ数は奇数であったりしても離散フーリエ変換は定義されますが、高速アルゴリズムでは 2 の冪乗の場合に一番効率がいいので、特に理由がなければデータ数は 2 の冪乗を選びましょう。

一番目の値は周波数 0 の直流成分を示しています。
6~8 番目の値は 5 番目のデータで折り返すように値が共役複素数になっていることが分かります。
この 5 番目 (n/2+1 番目) のインデックスがサンプリング定理のナイキスト周波数に対応します。
ただし、Python ではデータのインデックスは 0 から始まるので、ナイキスト周波数のインデックス番号は len(y)/2 ですので、間違えないように注意しましょう。

対称な共役複素数のデータはナイキスト周波数を中心とした対称な周波数に対応するデータであり、対称どうしの周波数の波は区別できません(エイリアシング現象)。
よって、離散信号を扱う場合はナイキスト周波数以降の周波数の信号は扱わないのが一般的です。
ナイキスト周波数は、サンプリング間隔 dt[sec] を用いて

15060701.png

となります。
前回の例では dt=0.01[sec] なのでナイキスト周波数は 50[Hz] です。
例えば計測データを 1000[Hz]まで周波数解析したい場合、データを計測する際にサンプリング間隔は 0.0005[sec]より小さく、サンプリング周波数では 2000[Hz] より大きく設定してデータを計測する必要があります。
まあ高いサンプリング周波数で計測しておき、必要に応じてデータ処理の段階でダウンサンプリングすればいいので、余裕があればサンプリング周波数は極力高く設定してデータは計測しましょう。

ただ、サンプリング周波数を高く設定しすぎて、データが大きくなって計測時間が短くなってしまうことがないように注意しましょう。
計測時間 t[sec]の逆数 1/t[Hz] が周波数領域のステップ幅になります。
例えば t=0.5[sec] なら周波数領域のステップ幅は 2[Hz]となってしまい、幅が大きすぎて特徴が抽出できなかったりします。
目的に合わせて適切なサンプリング周波数と計測時間を設定しましょう。

さて、以上を踏まえて前回の例を周波数スケールでプロットします。

# coding: utf-8
import numpy as np
from scipy.fftpack import fft, fftfreq
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)

# 周波数スケール作成
# t/dt/dt/N で作成しても良い
freq = fftfreq(N, dt)

# プロット
# 大きさ、位相
plt.figure(1)
plt.subplot(211)
plt.plot(freq[1:N/2], np.abs(yf)[1:N/2])
plt.axis('tight')
plt.ylabel("amplitude")
plt.subplot(212)
plt.plot(freq[1:N/2], np.degrees(np.angle(yf)[1:N/2]))
plt.axis('tight')
plt.ylim(-180, 180)
plt.xlabel("frequency[Hz]")
plt.ylabel("phase[deg]")

plt.show()

Sorry, your browser does not support SVG.

周波数スケールの作成は SciPy の fftfreq 関数を使うのが簡単です。
図から、ちゃんと 20[Hz] にピークがあるのが分かります。なお、データの表示には直流成分はプロットしないのが一般的です。

次回は fft 処理したデータと振幅の関係を説明していきます。

コメント

Comments powered by Disqus
書籍更新情報
2016-10-21
Pythonによる科学技術計算 基礎編
PDF版の販売を開始しました。
販売ページはこちら

2016-09-09
Pythonによる科学技術計算 基礎編
1.2版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
Kindle ストア、Leanpubで販売中です
Pythonによる科学技術計算 基礎編
PDF版の販売はこちら
同人誌
技術書典(2016/6/25)
Emacs/org-modeのPDF作成術
電子版をBOOTHで販売中です
Emacs/org-modeのPDF作成術
Share