ヒルベルト変換を用いた信号解析 (包絡線，瞬時周波数，瞬時位相) とその実装 (Matlab, Python)

本記事は，ヒルベルト変換を利用して，どのような信号処理ができるのか記載します．実装例として，MatlabとPythonのソースコードやそれの記載された他のページへのリンクを掲載します．ヒルベルト変換を利用することで，

**(1) 包絡線の検出**
**(2) 瞬時周波数の算出**
**(3) 瞬時位相の算出**

などを行うことができます．もちろん，ヒルベルト変換以外にも手法はありますが，ここでは割愛します．
まず，ヒルベルト変換について述べた後，(1)から(3) のそれぞれについてみていきます．

実装自体はとても簡単なので，具体的なソースコードや図はこのページ内には挙げず，使い方をまとめる形にしました．

# ヒルベルト変換とは
本節は，小野測器さんの資料 (https://www.onosokki.co.jp/HP-WK/eMM_back/emm180.pdf) から抜粋しています．最低限の内容しか書いていませんので，詳細はこちらをご覧ください．とてもわかりやすくまとまっています，すごい！

まず，実数の信号 $x(t)$ を考えます．
ヒルベルト変換 (Hilbert Transform) は，次のコーシー積分による変換を指します．

``` math
\hat{x}(t) = \frac{1}{\pi t}\int^{\infty}_{-\infty}\frac{x(\tau)}{t-\tau}d\tau
```

結論だけいうと，この $\hat{x}(t)$ は，$x(t)$ の各フーリエ成分 (すなわち，フーリエ展開して得られる各三角関数の項) の位相を $\frac{\pi}{2}$  ずらした信号を算出しています．

一般に，

``` math
z(t) = x(t) + i\hat{x}(t)
```

を解析信号 (Analytic Signal) と呼び，この値を用いて解析します．

MatlabおよびPythonでは，ヒルベルト変換はライブラリとして実装されており，1行で解析信号を算出できます，簡単！以下にそのリンクを掲載しました．

https://jp.mathworks.com/help/signal/ref/hilbert.html

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html


# 包絡線の算出
包絡線 (envelope; または瞬時振幅) の値は，次の式で求めることができます．具体的な実装は，リンクを掲載しました (といっても，hilbert(x)をして，そのabsをとるだけですが)．

``` math
envelope(t) = |z(t)| = \sqrt{x(t)^2 + \hat{x}(t)^2}
```

https://jp.mathworks.com/help/signal/ug/envelope-extraction-using-the-analytic-signal.html

https://org-technology.com/posts/Hilbert-transform.html


# 瞬時周波数と瞬時位相の算出
フーリエ変換の考え方からすれば，信号は三角関数の組み合わせによって表現されるので，「位相情報」を有しています．信号の持つ瞬時位相 (instantaneous phase) は，以下の式で求めることができます．

``` math
\theta(t) =  arg(z(t)) = tan^{-1}\frac{\hat{x}(t)}{x(t)}
```

また，瞬時周波数 (instantaneous frequency) は位相から求めることができます．

``` math
freq(t) = \frac{1}{2\pi}\frac{d}{dt}\theta(t)
```


まず，Matlabでのソースコードを見ていきます (以下のリンクより転載)．

``` matlab
z = hilbert(y);
instfrq = fs/(2*pi)*diff(unwrap(angle(z)));
```

https://jp.mathworks.com/help/signal/ug/hilbert-transform-and-instantaneous-frequency.html

解析信号 $z$ の位相を `angle(z)` で求め，それを`unwrap`によって成形します．こうして求まった瞬時位相から`diff`によって差をとって (微分に相当)，上記の計算式を適用しています．サンプリング周波数 $fs$ は，微分の幅になります．配列1つ差は，$\frac{1}{fs}$ に相当するので，これで割っているわけですね．

Pythonの場合もほぼ同様です (以下のリンクより転載).

``` python
import numpy as np
from scipy.signal import hilbert

z = hilbert(x)
inst_phase = np.unwrap(np.angle(z))
inst_freq = np.diff(inst_phase)/(2*np.pi)*fs
```

https://www.gaussianwaves.com/2017/04/extract-envelope-instantaneous-phase-frequency-hilbert-transform/


# おまけ
- 瞬時位相の算出時にunwrapしないとどうなるかは，次のサイトが図を示しています
    - http://retrofocus28.blogspot.com/2013/12/phase-unwrapping_26.html
    
- 信号処理においては，特定の周波数成分に着目することが多々あります．このような場合，特定の周波数帯域でのデジタルバンドパスフィルタの後に，ヒルベルト変換による処理をします．
- 分野によっては，ヒルベルト変換よりも適切な手法が存在する場合があります (私がかつて扱っていたデータでは，ヒルベルト変換よりもウェーブレット変換を使うべきだ，との論文をみかけました)
- 信号の周波数解析では，FFTによるパワー算出と，これら3つがポピュラーかな，と思います (個人の感想)．


# 参考文献
全て再掲になります．
- https://www.onosokki.co.jp/HP-WK/eMM_back/emm180.pdf
- https://jp.mathworks.com/help/signal/ref/hilbert.html
- https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html
- https://jp.mathworks.com/help/signal/ug/envelope-extraction-using-the-analytic-signal.html
- https://org-technology.com/posts/Hilbert-transform.html
- https://jp.mathworks.com/help/signal/ug/hilbert-transform-and-instantaneous-frequency.html
- https://www.gaussianwaves.com/2017/04/extract-envelope-instantaneous-phase-frequency-hilbert-transform/
- http://retrofocus28.blogspot.com/2013/12/phase-unwrapping_26.html