scipy.signal.hilbertでヒルベルト変換・包絡線検出・AM/FM復調・軸受診断をPython実装

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

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

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}")

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")

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

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

fc = 2000           # キャリア周波数 [Hz]
kf = 300            # 周波数偏移定数 [Hz/V]
f_msg = 5           # メッセージ信号周波数 [Hz]
message = np.sin(2 * np.pi * f_msg * t)

# 累積位相 = 2π * kf * ∫m(τ)dτ
phase_msg = 2 * np.pi * kf * np.cumsum(message) / fs
fm_signal = np.cos(2 * np.pi * fc * t + phase_msg)

# --- ヒルベルト変換で瞬時周波数を復元 ---
analytic = hilbert(fm_signal)
inst_phase = np.unwrap(np.angle(analytic))
inst_freq = np.diff(inst_phase) / (2 * np.pi) * fs

# キャリア周波数を差し引いてメッセージを復元
demodulated = (inst_freq - fc) / kf

# --- プロット ---
fig, axes = plt.subplots(2, 1, figsize=(10, 6))
axes[0].plot(t[:2000], fm_signal[:2000], alpha=0.7)
axes[0].set_title('FM Signal (carrier 2 kHz, message 5 Hz)')
axes[0].set_xlabel('Time [s]'); axes[0].grid(True, alpha=0.3)

axes[1].plot(t[:-1], message[:-1], 'k--', label='Original message')
axes[1].plot(t[:-1], demodulated, 'r-', alpha=0.8, label='Demodulated (Hilbert)')
axes[1].set_title('FM Demodulation via Instantaneous Frequency')
axes[1].set_xlabel('Time [s]'); axes[1].legend(); axes[1].grid(True, alpha=0.3)
plt.tight_layout(); plt.show()

rmse = np.sqrt(np.mean((message[:-1] - demodulated) ** 2))
print(f"FM復調RMSE: {rmse:.4f}")

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()

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程度の誤差が生じることがある

import numpy as np
from scipy.signal import hilbert as scipy_hilbert

def hilbert_manual(x):
    """scipy.signal.hilbert と等価な実装（教育用）"""
    N = len(x)
    # ステップ1: 実信号を FFT で周波数領域へ
    Xf = np.fft.fft(x)

    # ステップ2: 周波数領域フィルタ H(f) を構築
    # 正の周波数: 2倍, DC と Nyquist: そのまま, 負の周波数: 0
    h = np.zeros(N)
    if N % 2 == 0:
        h[0] = 1                 # DC
        h[N // 2] = 1            # Nyquist
        h[1:N // 2] = 2          # 正の周波数
        # h[N // 2 + 1:] = 0     # 負の周波数（既に0）
    else:
        h[0] = 1
        h[1:(N + 1) // 2] = 2

    # ステップ3: 周波数領域で乗算してから IFFT
    z = np.fft.ifft(Xf * h)
    return z

# scipy 実装との一致を確認
rng = np.random.default_rng(0)
x = rng.standard_normal(1024)
z_scipy = scipy_hilbert(x)
z_manual = hilbert_manual(x)
print(f"最大誤差: {np.max(np.abs(z_scipy - z_manual)):.2e}")
# 期待出力: 1e-14 程度（浮動小数点誤差）

import numpy as np
import pywt
from scipy.signal import hilbert, stft

fs = 1000
t = np.arange(0, 2, 1/fs)
msg = 0.5 * (1 + np.sin(2 * np.pi * 2 * t))   # 2 Hz 振幅変調
x = msg * np.cos(2 * np.pi * 100 * t)         # 100 Hz キャリア

# Hilbert: 直接包絡線
env_hilbert = np.abs(hilbert(x))

# STFT: 100Hz帯の振幅を時間軸にプロジェクション
f, tt, Z = stft(x, fs=fs, nperseg=128)
env_stft = np.abs(Z[np.argmin(np.abs(f - 100))])

# Wavelet: 100Hz相当スケールの絶対値
scales = pywt.scale2frequency('morl', np.arange(1, 64)) * fs
idx = np.argmin(np.abs(scales - 100))
coeffs, _ = pywt.cwt(x, np.arange(1, 64), 'morl', sampling_period=1/fs)
env_wavelet = np.abs(coeffs[idx])