窓関数の形状

以下の窓関数を考えてみる。
  • 矩形
  • han
  • hamming
  • blackman
  • blackman harris

各々の窓関数の形は以下のようになる。ただ、関数の形をみても性質はよくわからない。

win_all

窓関数の特徴をみるために、窓関数をフーリエ変換してみる。時間領域で対象関数に窓関数をかけるということは、周波数領域で畳み込みをするということであるため、窓関数をフーリエ変換した関数の形状をみることで、周波数の分解能やダイナミックレンジの性質を把握することができる。

窓関数の式

各々の窓関数の式は以下のようになる。

hann及びhammingは以下の式で表され

$a + b \cos(2\pi \frac{t}{N})$

hannは$a=b=1/2$、hamningは$a=0.54,b=0.46$としたもの。

blackmanは、

$a + b \cos(2\pi \frac{t}{N}) + \cos(4\pi \frac{t}{N})$

で$a=0.42,b=0.5,c=0.08$としたもの。

窓関数をフーリエ変換

上記の各々の窓関数をフーリエ変換してみる。

解析的に求めるための準備としてまず$\cos(2\pi n \frac{t}{N}) $のフーリエ変換を求めるておくとよい。

$\cos(2\pi n \frac{t}{N}) $のフーリエ変換は以下のようになる \[W_{n}(\omega) = (-1)^{n}\left(\frac{2 L^2 \omega}{L^2\omega^2 - 4 n^2 \pi^2}\right)\sin\left(\frac{N\omega}{2}\right) \]

展開式はこちらを参照。

例えば、hannやhammingのフーリエ変換は上記を利用して以下のように求められる。

\[a + b \cos(2\pi \frac{t}{N}) \xrightarrow{F} aW_{0} + bW_{1} \]

hannは$a=b=1/2$であるので

\[W(\omega) = \left(\frac{1}{\omega}-\frac{L^2 \omega}{L^2\omega^2 - 4 \pi^2}\right)\sin\left(\frac{N\omega}{2}\right) \]

となる。

グラフの形状は以下のようになり最初に0と交わる$\omega$は$\frac{4\pi}{N}$となる。$\sin$を見ると$\frac{2\pi}{N}$周期で0と交わるとはずだがそうななっていない。

$\frac{1}{\omega}-\frac{L^2 \omega}{L^2\omega^2 - 4 \pi^2}$の影響によりそうなっている。展開式はこちらの式(22)を参照。

各々の窓関数のフーリエ変換の形状をみてみる

矩形

フーリエ変換した式は$W_{0}=\frac{2}{\omega}\sin\left(\frac{N\omega}{2}\right)$となり、最初の節の位置は$\frac{2\pi}{N}$となる。

win_rec
hann

最初の節の位置は$\frac{4\pi}{N}$。矩形と比べると幅が広いので周波数分解能は落ちるがサイドローブが矩形と比べて小さい。 サイドロープがシャープに落ちているのでダイナミックレンジを広くとれる。 ただし、メインローブに近いサイドローブがあまり落ちていないので、 近接した周波数が含まれる信号に対しては向いていない。

win_hann
hamming

最初の節は$\frac{4\pi}{N}$。幅はhannと同じだが細い。hannと比較してメインローブのすぐ脇にでるサイドローブの大きさが 抑えられてるので周波数が近接した複数の信号に向いている。 ただ、サイドローブの減衰が鈍いので広いダイナミックレンジは望めない。

win_hamming
balckman

最初の節の位置は$\frac{6\pi}{N}$。幅がhannやhammingより広いがサイドローブがとても小さい。サイドローブの減衰がよいのでダイナミックレンジが広くとれ、 強い信号と弱い信号が周波数を隔てて存在する場合に有効。 ただ、メインローブの幅が広いので近接した信号の分離ができずに ひとつのピークになってしまう危険性がある。

win_blackman
balckman harris

blackmanと比べて少し太くなるがサイドローブがblackmanより小さくなる。

win_blackmanharris

上記の窓関数を重ね合わせて表示して、サイドローブの大きさやメインローブの幅の比較してみる。

win_all

窓関数の使い分け

窓関数は、サイドローブのレベルが低く、メインローブが細いほうが、望ましい。

  • サイドローブが低い=離れた周波数の信号の影響度が小さい → ダイナミックレンジが大きい
  • メインローブの帯域幅が狭い = 近い周波数に対する影響度が小さい → 周波数分解能が高い
ただ、両方は性質はトレードオフの関係にあるので、用途に合わせてどちらの特徴が重要であるかを見極めて窓関数を選ぶ必要がある。

こちらのサイトには以下のような指針が紹介されている。
  • 短時間の過渡信号 → 矩形
  • ピークの値を正確にしたい → フラット・トップ
  • 非常に近い帯域の信号をはっきり分離したい → カイザー
  • 汎用的、どういった信号か不明 → ハニング、最小 4 項ブラックマン・ハリス
  • とにかくダイナミックレンジを広くしたい → ドルフ・チェビシェフ、ガウス、DPSS でパラメータ調整
まだ、大きなポイント数でFFTでき周波数分解能を確保できる場合は、サイドローブを低く抑える性質をもった窓関数を選ぶ傾向にあると。 参考資料としてNIのこちらのホワイトペーパーがあげられている。
"""
Compare window function
"""
# -*- coding: utf-8 -*-
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import sys, os


N = 51
#window function
w_rec = np.ones(N)
w_hann = signal.hann(N)
w_hamming = signal.hamming(N)
w_blackman = signal.blackman(N)
w_blackmanharris = signal.blackmanharris(N)
w_dic = {'rectangle':w_rec, 'hann':w_hann, 
    'hamming':w_hamming, 'blackman':w_blackman, 'blackman_harris':w_blackmanharris}
#display window function
w_list = []
for k in w_dic:
    if k == 'rectangle':
        continue
    w_list.append(k)
    plt.plot(w_dic[k])
plt.legend(w_list)
plt.ylabel("Amplitude")
plt.xlabel("sample")
plt.title('compare window function shape')
plt.show()
plt.close('all')
#display frequency response of window function
fftpt=8192
f = lambda x : np.fft.fftshift(np.fft.fft(x, fftpt)/np.sqrt(fftpt))
freq_list = 2 * np.pi * np.linspace(-0.5, 0.5, fftpt)
#each 
for k in w_dic:
    fftdt_lin = f(w_dic[k])
    fftdt_log = 20*np.log10(np.abs(fftdt_lin))
    #fftdt = np.abs(fftdt_lin)
    plt.plot(freq_list, fftdt_log)
    plt.legend([k])
    plt.ylabel("Amplitude[dB]")
    plt.xlabel("bin")
    plt.title('frequency response of window function')
    plt.xlim([-1,1])
    plt.ylim([-170,0])
    plt.show()
#all
w_list = []
for k in w_dic:
    w_list.append(k)
    fftdt_lin = f(w_dic[k])
    fftdt_log = 20*np.log10(np.abs(fftdt_lin))
    plt.plot(freq_list, fftdt_log)
plt.legend(w_list)
plt.ylabel("Amplitude[dB]")
plt.xlabel("bin")
plt.title('compare frequency response of window function')
plt.xlim([-1,1])
plt.show()
plt.close('all')