Python NumPy SciPy : FFT 処理による波形整形(スムーザ)

前回 はデジタルフィルタによる波形整形を紹介しました。
デジタルフィルタはリアルタイム処理できるのが利点ですが、位相ずれがあったり、慣れるまで設計が難しいなどの弱点があります。

リアルタイム処理が不要ならば、スムーザと呼ばれる波形整形の方法が使え、位相ずれもなく、簡単に任意の周波数帯域をカットすることができます。
スムーザでは FFT 処理による方法が簡単で、FFT 処理したデータの内カットしたい周波数に対応する値を 0 にし、逆 FFT 処理を行ないます。
以下にスムーザによるローパスフィルタ処理の例を示します。

# coding: utf-8
import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt

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

# FFT 処理と周波数スケールの作成
yf = fftpack.fft(y)/(n/2)
freq = fftpack.fftfreq(n, dt)

# フィルタ処理
# ここではカットオフ周波数以上に対応するデータを 0 にしている
fs = 2                          # カットオフ周波数[Hz]
yf2 = np.copy(yf)
yf2[(freq > fs)] = 0
yf2[(freq < 0)] = 0

# 逆 FFT 処理
# FFT によるフィルタ処理では虚数部が計算されることがあるため
# real 関数が必要(普段は必要ない)
y2 = np.real(fftpack.ifft(yf2)*n)

# rfft 関数を使う場合には yf, y2 は以下のように求まる
# rfft はナイキスト周波数以上の帯域を計算しないため負荷が軽い
# yf = fftpack.rfft(y)/(n/2)
# y2 = fftpack.irfft(yf2)*(n/2)

# プロット
plt.figure(1)
plt.subplot(211)
plt.plot(freq[1:n/2], np.abs(yf[1:n/2]))
plt.ylabel("Amplitude")
plt.axis("tight")
plt.subplot(212)
plt.plot(freq[1:n/2], np.abs(yf2[1:n/2]))
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude")
plt.axis("tight")
plt.savefig('../files/150611Smoother01.svg')

plt.figure(2)
plt.plot(t, y, "b", label="original")
plt.plot(t, y2, "r", linewidth=2, label="filtered")
plt.axis("tight")
plt.legend(loc="upper right")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")

plt.savefig('../files/150611Smoother02.svg')
# plt.show()

Sorry, your browser does not support SVG.

Figure 1: 周波数スペクトル[上段:元データ, 下段:フィルタ処理後]

Sorry, your browser does not support SVG.

Figure 2: フィルタ処理の結果

コメント

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