HELLO CYBERNETICS

深層学習、機械学習、強化学習、信号処理、制御工学、量子計算などをテーマに扱っていきます

ウェーブレット変換の基本

 

 

follow us in feedly

 

 はじめに

ウェーブレット変換は、生まれながらにしての時間周波数解析手法であり、短時間フーリエ変換などに変わる解析手法として使われています。今回はウェーブレット解析の基本を見ながら、短時間フーリエ変換に対するメリットなどを見て行きたいと思います。

 

フーリエ変換や短時間フーリエ変換に関しては以下の記事で解説しています。

s0sem0y.hatenablog.com

s0sem0y.hatenablog.com

 

ウェーブレット変換を考える動機

周波数解析

まず、フーリエ変換より、ある波形を三角関数の和に展開することができます。様々な周波数の三角関数の線型結合として波形を表現することによって、その波形に、とある周波数がどれだけ含まれているのかを定量的に分析することが可能となりました。

f:id:s0sem0y:20161208002121p:plain

f:id:s0sem0y:20170106163055p:plain

上の図の波形は、フーリエ解析をすると下の図のようにいくつかの波形が重なりあってできたものだと分かる。

 

例えば、普段からとある波形を周波数解析をしておいて、その構成を把握しておけば、何か異常が生じた時にその波形を周波数解析して見比べることで、どの周波数が異常な波形として現れているのかを知ることができます。

その異常な周波数を把握することができれば、バンドカットフィルターなどを通すという単純な方法で異常な波形をいつでも削除することができるようになるかもしれません。

 

時間周波数解析

周波数解析の限界

ある波形が時間的に、含んでいる三角関数を変化させていく場合は周波数解析は無力です。

なぜなら周波数解析は、ある時間計測した波形に、どれくらいの三角関数が平均的に含まれているのかを知るだけであり、いつそのような構成であったのかを知ることができないからです。

 

例えば、単純に時刻0〜10秒までが10Hzの波形であり、10〜20秒までが20Hzの波形で構成されている場合に、0〜20秒まで観測した波形を周波数解析すれば、出てくる情報はあくまで、10Hzと20Hzの信号が含まれていたということだけなのです。

このような波形は、波形を見た時点で周波数解析をするまでもありませんが、普通に自然界に存在する波形はもっと複雑で、見ただけでは分かりません。どこのタイミングで波形の構成が変わったのかなど分からないため、解析する時間区間を区切るタイミングなど知るよしも無いのです。

 

短時間フーリエ変換

そこで、時間間隔を非常に短くして、例えば1秒毎に周波数解析を行っていくという方法を考えましょう。この場合、1秒間という短い間隔ならば波形の構成は大きく変化していないだろうと目論むことになります。そして次の1秒で波形の構成が変化していたとすれば、1秒前に比べ周波数解析で得られる波形の構成が変わり、時間的に変化する構成の変化を捉えることができます。

 

これが短時間フーリエ変換の素朴な発想です。

 

0〜1秒のデータに周波数解析

1〜2秒のデータに周波数解析

 

としてもいいのですが、1.5〜2.5秒の間に大きな変化があるかもしれません。

その場合は

 

0〜1秒のデータに周波数解析

0.1〜1.1秒のデータに周波数解析

 

と、解析する時間幅を1秒のまま、これをスライド移動して解析するという方法が取られます。

通常こちらを短時間フーリエ変換と言います。

 

不確定性関係

フーリエ変換には不確定性関係という厄介な関係があります。

時間周波数解析においては、時間と周波数の間にこの不確定性関係が出てきます。

これは簡単に言えば、時間幅を短くすると周波数の情報がハッキリと分からなくなり、逆に周波数の情報をきめ細かく知ろうとすると、時間幅を大きく取らなければならないという関係です。

 

時間幅を1秒と短くとった場合には、高々1秒の観測で得られる波形から周波数の構成を知ろうというので無理が生じてきます。

かと言って、周波数の構成を正しく知るために時間幅を5秒などにしてしまえば、周波数の構成の時間的変化を細かく知ることができなくなってしまいます。

 

しかし、これは数学的に内在するものであり、解決の方法がありません。これを認めて解析をするしか無いのです。

 

時間幅をΔt秒で解析すると決めたのならば、周波数の構成を知る際に、どれほどの誤差Δωが生じるのかが決定されます。要するにΔt程度の波形の情報では、ωの波形がただ1つ含まれているだけだとしても、解析するとω±Δωの波形が含まれているように見えてしまうのです。

 

当然実際にはωという波形がただ1つ含まれているわけではなく、ω_1ω_2などとたくさんの波形が含まれているはずです。その全てに対して誤差が乗っかってくると、もはやどこが真の周波数であったのかなどわからなくなってしまうのです。

 

時間と周波数のどちらを立てるのか

今から解析する波形は1〜60Hzまでの波形が含まれていると見込まれているとしましょう。このときにその波形を時間幅Δtで時間周波数解析した場合には、すべての周波数に対してΔωの誤差が生じうるということになります。

 

これには直感的に問題があるように感じます。

例えばΔωが2Hzくらいの誤差だとしましょう。

そうした場合

 

実際に4Hzの波形が含まれていたとして、解析すると

 

4±2 Hz

 

の信号が含まれているように見えるのです。

一方で50Hzの信号が含まれていた場合は

 

50±2 Hz

 

の信号が含まれているように見えます。

他にもいろいろな周波数が含まれていることが見込まれますが、単純に考えて、上の2つでは誤差の割合が全く違っています。

 

4を6と見誤るのと、50を52と見誤るのでは全く状況が異なるのです。

 

しかし短時間フーリエ変換ではこういうことが起きてしまうのです。

 

やはり、4Hzなどの低周波領域では、時間と周波数に関して周波数の精度を高めたいですし、高周波領域では周波数の精度は荒っぽくする代わりに、時間的な変化を捉えたいなどと使い分けが出来たほうが嬉しいはずです。

 

これを考え始めたのがウェーブレット変換です。

 

ウェーブレット変換

ウェーブレット変換のモチベーション1

フーリエ変換の考えでは、どんな波形もいろいろな周波数の\sin\cosを重ねあわせることで表現できるという観点に基づいて変換を行っていました。しかしフーリエ変換には不確定性関係によってピンぼけするという問題があります。

 

これは時間と周波数のどちらを立てるべきかという問題として捉えることができ、フーリエ変換は時間幅を固定するために、周波数のピンぼけ度合いが調整できないという問題があると考えることができます。

 

これに対しウェーブレット変換では、低周波領域では周波数の分解能を優先し、高周波領域では時間の分解能を優先するように変換を行います。

 

ウェーブレット変換のモチベーション2

フーリエ変換では色々な周波数の正弦波を用いましたが、ウェーブレット変換では好きな形の波形を使うことできます。任意の波形を正弦波の和に分解するのではなく、任意の波形を適当な波形の和に分解するのです。

 

ここで適当な波形Ψ(t)というのはただ1つだけ選ばれて、それをマザーウェーブレットΨ(t)と呼びます。仮にある波形に特徴的な波形が混在しているとして、その特徴的な波形の形にある程度目処が付いているのならば、その波形をマザーウェーブレットΨ(t)として準備しておくことで、いとも簡単にその特徴的な波形が生じている時間を特定することができます。

 

またマザーウェーブレットの形は、Ψ \left( \frac{t-b}{a} \right)とすることで、a,bを調節すれば伸縮することができます。

 

ウェーブレット変換とは、任意(ある条件を満たしている必要があるが)の波形を色々な値のaを使ったウェーブレット基底Ψ \left( \frac{t-b}{a} \right)の和で表現する手法なのです。

 

 

ここでbの方は時間的な移動を表しています。

 

ハールウェーブレット

ハールウェーブレットはウェーブレット変換がなにを達成しようとしているかを知るのにはとても良い題材です。ハールウェーブレットは以下の関数です。

 

f:id:s0sem0y:20170617015042p:plain

 

形としては以下のようになっています。(線が無い部分は0と見なせばいい)

f:id:s0sem0y:20170617015350p:plain

マザーウェーブレットから生まれる基底たち

マザーウェーブレットΨ(t)について

 

\displaystyle \frac{1}{\sqrt{a}}Ψ(\frac{t-b}{a})

 

a,bを用いた平行移動と拡大縮小を行うと以下のようになります。

 

 

a = 2,b = 0

f:id:s0sem0y:20170617015848p:plain

 

a = 2,b = 1

 

f:id:s0sem0y:20170617020301p:plain

 

aの値が増えると、元々の波形は横方向にも縦方向にも圧縮されます。より細かい波形を表現する素となるのです(このaが周波数の調整に相当する)。そして、bを変更することで波形が値を持つ時間区間を移動できます(つまり時間の表現をbで行う)。

 

基底を足し算していろんな波形を作る 

マザーウェーブレットはa=1,b=0に相当しており、今、a=2,b=0と、a=2,b=1のパターンも見てきました。これら3つの波形を以下のように係数を付けて足し算すれば色々な波形x(t)を表現できるということです。

 

\displaystyle x(t) = w_{1,0}\frac{1}{\sqrt{1}}Ψ(\frac{t-0}{1}) + w_{2,0}\frac{1}{\sqrt{2}}Ψ(\frac{t-0}{2}) + w_{2,1}\frac{1}{\sqrt{2}}Ψ(\frac{t-1}{2})

 

仮にw_{a,b}が全て1だとして足し算すれば

 

f:id:s0sem0y:20170617022126p:plain

 

のような形になります。これを係数w_{a,b}をいろいろ調整することで、異なる波形にできるのは言うまでもありません。更にaのパターンを更に増やせば、もっと細かい波形も表現できるようになるでしょう。bをずらしていくことで、時間的に細かく波形の形を調整できます。

 

逆ウェーブレット変換から見る

波形x(t)があったときに、ウェーブレット基底を足しあわせて線型結合で表現できそうだということを見ました。フーリエ変換の記事でもそうしたように、まずはこの基本的なアイデアを見るために逆ウェーブレット変換の方から見ていきましょう。

 

先ほど見た例の


\displaystyle x(t) = w_{1,0}\frac{1}{\sqrt{1}}Ψ(\frac{t-0}{1}) + w_{2,0}\frac{1}{\sqrt{2}}Ψ(\frac{t-0}{2}) + w_{2,1}\frac{1}{\sqrt{2}}Ψ(\frac{t-1}{2})

 

a,bを無数のパターンで試して

 

\displaystyle x(t) = \sum_{a}\sum_{b}w_{a,b}\frac{1}{\sqrt{a}}Ψ(\frac{t-b}{a})

 

というように書くことができます。a,bの値を連続的にずらすのであれば(たとえば少数点を非常に細かく刻んでいく)、積分の形になり

 

\displaystyle x(t) = \int_{a}\int_{b}w_{a,b}\frac{1}{\sqrt{a}}Ψ(\frac{t-b}{a})dadb

 

と表現されることになります。

 

マザーウェーブレットを1つ決めて、それを拡大縮小a、平行移動bしたものを基底として準備して、w_{a,b}をいろいろ調整すればいろんな波形が表現できます、ということを言っています。

 

w_{a,b}とは

aは波形を横方向に縮小する効果があるため、周波数に相当する概念を持ってくれます。そしてbはウェーブレットが値を持つ時間的な位置を調整するため、時間に相当する概念を持っています。

 

従って、w_{a,b} というのはaから見積もられる周波数の波形が、bという時間にどれくらいの振幅で入っているのかを表すことになります。こうして、自然と時間周波数解析の概念を持つことができるのです。

 

当然、|w_{a,b}|^2を取ることによってパワーに相当する概念を得ることができ、スペクトログラムに相当するものも獲得できます。(ウェーブレット変換ではスカログラムという)

 

ウェーブレット変換

さて、マザーウェーブレットΨ(t)を決めることで、a,bを調整したウェーブレット基底を量産し、w_{a,b}で足し合わせることを見ました。

 

しかし実際の問題では、x(t)が得られた時に適当なマザーウェーブレットΨ(t)を準備し、x(t)に対してどのような係数w_{a,b}を使えば足し合わせが上手くいくかを知りたいということになります。

 

つまり、波形x(t)とマザーウェーブレットΨ(t)から、波形を表現できるようなw_{a,b}を上手く決定し、スカログラムを得るなどして時間周波数解析がしたいというのが本来の使い方になります。ここではw_{a,b}a,bの関数として明示しw(a,b)と書きます(意味は何も変わらない)。

 

この目的を達成するのがウェーブレット変換であり、以下の式で表されます。

 

\displaystyle w(a,b)= \int x(t) \frac{1}{\sqrt{a}} \overline{Ψ\left( \frac{x-b}{a} \right)} dt

 

Ψの上に\overline Ψがついていますが、これは共役複素数を取るという意味で、数学的な都合です。マザーウェーブレットが実数である限りは何も変わりはありません。

(フーリエ変換は\exp(-iωt)を基底に使う。その逆変換は\exp(iωt)となり、共役複素数をとっていることになる。この形式に合わせておくと、ウェーブレット変換の計算式にフーリエ変換の色々な定理が使えて便利というだけの話)

 

 

 

ウェーブレット変換を更に詳しく

基底について

基底の選び方で結果が変わる

ウェーブレット変換ではΨ(t)を勝手に決めて、これを拡大縮小、平行移動した基底を量産することでw(a,b)を求めて時間周波数解析を行います。

 

\displaystyle w(a,b)= \int x(t) \frac{1}{\sqrt{a}} \overline{Ψ\left( \frac{x-b}{a} \right)} dx

 

当然、このw(a,b)は使うΨ(t)によって得られる値が異なってきます。

 

短時間フーリエ変換の場合は、あくまでフーリエ変換を時間をずらしながら行っていく手法であるため、各時刻の周波数に対する解析は、通常のフーリエ変換と基本的には変わりません(精度の問題はもちろんあるが)。

 

一方で、勝手に選んだマザーウェーブレットΨ(t)によって結果が変わってくるウェーブレット変換に関しては、ある程度の試行錯誤も必要になってくることになります。得られる数値が違ったとしても、観測した波形の物理的現象が変わるわけではありません(人間の解析の仕方で、対象の本質が揺らぐなんてことはおかしい)。

 

ですから、時間周波数解析を上手く行える手法であったとしても、結果に対する解釈は慎重にならざるを得ません。

 

基底の性質の差異

フーリエ変換で用いられる基底\exp(-iωt)、あるいは\cos(ωt)\sin(ωt)は、数学的に素晴らしい性質を有しています。

 

まず\sin(ωt)\sin(2ωt)の内積は0です。周波数が一致していない限りは、すなわちm,nを整数とした時に

 

(sin(mωt),sin(nωt)) = 0  if m≠n

 

となります。一致している時に限って\sqrt(π)(どんな周波数でも)となります。

 

言い換えれば、\sinのノルムは\sqrt(π)であり、かつ周波数が異なる場合には直交していると言えます。これは\cosにも同じことが言え、更に\sin\cos間は完全に直行しています。

 

波形x(t)をフーリエ吸収展開するとは、優れた直交基底によって展開するということです。基底は別に直行していなくてもいいです。別の方向を向いてさえいれば使うことができます。しかし直交しているというのは、その座標系において、その基底でしか絶対に表現できない成分があるということです。

 

ウェーブレット変換ではこのような性質は一般的に得られません。\frac{1}{\sqrt(a)}という値を掛けることで、ウェーブレット基底のノルムを保つことはできますが、直交性は難しいのです。(もちろん有しているウェーブレット基底もある)

 

利用方法や得意領域

負の側面を強調してしまいましたが、利用する上では大事なことなので大げさに述べました。実際にはウェーブレット変換はとても強力な解析を可能にしてくれます。

 

非定常な波形の時間周波数解析

フーリエ変換は本来、一定の周期の波形に対して用いるものです。波形の周波数が時間的に変化するようなケースには元来使うものではありません。しかし、解析区間を非常に短くする短時間フーリエ変換によって、「その区間では波形の特性は変化していない」という仮定のもとで解析を行います。

 

この仮定をクリアするためには極力時間幅を狭めて、細かく解析が必要(時間分解能を向上)ですが、そうすると行いたかった解析に、多くの場合で支障を来たします。

 

ウェーブレット変換は基本的なアイデアは、勝手に作ったΨ(\frac{t-b}{a})を足し合わせていくというものです。波形が時間的に特性を変えていったとしても、Ψbによって時間を調整され、そのbの地点に、スケールaの波形がどの程度含まれているかを計算するだけです。しかもその波形は、時間的に局在しており(定義域外で値が0)、係数を大きく取ろうがその場所以外に全く影響を及ぼしません。

 

フーリエ変換は時間的に同じ波がずっと続いているという仮定(要するに定常性)があるため厳しかったのですが、ここを生まれながらにしてクリアしているのです。

 

不確定性関係のバランス取り

スケールを変えることはウェーブレット基底の横幅を変えることになります。そして、それにより事実上時間窓を調整していることとなり、「低周波領域では時間窓を大きく取り(時間分解能は悪くなる)、周波数分解能を高め」、「高周波領域では時間窓を小さく取り、周波数分解能を低くする」ということが自然と行えます。

 

低周波領域では、6Hzなのか10Hzなのかは大問題です、一方で高周波の1000Hzなのか1004Hzなのかはさほど問題ではないでしょう。

 

逆に高周波というのはノイズの典型例ですし、あるいは着目すべき高エネルギーの反応であるかもしれませんが、そのような波形が「いつどのくらい生じたか」に関して時間分解能を優先するというのは自然な発想です。

 

ある程度性質が分かっている場合の異常検知

機械に振動が生じる場合に、その波形をセンサーでとっているとしましょう。

異常時には特有の波形が現れるということが分かっていれば、その波形の形に似たマザーウェーブレットを使って、ウェーブレット解析を行うことで、異常を鮮明に捉えることができます。

 

その特有の波形はもしかしたら、複数の周波数が混じった波形をしているかもしれません。フーリエ変換では、複数の周波数の反応を同時に見なければならない場面でも、ウェーブレット変換では特有のマザーウェーブレットを使ってしまうことで、その波形をドンピシャで表現できるスケールaについて、w(a,b)が大きな値を持つことが期待できます。

 

ノイズの除去

ウェーブレット変換には紹介したとおり逆変換が存在しますから、これを利用しない手はありません。

 

通常であればバンドパスフィルタを用いる場合であったとしても、ノイズに相当する成分が複数の飛び飛びの周波数に跨っている可能性もありますし、また時間遅れという作用が無視できない場合もあります。

 

ノイズの形がある程度想定できるならば、その形に近いマザーウェーブレットを用いてw(a,b)を解析し、ノイズ成分を0に置き換えてから逆変換を行えば、ノイズが無い波形が復元されます。

 

最後に

他にも時間周波数解析としては、Stockwell変換やHilbert Huang変換(経験的モード分解)などが存在します。これらについても解説できるように勉強していきたいと思います。

 

今回は時間周波数解析として紹介しましたが、変数を時間ではなく、位置などで扱ってもいいです。例えば画像処理でもウェーブレット変換(あるいは多重解像度解析)は行われており、ピクセルが位置的にどのように変化しているかを捉えるのに用いられます。

 

明示的にスカログラム|w(a,b)|^2を獲得しなくとも、ウェーブレットが時間周波数の特徴量を獲得できることを利用し、ウェーブレットで得られた係数を機械学習の特徴量に用いる手法も数多くあります。また、特徴量として抽出することも明示的に行わず、マザーウェーブレットが活性化関数として用いられたウェーブレットニューラルネットワークなどもあります。

 

 

 

 

s0sem0y.hatenablog.com

 

s0sem0y.hatenablog.com

s0sem0y.hatenablog.com

s0sem0y.hatenablog.com