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から除算することで、単位減の平均パワーを保つようにしている。

 

窓関数の応用

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

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

そもそもウェーブレット変換とは？〜その２〜

１．前回の内容(この記事は、「そもそもウェーブレット変換とは？〜その１〜」を受けて書いています。)

「そもそもウェーブレット変換とは？〜その１〜」では、ウェーブレット変換を用いようと思った理由と、連続と離散の応用のされ方のおおまかな違いについて説明した。今回から、連続ウェーブレット変換について説明していく。

 

２．連続ウェーブレット変換

 マザーウェーブレットと呼ばれる基本となる関数を拡大、縮小、シフトすることで、信号の周波数と時間の解析を行う方法である。

このマザーウェーブレットには特徴がある。

 

２－１．マザーウェーブレット

ここで、フーリエ変換について少し復習する。フーリエ変換は、信号をサイン波とコサイン波に分解して、解析する方法であった。ある周波数のサイン波を、どれほどの含んでいるかを考えていた。

一方で、ウェーブレット変換は、三角関数ではなく、別の関数(マザーウェーブレット)である。

マザーウェーブレットとそのフーリエ変換を

 

ψ(t)

 

\hat{ψ}(ω)

 

としたとき、以下の条件を満たす。

  

\int_{-∞}^∞ |{ψ(t)}|^2 dt ＜ ∞\tag{1}

 

\int_{-∞}^∞ \frac{|\hat{ψ}(ω)|^2}{ω} dω  ＜ 0\tag{2}

 

\int_{-∞}^∞ ψ(t) dt  = 0\tag{2’}

 

式1は、「ψ(t)は、有限のエネルギーを持つ」ということである。

また、式2は、アドミッシブル条件と言われるものであり、式2’はその代用として使用される場合がある。式2が示していることは、

 

\hat{ψ}(0) = 0

 

である。つまり、ψ(t)はゼロの周波数成分を持たない。

 

 

これはすなわち、

・遠くでは０に近づく　＝　局在している　(∵ 式1の有限値であるという条件)

・おおざっぱにみれば波のような形をしている　(∵式2、2’から、正と負の面積が等しいという条件)

具体的に見ていく。マザーウェーブレットの代表的な例として、Morlet関数を用いる。式は

 

ψ(t) = \frac{1}{2\sqrt{π}σ}\exp(-\frac{t^2}{σ^2})\exp(-it) \tag{3}

 

f:id:nitsuta_s:20181109183633p:plain

 

である。

 

次回は、このマザーウェーブレットを拡大、縮小、シフトするパラメータについて説明する。

そもそもウェーブレット変換とは？〜その1〜

1．はじめに

マウスの脳から取ってきた神経活動について、解析を行っています。

これまでは、フーリエ変換を用いた解析を行っていたが、手詰まり感が出てきた。

そこで、時間的変化にも対応したウェーブレット変換を用いようと思った。

 

ウェーブレット変換は、連続ウェーブレット変換と離散ウェーブレット変換に分けられるが、前者は信号のパターンや相違性の解析など、後者は、データの圧縮やエネルギー解析に用いられる。

ここでは、連続ウェーブレット変換について書こうと思う。

そのうち、離散ウェーブレット変換についても書きたいと思う。

 

 

次回から、連続ウェーブレット変換について説明をしていこうと思う。

pythonでウェーブレット変換(つまづき.1)

・連続ウェーブレット変換とは？

 

・なんで cwt などのプログラミング(離散値)で連続ウェーブレット変換の関数が存在するの？

answer : np.sin とかは三角関数だから(？)

 

・pywt.cwt の戻り値は何を表している？

 

・離散ウェーブレット変換とは？

answer : デジタル信号など離散信号に用いる

 

・pywt.dwt の戻り値は何を表している？

 

pythonでウェーブレット変換

  • ウェーブレット変換

式1．https://wikimedia.org/api/rest_v1/media/math/render/svg/ab86833f90147f112543224b11006d583047f49e(wikipediaから引用)

 

ψ(x)は以下の条件を満たせばよい

 

・　f:id:nitsuta_s:20181107160904p:plain

→局在した関数(xが大きいときは、0である)

→無限の足し合わせが0だから、おおざっぱに見れば、波に見える

 

  • ウェーブレット変換のa,bの意味

ψ(x)が x=b0 に局在している関数であるとすると、ψ((x-b)/a)は x=b0 + b に局在している関数となる。

だから式1．は、x=b0 + b 付近の f(x) の特徴を表している。

一方aは、ψ(x)をフーリエ変換すると、周波数成分がわかる。ψ(x)は三角関数ではないので、一意の周波数成分とはならないが、おおざっぱに波の関数であったので、周波数f0 の成分を持つと考えられる。よって、ψ((x-b)/a)は周波数f0/a の成分を持つ関数であることが考えられる。

 

まとめると、式1.は

・x=b0 + b 付近

・周波数f0/a 付近の成分

を表していることとなる。

 

引用：

d.hatena.ne.jp