位相スペクトルと群遅延の理論とPython実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, sosfilt, freqz, firwin, group_delay, lfilter, iirfilter

# --- フィルタ設計 ---
fs = 1000          # サンプリング周波数 [Hz]
fc = 100           # カットオフ周波数 [Hz]
N_iir = 8          # IIRフィルタ次数
N_fir = 64         # FIRフィルタタップ数

# バターワースIIRフィルタ（ba形式）
b_iir, a_iir = butter(N_iir, fc, fs=fs, output='ba')

# 線形位相FIRフィルタ（hamming窓）
b_fir = firwin(N_fir, fc, fs=fs, window='hamming')
a_fir = np.array([1.0])

# --- 周波数応答の計算 ---
worN = 4096
w, H_iir = freqz(b_iir, a_iir, worN=worN, fs=fs)
_, H_fir = freqz(b_fir, a_fir, worN=worN, fs=fs)

# --- 位相応答の計算 ---
phase_iir_wrapped = np.angle(H_iir)              # 折り返し位相
phase_iir = np.unwrap(np.angle(H_iir))           # アンラップ位相
phase_fir = np.unwrap(np.angle(H_fir))           # FIR（線形位相）

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

# 振幅応答
axes[0].plot(w, 20 * np.log10(np.abs(H_iir) + 1e-12), label=f'IIR Butterworth N={N_iir}')
axes[0].plot(w, 20 * np.log10(np.abs(H_fir) + 1e-12), label=f'FIR Hamming {N_fir} taps', linestyle='--')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title('Magnitude Response')
axes[0].set_xlim(0, fs / 2)
axes[0].set_ylim(-80, 5)
axes[0].axvline(fc, color='gray', linestyle=':', alpha=0.5, label='Cutoff')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 折り返し位相（IIRのみ）
axes[1].plot(w, np.degrees(phase_iir_wrapped), label='IIR (wrapped)', color='orange')
axes[1].plot(w, np.degrees(phase_iir), label='IIR (unwrapped)', color='blue')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Phase [degrees]')
axes[1].set_title('Phase Response: Wrapped vs Unwrapped (IIR)')
axes[1].set_xlim(0, fs / 2)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# IIR vs FIR 位相比較
axes[2].plot(w, np.degrees(phase_iir), label=f'IIR Butterworth N={N_iir} (nonlinear)')
axes[2].plot(w, np.degrees(phase_fir), label=f'FIR Hamming {N_fir} taps (linear)', linestyle='--')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Phase [degrees]')
axes[2].set_title('Phase Response: IIR vs FIR')
axes[2].set_xlim(0, fs / 2)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('phase_response_comparison.png')

# --- 折り返しの検証：|Δφ| > 0.9π となる不連続点の個数 ---
wrap_jumps = np.sum(np.abs(np.diff(phase_iir_wrapped)) > np.pi * 0.9)
print(f"折り返し位相の不連続ジャンプ数: {wrap_jumps}")
print(f"アンラップ位相の範囲: {np.degrees(phase_iir.max()):.2f}度 〜 {np.degrees(phase_iir.min()):.2f}度")

折り返し位相の不連続ジャンプ数: 13
アンラップ位相の範囲: 0.00度 〜 -905.56度

# --- scipy.signal.group_delay による群遅延計算 ---
worN = 4096
w_gd, gd_iir = group_delay((b_iir, a_iir), w=worN, fs=fs)
w_gd, gd_fir = group_delay((b_fir, a_fir), w=worN, fs=fs)

# --- 差分近似による群遅延計算（教育的目的） ---
w_arr, H_iir_arr = freqz(b_iir, a_iir, worN=worN * 2, fs=fs)
omega_arr = w_arr * 2 * np.pi / fs  # Hz → rad/sample
phase_arr = np.unwrap(np.angle(H_iir_arr))

# τ_g(ω) ≈ -Δφ/Δω
gd_diff = -np.diff(phase_arr) / np.diff(omega_arr)
w_diff = (w_arr[:-1] + w_arr[1:]) / 2  # 中点

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

# IIR vs FIR 群遅延比較
axes[0].plot(w_gd, gd_iir, label=f'IIR Butterworth N={N_iir}', color='blue')
axes[0].plot(w_gd, gd_fir, label=f'FIR Hamming {N_fir} taps (constant = {(N_fir-1)/2:.1f})',
             linestyle='--', color='red')
axes[0].axhline((N_fir - 1) / 2, color='red', linestyle=':', alpha=0.5)
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Group Delay [samples]')
axes[0].set_title('Group Delay Comparison: IIR vs FIR')
axes[0].set_xlim(0, fs / 2)
axes[0].set_ylim(-5, 50)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 差分近似による群遅延（IIR）
axes[1].plot(w_diff, gd_diff, label='IIR (finite difference approximation)', color='purple', alpha=0.7)
axes[1].plot(w_gd, gd_iir, label='IIR (scipy.signal.group_delay)', color='blue', linestyle='--')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Group Delay [samples]')
axes[1].set_title('Group Delay: Finite Difference vs scipy.signal.group_delay')
axes[1].set_xlim(0, fs / 2)
axes[1].set_ylim(-5, 50)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('group_delay_comparison.png')

print(f"FIR群遅延: 最小={gd_fir.min():.6f}, 最大={gd_fir.max():.6f}, 標準偏差={gd_fir.std():.2e}")
print(f"理論値 (N_fir-1)/2 = {(N_fir - 1) / 2}")
mask_pb = (w_gd >= 0) & (w_gd <= 90)
print(f"IIR群遅延(通過域 0-90Hz): 最小={gd_iir[mask_pb].min():.3f}, 最大={gd_iir[mask_pb].max():.3f} サンプル")

gd_iir_interp = np.interp(w_diff, w_gd, gd_iir)
mask_valid = (w_diff > 5) & (w_diff < fs / 2 - 5)
diff_err = np.abs(gd_diff[mask_valid] - gd_iir_interp[mask_valid])
print(f"差分近似とscipy.group_delayの誤差: 平均={diff_err.mean():.4f}, 最大={diff_err.max():.4f} サンプル")

FIR群遅延: 最小=31.499999, 最大=31.500002, 標準偏差=3.59e-08
理論値 (N_fir-1)/2 = 31.5
IIR群遅延(通過域 0-90Hz): 最小=7.888, 最大=14.584 サンプル
差分近似とscipy.group_delayの誤差: 平均=0.0009, 最大=0.9531 サンプル

# --- テスト信号の生成 ---
fs = 1000
t = np.arange(0, 0.5, 1/fs)

# 複数の周波数成分からなる矩形パルス的な信号
# （10Hz + 30Hz + 50Hz + 70Hz + 90Hz の重ね合わせ）
freqs_sig = [10, 30, 50, 70, 90]
x = sum(np.sin(2 * np.pi * f * t) / f for f in freqs_sig)

# --- フィルタリング ---
fc_test = 200  # カットオフを高くして通過させる

# IIRフィルタ（非線形位相）
b_iir_test, a_iir_test = butter(6, fc_test, fs=fs, output='ba')
y_iir = lfilter(b_iir_test, a_iir_test, x)

# FIRフィルタ（線形位相）
b_fir_test = firwin(101, fc_test, fs=fs, window='hamming')
y_fir = lfilter(b_fir_test, [1.0], x)

# FIRの群遅延補正（位相を合わせるためのサンプルシフト）
fir_delay = 50  # (101 - 1) / 2 = 50 サンプル
y_fir_aligned = np.concatenate([y_fir[fir_delay:], np.zeros(fir_delay)])

# --- プロット ---
fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)

axes[0].plot(t, x, 'k', linewidth=1.5)
axes[0].set_title('Original Signal (sum of 10, 30, 50, 70, 90 Hz)')
axes[0].set_ylabel('Amplitude')
axes[0].grid(True, alpha=0.3)

axes[1].plot(t, y_iir, color='blue')
axes[1].plot(t, x, 'k--', alpha=0.3, linewidth=1, label='Original')
axes[1].set_title('After IIR Butterworth (nonlinear phase — waveform distortion)')
axes[1].set_ylabel('Amplitude')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].plot(t, y_fir_aligned, color='red')
axes[2].plot(t, x, 'k--', alpha=0.3, linewidth=1, label='Original')
axes[2].set_title('After FIR (linear phase, delay-corrected — shape preserved)')
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Amplitude')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('waveform_distortion.png')

# --- 定量評価：過渡部を除く中央区間でのRMSE ---
mask_mid = (t >= 0.15) & (t <= 0.35)
rmse_iir = np.sqrt(np.mean((y_iir[mask_mid] - x[mask_mid]) ** 2))
rmse_fir = np.sqrt(np.mean((y_fir_aligned[mask_mid] - x[mask_mid]) ** 2))
print(f"元信号との一致度(中央区間0.15-0.35秒, RMSE): IIR={rmse_iir:.5f}, FIR(遅延補正後)={rmse_fir:.5f}")
print(f"信号のRMS振幅: {np.sqrt(np.mean(x[mask_mid] ** 2)):.5f}")

元信号との一致度(中央区間0.15-0.35秒, RMSE): IIR=0.02584, FIR(遅延補正後)=0.00011
信号のRMS振幅: 0.07694

def allpass_filter_2nd(r, theta):
    """
    2次全域通過フィルタ係数の計算。

    H_AP(z) = (r^2 - 2r*cos(θ)*z^{-1} + z^{-2}) / (1 - 2r*cos(θ)*z^{-1} + r^2*z^{-2})

    Parameters
    ----------
    r : float
        極の半径（0 < r < 1）
    theta : float
        極の角度 [rad]

    Returns
    -------
    b, a : ndarray
        フィルタ係数（b = 分子, a = 分母）
    """
    c = 2 * r * np.cos(theta)
    r2 = r ** 2
    b = np.array([r2, -c, 1.0])
    a = np.array([1.0, -c, r2])
    return b, a


# --- 全域通過フィルタの位相・群遅延特性の可視化 ---
worN = 4096
fig, axes = plt.subplots(3, 1, figsize=(10, 10))

r_values = [0.5, 0.7, 0.9]
theta_ap = np.pi / 4  # 極の角度（ω = π/4 付近に群遅延ピーク）

peak_info = []
for r in r_values:
    b_ap, a_ap = allpass_filter_2nd(r, theta_ap)
    w_ap, H_ap = freqz(b_ap, a_ap, worN=worN, fs=fs)
    w_ap_gd, gd_ap = group_delay((b_ap, a_ap), w=worN, fs=fs)

    # 振幅応答（常に 0 dB）
    mag_db = 20 * np.log10(np.abs(H_ap) + 1e-12)
    axes[0].plot(w_ap, mag_db, label=f'r={r}')

    # 位相応答
    axes[1].plot(w_ap, np.degrees(np.unwrap(np.angle(H_ap))), label=f'r={r}')

    # 群遅延
    axes[2].plot(w_ap_gd, gd_ap, label=f'r={r}')

    peak_freq = w_ap_gd[np.argmax(gd_ap)]
    peak_info.append((r, mag_db.max(), mag_db.min(), gd_ap.max(), peak_freq))

axes[0].set_title('All-Pass Filter: Magnitude Response (should be 0 dB)')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_ylim(-3, 3)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].set_title('All-Pass Filter: Phase Response')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Phase [degrees]')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].set_title('All-Pass Filter: Group Delay (controllable peak)')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Group Delay [samples]')
axes[2].set_ylim(0, 30)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('allpass_response.png')

for r in r_values:
    idx = r_values.index(r)
    r_v, mag_max, mag_min, gd_peak, peak_freq = peak_info[idx]
    print(f"r={r_v}: 振幅範囲=[{mag_min:.2e}, {mag_max:.2e}] dB, "
          f"群遅延ピーク={gd_peak:.3f}サンプル @ {peak_freq:.2f}Hz")

r=0.5: 振幅範囲=[8.68e-12, 8.69e-12] dB, 群遅延ピーク=3.610サンプル @ 118.16Hz
r=0.7: 振幅範囲=[8.68e-12, 8.69e-12] dB, 群遅延ピーク=6.010サンプル @ 124.39Hz
r=0.9: 振幅範囲=[8.67e-12, 8.70e-12] dB, 群遅延ピーク=19.105サンプル @ 125.00Hz

# --- バターワース + 全域通過フィルタによる群遅延等化 ---
fs = 1000
fc = 100
N_butter = 6

b_butter, a_butter = butter(N_butter, fc, fs=fs, output='ba')
w_hz, gd_butter = group_delay((b_butter, a_butter), w=4096, fs=fs)

# 全域通過フィルタを設計して群遅延を平坦化
# バターワースの群遅延が低い低周波数帯にAPFの群遅延ピークを配置し、底上げする
# 複数のAPFを直列に接続
r_ap = 0.5
theta_ap1 = 2 * np.pi * 5 / fs   # 5Hz付近に群遅延ピークを配置
theta_ap2 = 2 * np.pi * 10 / fs  # 10Hz付近に群遅延ピークを配置

b_ap1, a_ap1 = allpass_filter_2nd(r_ap, theta_ap1)
b_ap2, a_ap2 = allpass_filter_2nd(r_ap, theta_ap2)

# 直列接続の係数（畳み込みで合成）
b_total = np.convolve(np.convolve(b_butter, b_ap1), b_ap2)
a_total = np.convolve(np.convolve(a_butter, a_ap1), a_ap2)

w_hz, gd_total = group_delay((b_total, a_total), w=4096, fs=fs)
w_hz, gd_ap1 = group_delay((b_ap1, a_ap1), w=4096, fs=fs)
w_hz, gd_ap2 = group_delay((b_ap2, a_ap2), w=4096, fs=fs)

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

# 振幅応答（APFは振幅に影響しない）
w_hz_mag, H_butter = freqz(b_butter, a_butter, worN=4096, fs=fs)
_, H_total = freqz(b_total, a_total, worN=4096, fs=fs)
axes[0].plot(w_hz_mag, 20 * np.log10(np.abs(H_butter) + 1e-12),
             label='Butterworth only', color='blue')
axes[0].plot(w_hz_mag, 20 * np.log10(np.abs(H_total) + 1e-12),
             label='Butterworth + APF (equalized)', color='red', linestyle='--')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title('Magnitude Response (APF does not change amplitude)')
axes[0].set_xlim(0, fs / 2)
axes[0].set_ylim(-80, 5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 群遅延の比較
axes[1].plot(w_hz, gd_butter, label=f'Butterworth N={N_butter}', color='blue')
axes[1].plot(w_hz, gd_ap1 + gd_ap2, label='APF contribution (ap1+ap2)', color='green', linestyle='--')
axes[1].plot(w_hz, gd_total, label='Butterworth + APF (equalized)', color='red', linewidth=2)
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Group Delay [samples]')
axes[1].set_title('Group Delay Equalization with All-Pass Filter')
axes[1].set_xlim(0, fc * 2)
axes[1].set_ylim(0, 30)
axes[1].axvline(fc, color='gray', linestyle=':', alpha=0.5, label='Cutoff')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('group_delay_equalization.png')

mask_pb = (w_hz >= 0) & (w_hz <= fc)
gd_butter_range = gd_butter[mask_pb].max() - gd_butter[mask_pb].min()
gd_total_range = gd_total[mask_pb].max() - gd_total[mask_pb].min()
mag_diff = np.abs(20 * np.log10(np.abs(H_butter) + 1e-12) - 20 * np.log10(np.abs(H_total) + 1e-12))
print(f"通過域(0-{fc}Hz)の群遅延変動幅: Butterworthのみ={gd_butter_range:.3f}サンプル, "
      f"Butterworth+APF={gd_total_range:.3f}サンプル")
print(f"振幅応答差の最大値(全帯域): {mag_diff.max():.6f} dB")

通過域(0-100Hz)の群遅延変動幅: Butterworthのみ=4.994サンプル, Butterworth+APF=1.473サンプル
振幅応答差の最大値(全帯域): 0.000000 dB

from scipy.signal import cheby1, ellip, bessel

fs = 1000
fc = 100
N = 6

# 各IIRフィルタの設計
filters = {
    'Butterworth': butter(N, fc, fs=fs, output='ba'),
    'Chebyshev I (1dB)': cheby1(N, 1, fc, fs=fs, output='ba'),
    'Elliptic (1dB, 60dB)': ellip(N, 1, 60, fc, fs=fs, output='ba'),
    'Bessel': bessel(N, fc, norm='mag', fs=fs, output='ba'),
}

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

gd_stats = {}
for name, (b, a) in filters.items():
    w_hz, H = freqz(b, a, worN=4096, fs=fs)
    w_gd, gd = group_delay((b, a), w=4096, fs=fs)

    axes[0].plot(w_hz, 20 * np.log10(np.abs(H) + 1e-12), label=name)
    axes[1].plot(w_gd, gd, label=name)

    mask_pb = (w_gd >= 0) & (w_gd <= fc)
    gd_stats[name] = (gd[mask_pb].min(), gd[mask_pb].max(), gd[mask_pb].std())

axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title(f'Magnitude Response: IIR Filter Comparison (N={N}, fc={fc} Hz)')
axes[0].set_xlim(0, 300)
axes[0].set_ylim(-80, 5)
axes[0].axvline(fc, color='gray', linestyle=':', alpha=0.5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Group Delay [samples]')
axes[1].set_title('Group Delay Comparison: IIR Filters')
axes[1].set_xlim(0, 300)
axes[1].set_ylim(0, 50)
axes[1].axvline(fc, color='gray', linestyle=':', alpha=0.5)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('iir_filter_comparison.png')

print("通過域(0-100Hz)群遅延統計（最小, 最大, 標準偏差）[サンプル]:")
for name, (mn, mx, sd) in gd_stats.items():
    print(f"  {name}: 最小={mn:.3f}, 最大={mx:.3f}, 標準偏差={sd:.3f}")

通過域(0-100Hz)群遅延統計（最小, 最大, 標準偏差）[サンプル]:
  Butterworth: 最小=5.946, 最大=10.940, 標準偏差=1.676
  Chebyshev I (1dB): 最小=6.858, 最大=31.263, 標準偏差=4.781
  Elliptic (1dB, 60dB): 最小=5.412, 最大=46.929, 標準偏差=7.203
  Bessel: 最小=4.160, 最大=4.596, 標準偏差=0.129