位相スペクトルの時間方向の差分と arctan経由の時間微分（瞬時周波数）の比較を行うPythonコード

phase-frequency-comparison.py

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

# パラメータ設定
fs = 1000  # サンプリング周波数
t = np.arange(0, 1, 1 / fs)  # 時間軸
f0 = 5  # 基本周波数
f1 = 10  # 周波数変調の上限

# テスト信号の生成（周波数変調信号）
freq = f0 + f1 * np.sin(2 * np.pi * 2 * t)  # 時間変化する周波数
phase = 2 * np.pi * np.cumsum(freq) / fs  # 位相の積分
x = np.sin(phase)  # テスト信号

# 解析信号の生成
z = signal.hilbert(x)
phase_spectrum = np.unwrap(np.angle(z))  # 位相スペクトル

# 方法1: 位相スペクトルの時間差分による瞬時周波数の計算
inst_freq_diff = np.diff(phase_spectrum) * fs / (2 * np.pi)
inst_freq_diff = np.append(inst_freq_diff, inst_freq_diff[-1])  # 長さを合わせる

# 方法2: arctanの微分による瞬時周波数の計算
y = signal.hilbert(x)
Q = y.imag
I = y.real
dQ = np.gradient(Q, 1 / fs)
dI = np.gradient(I, 1 / fs)
inst_freq_arctan = (I * dQ - Q * dI) / (2 * np.pi * (I**2 + Q**2))

# 誤差の計算
error = inst_freq_diff - inst_freq_arctan
rms_error = np.sqrt(np.mean(error**2))

# プロット
plt.figure(figsize=(12, 10))

# 元信号
plt.subplot(411)
plt.plot(t, x)
plt.title("元信号")
plt.xlabel("時間 [s]")
plt.ylabel("振幅")
plt.grid(True)

# 位相スペクトルの差分による瞬時周波数
plt.subplot(412)
plt.plot(t, inst_freq_diff)
plt.title("位相スペクトルの差分による瞬時周波数")
plt.xlabel("時間 [s]")
plt.ylabel("周波数 [Hz]")
plt.grid(True)

# arctanの微分による瞬時周波数
plt.subplot(413)
plt.plot(t, inst_freq_arctan)
plt.title("arctanの微分による瞬時周波数")
plt.xlabel("時間 [s]")
plt.ylabel("周波数 [Hz]")
plt.grid(True)

# 誤差
plt.subplot(414)
plt.plot(t, error)
plt.title(f"誤差 (RMS = {rms_error:.6f} Hz)")
plt.xlabel("時間 [s]")
plt.ylabel("周波数差 [Hz]")
plt.grid(True)

plt.tight_layout()
plt.show()

print(f"RMS誤差: {rms_error:.6f} Hz")
print(f"最大誤差: {np.max(np.abs(error)):.6f} Hz")
print(f"平均誤差: {np.mean(error):.6f} Hz")