はじめに
ヒルベルト変換は非常に便利な変換で、qiita だと、
など素晴らしい記事があるので、こちらを参考ください。
その他、ヒルベルトファン変換という、時系列データに対して、Emperical Mode Decomposition (EMD) を用いて、固有振動モード(Independet Mode Function) に分割してから、ヒルベルト変換を用いる方法など、様々な応用が広がっています。
-
Hilbert-Huang Transform の wikipedia の説明
など。
ただし、これらの情報を読んでも、ヒルベルト変換の説明の部分の、「位相が90度回る成分を取り出す、1/x との畳み込み... 」 あたりが、なんだかもやもやすることないでしょうか。私もちょっと混乱してしまったので、そういう人向けにまとめておきたいと思います。
ヒルベルト変換を理解する手順
数学的な準備
ヒルベルト変換を理解するためには、いくつか数学的な準備をしておきます。
任意の関数は、偶関数と奇関数に分割できる
まずは、任意の関数 $g(x)$ は、次のように２つの関数の足し算、
g(x)=\frac{g(x)+g(-x)}{2}+\frac{g(x)-g(-x)}{2}
に分割できます。それぞれ、x --> -x に変換すると、符号が変わらない偶関数と、-1 倍される奇関数、の和に分割できます。
これはヒルベルト変換に限らず、複雑な関数を、より単純な２つの関数の和に分割する、という意味で、極めて重要な性質になりますので、暗記しておいて損はないです。
フーリエ変換は関数の偶奇性を保存する
偶関数をフーリエ変換すると、フーリエ空間でも偶関数になります。奇関数も同様で、関数の偶奇性はフーリエ変換で変わりません。
言葉で説明してみるとこんな感じ。フーリエ変換は、偶関数(cos)と奇関数(sin)を掛け合わせて、無限区間で積分(有限区間の場合はその波の整数倍の範囲)するので、積分[偶関数 x 偶関数(cos)] は有限の値になりますが、積分[偶関数 x 奇関数(sin)]はゼロになります。逆に、 積分[奇関数 x 偶関数(cos)] はゼロになりますが、積分[奇関数 x 奇関数(sin)]は有限の値になります。周波数空間での偶奇性は、波数 k の偶奇を変えた場合の符号の変化なので、cos(-kx) = cos(kx), sin(-kx) = -sin(kx) のように、 sin, cos の偶奇性だけで決まります。つまり、時間空間で偶関数の場合は 偶関数の cos との積分だけ残り、時間空間で奇関数の場合は 奇関数の sin との積分だけ残るため、フーリエ変換後も偶奇性が変わらない、ということになります。
式で見た方がわかりやすい話でもあるので、
が素晴らしくよくまとまっておりますので、こちらをご参照ください。
複素フーリエ変換で、偶関数は実部、奇関数は虚部に移る
これは、関数の偶奇性が保存する話と同じことです。複素フーリエ変換の実数部は cos 関数(偶関数)との積の和、虚数部は sin 関数(奇関数)との積の和、で表されるので、偶関数は実部のみ、奇関数は虚部のみを持つことになります。
実数関数の複素フーリエ変換は、実部が偶関数、虚部が奇関数になる
実数関数は、複素共役をとっても、何も変わらない。一方、複素フーリエ変換するときにかかる exp(iωt) の部分は、複素共役を取ると、 exp(-iωt) に変換される。これは、exp(i x (-ω) x t) と同じなので、 ω --> -ω という操作と、複素共役を取ることは同じことになる。
これは頭で考えるよりも、ちょこっと式を見た方が良いので、実関数のフーリエ変換が G(f) をかけるとき、 G(f) の複素共役は、積分[実関数 x exp(iωt)] なので、 exp(iωt) だけが複素数で虚部の符号が変わり、かつそれは周波数を負にすることと同値なので、
$ G^*(f)=G(-f) $ とかける。これをそのまま、実部と虚部の関係式に書き下すと、
\begin{cases}\operatorname{Re} G(f)=\operatorname{Re} G(-f) & (G \text { の実部が偶関数 }) \\ -\operatorname{Im} G(f)=\operatorname{Im} G(-f) & (G \text { 虚部が奇関数 })\end{cases}
となる。つまり、実数関数を複素フーリエ変換すると、かならず実部が偶関数、虚部が奇関数になる。
なども参照ください。
純虚数関数は複素フーリエ変換で、実部が奇関数、虚部が隅関数になる
上と同様だが、純虚数関数は、複素共役をとると、-1 がかかるため、偶奇性が逆になる。
負の周波数とはなにか？
ヒルベルト変換の説明には、よく負の周波数の絵が出てきますので、そのイメージを確認しておくのがよいです。cos(ft) の f の符号が変わっても、cos(ft)の符号は何も変わらないので、実関数だけ考えている場合はあまり意味のある定義にはなりません。
負ということは、複素フーリエ空間で意味が出てきます。exp(-iωt)とexp(iωt)で何が違うか、、というと、時間が進むときに、位相 exp(iφ) の φ = -ωt と φ = ωt の違いと思えば、時間が進むと、位相角が増えるのか、減るのか、という違いになります。(それ以上の意味はあまりないような。。)
ヒルベルト変換の理解
ヒルベルト変換に直結する部分になります。
解析信号 Analysic Signal の理解
解析信号 「Analysic Signal」 というのを理解する必要があります。
に、説明があるのですが、日本語では Analysic Signal という言葉がそもそも浸透してない気がしますし、その直訳の"解析信号"というのが固有名詞っぽくない名前なのでこれも浸透してない気がします。
ここで、フーリエ変換を $ \mathcal{F} $ で表すことにして、時系列データ s(t) をフーリエ変換して、S(f)となる状況を考えます。
S(f) =\mathcal{F}[s(t)]
と掛けます。
解析信号 $s_{\mathrm{a}}(t)$ は、
\begin{aligned}
s_{\mathrm{a}}(t) & \overset{def}{=} \mathcal{F}^{-1}[S(f)+\operatorname{sgn}(f) \cdot S(f)]
\end{aligned}
で定義されます。
sgn関数 $\operatorname{sgn}(f)$ は、中身が正なら何もしない、負なら -1 を掛ける、という操作をするものです。この部分がヒルベルト変換にとって、最重要なパートです。
$\operatorname{sgn}(f)$ がなければ、ただのフーリエ変換して、逆変換するだけなので、何も起こらないことになります。
$\operatorname{sgn}(f)$ は、一体なにをやってるのでしょうか？ 中身が正ならそのまま、負なら -1 を掛ける、という操作は、「偶関数を奇関数に変換」する、あるいは、「奇関数を偶関数に変換」する、というオペレーションに相当します。つまり、関数の偶奇性を入れ替える操作をしています。
これが、どういう意味を持つか？ですが、数学的な準備で考えたように、偶奇性を変えることは、フーリエ変換で保存するはずの偶奇性を破ってしまうので、逆変換しても元には戻りません。つまり、$\operatorname{sgn}(f)$ があることが極めて重要で、偶奇性の保存を破ることで、解析信号 s_a(t) と 元の関数 s(t) に違いが生まれるのです。
では、生まれた違い、とは何なのでしょうか？ 元の関数 s(t) から情報が増えることはありません。生まれた違いは、 s(t) は一つの側面にすぎません。 ここで、s(t) を複素フーリエ級数展開することを考えましょう。
s(t) = \sum_k a_k e^{-i \omega t} = \sum_k a_k ( \cos(\omega t) + i \sin(\omega t) )
フーリエ空間で、偶奇性を変換する、ということは、cos --> sin, sin --> cos に変換する操作に対応します。これは、複素数空間で幾何学的に考えると、位相が90度くるっと回す操作に対応します。(+90, -90度とかはとりあえず気にしません)。
つまり、解析信号 $s_{\mathrm{a}}(t)$ は、フーリエ空間で偶奇性を交換することで、"s(t) の位相を90度ずらした関数" を生成し、
\begin{aligned}
s_{\mathrm{a}}(t) & \overset{def}{=} s(t) + s(t) の位相を90度ずらした純虚数関数
\end{aligned}
の２つの和に分解する操作に対応します。
ヒルベルト変換と Analysitc Signal の関係
Analysitc Signal s_a(t) は、ヒルベルト変換を用いると、
\begin{aligned}
s_{\mathrm{a}}(t) & \overset{def}{=} s(t) + s(t) の位相を90度ずらした純虚数関数 \\
& \overset{def}{=} s(t) + i ~ s(t)のヒルベルト変換した実関数 \\
\end{aligned}
と掛ける。
ここで、ヒルベルト変換 $\mathcal{H}$ を、
\mathcal{H}[s(t)] = \left[\frac{1}{\pi t} * s(t)\right]
と定義される。 * は畳み込みです。1/t も s(t) もどちらも実関数で、それらの時間空間での畳み込みも実関数にしかならないので、実関数のヒルベルト変換は、実関数になります。
Analysic signal は、s(t)をヒルベルト変換した関数を、$\hat{s}$で表すと、
\begin{aligned}
s_{\mathrm{a}}(t) & \overset{def}{=}
\mathcal{F}^{-1}[S(f)+\operatorname{sgn}(f) \cdot S(f)] \\
& =s(t)+ i \underbrace{\left[\frac{1}{\pi t} * s(t)\right]}_{\mathcal{H}[s(t)]} \\
& = s(t)+ i \hat{s}(t),
\end{aligned}
となります。(上の日本語の式と同じです)
ヒルベルト変換やAnalysic signalのご利益
大元の情報は s(t) で尽きているわけで、ごにょごにょ変換して、何か意味があるのだろうか？と思う人もいるかもしれません。でも、フーリエ変換も情報は増えてませんが、よく使いますよね。それはなぜかというと、より素性の良い関数の和に分割できるためです。(厳密には、どんな関数も無限個の三角関数で分割できるかは自明ではないですが、分割できるとすると、それが良い分割であれば物事が見えやすくなります。)
ヒルベルト変換やAnalysic signalも同様に考えると、元の情報に対して、フーリエ変換を施し、フーリエ空間で位相を90度回した信号を生成することで、時事刻々と波の振幅や位相を一回の演算で出せる、という点がメリットになります。素性のよい分割をすることで、現象を記述する良い情報が簡単に得られる、そのための変換という感じですね。
ヒルベルト変換の計算方法(scipy.signal.hilbert)
信号を $x(t)$ として、解析信号(analytic signal)を $x_a(t)$ と置くと、$x(t)$のフーリエ変換(F)とステップ関数(U)を用いて、
$$
x_a(t) =F^{-1}(F(x(t)) 2 U) = x(t) + i y(t)
$$
で一発で計算できて、scipy.signal.hilbert を用いて、
計算できます。
ヒルベルト・ファン (Hilbert Huang) 変換 の計算方法
ヒルベルト・ファン変換は、pyhtt を使えば一発で計算してくれます。
非定常的な時系列データや、衝撃データなどの解析には Hilbert Huang 変換のような、経験的な方法で積極的にデータを分割する方法が良い場合があります。(数学的な一意性や擬似信号の扱いなどには注意が必要です。)