【自分用メモ】時間周波数解析手法のまとめ

f:id:s0sem0y:20170608011337p:plain

波形の一部を切り取りフーリエ変換を行います。

切り取る場所を少しずつずらしながらフーリエ変換することで、スペクトルの時間的な変化を追うことが可能になる時間周波数解析の最も古典的な方法です。

英語ではShort-Time Fourier Transform(STFT)と略記します。

変換は波形 $x(t)$ 時間窓 $w(t)$ として

$\displaystyle STFT(t,ω)=\int_{-∞}^{∞}x(τ)w(τ-t)e^{-iωτ}dτ$

と、 $x(t)e^{-iωt}$ に $w(t)$ を畳み込み積分した形になります。

フーリエ変換では時間と周波数の間にトレードオフが存在し、 $w(t)$ の定義域に応じて $ω$ の方の誤差が決定されます。

短時間フーリエ変換では１つの $w(t)$ を使うため、基本的にどの周波数帯域に置いても誤差は一定となり、 $ω=1$ に対する誤差 $δω$ も $ω=100$ に対する誤差 $δω$ も同じ大きさになっており、低周波域に対する周波数分解能が相対的に悪くなる傾向にあります。

短時間フーリエ変換では時間窓が固定されていましたが、ウェーブレット変換では高周波領域は時間分解能を大きく（周波数分解能を小さく）して、低周波領域では時間分解能を小さく（周波数分解能を大きく）することが可能になります。

$\displaystyle T_{wave}f(a,b)=\frac{1}{\sqrt{a}}\int x(t)Ψ \left( \frac{t-b}{a} \right)$

係数 $\frac{1}{\sqrt{a}}$ はエネルギー（振幅の二乗の総和）を保存するための係数で、周波数という概念からすれば特に関係のないものです。

$Ψ(t)$ がマザーウェーブレットであり、 $Ψ \left( \frac{t-b}{a} \right)$ はウェーブレット基底と呼ばれ、色々な $a$ の値を使うことで、ウェーブレット基底はマザーウェーブレットを時間方向に伸縮した形をすることになります。

要するに周波数を高くする代わりに、波形を時間方向に押しつぶして「小さな波形」として扱うのです。これによって波形の定義域も $a$ の増加によって狭まるため、結果的に時間分解能と周波数分解能の関係を、周波数帯域毎に変更する役割を担うことができます。

これは短時間フーリエ変換の問題を、短時間フーリエ変換の時間窓を直接いじることで解決する方法です。英語では「S transform」と言います。

$\displaystyle S(t,f)=\int_{-∞}^{∞}|f|x(τ)e^{-π(t-τ)^2f^2}e^{-i2πfτ}dτ$

$\displaystyle f=\frac{ω}{2π}$

短時間フーリエ変換において、時間窓をガウス関数にして、ガウス関数の標準偏差を $f$ にしていることに相当します。すなわち周波数 $f$ に応じて時間窓の幅が変更され、結果として周波数分解能と時間分解能を調整することになるのは同じでしょう。

これまでの時間周波数解析は、例えば正弦波、あるいはマザーウェーブレットなどを基底として、その基底の和という形式で波形を表現していました。どの基底がどれくらい含まれているかを知ることで解析をするという方法です。

（Empirical Mode Decomposition:EMD）とは、「基底」をヒューリスティックな方法によって獲得して、その基底により分解を行うことです。すなわち、基底を予め決めておくのではなく、基底すらも波形から見つけてしまおうと言う方法です。

波形 $x(t)$ に対して以下の操作によってIntrinsic mode functinon(IMF)（こいつが基底に相当する）を複数抽出していきます。

入力信号 $x(t)$ から1つのIMF $y(t)$ を取り出すアルゴリズムは以下です。

入力信号 $x(t)$ の全ての極値を検出
極大点と極小点をそれぞれ補間し,上側包絡線 $emax(t)$ と下側包絡線 $emin(t)$ を獲得
$emax(t)$ と $emin(t)$ の局所平均 $m(t) = \frac{emin(t) + emax(t) }{2}$ を算出
入力と局所平均の差分 $y(t) = x(t) - m(t)$ を入力とみなして 1~4 を繰り返す
$y(t)$ が閾値以下になった時点で $y(t)$ をIMFと見なす

これで1つのIMF $y(t)$ を取り出せたら、更にこの結果のもとでIMFを取り出していきます。
$x(t)$ と $y(t)$ の残差 $x_1(t) = x(t) - y(t)$ を新たな入力信号として 1~6 を繰り返す
上記を繰り返し、極値が１つだけの $x_n(t)$ が得られたら終了（全てのIMFが抽出終了）