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

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

FFT 処理したデータと振幅の関係

前回はサンプリング定理との関係から、fft 関数から出力されたデータのナイキスト周波数以降のデータは無視することを説明しました。
しかし、正しい周波数解析を行なうにはもう少しデータを処理してやる必要があります。

FFT 処理したデータの大きさを見てみると、元の信号の振幅と全く異っていることがわかります。FFT 処理したデータをちゃんと元の信号と対応させないといけません。
ではどうするのかと言うと、FFT 処理したデータに 1/N を掛け、交流成分については更に 2 倍してやります。
離散フーリエ逆変換の定義から、正規化係数 1/N を掛けることはすぐにわかります。

交流成分について 2 倍するというのは、前回のエイリアシング現象の話しに関連します。
ナイキスト周波数を中心とした対称な周波数の波の区別ができず、それぞれにピークが立っていたように、波の大きさが等分されているので、交流成分については 2 倍する必要があるのです(1 周期がデータ数 N/2 に対応しているからと言った方が正しいかもしれません)。

さて、以上を踏まえて前回の例をちゃんと元の信号の振幅に正規化してやります(前回から周波数を変えています)。
周波数解析では一般的には log スケールで表示することが多いので、ついでに loglog メソッド、semilogx メソッドでのプロットも示します。

# 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 = 15                          # 周波数
t = np.linspace(1, N, N)*dt-dt
y = np.sin(2*np.pi*f*t)

# 離散フーリエ変換+正規化
yf = fft(y)/(N/2)
# 直流成分の振幅を揃える(実用上は不要)
yf[0] = yf[0]/2

# 周波数スケール作成
# 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.angle(yf)[1:N/2]*180/np.pi)
plt.axis('tight')
plt.ylim(-180, 180)
plt.xlabel("frequency[Hz]")
plt.ylabel("phase[deg]")

plt.figure(2)
plt.subplot(211)
plt.loglog(freq[1:N/2], np.abs(yf)[1:N/2])
plt.axis('tight')
plt.ylabel("amplitude")
plt.subplot(212)
plt.semilogx(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.

Figure 1: 振幅の正規化

Sorry, your browser does not support SVG.

Figure 2: log スケールのプロット

なお、そもそも信号から直流成分は大抵カットしてしまいますし、対数表示する場合は使わないので、実用上は直流成分の振幅を揃えることにあまり気を使うことはありません。
今回までで一応 Python での 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