HELLO CYBERNETICS

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

時間周波数解析:短時間フーリエ変換

 

 

follow us in feedly

短時間フーリエ変換

時間周波数解析とは、最も単純に言えば信号f(t)のフーリエ変換F(ω)を見る際に、これをF(ω,t)として考えることで、時間変化する周波数情報を解析することです。

 

時間周波数解析と普通に言った時には、必ずしもフーリエ変換を用いるとは限りませんが、時間と周波数の間で変換を行うフーリエ変換から入るのが最も理解がしやすいかと思われます。

 

フーリエ変換

フーリへ解析に関する基本事項は以下の記事でまとめられています。

s0sem0y.hatenablog.com

 

 

 

フーリエ変換おさらい

ここでは時間周波数解析を見ていく上でのフーリエ変換についておさらいします。

まず、フーリエ逆変換(上の式)とフーリエ変換(下の式)を見ましょう。

 

f(t)=\frac{1}{2π}\int_{-∞}^{∞}F(ω)e^{jωt}dω

 

F(ω)=\int_{-∞}^{∞}f(t)e^{-jωt}dt

 

なぜフーリエ逆変換から見るのかというと、これが波形f(t)を周期関数e^{jωt}に重みF(ω)を掛けて線形結合したものであることを表しており、物理的な意味が理解しやすいからです。実際には線型結合ではなく、ωを連続的に動かすことで積分計算になっていますが、心としては線型結合をしているというふうで良いでしょう。

 

波形f(t)は角周波数ωを持つ周期関数によって表現できるということです。

 

そこで、f(t)には角周波数ωの周期関数がどれくらいの重みF(ω)で含まれているのかを求めたいということになります。それがフーリエ変換です。通常この重みF(ω)のことをスペクトルと言います。

 

フーリエ変換の落とし穴的な仮定

ところでフーリエ変換とは、周期的波形をフーリエ級数展開するというところから始まり、任意の波形を無限の大きさの周期を持つ周期関数とみなしたものです。

いわば、いまから波形f(t)をフーリエ変換すると言った時には、気持ちの面では、ある周期を持つ波形をしっかりと一周期分観測できており、その周期は無限大に大きいということを認めていることになります。言い換えれば、波形の観測時間は無限大であると言っているのです。

実用上は、解析したい波形を十分に長いと思われる時間観測する必要があります。

 

もしも観測した波形が、観測を終えた次の瞬間に全く違う波形に変化するとしたならば、それは無限大の周期を持つ波形の一周期分を観測できたとは言えません。なぜならば、周期的であるとは同じ波形が繰り返されることを指すからです。無限大の周期を観測できたというならば、観測の外で波形が大きく変わっているような状況は起こってはいけないのです。

 

もちろん、その変化自体の原因が分かっており、それは全く持って解析対象ではないというのであれば、どうでもいいことでしょう。例えば回路に流れる交流電圧を掛けて、振動的な電流が流れているところをスペクトル解析したとしましょう。電源を切って観測をやめた瞬間に、電源の影響で回路の振る舞いが変わるのは当然ですが、これは興味のない現象です。

 

そうだとしても一周期分をぴったりと観測できているケースは少ないでしょう。

周期的データであるならば、観測データの始まりと終わりが、全く同じ値になっていなければなりません。そんなことは稀です。フーリエ級数展開ならば「周期的な波形であること」を、フーリエ変換であるならば、「無限大の周期が一周期分であること」が要求されますが、観測データはそうではありません。

 

このような場合には、観測データの始まりと終わりが、全く同じ値になるように修正を施します。単純には、観測データの端の方で波形f(t)が減衰するような関数w(t)を掛けてしまうのです。このようなw(t)を窓関数と言います。

波形の元データは少し変更を受けますが、それでもフーリエ変換自体の重大な仮定を破るよりはマシです。実際窓関数を準備しなければ、波形が周期的であることを無理やり満たしているようなスペクトルが計算されてしまいます。(ギブス現象)

 

窓関数利用の例

f:id:s0sem0y:20170108164853p:plain

上記の画像は

小野測器-FFTアナライザについて (page11)

から引用しています。実用的な面で非常に詳しく分かりやすい開設があるので一読してみてください。

窓関数は、フーリエ解析の理論ではあまり扱われませんが、実解析をする上では必須です。

 

時間周波数解析への誘い

フーリエ変換に関してはここまでで、時間周波数解析に入っていきます。

 

波形f(t)のスペクトルF(ω)とは、真の意味で波形の一周期分を観測してあるのならば、未来永劫その波形のスペクトルF(ω)には変化が無いはずなのです。

 

しかし無限に長い周期を持つ波形を完全に観測するのは不可能です。従って次に全く同じ波形を解析したつもりでも、結果がわずかに変わる場合もあるでしょう。僅かな変化ならば些細な問題です。しかし解析するたびに全然違う結果が出るような、全く通用しないケースもあります。

 

そのようなケースは、そもそも波形f(t)の性質が時間と共に著しく変化していると考えるべきです。すなわちスペクトルF(ω)自体が刻一刻と変化しているだろうと考えるべきなのです。つまりスペクトル自体が時間の関数でありF(ω,t)と表記すべきということです。

 

時間周波数解析の出発点

波形f(t)の様相が時間と共に変化していくのであれば、波形f(t)を一部の時間だけ繰り抜いて、その部分だけフーリエ変換をするということが考えられます。

 

例えば、f(t)[t_0,t_1]の区間をフーリエ変換してF_1(ω)を作り、その次に[t_1,t_2]の区間に対してF_2(ω)を作って行けば、繰り抜いた時間区間の関数F(ω,t)を構成することができます。

 

このように時間区間を繰り抜くような関数を時間窓と呼び、このような手法でF(ω,t)を構成する手法を短時間フーリエ変換と言います。

 

コレは結局、ある区間にフーリエ変換をしているだけの単純な話です。

もしも十分に観測された、性質の変わらない波形f(t)に対して短時間フーリエ変換をしたとすると、そのスペクトルF(ω,t)は時間tによらないはずです。どこを区切って解析してみても、その区切った時間窓自身も十分な大きさならばフーリエ変換の結果に影響は及ぼされないはずだからです。

 

 

いや、実を言うと、実世界のフーリエ変換自体がある種の短時間フーリエ変換の一部を使っていると言えます。

 

要するに本当は時間区間[0,∞]まで観測しなければならないものを、[0,t]まででやめてフーリエ変換しているのです。そして後日、更に解析を続け(前日の時間を基準に0と考えて)時間[t',t'']まで観測したデータに対してフーリエ変換を施します。そして前日の観測データに対するスペクトルに比べて、やっぱり今回も変化がないので、波形f(t)の性質は変化しない、普遍的なものだったんだな、とするわけです。

 

そういう意味で短時間フーリエ変換は全く特別なことはしていません。

ただその時間窓を著しく狭く取るというだけの話です。

 

時間周波数解析は、f(t)に対する時間窓の位置に応じて、スペクトルF(ω)が変化するような場合に有効になります。

 

短時間フーリエ変換

短時間フーリエ変換を用いる場合には、時間窓で波形を切り抜いていくと言いました。

しかしその切り抜き方は先程のように[t_0,t_1][t_1,t_2]、などと隣り合わせに作る必要はありません。例えば[0,2][0.2,2.2][0.4,2.4]などのように、少しずつずらしても良いのです。このずらし具合を十分に小さくすれば、スペクトルF(ω,t)という連続的な関数だと思ってしまっても差し支えないでしょう。

 

この方がよりキメの細かい時間変化を追うことができます。当然、これは各区間に単なるフーリエ変換を施しているに過ぎないので、実用上は窓関数をさらに準備する必要があります。

そのような場合は、そもそも窓関数自体を時間窓に使ってしまえばいいです。

なぜなら、窓関数w(t)は観測外で値が0になる関数であるので、これを切り取りたい波形の幅で準備すれば、周期を整えつつ解析したい部分だけ抜き取れるからです。

 

つまり元々のフーリエ変換

 

F(ω)=\int_{-∞}^{∞}f(t)e^{-jωt}dt

 

に対して以下のような変更を施します。

 

F(ω,t)=\int_{-∞}^{∞}f(τ)w(τ-t)e^{-jωτ}dτ

 

このとき積分変数はτであり、積分の中身で移動するのはτの方です。

従ってフーリエ変換に対して変更があるのは、w(τ-t)という値だけです。

この計算によって、例えば時間t=0に対してはw(τ)という窓関数を掛けた際のスペクトルF(ω,0)を得ることができます。これは単に窓関数を掛けただけのフーリエ変換です。

次に時間t=0.1にはw(τ-0.1)を窓関数としたフーリエ変換をするわけですから、これはw(τ)0.1だけ右にずらした窓となっていますから、しっかりと繰り抜く場所を変動させられているのがわかります。

 

これで先ほど言った、繰り抜く場所を変えながらフーリエ変換をしていくことが達成出来ました。

 

短時間フーリエ変換の欠点

こうしてフーリエ変換を工夫することで短時間フーリエ変換という、スペクトルが時間変化する波形にも対応できる術を手に入れました。しかし、その本質はやはりフーリエ変換に宿っており、観測区間を細かく区切ってしまうだけという単純な方法です。

 

これほどまでに単純で、かつ、時間変化にも対応できるのならば、フーリエ変換の完全なる上位互換のように感じるかもしれません。しかも、時間幅を狭めれば狭めるほど、窓のスライドの幅も細かくすればするほど、たくさん情報を得られるではないかと考えるかもしれません。

 

やはりそうは上手くいかず、実際には不確定性関係によってその利用は著しく制限されます。

 

不確定性関係

(時間周波数解析で)不確定性関係とは時間tと角周波数ωに関する

 

ΔtΔω \geq \frac{1}{2}

 

という関係です。Δというのは簡単に言えば精度のことです。

非常に大雑把ですが、もしも時間幅[0,0.5]という0.5秒間でのデータを使った場合には、この時間幅での波形f(t)のスペクトルF(ω)に関しては、不確定性関係が

 

0.5 * Δω \geq \frac{1}{2}

 

Δω \geq 1

 

と言った具合に計算され、約1(rad/s)もズレが生じます。スペクトルF(ω)でのωの位置にブレがあるということです。

例えば、sin(2t)をフーリエ変換すれば、F(ω)ω=2にのみ値を持つはずなのに、実際には不確定性関係のΔωだけズレた部分にも値を持ってしまい、完全な真の情報は得られないのです。

 

従って正しいスペクトルF(ω)を獲得したければ、波形が性質の一定のものであろうとなかろうと、なるべく時間幅を大きくする必要があります。

 

短時間フーリエ変換に至っては、大抵今から解析しようという波形は、波形の性質が時間と共に変化しているような扱いづらいものです。そのような波形を短い時間幅で繰り抜いて、フーリエ変換していくのですから、そもそも1回1回のフーリエ変換の結果自体が怪しいものだと言えます。

 

とは言っても、詳細なスペクトルの情報は分からなくとも、スペクトルのピーク(つまりどのωが最も含まれているか)を知る分にはそれなりの手がかりになります。

 

例えばsin(t・t)という波形を短時間フーリエ変換しましょう。

これは、角周波数ω=tであって、要するに時間と共に角周波数が大きくなっていくような波形です。本来t=2という時点において、スペクトルF(ω)ω=2でしか値を持たないはずなのですが、不確定性関係によって、その周辺でも値を持ってしまうでしょう。しかし、この程度ならば、ピークの位置が正しくω=2の位置に来るはずです。

時間周波数解析によってF(ω,t)を得た場合には、当然ピークの位置が線形に増加していくのが確認できるでしょう。

 

精度が必要な場所とそうでない場所

不確定性関係によって、時間窓の取り方でスペクトルの精度に影響を及ぼすことが分かりました。

ところで、時間窓を決めたら、そのスペクトルのぼやけ具合Δωは決まります。さて、いまぼやけ具合が適当に5だったとしましょう。そして今から解析する波形f(t)は非常に複雑であり、含まれている波形e^{jωt}の種類も膨大であるとします。つまりたくさんのωF(ω)が値を持つということです。そしてそれが、5だけぼやけるという状況です。

 

果たしてω=1でスペクトルが値を持つ場合に、これが5の分だけぼやけるのと、ω=100でスペクトルが値を持つ場合に、これが5の分だけぼやけるのは同じでしょうか?

 

ぼやけてる値は同じにしても、割合が全く違います。

 

短時間フーリエ変換で解析をする限りは、波形f(t)に含まれるe^{jωt}は、ωの値によらずすべて均等にぼやけることになります。

 

これは如何なものでしょうか。

 

ウェーブレット変換

基本的には低周波成分に関しては、周波数側になるべく精度を与えたくて、高周波側には周波数側には精度を与えなくてもいいと考えられるでしょう。言い換えれば、低周波成分の波形に関しては十分に時間窓を大きくとってやって、高周波側には幅の狭い時間窓を準備したいということです。

 

これってよく考えれば当たり前です。

低周波の波というのは、1秒間に少ししか揺れないわけですから、大きく時間窓を取ってやらなければ、波の形を観測できません。一方高周波は少ない観測時間でも、何周期も波を観測できてしまいます。

 

ですから、解析したい周波数に応じて時間窓の幅を変えるのは当たり前といえば当たり前です。

 

しかし通常、解析しなければいけないようなたちの悪い波形というのは、どんな周波数が含まれており、どんな周波数が重要なのかもわからないケースがあります。

短時間フーリエ変換であれこれ時間窓の幅を変えながら何度も解析することも考えられますが、世の中には、周波数帯域に応じて時間窓を変えるような便利な時間周波数解析があります。

 

それがウェーブレット変換です。

 

ここでは詳細な説明は割愛しますが、ウェーブレット変換の出現の意義はこんなところにあったのです。それを知っているだけでも、ウェーブレット変換というものの活用方法がなんとなく分かるでしょう。

 

また、フーリエ変換ではe^{jωt}を使って波形を表現する一方で、ウェーブレット変換はいろんな関数の形を自分で設計して、それを用いて元の波形を表現することができます。つまり線型結合という概念に戻った時、その基底関数自体を選べるのです。

そのようなことができるアルゴリズムが出現したのは驚くべきことです。(ただしフーリエ変換にはe^{jωt}が完全直交基底になるという数学的な強みがあり、ウェーレット変換は一般にそうではないですが)

 

このことによって、例えば工作機械の異常な波形というものが既に調べつくされている時に、その波形をウェーブレット変換の基底(マザーウェーブレット)として設計しておけば、工作機械が異常な波形を出してしまったときに、ウェーブレット変換の係数(フーリエ変換で言うところのスペクトル)が直ちにそれを知らせてくれるようになります。

 

ウェーブレット変換は生まれながらにしての時間周波数解析ですから、元々F(ω,t)に相当する概念しか持ちあわせておらず、スペクトルF(ω)に相当するものはありません。

時間変化するスペクトルF(ω,t)の二乗|F(ω,t)|^2のことを普通はスペクトログラムと言い、ウェーブレット変換でのF_w(ω,t)の二乗|F_w(ω,t)|^2のことはスカログラムと言います。