Python NumPy SciPy : 窓関数による前処理

前回 までで fft 関数の基本的な使い方を説明しました。 しかし周波数解析を行うには、窓処理と呼ばれる前処理が大抵必要となります。

測定データは N 点の長さの有限区間 \((0, N)\) で定義されていますが、フーリエ変換ではこれを無限区間 \((-\infty, \infty)\) で定義される三角関数の重なりで近似しようというわけです。 そのため、周波数解析では図 2 のように測定データを周期 N の周期信号の 1 周期と仮定しています。 しかし、図 2 のような 0 番目と N 番目の値と傾きが同じでないデータをフーリエ変換すると、本来含まれていない波長の波が解析結果に出てきて本来の周波数特性がぼやけてしまいます。

この問題に対処するため、測定データに窓関数を掛ける窓処理を行ないます。 この処理により図 3 のような滑かに繋る周期関数に波形が変換され、より周波数特性がはっきりと確認できるようになります。 当然このような波形整形をしているため、窓処理をしても多少は周波数特性に歪みが生じます(周波数解析ではどのようにしても特性に歪みが出るのは避けられません)。 窓関数の形はいくつも考案されており、注目する波の帯域やどの程度高い周波数分解能が必要かなど、目的によって適宜窓関数を選択します。

NumPy、SciPy には多くの窓関数が用意されていますが、SciPy の方が種類が豊富です。 以下のように、所定の関数にデータ点数を入力すれば、そのデータ点数で窓を作成することができます。

from scipy import signal
import matplotlib.pyplot as plt

# 時系列のサンプルデータ作成
N = 512                         # データ数

# 窓関数の一例
w1 = signal.hann(N)             # ハニング窓
w2 = signal.hamming(N)          # ハミング窓
w3 = signal.blackman(N)         # ブラックマン窓

# プロット
plt.figure()
plt.plot(w1, "r", w2, "b", w3, "g")
plt.axis("tight")
plt.ylabel("Amplitude")
plt.xlabel("data number")
plt.legend(["hanning", "hamming", "blackman"])
# plt.show()