時系列データの周波数解析を実施する際のフーリエ変換で使用する窓関数の処理についてまとめます。
開発環境について
$wmic os get caption >Microsoft Windows 10 Pro $wmic os get osarchitecture >64ビット $python -V Python 3.7.6 $pip list -V scipy 1.4.1
目次
窓関数の役割
測定データはある程度のデータ長を持っており、フーリエ変換するデータを切り出す役割を示します。
また、以下に示すような計測データを周期拡張(コピーして繰り返す)した場合に、本来含まれない波長の波が結果に含まれないようにする役割も示します。
計測データと切り出し区間(青枠)
周期拡張した計測データ
切り出した区間によっては、上記のように信号がつながらなくなり、本来含まれない波長の波が結果に乗ってきてしまう。
サンプルプログラム
SciPyライブラリには多くの窓関数が用意されているので、指定の窓関数を作成できます。
以下サンプルプログラムは、代表的な窓関数であるハニング窓, ハミング窓, ブラックマン窓を示します。
from scipy import signal
from scipy.fftpack import fft, fftshift
import matplotlib.pyplot as plt
import numpy as np
N = 256 # サンプル数
dt = 0.01 # サンプリング間隔
t = np.arange(0, N * dt, dt) # 時間軸
# 窓関数の一例
fw1 = signal.hann(N) # ハニング窓
fw2 = signal.hamming(N) # ハミング窓
fw3 = signal.blackman(N) # ブラックマン窓
# 時間信号
fig = plt.figure(figsize=(12.0, 5.0))
plt.plot(t, fw1, label='hann')
plt.plot(t, fw2, label='hamming')
plt.plot(t, fw3, label='blackman')
plt.xlabel("data number", fontsize=12)
plt.ylabel("Amplitude", fontsize=12)
plt.grid()
leg = plt.legend(loc=1, fontsize=15)
#周波数特性
num = 0
window = [fw1, fw2, fw3]
fig = plt.figure(figsize=(12.0, 5.0))
for i in window:
A = fft(i, 2048) / (len(i)/2.0)
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
if num == 0:
plt.plot(freq, response, label='hann', linewidth = 0.5)
elif num == 1:
plt.plot(freq, response, label='hamming', linewidth = 0.5)
else:
plt.plot(freq, response, label='blackman', linewidth = 0.5)
num += 1
plt.axis([-0.5, 0.5, -120, 0])
plt.title("Frequency response of window")
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.grid()
leg = plt.legend(loc=1, fontsize=15)
実行結果
グラフ1つ目が窓関数、
グラフ2つ目が窓関数の周波数特性を示しています。
どの窓関数とも、データ軍の両端は0に近いため、窓関数を掛けたデータの両端が0近くなり、周期拡張した時に信号がつながらなくなるという問題がなくなります。
ちなみにですが、フーリエ変換で使用するScipyモジュール (scipy.signal.stft)はデフォルトの窓関数がハニング窓になっているようです。
まとめ
今回はScipyを使った窓関数作成方法を紹介しました。
この辺は特に気にせずフーリエ変換、周波数特性を確認していたので、役割含めて良い勉強になりました。
