ヒルベルト変換と解析信号：瞬時振幅・位相・周波数のPython実装

はじめに

音声・振動・生体信号などを扱う場面では、信号の「瞬間的な強さ（振幅）」や「瞬間的な周波数」を時刻ごとに知りたい場合があります。しかし通常のフーリエ変換では信号全体の周波数成分しか得られず、「いつ・どの周波数が強いか」という情報が失われます。

ヒルベルト変換（Hilbert Transform）はこの問題を解決する強力な手法です。ヒルベルト変換を用いると解析信号（Analytic Signal）を構成でき、そこから瞬時振幅・瞬時位相・瞬時周波数を直接算出できます。

本記事では、ヒルベルト変換の数学的定義から始め、scipy.signal.hilbert を使った実践的なPython実装まで解説します。

ヒルベルト変換の定義

連続時間ヒルベルト変換

実信号 \(x(t)\) のヒルベルト変換 \(\hat{x}(t)\) は次のように定義されます。

\[\hat{x}(t) = \mathcal{H}\{x(t)\} = \frac{1}{\pi} \text{P.V.} \int_{-\infty}^{\infty} \frac{x(\tau)}{t - \tau} d\tau \tag{1}\]

ここで P.V. はコーシーの主値積分を表します。この積分は \(t = \tau\) の特異点を対称的に回避することで定義されます。

式 \((1)\) は \(x(t)\) と \(h(t) = \frac{1}{\pi t}\) の畳み込みとみなせます。

\[\hat{x}(t) = x(t) * \frac{1}{\pi t} \tag{2}\]

周波数領域での解釈

畳み込みの定理より、ヒルベルト変換は周波数領域で次のように表現できます。

\[\hat{X}(f) = H(f) \cdot X(f) \tag{3}\]

ここでフィルタ \(H(f)\) は次の形をとります。

\[H(f) = \begin{cases} -j & (f > 0) \\ 0 & (f = 0) \\ +j & (f < 0) \end{cases} \tag{4}\]

つまり、ヒルベルト変換は正の周波数成分を \(-90°\) 位相シフトし、負の周波数成分を \(+90°\) 位相シフトするオールパスフィルタです。振幅スペクトルは変化せず、位相のみが変化します。

この性質は直感的に理解できます。例えば \(x(t) = \cos(2\pi f_0 t)\) に対してヒルベルト変換を適用すると、\(-90°\) 位相シフトにより \(\hat{x}(t) = \sin(2\pi f_0 t)\) が得られます。

\[\mathcal{H}\{\cos(2\pi f_0 t)\} = \sin(2\pi f_0 t) \tag{5}\]

解析信号

解析信号の定義

実信号 \(x(t)\) から解析信号（Analytic Signal） \(z(t)\) を次のように構成します。

\[z(t) = x(t) + j\hat{x}(t) \tag{6}\]

解析信号は複素数値をとります。実部が元の信号 \(x(t)\)、虚部がそのヒルベルト変換 \(\hat{x}(t)\) です。

解析信号の周波数スペクトルは次のようになります。

\[Z(f) = \begin{cases} 2X(f) & (f > 0) \\ X(0) & (f = 0) \\ 0 & (f < 0) \end{cases} \tag{7}\]

すなわち、解析信号は元の信号の負の周波数成分を除去し、正の周波数成分を2倍にした複素信号です。この片側スペクトル表現が瞬時特性の抽出を可能にします。

瞬時振幅・瞬時位相・瞬時周波数

解析信号 \(z(t)\) を極形式で表すと次のようになります。

\[z(t) = A(t) \cdot e^{j\phi(t)} \tag{8}\]

ここから以下の瞬時特性を直接算出できます。

瞬時振幅（包絡線）:

\[A(t) = |z(t)| = \sqrt{x(t)^2 + \hat{x}(t)^2} \tag{9}\]

瞬時位相:

\[\phi(t) = \angle z(t) = \arctan\!\left(\frac{\hat{x}(t)}{x(t)}\right) \tag{10}\]

瞬時周波数:

\[f_i(t) = \frac{1}{2\pi} \frac{d\phi(t)}{dt} \tag{11}\]

瞬時周波数は位相の時間微分であり、信号の周波数が時間とともに変化する様子を追跡できます。

Pythonによる実装

scipy.signal.hilbert の使い方

SciPyには scipy.signal.hilbert 関数が用意されており、解析信号を簡単に計算できます。

from scipy.signal import hilbert
import numpy as np

# 入力信号
x = np.array([...])  # 実信号

# 解析信号の計算（離散ヒルベルト変換を使用）
z = hilbert(x)

# 瞬時振幅（包絡線）
amplitude = np.abs(z)

# 瞬時位相（ラジアン）
phase = np.angle(z)

# 瞬時周波数（Hz）
fs = 1000  # サンプリング周波数
inst_freq = np.diff(np.unwrap(phase)) / (2 * np.pi) * fs

scipy.signal.hilbert は FFT ベースで実装されており、計算量は \(O(N \log N)\) です。内部的には式 \((4)\) の周波数領域フィルタを適用し、解析信号を返します。

注意: scipy.signal.hilbert の戻り値は（実部が元の信号、虚部がヒルベルト変換）という解析信号 \(z(t)\) です。ヒルベルト変換 \(\hat{x}(t)\) のみ必要な場合は np.imag(hilbert(x)) で取り出せます。

基本的な動作確認

正弦波・余弦波に対するヒルベルト変換を確認します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert

# --- 信号の生成 ---
fs = 1000          # サンプリング周波数 [Hz]
f0 = 10            # 信号周波数 [Hz]
t = np.arange(0, 1, 1/fs)

cos_signal = np.cos(2 * np.pi * f0 * t)
sin_signal = np.sin(2 * np.pi * f0 * t)

# --- ヒルベルト変換 ---
cos_hilbert = np.imag(hilbert(cos_signal))  # cos → sin になるはず
sin_hilbert = np.imag(hilbert(sin_signal))  # sin → -cos になるはず

# --- プロット ---
fig, axes = plt.subplots(2, 1, figsize=(10, 6))

axes[0].plot(t[:200], cos_signal[:200], label='cos(t)')
axes[0].plot(t[:200], cos_hilbert[:200], '--', label='H{cos(t)} ≈ sin(t)')
axes[0].set_title('Hilbert Transform of cos(t)')
axes[0].set_xlabel('Time [s]')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(t[:200], sin_signal[:200], label='sin(t)')
axes[1].plot(t[:200], sin_hilbert[:200], '--', label='H{sin(t)} ≈ -cos(t)')
axes[1].set_title('Hilbert Transform of sin(t)')
axes[1].set_xlabel('Time [s]')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 理論値との誤差確認
print(f"H{{cos}} と sin の最大誤差: {np.max(np.abs(cos_hilbert - sin_signal)):.2e}")
print(f"H{{sin}} と -cos の最大誤差: {np.max(np.abs(sin_hilbert + cos_signal)):.2e}")
# 出力例:
# H{cos} と sin の最大誤差: 1.33e-15
# H{sin} と -cos の最大誤差: 1.33e-15

誤差は数値精度の範囲（\(10^{-15}\) 程度）であり、\(\mathcal{H}\{\cos\} = \sin\)、\(\mathcal{H}\{\sin\} = -\cos\) が正確に成立することが確認できます。

実践1：包絡線検出

AM変調信号の包絡線抽出

AM（振幅変調）信号では搬送波（キャリア）に情報信号が乗っています。ヒルベルト変換を用いた包絡線検出により、元の情報信号を復元できます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert

# --- AM変調信号の生成 ---
fs = 10000          # サンプリング周波数 [Hz]
t = np.arange(0, 1, 1/fs)

# メッセージ信号（情報信号、5Hz）
f_msg = 5
message = 0.5 * np.sin(2 * np.pi * f_msg * t)

# 搬送波（キャリア、500Hz）
f_carrier = 500
carrier = np.cos(2 * np.pi * f_carrier * t)

# AM変調: x(t) = (1 + m(t)) * carrier
am_signal = (1 + message) * carrier

# --- ヒルベルト変換による包絡線検出 ---
analytic = hilbert(am_signal)
envelope = np.abs(analytic)  # 瞬時振幅 = 包絡線

# --- プロット ---
fig, axes = plt.subplots(3, 1, figsize=(10, 9))

# AM変調信号と包絡線
axes[0].plot(t, am_signal, alpha=0.5, label='AM Signal')
axes[0].plot(t, envelope, 'r-', linewidth=2, label='Envelope (Hilbert)')
axes[0].plot(t, -envelope, 'r-', linewidth=2)
axes[0].set_xlim(0, 0.4)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('AM Signal and Envelope')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 元のメッセージ信号と検出された包絡線の比較
axes[1].plot(t, 1 + message, 'k--', label='True Envelope (1 + message)')
axes[1].plot(t, envelope, 'r-', alpha=0.8, label='Detected Envelope')
axes[1].set_xlim(0, 0.4)
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('True Envelope vs Detected Envelope')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# 復調されたメッセージ信号
demodulated = envelope - 1  # DC成分（1）を除去
axes[2].plot(t, message, 'k--', label='Original Message')
axes[2].plot(t, demodulated, 'r-', alpha=0.8, label='Demodulated (Hilbert)')
axes[2].set_xlim(0, 0.4)
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Amplitude')
axes[2].set_title('Demodulation Result')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 復調精度の確認
rmse = np.sqrt(np.mean((message - demodulated)**2))
print(f"復調RMSE: {rmse:.6f}")

包絡線として検出された波形が元のメッセージ信号（5Hz）の形状に正確に対応していることが確認できます。

実践2：瞬時周波数の推定

チャープ信号への適用

周波数が時間とともに変化するチャープ信号に対して、瞬時周波数の推定を行います。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert, chirp

# --- チャープ信号の生成（50Hz → 200Hz、1秒間） ---
fs = 4000  # サンプリング周波数 [Hz]
t = np.arange(0, 1, 1/fs)
f_start, f_end = 50, 200
signal = chirp(t, f0=f_start, f1=f_end, t1=1, method='linear')

# --- ヒルベルト変換による瞬時周波数の推定 ---
analytic = hilbert(signal)
phase = np.unwrap(np.angle(analytic))

# 位相の数値微分で瞬時周波数を計算
inst_freq = np.diff(phase) / (2 * np.pi) * fs

# 理論値（線形に変化）
f_theory = f_start + (f_end - f_start) * t[:-1]

# --- プロット ---
fig, axes = plt.subplots(3, 1, figsize=(10, 9))

# チャープ信号の時間波形
axes[0].plot(t, signal, alpha=0.7)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title(f'Chirp Signal ({f_start} Hz → {f_end} Hz)')
axes[0].grid(True, alpha=0.3)

# 瞬時振幅（包絡線）
axes[1].plot(t, np.abs(analytic), 'r-')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Instantaneous Amplitude (Envelope)')
axes[1].set_ylim(0, 1.5)
axes[1].grid(True, alpha=0.3)

# 瞬時周波数と理論値の比較
axes[2].plot(t[:-1], inst_freq, label='Instantaneous Frequency (Hilbert)')
axes[2].plot(t[:-1], f_theory, 'k--', label='Theoretical Frequency')
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Frequency [Hz]')
axes[2].set_title('Instantaneous Frequency vs Theoretical')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 推定精度
mae = np.mean(np.abs(inst_freq - f_theory))
print(f"瞬時周波数の平均絶対誤差: {mae:.4f} Hz")

チャープ信号の瞬時周波数が理論値（50Hz→200Hzの線形変化）と非常に近い精度で推定されていることが確認できます。

実践3：振動信号の包絡線スペクトル解析（HFRT法）

ヒルベルト変換を応用した**包絡線スペクトル解析（HFRT: Hilbert-based Fault-Related Technique）**は、回転機械の故障診断に広く用いられます。転がり軸受の疲労剥離はインパルス状の振動を生じさせ、その包絡線スペクトルに特徴的な周波数ピークが現れます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert, butter, sosfiltfilt

# --- 軸受故障を模した振動信号の生成 ---
np.random.seed(42)
fs = 20000          # サンプリング周波数 [Hz]
t = np.arange(0, 1, 1/fs)

# 高周波共振成分（軸受の固有振動数: 5kHz）
f_resonance = 5000
resonance_carrier = np.cos(2 * np.pi * f_resonance * t)

# 故障周波数（インパルスの繰り返し周期: 100Hz）
f_fault = 100
# インパルス列（故障によるパルス）
impulse_train = np.zeros_like(t)
impulse_indices = (np.arange(0, 1, 1/f_fault) * fs).astype(int)
impulse_indices = impulse_indices[impulse_indices < len(t)]
impulse_train[impulse_indices] = 1.0

# 故障信号 = インパルス列 × 減衰する共振
from scipy.signal import lfilter
alpha = 500  # 減衰率
decay = np.exp(-alpha * np.mod(t, 1/f_fault))
fault_component = np.convolve(impulse_train, decay[:len(t)//f_fault], mode='same')
fault_signal = fault_component * 0.3 * resonance_carrier

# ガウスノイズを加えた観測信号
observed = fault_signal + 0.5 * np.random.randn(len(t))

# --- ステップ1: バンドパスフィルタで共振周波数帯域を抽出 ---
sos = butter(4, [4000, 6000], btype='bandpass', fs=fs, output='sos')
filtered = sosfiltfilt(sos, observed)

# --- ステップ2: ヒルベルト変換で包絡線を計算 ---
analytic = hilbert(filtered)
envelope = np.abs(analytic)

# --- ステップ3: 包絡線のFFTでスペクトル解析 ---
N = len(envelope)
envelope_spectrum = np.abs(np.fft.rfft(envelope - np.mean(envelope))) / N
freqs = np.fft.rfftfreq(N, 1/fs)

# --- プロット ---
fig, axes = plt.subplots(3, 1, figsize=(10, 9))

# 観測信号
axes[0].plot(t[:3000], observed[:3000], alpha=0.7)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Observed Vibration Signal')
axes[0].grid(True, alpha=0.3)

# バンドパス後の包絡線
axes[1].plot(t[:3000], filtered[:3000], alpha=0.5, label='Bandpass Filtered')
axes[1].plot(t[:3000], envelope[:3000], 'r-', linewidth=1.5, label='Envelope')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Bandpass Filtered Signal and Envelope')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# 包絡線スペクトル（0〜500Hz）
mask = freqs <= 500
axes[2].plot(freqs[mask], envelope_spectrum[mask])
axes[2].axvline(f_fault, color='r', linestyle='--', alpha=0.7,
               label=f'Fault Freq = {f_fault} Hz')
for harmonic in range(2, 6):
    axes[2].axvline(f_fault * harmonic, color='r', linestyle=':', alpha=0.4)
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Amplitude')
axes[2].set_title('Envelope Spectrum (Fault Diagnosis)')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

包絡線スペクトルには故障周波数（100Hz）とその高調波（200Hz、300Hz…）にピークが現れます。この手法は軸受やギアの早期故障検出に実用されています。

離散ヒルベルト変換の注意点

端点効果（Edge Effect）

FFTベースのヒルベルト変換は信号が周期的であることを前提とします。実際の信号では両端で非連続性が生じるため、端点付近で精度が低下する端点効果が生じます。

import numpy as np
from scipy.signal import hilbert

# 端点効果の確認
fs = 1000
t = np.arange(0, 1, 1/fs)
f0 = 50
signal = np.cos(2 * np.pi * f0 * t)  # 両端が非連続でない（完全な周期）

# 信号長をずらして非連続にする
t2 = np.arange(0, 1 + 1/(3*f0), 1/fs)  # 1.33サイクル（非整数周期）
signal2 = np.cos(2 * np.pi * f0 * t2)

analytic2 = hilbert(signal2)
envelope2 = np.abs(analytic2)

# 理論包絡線は1.0（一定）なので、端点の誤差を確認
edge_error = np.max(np.abs(envelope2[:20] - 1.0))
print(f"端点付近の最大誤差: {edge_error:.4f}")
# → 端点で0.1〜0.3程度の誤差が生じることがある

対策:

信号の両端にゼロパディングを追加し、端点の影響が有効区間に入らないようにする
解析対象の前後に余分なデータを取得し、処理後に端点を除去する（バッファリング）

ナロウバンド信号の条件

瞬時周波数の意味ある解釈には、対象信号がナロウバンド（狭帯域）信号であることが前提となります。広帯域信号や複数の周波数成分が混在する信号では、瞬時周波数が負になったり物理的に意味のない値をとったりします。

条件	推奨する前処理
広帯域信号	バンドパスフィルタで関心帯域を抽出
複数成分が混在	EMD（経験的モード分解）でIMFに分解
非定常信号	STFTやウェーブレット変換との併用を検討

まとめ

ヒルベルト変換は信号を \(-90°\) 位相シフトする演算で、周波数領域では \(H(f) = -j\,\text{sgn}(f)\) と表現できる
元の信号とそのヒルベルト変換を組み合わせた解析信号は負の周波数成分を持たず、瞬時特性の抽出に適した複素信号である
解析信号の絶対値・偏角・位相の時間微分から瞬時振幅・瞬時位相・瞬時周波数を直接計算できる
scipy.signal.hilbert は FFT ベースで \(O(N \log N)\) の効率的な実装を提供しており、包絡線検出・AMデモジュレーション・故障診断など幅広い用途に使用できる
端点効果とナロウバンド条件の注意点を理解した上で適切に前処理を施すことが重要である

高速フーリエ変換（FFT）の仕組みとPython実装 - ヒルベルト変換がFFTを内部で利用していることを理解するための基礎です。
指数移動平均（EMA）フィルタの周波数特性 - Z変換による伝達関数の導出など、フィルタ解析の基礎を解説しています。
窓関数とパワースペクトル密度（PSD）の理論とPython実装 - ヒルベルト変換と組み合わせると有用なスペクトル解析手法を解説しています。
ノッチフィルタの設計とPython実装 - 電源ノイズ除去後に包絡線解析を適用するユースケースで関連します。
バンドパスフィルタの設計とPython実装 - ヒルベルト変換の前処理として用いるバンドパスフィルタの設計を解説しています。
ウェーブレット変換の理論とPython実装 - 非定常信号の時間-周波数解析においてヒルベルト-黄変換（HHT）と比較されることが多い手法です。

参考文献

Gabor, D. (1946). “Theory of communication.” Journal of the Institution of Electrical Engineers, 93(3), 429-457.
Huang, N. E., et al. (1998). “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis.” Proceedings of the Royal Society of London A, 454, 903-995.
Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
SciPy signal.hilbert documentation