位相スペクトルと群遅延の理論と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.show()

# --- 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.show()

# --- テスト信号の生成 ---
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.show()

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 付近に群遅延ピーク）

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）
    axes[0].plot(w_ap, 20 * np.log10(np.abs(H_ap) + 1e-12), 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}')

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

# --- バターワース + 全域通過フィルタによる群遅延等化 ---
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.65
theta_ap1 = 2 * np.pi * (fc * 0.9) / fs  # カットオフより少し低い周波数
theta_ap2 = 2 * np.pi * (fc * 0.7) / fs

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

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

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)

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