FFT用のちぇびしぇふ窓を作ってみよう
（2011/6/19）
自分はあんまり見たことないですけれど、
FFT窓関数として「チェビシェフ窓」なる関数があるらしいので、
チョロリとプログラムを作って、パラメーターを振って、周波数領域/時間領域の特性を見てみよう、
・・・という趣向です。
窓関数というのは、FFT長毎にサンプルを切り出してスペクトル分析をする際に、切り出した両端が不連続になっちゃって、
本来の信号には存在しないはずの不要なスペクトルが発生してしまうのを抑制するための処理です。
具体的には切り出したサンプルに対して窓関数を掛けて両端の振幅を抑え込みます。
【１．参考にしたサイト】
[1]Windowing Functions Improve FFT Results, Part II
Flexible windowing functions let you adjust frequency and leakage effects.
Part I of this article appeared in our June 1998 issue.
(Test and Measurement World）
http://www.tmworld.com/article/325630-Windowing_Functions_Improve_FFT_Results_Part_II.php
[2]Computing Chebyshev Window Sequences (DSP Related)
http://www.dsprelated.com/showarticle/42.php
【２．チェビシェフ窓の設計に必要なパラメータとはなんぞや？】
チェビチェフ窓の設計は、２つのパラメーターからなるようです。
①窓サイズ（N）：
Nサンプル長の窓を作る、という指定です。こりゃ必要ですわね。
②サイドローブの抑制具合（Gとかγとかg色々記載）：
20dBとこのGを掛けた値がサイドローブの抑制値になります。
例えばG=2を指定したら、サイドローブは-40dBに抑制される・・・って具合。
【作りかた】
[1][2]のサイトに記述がある通りですが（特に[2]はとても親切で詳細！）大体こんな感じです。↓：
①まず周波数領域の応答を求める。
W(m) = (-1)m•cosh{N•cosh-1[A(m)]}, when |A(m)| > 1,
or
W(m) = (-1)m•cos{N•cos-1[A(m)]}, when |A(m)| ≤ 1.
※ここで
A(m) = |α•cos(πm/N)|, (0 ≤ m ≤ (N-1).
α = cosh[cosh-1(10γ)/N].
②次に上で求めた周波数領域の応答W(m)に対して逆DFTをかけて、時間領域にしてやる。
これを小文字のw(m)とすると、w(m)がチェビチェフ窓関数です。（まだ正規化されてないけど。）
③上で求めたw(m)を正規化（なんのことはない、maxの値が1となるように割り算してるだけです。）
これにてチェビシェフ窓の出来上がり！・・・のはず。
【はたして上手くできるでしょうか】
VBAで上の処理を作ってみました。
コアの部分（上の②の処理）はこんだけ。
③のDFTの処理は窓関数は使わないで実行します。
まあその・・やっつけなんで許したってや～。
VBの関数でとかcoshとかacosh双極関数が使えず、Excelのシート関数を参照するという謎の挙動もあるヨ。
ドキドキしつつ動かした結果！
| 周波数ドメインの表示。 |
左をDFTして時間領域にヘンカンした結果。 |
ふむふむ。ちゃんとパラメータで与えたとおり、サイドローブのリップル（波うち）は抑制されているぜ！
・・・あたりまえか。
サイドローブを抑制しようとすると、メインローブの幅が広がって周波数分解能が落ちちゃいます。
世の中そうウマイ話はありません。 |
G=2（-40dBに抑制）したときなんかはピョコンと両端上がって塩梅悪そうです。
他はそれなりにいけてるように見えます。
|
[2]のサイトのグラフなんかとも比較してみてください。
http://www.tmworld.com/article/325630-Windowing_Functions_Improve_FFT_Results_Part_II.php
＜感想＞
ということで、FFTのチェビシェフ窓は以上で求まるんでないのか？
上の処理以上のなにかあるんだらうか？・・っていうのが自分の疑問です。
例えば、周波数領域W(m)の状態なんかで、リップルの手前でDFTを打ち切ったり、
・・・などといった特別な処理ってないんじゃあないか、と思うのですが、良く分からないのでこれにて！
（詳しい方教えてくださいね。）
以上
TOPへ戻る