FFT ~ピリオドグラム_パワースペクトル~

ピリオドグラムとパワースペクトルについて

ピリオドグラム法とは、ノンパラメトリック推定法のひとつである。

変動過程における、パワースペクトル密度関数(PSD)の推定方法。

パワースペクトルとの違い(イメージ)は、

無限長のデータXから、矩形波の窓関数W(n)で切り取ったデータX1(n),X2(n),,,,,について、FFTしたものをx1(f),x(f),,,,,とすると、

パワースペクトル：xi(f)^2 の、期待値

$E[xi^2(f)$ ]

ピリオドグラム：Xi(n)の平均のFFT

$E[\bar{xi}(f)$ ]

長さLの信号X(n)のPSDのピリオドグラム推定は、

$P(f) = \frac{1}{LFs} |\sum_{n=0}^{L-1} X(n)exp(-j2 \pi fn/Fs)|^2 \tag{1}$

$P(f)　：　PSDの推定値$

$Fs　：　サンプリング周波数$

$f　：　周波数$

$n = 0,1,.....,L-1$

また、ピリオドグラム法を実装する場合、次のL点での周波数で計算する。

$f = \frac{kFs}{L} \tag{2}$

$k = 0,1,.....,L-1$

アルゴリズム

１．解析したい信号にFFTをかける。この時、信号が実数値ならば、FFTの振幅は、ゼロ周波数に対して対照的になるので、偶数長のでは、最初の(1+(FFTの長さ)/2)点だけが一意である。

２．FFTの振幅を2乗する。ゼロ周波数を除いて、振幅の2乗を 2/(FsL) でスケーリングする。ゼロ周波数は、 1/(FsL) でスケーリングする。

３．「１．で作成した、一意な点の数」と「サンプリング周波数」から周波数の配列を作成(式2)

４．1．と3．をプロット

FFTの性能について

１．スペクトル漏れ

理論上では、無限長の信号XについてPSDを計算するが、

実際には、Xから矩形波の窓関数、W(n)を用いて、有限長の一部分X(n)を切り出し、それのPSDを計算する。

XとW(n)が適当な場合は、X(n)にFFTで計算したPSDは正しくなる。

しかしそうではない場合は、窓関数の端で不連続となり信号にゆがみが生じる。このゆがみをスペクトル漏れという。これは、Xのa番目の周波数成分faが、faを中心に、窓関数をFFTした波形で漏れ出る。

W(n)の時間領域と周波数領域は図１である。

　　　　　　　　　　　図１、引用：

https://ja.wikipedia.org/wiki/%E7%AA%93%E9%96%A2%E6%95%B0

解像度

信号X(n)のスペクトルを区別する能力である。

X(n)が周波数faとfbのサイン波を持っているとき、「スペクトル漏れ」により、互いが互いに影響しあい、ピークが低くなってしまう。

明確にピークを分けるためには、解像度を上げればよい。解像度と周波数fa、fbの関係は、

「メインローブの幅よりも、faとfbの差が大きいことである。」

メインローブの幅は、近似的に、

$Fs/L$

とも表される。よって、

$|fa-fb| > Fs/L$

FFTのバイアス、分散

データ長が無限大に近づくにつれ、矩形波の周波数領域は、デルタ関数に近づく。

詳しくは引用。

参考：

ノンパラメトリック法 - MATLAB & Simulink - MathWorks 日本

修正方法

スペクトル漏れの改善

窓関数を別のものにすることで、サイドローブやスペクトル漏れが変えられる。一方で、矩形波窓関数は、メインローブが鋭く、解像度は理論上最適であったので、解像度に関しては低くなってしまう。

また、窓関数は、X(n)の端を0に収束させることで、滑らかにしている。このことにより、X(n)の平均パワースペクトルに影響を与える。

そこで、窓関数を正規化し、PSDから除算することで、単位減の平均パワーを保つようにしている。

窓関数の応用

窓関数を変えることで、スペクトル漏れが軽減されることが分かったが、同時に、解像度が低下してしまうことが問題である。

次回にマルチテーパー法について書く