11. スペクトル解析と窓関数

11.1 離散フーリエ変換によるスペクトル解析

やる夫
離散フーリエ変換のおかげで,時間領域から周波数領域への変換が有限の数列から有限の数列への変換として扱えるようになったわけだお.連続とか無限とかを扱わなくて済むので,実際の信号をコンピュータで解析できるわけだお.

やらない夫
そうだな.まあ解析といってもいろいろあるが,信号がどんな周波数成分を持っているか,同じことだが別の言い方をすると,どのようなスペクトルで構成されているかを調べることができるようになるわけだ.スペクトル解析とか,周波数解析とか呼ぶ.解析の代わりに分析でもいいが,まあどれも同じようなことを指している.

やる夫
実際に計算するには高速フーリエ変換のアルゴリズムを使うのがいいんだお.単純に入力信号を高速フーリエ変換の処理にぶち込めばいいのかお?

やらない夫
うーん,実はそこはいろいろと厄介な問題をはらんでいる.どんな問題かを話す前に,典型的なスペクトル解析の流れを説明しておこうと思う.

まず,解析対象の信号が連続時間信号だった場合は,サンプリングして離散時間信号に変換する必要がある.もちろん最初から離散時間信号を扱う場合もあるだろうからそういうときはそのままでいいんだが,物理的な信号の多くは,例えば音声信号にしろ電気信号にしろ連続時間信号だからな.離散化が必要だ.

やる夫
まあそりゃそうだお.

やらない夫
次に,その離散時間信号を離散フーリエ変換にかけるわけだが,このとき信号の一部の区間を切り出して処理する.$ N$点離散フーリエ変換をする場合は, $ N$時刻分の信号を切り出すわけだ.


\includegraphics[scale=0.5]{fig_win/window.eps}

やる夫
で,計算結果として$ N$点の周波数スペクトルが得られるわけだお.でも,切り出されなかった他の区間はどうするのかお?

やらない夫
普通は,一連のプロセスをえんえんと繰り返す.つまり,入力信号からまた新たに $ N$時刻分の区間を切り出してきてスペクトルを計算する.これを入力が続く限り繰り返すわけだ.結果として「周波数スペクトルのグラフが時間とともに変化する様子」みたいのが得られる.

やる夫
てことは,入力信号が $ N$ 時刻分たまるごとに離散フーリエ変換を走らせるってことかお.入力信号を $ N$ 時刻ずつに区切るってことになるお.


\includegraphics[scale=0.5]{fig_win/multi_window.eps}

やらない夫
そうとは限らないな.例えば,直前に処理した区間とオーバラップした区間を切り出してきて処理したってよいわけだ.


\includegraphics[scale=0.5]{fig_win/multi_window_overlap.eps}

やる夫
なるほど,そうすれば,極端な話 1 時刻ごとにスペクトルを計算することもできるわけだお.

やらない夫
まあ実際にどうするかは,そんなに細かい時間間隔でスペクトルを計算する必要があるかとか,計算が間に合うかとか,そういう辺りを考慮して決めるわけだな.

11.2 アンチエイリアスフィルタ

やる夫
で,問題をはらんでるってのはどういうことだお?

やらない夫
ああ,まず一点目は最初のサンプリングのところだ.前回学んだ通り, サンプリング定理の条件を満たさないとまずいことになる.

やる夫
えーと,エイリアシングが起きて偽の周波数成分が出てきてしまうんだお.それを防ぐには,信号の帯域の2倍以上の周波数でサンプリングしなきゃならないお.

やらない夫
そうだったな.じゃあ,もしサンプリング周波数が固定されてしまっている場合はどうする?

やる夫
ええー,サンプリング周波数の半分の帯域までしか入力信号は持てないんだお.でも現実の音声信号とか電気信号とかがそうなっていることを保証するのは難しいと思うお.

やらない夫
なので,普通はサンプリングの前に低域通過フィルタを入れてやる必要がある.ここでナイキスト周波数以上の成分を除去してしまうんだな.

やる夫
んー,でもサンプリングする前にかけるフィルタってことは…ディジタル処理はできないんだお?

やらない夫
そうだな.ここはアナログのシステムで処理してやる必要がある.この役割のフィルタのことを,アンチエイリアスフィルタと呼ぶ.エイリアシングに抗するフィルタってことだな.元の信号がナイキスト周波数以下に収まっていると確信できる場合以外は,サンプリングの前にこれを置いてやらなくてはならない.


\includegraphics[scale=0.5]{fig_win/block_antialias.eps}

やる夫
でも,まあ,要はローパスフィルタを置けばいいんだお.そんなに問題とも思わないお.

やらない夫
ああ,本当に厄介なのはサンプリング後の話だ.

やる夫
まだあるんかお.

11.3 有限区間の切り出し

やらない夫
離散時間信号から $ N$時刻分を切り出してきて離散フーリエ変換にかけるわけだが,この切り出し方が問題だ.

まず復習だが,離散フーリエ変換ってのはどういう信号に適用できるんだった?

やる夫
えっと,まず当然離散時間信号だお.それから…,あ,そうか,周期的であることも前提だお.時間領域で離散的で周期的だから,周波数領域も離散的で周期的になるんだったお.

やらない夫
そうだな.さて今何らかの入力信号から $ N$ 時刻分の区間を切り出してきた.これを離散フーリエ変換にかけるということは,どういう信号を扱っていることになる?

やる夫
その区間が,全体のうちの1周期分になるような信号だってことだお.だから,その区間の信号をそのままコピペして繰り返した信号を解析しているってことになるお.

やらない夫
その「コピペして繰り返す」ことを周期拡張するという.具体例を考えてみよう.極端に簡単な例として,単一周波数のサイン波を考える.本当は離散時間信号だが,面倒なので図は連続時間信号のようにかくことにしよう.まあ,サンプリング周期がが十分に短いと思ってくれ.

これを,こういう区間で切り出してくれば,何の問題もないわけだ.


\includegraphics[scale=0.5]{fig_win/window_ok.eps}

やる夫
この区間を周期拡張しても元の信号と同じになるわけだお.

やらない夫
じゃあ,こういう区間だったらどうなるか.


\includegraphics[scale=0.5]{fig_win/window_ng.eps}

やる夫
んー,この場合は周期拡張するとこんな風になるお.元の信号とはえらい違いだお.


\includegraphics[scale=0.5]{fig_win/window_ng_repeated.eps}

やらない夫
そう,こんな風に,元の信号と切り出す区間の関係によっては,周期拡張したときにうまくつながらなくて,本来スペクトル解析したい信号とは別のものになってしまう.

やる夫
そっか.切り出す区間の長さ $ N$ は元の信号の周波数に合わせてうまく決めなくちゃならないんだお.

やらない夫
おいおい,そういうことじゃないぞ.

やる夫
えっ,違うのかお?

やらない夫
よく考えてみろ.スペクトル解析をしたいときってのは,元の信号がどんな周波数かわからないから解析したいんだ.わからないのにうまく合わせようがないだろう.

やる夫
あ…そうか,その通りだお.ていうか,そもそも複数の周波数成分が含まれていたら,すべての周波数にうまく合わせるなんてこと自体が無理だお.

やらない夫
というわけで,こういう問題が起きないようなうまい $ N$ を選ぶなんてことは不可能だ.そのことは認めたうえで,周期拡張するときにつなぎ目がおかしくなるのを回避するための手段が必要になる.

やる夫
そんなことできるのかお.

やらない夫
完璧な方法とは言えないが,広く用いられているのが,窓関数を使って切り出しを行い,それから離散フーリエ変換につっこむという方法だ.


\includegraphics[scale=0.5]{fig_win/block_window.eps}

11.4 窓関数とその特性

やらない夫
まず,そもそも一定の長さの区間を切り出してくるっていう操作は,数式としてはこう書ける.

$\displaystyle x_\textnormal{w}[n]$ $\displaystyle = w_\textnormal{r}[n] x[n]$ (11.1)
  where$\displaystyle   w_\textnormal{r}[n] = \left\{ \begin{array}{ll} 1, & n = 0, 1, \cdots, N-1 0, & \text{otherwise} \end{array} \right.$ (11.2)


\includegraphics[scale=0.5]{fig_win/window_rect.eps}

やる夫
$ n = 0, 1, \cdots N-1$ の間だけ 1 であと全部 0 になるような矩形関数をかけるんだお.

やらない夫
これがどうしてまずいのかというと,$ n = 0$$ n = N-1$ のところで $ w_\textnormal{r}[n]$ の値が 0 から 1,1 から 0 に急激に変化するからだ.だから $ x[n]$ を周期拡張したときにつながらなくなってしまう.じゃあどうするか.

やる夫
…なめらかに変化させる,ってことかお.

やらない夫
そういうことになるな.例えば

$\displaystyle w_\textnormal{h}[n] = \left\{ \begin{array}{ll} \displaystyle 0.5...
...i n}{N}}, & n = 0, 1, \cdots, N-1 0, & \text{otherwise} \end{array} \right.$ (11.3)

なんていう形のものがよく使われる.


\includegraphics[scale=0.5]{fig_win/window_hamming.eps}

やる夫
両端がなめらかに小さくなっているお.

やらない夫
こうやって信号の一部を切り出すためにかける関数のことを,一般に窓関数と呼ぶ.最初の $ w_\textnormal{r}[n]$ は矩形窓とか方形窓とか呼ばれる.対して $ w_\textnormal{h}[n]$ の方はハミング窓と呼ばれるものだ.スペクトル解析をするときは,こんな風に両端がなめらかに絞られた窓関数を使って切り出すことで,周期拡張の際につながりがおかしくなることを避けるのが重要だ.

やる夫
うーん,でも,こんなのかけたら波形の形状が変わってしまうお.それは問題にならないのかお?

やらない夫
ああ,大問題だな.

やる夫
あっさり認めたお.

やらない夫
問題だが,矩形窓を使うよりはまだましということだな.もちろんハミング窓を使うのが常に最善というわけじゃなく,他にもいろんな窓関数が提案されている.結局どうやっても何らかの形で波形は歪んでしまって,元の信号のスペクトルとは大なり小なり異なる結果が得られてしまう.そのときどきの目的に応じて都合のよい窓関数を選ぶことが重要だ.

やる夫
んー,一般論過ぎてピンと来ないお.具体的にはどうすればいいのかお.

やらない夫
そうだな.元の信号のスペクトルが,窓関数をかけることでどんな風に変化するかを考えるといい.

まず矩形窓から考えるか. $ w_\textnormal{r}[n]$ の離散時間フーリエ変換 $ W_\textnormal{r}(\omega)$ はどんな形になる?

やる夫
えーと,

$\displaystyle W_\textnormal{r}(\omega)$ $\displaystyle = \sum_{n = -\infty}^{\infty} w_\textnormal{r}[n] e^{-j\omega n}$ (11.4)
  $\displaystyle = \sum_{n = 0}^{N-1} e^{-j\omega n}$ (11.5)
  $\displaystyle = \frac{1 - e^{-j\omega N}}{1 - e^{-j\omega}}$ (11.6)

になるお.でもこれじゃ形はわからないお….

やらない夫
連続時間の矩形関数をフーリエ変換したときのことを思い出すといいぞ.

やる夫
んー,あ,そうか,オイラーの公式で sin に持ち込むんだお.

$\displaystyle W_\textnormal{r}(\omega)$ $\displaystyle = \frac{e^{-j\omega N/2}}{e^{-j\omega/2}}\cdot \frac{e^{j\omega N/2} - e^{-j\omega N/2}}{e^{j\omega/2} - e^{-j\omega/2}}$ (11.7)
  $\displaystyle = e^{-j\frac{\omega(N-1)}{2}}  \frac{\sin \frac{\omega N}{2}}{\sin \frac{\omega}{2}}$ (11.8)

ふーん,同じく矩形関数なのに,離散時間の場合は sinc 関数とは違うものになるんだお.

やらない夫
そうだな.振幅スペクトルの形状はこんな感じで sinc 関数に似てはいるんだが,ちょっと違う.


\includegraphics[scale=0.5]{fig_win/spectre_rect.eps}

sinc 関数は sin の振幅が原点から離れるに従って反比例で減っていくんだったが,これは $ 1 /(\sin{\omega/2})$ で減っていくんだな.しかもちゃんと周期 $ 2\pi$ で周期的になっているのに注意しよう.

やる夫
あ,そうか,離散時間フーリエ変換だから,スペクトルは周期 $ 2\pi$ のはずなんだったお.辻褄が合ってるお.

やらない夫
時間領域で $ w_\textnormal{r}[n]$ をかけるということは,周波数領域では元のスペクトルに今の形の関数 $ W_\textnormal{r}(\omega)$ をたたみこむことになるわけだ.イメージできるか?

やる夫
そろそろ慣れてきたお.元のスペクトルの各点から,その点の高さに合わせた $ W_\textnormal{r}(\omega)$ をぶら下げて,全部重ね合わせるんだお.

やらない夫
その結果どうなるかは, $ W_\textnormal{r}(\omega)$ の形状の各部ごとに考えることができる.まず真ん中の大きな山の部分だな.この部分をよくメインローブと呼ぶ.たたみこまれたときに,メインローブは近傍の値同士を互いと平均化する効果を持つのがイメージできるか?


\includegraphics[scale=0.5]{fig_win/mainlobe.eps}
やる夫
ある点に $ W_\textnormal{r}(\omega)$ をぶら下げたときに,メインローブの分だけその点の値が広がるんだお.全部の点についてそう作用するので,確かに近傍での平均化のように作用しそうだお.

やらない夫
なので,例えばスペクトルに非常に近接した2本のピークがあったとすると,メインローブの影響で両者がくっついてしまったりする.本当は2本のピークがあることを見逃してしまったりしがちなわけだ.

やる夫
なるほど,てことはメインローブはできるだけ細い方がいいんだお.

やらない夫
次にメインローブ以外の背の低い山の部分,これらはサイドローブと呼ばれる.サイドローブは,ある点の値の影響を,ずっと遠くまで引きずる効果を持つことになる.


\includegraphics[scale=0.5]{fig_win/sidelobe.eps}

やる夫
てことはどこかに大きなピークがあると,関係ないところまでちょっとずつ底上げされることになるお.

やらない夫
その結果,スペクトルの小さなピークが,無関係な大きなピークからの影響に埋もれてしまったりする恐れがある.だからサイドローブはできるだけ低いことが望ましい.

やる夫
結局,メインローブが限りなく細くて,サイドローブが限りなく存在しない,単位インパルスみたいな形状のものをたたみこむのが一番ってことだお?

やらない夫
いや,そりゃそうなんだが,周波数領域で単位インパルスになるってことは時間領域では定数 1 だからな.窓をかけることにならんだろ.

やる夫
あ,そっか,ダメだお.

やらない夫
一定の長さの窓関数として機能するものを考えている限り,すべての要件を同時に満たすのは不可能で,トレードオフの関係にある.だから,目的に応じて選択する必要があるわけだな.

で,具体例としてさっきのハミング窓のスペクトルを見てみよう.


\includegraphics[scale=0.5]{fig_win/spectre_hamming.eps}

やる夫
メインローブはむしろやや太くなってるお.でもサイドローブはかなり小さくなってるお.

やらない夫
矩形窓は,実はメインローブの細さについてはかなりよい性能を持っている.ただしサイドローブが高すぎてまずいわけだ.時間領域との対応で考えると,矩形窓は周期拡張の境目で信号が不連続になるのが問題だった.不連続な関数を作るためには,無限に高い周波数成分まで重ね合わせてやらなくてはならない.低周波数のサイン波では滑らかな信号にしかならないからな.つまり広い周波数に渡って偽の成分が発生してしまうことになる.まさにサイドローブの効果なわけだ.

やる夫
ハミング窓みたいに両端を絞ってやるとその効果は抑えられるけど,波形が歪んでスペクトルの形状が変わっちゃうわけだお.メインローブの太さとか,残ったサイドローブがその効果をもたらしていることになるお.

やらない夫
他の代表的な窓関数としては,ハニング窓,ブラックマン窓などがある.興味があれば調べてみるといい.

それから,ちょっと注意点だ.今回ハミング窓は周期 $ N$ の cos と矩形関数を組み合わせたような関数で表した.教科書によっては代わりに周期 $ N-1$ の cos を使っている場合もある.

やる夫
んーと,それは別の窓関数にならないかお.

やらない夫
もちろん厳密には別の関数だし,特性も全く同じじゃないが,まあ流儀の違い程度のものだと思ってもらうのがいいかな.いずれにせよ長さ $ N$ の区間の両端を絞るという目的は同じだ.

swk(at)ic.is.tohoku.ac.jp
2016.01.08