正弦波信号のDFTによるスペクトル:窓関数の形状と裾野


正弦波 (サンプリング周波数:1 kHz) 周波数揃え
周波数 250 Hz 250 Hz
Type 1
Type 2
周波数微調整 0 Hz
窓関数
長さ 8 bits (\(N=\) 256 points)
種類 矩形 hann hamming Blackman
離散化方法 \(x=n/(N-1)\) \(x=(n+1/2)/N\)

画面の説明

離散フーリエ変換 (DFT) およびその高速計算アルゴリズムFFTを使った周波数解析の最も単純な例として,正弦波のスペクトルを示します.

ディジタル信号処理の重要な用途の一つは,信号の周波数解析です. 単一の周波数成分しか持たない正弦波を周波数解析したら,きれいな線スペクトルが得られそうな気がしますが,実は,DFT/FFTにより観察されるスペクトルは,そうなりません. もちろんある程度の長さの信号を蓄積してからDFTにかけるときれいなピークは出るのですが,特別な周波数である場合を除いて「裾野」を伴います. この裾野は,信号を有限長で切り出す際の窓関数(ウィンドウ)に由来します.

このページを最初に読み込んだときの初期値は,矩形窓,250 Hz(サンプリング周波数のちょうど1/4)です. このときのスペクトルは,250 Hzのところに鋭いピークが出ていますが,周波数をずらすとどうなるか,確認してみましょう. 初期状態のような鋭いピークではなくなり,ピークの両脇にだらだらと裾野が見えます. メインの周波数スライダを動かすと1 Hz単位で,微調整のスライダを動かすとより細かいステップで,正弦波の周波数が変化します.微調整スライダーで動かせる周波数範囲は,プラスマイナスそれぞれでDFTスペクトルの周波数間隔2つ分です. 各スライダは選択した後矢印キーを使うと,1ステップずつ動せます. この裾野が広がっている状態で,ウィンドウをhann窓,hamming窓,Blackman窓に切り替えると,裾野の高さが下がる様子が確認できます(その代わり,DFTスペクトルのピークの幅は広がります). また,窓関数を長くすると,山の幅が小さくなり,周波数分解能が向上することが確認できます.

解説

一般的な信号の理論上の周波数解析手法は離散時間フーリエ変換 (DTFT) ですが,DTFTは無限長の信号を解析するものですので,現実には実行不可能です. そこで,無限長の信号から有限長の区間を切り出して解析することになります. 有限長区間の切り出しは,時間領域で信号に窓関数をかける操作であり,周波数領域では信号のDTFT(本来のスペクトル)と窓関数スペクトルとの畳み込みになります. 信号が正弦波の場合,正弦波のDTFTは周波数領域におけるインパルスですから,有限長で切り出した正弦波のスペクトルは,中心周波数を正弦波の周波数にシフトした窓関数スペクトル(正確には,プラスマイナスの周波数にシフトした窓関数スペクトルの重ね合わせ)になります. このため,有限長で切り出された正弦波のスペクトルは,周波数シフトされた窓関数のスペクトルと等価です.

次に,このスペクトルをDFT/FFTにより計算することの意味について考えます. DFTは,上記のように有限長で切り出した信号のDTFT(連続周波数)を,周波数領域でサンプリングしたものです. このため,DFT/FFTにより観察される正弦波信号のスペクトルは,周波数シフトした窓関数のスペクトルを周波数領域でサンプリングしたものです. 時間領域で見れば,切り出した長さで信号を周期化し,周期信号と見なして離散時間フーリエ級数 (DTFS) を計算することと等価です. 窓関数は,その振幅がゼロに落ち込む周波数が周期的に現れる構造,つまりサイドローブの山と谷が周期的に現れる性質を持ちますが,周波数軸上でのサンプル点と,窓関数のサイドローブの山と谷の位置関係により,観察される正弦波のDFTスペクトルの裾野の高さが変わってきます. 周波数軸上で窓関数のスペクトルがサンプリングされる様子は,ウィンドウ長を短く(ビット数を小さく)して正弦波の周波数を変化させると,わかりやすいと思います.

発振器から出力した正弦波をスペクトラムアナライザに入力してスペクトルを観察した経験がある人は,ウィンドウの形状を変更するだけでなく,正弦波の周波数をほんの少し変化させるだけで,表示される振幅スペクトルの裾野が大きく変化する様子を見たことがあると思います. 観察されるスペクトルの裾野の変化は矩形窓の場合に最も顕著で,正弦波周波数がサンプリング周波数をウィンドウ長で割った値の整数倍になるとき,サイドローブの谷をサンプリングすることになりますので裾野が低くなります. それらの中間ぐらいの周波数だと,サイドローブの山をサンプリングすることになり,裾野が高くなります.

「周波数揃え」オプションは,2つのスライダで指定した周波数の正弦波をそのまま表示するのではなく,(矩形窓で)裾野のないDFTスペクトルがちょうど出現するよう正弦波周波数にオフセットを加えて表示する機能です. 「Type 1」では,「周波数微調整」のスライダを動かすとどこかで裾野のないDFTスペクトルが出現するように微調整します. 「Type 2」では,DTFTスペクトルの山の幅単位で正弦波周波数を調整し,「周波数微調整」スライダが中央の時に裾野のないDFTスペクトルが観察されるよう,正弦波の周波数にオフセットを加えています. このとき,メインの周波数スライダを1 Hz刻みで動かしても,周波数自動調整により表示される周波数が変化しない場合があります.さらにスライダを動かしていくと,サイドローブの谷がとDFTの周波数サンプリングの位置関係を保ったまま,隣の山が高くなります(このとき,DFTスペクトルの裾野の高さはほぼ維持されます). 「なし」の場合は,このような調整を行いませんので,「周波数微調整」スライダを動かしても裾野のないDFTスペクトルは必ずしも現れません.

窓関数の「離散化方法」のオプションは,連続時間の窓関数を離散時間化する際のサンプリング方法に2つの流儀が存在することに対応するものです. 方法 (a) は,\(0 \le x < 1\) の範囲に対して \(x = n / (N-1)\) という対応関係でサンプリングし,方法 (b) は,\(x = (n + 1/2) / N\) という対応関係でサンプリングします. 窓関数の振る舞いとして両者の間にはほとんど差がなく,特に矩形窓の場合は全く同一なのですが,コサイン窓やBlackman窓においては,裾野のないDFTスペクトルが観察されるかどうかという点で違いが見られます. これは,DFTにおける周波数軸でのサンプリング位置が,DTFTスペクトルのサイドローブの谷にちょうど一致するかどうかという話です. もっとも,これらの窓関数ではもともとサイドローブの値は小さく,かつその谷の正確な周波数でサンプリングするか少しだけずれた周波数でサンプリングするかだけの違いですので,実用上の違いはほぼないと言って差し支えないでしょう.

ちなみに,このグラフの描画に,DFT/FFTは使用していません.窓関数の振幅スペクトルを周波数の関数として数式で導出し,プラスマイナスの正弦波周波数だけシフトした窓関数振幅スペクトルの合計をプロットしています.また,-100 dB以下の値は-100 dBとして表示しています.(クリップすると表示がワンテンポ遅くなるので...)

参考:窓関数とそのスペクトルの形状

窓関数

名称 \(0 \le x < 1\) \(x < 0, 1 \le x\)
矩形窓 \(w(x) = 1\) \(w(x)=0\)
コサイン窓
(hann窓,hamming窓)
\(w(x) = A - (1-A) \cos(2\pi x)\)
hann窓: \(A = 0.5\)
hamming窓: \(A = 0.54\)
Blackman窓 \(w(x) = 0.42 - 0.5 \cos(2\pi x) + 0.08 \cos(4\pi x)\)

窓関数サンプリング

方法 変数対応関係 コサイン窓・Blackman窓スペクトル
周波数シフト量(下記参照)
(a) \(\displaystyle x = \frac{n}{N-1}\) \(\displaystyle \alpha = \frac{2\pi}{N-1}\)
(b) \(\displaystyle x = \frac{n + 1/2}{N}\) \(\displaystyle \alpha = \frac{2\pi}{N}\)

窓関数スペクトル

\(\displaystyle S(\omega) = \frac{\sin(\omega N/2)}{\sin(\omega/2)} \)
として,
矩形窓
\(\displaystyle W_{\rm squ}(e^{j\omega}) = S(\omega) e^{-j\omega(N-1)/2} \)
コサイン窓 (hann窓,hamming窓)
\(\displaystyle W_{\rm cos}(e^{j\omega}) = \Bigl( A S(\omega) + \frac{1-A}{2}(S(\omega+\alpha) + S(\omega-\alpha) \Bigr) e^{-j\omega(N-1)/2} \)
Blackman窓
\(\displaystyle W_{\rm blk}(e^{j\omega}) = \Bigl( 0.42 S(\omega) + 0.25 \bigl(S(\omega+\alpha) + S(\omega-\alpha)\bigr) + 0.04 \bigl(S(\omega+2\alpha) + S(\omega-2\alpha)\bigr) \Bigr) e^{-j\omega(N-1)/2} \)


変更履歴:
2020/06/21 初版
2020/08/02 htmlとJavaScriptをファイル分離
2021/06/16 正弦波周波数範囲を拡大,微調整スライダー追加
2021/07/19 周波数揃え,窓関数サンプリング方式設定追加
2022/08/04 細かい文言の修正
2023/03/23 URL変更

作成:2020年6月21日