SciPy で使用可能な窓関数

SciPy で使える窓関数についてまとめていきます。 SciPy には表 1 の窓関数が用意されています。 これらの関数は get_window 関数から呼び出しても使えます。 表の数値は参考まで(窓関数のサンプル点数 51、FFT の長さ 65536 で計算)。

Table 1: SciPy の窓関数
関数 サイドロープ メインローブ メインローブ スカロップロス NENBW
    レベル[dB] -3dB 帯域幅[bin] -6dB 帯域幅[bin] 最大値[dB] [bin]
barthann 修正バートレット・ハン -35.89 1.42 1.98 -1.46 1.49
bartlett バートレット -26.43 1.30 1.80 -1.75 1.36
blackman ブラックマン -58.11 1.67 2.34 -1.05 1.76
blackmanharris 最小 4 項ブラックマン・ハリス -92.10 1.93 2.72 -0.79 2.04
bohman ボーマン -46.00 1.73 2.42 -0.98 1.82
boxcar 矩形 -13.25 0.88 1.20 -3.91 1.00
chebwin(N=51, at=100) ドルフ・チェビシェフ -100.0 1.85 2.60 -0.86 1.96
cosine コサイン -23.02 1.19 1.64 -2.10 1.23
flattop フラット・トップ -67.74 3.80 4.68 -0.02 3.84
gaussian(N=51, sig=7) ガウス -75.49 1.93 2.73 -0.80 2.06
general_gaussian(N=51, p=1.5 , sig=7) 一般化ガウス -24.89 2.46 3.41 -0.48 2.57
hamming ハミング -42.31 1.32 1.84 -1.70 1.38
hann ハン(ハニング) -31.47 1.47 2.04 -1.37 1.53
kaiser(N=51, beta=14) カイザー -105.73 2.08 2.93 -0.68 2.20
nuttall ナットールの定義による最小 4 項ブラックマン・ハリス -92.28 1.91 2.68 -0.82 2.02
parzen パルザン -53.04 1.82 2.55 -0.90 1.92
stepian(N=51, width=0.3) 離散扁長回転楕円体配列(DPSS) -91.69 1.87 2.63 -0.84 1.98
triang 三角 -26.44 1.25 1.73 -1.89 1.31

窓関数はサイドローブのレベルが低く、メインローブの帯域幅が狭いことが望ましいです。

  • サイドローブのレベルが低い = ダイナミックレンジが広い = 遠く離れた周波数に対する影響が小さい
  • メインローブの帯域幅が狭い = 周波数分解能が良い = 近接する周波数に対する影響が小さい

しかし、この 2 つはトレード・オフの関係にあり、場合によって適宜窓関数を選択することが必要です。

個人的には以下の指針で窓関数を選んでいます(あくまで参考まで)。 選定については ここ が非常に参考になります。

  • 短時間の過渡信号 → 矩形
  • ピークの値を正確にしたい → フラット・トップ
  • 非常に近い帯域の信号をはっきり分離したい → カイザー
  • 汎用的、どういった信号か不明 → ハニング、最小 4 項ブラックマン・ハリス
  • とにかくダイナミックレンジを広くしたい → ドルフ・チェビシェフ、ガウス、DPSS でパラメータ調整

近年はハードウェアが進歩して大きいデータの FFT も可能になり、高いサンプリング周波数でも周波数ステップ幅を小さくできるようになったため、多少メインローブの帯域幅が広くてもサイドローブレベルの低い窓関数が使われるようになって来ています。

機械系だと構造物の振動なんて狭帯域の信号しか扱わないですし、近接するモードが分離できれば良かったりするので、ハニングやカイザーあたりをよく使いますかね。 電気系で GHz とかの帯域の信号を扱う人はダイナミックレンジが広い窓関数やフラット・トップをよく使うそうです。 まあ、どんな場合でも窓関数はいろいろ試して結果を比較することが必要です。

SciPy の窓関数とその周波数特性をプロットすると以下のようになります。

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

N = 51
w_barthann = signal.barthann(N)
w_bartlett = signal.bartlett(N)
w_blackman = signal.blackman(N)
w_blackmanharris = signal.blackmanharris(N)
w_bohman = signal.bohman(N)
w_boxcar = signal.boxcar(N)
w_chebwin = signal.chebwin(N, at=100)
w_cosine = signal.cosine(N)
w_flattop = signal.flattop(N)
w_gaussian = signal.gaussian(N, std=7)
w_general_gaussian = signal.general_gaussian(N, p=1.5, sig=7)
w_hamming = signal.hamming(N)
w_hann = signal.hann(N)
w_kaiser = signal.kaiser(N, beta=14)
w_nuttall = signal.nuttall(N)
w_parzen = signal.parzen(N)
w_slepian = signal.slepian(N, width=0.3)
w_triang = signal.triang(N)

plt.figure(figsize=(6, 18))
plt.subplot(411)
plt.plot(w_barthann)
plt.plot(w_bartlett)
plt.plot(w_blackman)
plt.plot(w_blackmanharris)
plt.plot(w_bohman)
plt.ylim(-0.2, 1.2)
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.legend(["barthann", "bartlett", "blackman", "blackmanharris", "bohman"],
	   bbox_to_anchor=(1.05, 1),
	   loc='upper left', borderaxespad=0)

plt.subplot(412)
plt.plot(w_boxcar)
plt.plot(w_chebwin)
plt.plot(w_cosine)
plt.plot(w_flattop)
plt.ylim(-0.2, 1.2)
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.legend(["boxcar", "chebwin", "cosine", "flattop"],
	   bbox_to_anchor=(1.05, 1),
	   loc='upper left', borderaxespad=0)

plt.subplot(413)
plt.plot(w_gaussian)
plt.plot(w_general_gaussian)
plt.plot(w_hamming)
plt.plot(w_hann)
plt.plot(w_kaiser)
plt.ylim(-0.2, 1.2)
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.legend(["gaussian", "general_gaussian", "hamming", "hann", "kaiser"],
	   bbox_to_anchor=(1.05, 1),
	   loc='upper left', borderaxespad=0)

plt.subplot(414)
plt.plot(w_nuttall)
plt.plot(w_parzen)
plt.plot(w_slepian)
plt.plot(w_triang)
plt.ylim(-0.2, 1.2)
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.legend(["nuttall", "parzen", "slepian", "triang"],
	   bbox_to_anchor=(1.05, 1),
	   loc='upper left', borderaxespad=0)

15071301.png

# coding: utf-8
def plotFR(w):
    A = fft(w, 2048) / (len(w)/2.0)
    freq = np.linspace(-0.5, 0.5, len(A))
    response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
    plt.plot(freq, response)

plt.figure(figsize=(6, 18))    
plt.subplot(411)  
plotFR(w_barthann)
plotFR(w_bartlett)
plotFR(w_blackman)
plotFR(w_blackmanharris)
plotFR(w_bohman)
plt.axis([-0.5, 0.5, -120, 0])
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.legend(["barthann", "bartlett", "blackman", "blackmanharris", "bohman"],
	   bbox_to_anchor=(1.05, 1),
	   loc='upper left', borderaxespad=0)

plt.subplot(412)
plotFR(w_boxcar)
plotFR(w_chebwin)
plotFR(w_cosine)
plotFR(w_flattop)
plt.axis([-0.5, 0.5, -120, 0])
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.legend(["boxcar", "chebwin", "cosine", "flattop"],
	   bbox_to_anchor=(1.05, 1),
	   loc='upper left', borderaxespad=0)

plt.subplot(413)
plotFR(w_gaussian)
plotFR(w_general_gaussian)
plotFR(w_hamming)
plotFR(w_hann)
plotFR(w_kaiser)
plt.axis([-0.5, 0.5, -120, 0])
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.legend(["gaussian", "general_gaussian", "hamming", "hann", "kaiser"],
	   bbox_to_anchor=(1.05, 1),
	   loc='upper left', borderaxespad=0)

plt.subplot(414)
plotFR(w_nuttall)
plotFR(w_parzen)
plotFR(w_slepian)
plotFR(w_triang)
plt.axis([-0.5, 0.5, -120, 0])
plt.ylabel("Normalized magnitude [dB]")
plt.xlabel("Normalized frequency [cycles per sample]")
plt.legend(["nuttall", "parzen", "slepian", "triang"],
	   bbox_to_anchor=(1.05, 1),
	   loc='upper left', borderaxespad=0)

15071302.png

コメント

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