More than 3 years have passed since last update.

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

Last updated at 2021-03-03Posted at 2021-03-02

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

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

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

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

ヒルベルト変換とは

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

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

\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}$ ずらした信号を算出しています．

一般に，

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

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

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

包絡線の算出

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

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

瞬時周波数と瞬時位相の算出

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

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

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

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

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

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

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

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

おまけ

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

参考文献

全て再掲になります．

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up