Z変換とデジタルフィルタの伝達関数：極零点解析とPython実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, freqz, tf2zpk, zpk2tf, sosfilt

def zplane(b, a, ax=None, title="Pole-Zero Plot"):
    """
    デジタルフィルタの極零点プロットを描画する。

    Parameters
    ----------
    b : array_like
        分子係数（b[0] + b[1]*z^-1 + ...）
    a : array_like
        分母係数（a[0] + a[1]*z^-1 + ...）
    ax : matplotlib Axes, optional
        プロット先の Axes。None の場合は新規作成。
    title : str
        プロットタイトル

    Returns
    -------
    zeros : ndarray
        零点の配列
    poles : ndarray
        極の配列
    """
    if ax is None:
        fig, ax = plt.subplots(figsize=(6, 6))

    # 極・零点を計算
    zeros, poles, gain = tf2zpk(b, a)

    # 単位円を描画
    theta = np.linspace(0, 2 * np.pi, 512)
    ax.plot(np.cos(theta), np.sin(theta), "k--", linewidth=0.8, alpha=0.5, label="Unit Circle")

    # 実軸・虚軸
    ax.axhline(0, color="k", linewidth=0.5, alpha=0.4)
    ax.axvline(0, color="k", linewidth=0.5, alpha=0.4)

    # 零点（○）
    ax.scatter(zeros.real, zeros.imag, s=80, marker="o", facecolors="none",
               edgecolors="blue", linewidths=1.5, zorder=5, label=f"Zeros ({len(zeros)})")

    # 極（×）
    ax.scatter(poles.real, poles.imag, s=80, marker="x",
               color="red", linewidths=1.5, zorder=5, label=f"Poles ({len(poles)})")

    # 軸設定
    ax.set_xlabel("Real Part")
    ax.set_ylabel("Imaginary Part")
    ax.set_title(title)
    ax.set_aspect("equal")
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

    # 表示範囲を自動調整（単位円が収まるよう余白を確保）
    lim = max(1.5, np.max(np.abs(poles)) * 1.2, np.max(np.abs(zeros)) * 1.2)
    ax.set_xlim(-lim, lim)
    ax.set_ylim(-lim, lim)

    return zeros, poles

# --- 安定システム（全極が単位円内） ---
b_stable = np.array([1.0])
a_stable = np.array([1.0, -0.5, 0.3])  # 極: |p| < 1

# --- 不安定システム（極が単位円外に出る） ---
b_unstable = np.array([1.0])
a_unstable = np.array([1.0, -1.5, 0.9])  # 極: |p| > 1

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
z_s, p_s = zplane(b_stable, a_stable, ax=axes[0], title="Stable System\n(All Poles Inside Unit Circle)")
z_u, p_u = zplane(b_unstable, a_unstable, ax=axes[1], title="Unstable System\n(Poles Outside Unit Circle)")

print(f"Stable system poles:   {np.roots(a_stable)}")
print(f"  Magnitudes: {np.abs(np.roots(a_stable))}")
print(f"Unstable system poles: {np.roots(a_unstable)}")
print(f"  Magnitudes: {np.abs(np.roots(a_unstable))}")

plt.tight_layout()
plt.show()

Stable system poles:   [0.25+0.479j  0.25-0.479j]
  Magnitudes: [0.541  0.541]
Unstable system poles: [0.75+0.661j  0.75-0.661j]
  Magnitudes: [1.0  1.0]  ← 単位円上（限界安定）

# バターワースフィルタ（ローパス、カットオフ 100 Hz、サンプリング 1000 Hz）
fs = 1000   # Hz
fc = 100    # Hz
orders = [2, 4, 6]

fig, axes = plt.subplots(len(orders), 2, figsize=(14, 4 * len(orders)))

for i, N in enumerate(orders):
    # SOS形式で設計
    sos = butter(N, fc, fs=fs, output='sos')
    # 伝達関数形式に変換
    b, a = butter(N, fc, fs=fs, output='ba')

    # 左列: 極零点プロット
    zplane(b, a, ax=axes[i, 0], title=f"Butterworth N={N} Pole-Zero Plot")

    # 右列: 周波数応答
    w, h = freqz(b, a, worN=4096, fs=fs)
    axes[i, 1].plot(w, 20 * np.log10(np.abs(h) + 1e-12), "b", linewidth=1.5)
    axes[i, 1].axvline(fc, color="gray", linestyle="--", alpha=0.7, label=f"fc={fc}Hz")
    axes[i, 1].axhline(-3, color="red", linestyle=":", alpha=0.7, label="-3 dB")
    axes[i, 1].set_xlim(0, fs / 2)
    axes[i, 1].set_ylim(-80, 5)
    axes[i, 1].set_xlabel("Frequency [Hz]")
    axes[i, 1].set_ylabel("Magnitude [dB]")
    axes[i, 1].set_title(f"Butterworth N={N} Frequency Response")
    axes[i, 1].legend()
    axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

def freq_response_from_zpk(zeros, poles, gain, omega_range=None):
    """
    極零点から周波数応答の大きさと位相を幾何学的に計算する。

    Parameters
    ----------
    zeros : array_like
        零点の配列（複素数）
    poles : array_like
        極の配列（複素数）
    gain : float
        ゲイン定数
    omega_range : array_like, optional
        評価する角周波数の配列（rad/sample）

    Returns
    -------
    omega : ndarray
        角周波数
    H_mag : ndarray
        振幅応答
    H_phase : ndarray
        位相応答 [rad]
    """
    if omega_range is None:
        omega_range = np.linspace(0, np.pi, 1024)

    H_mag = np.ones(len(omega_range)) * np.abs(gain)
    H_phase = np.zeros(len(omega_range))

    for omega in range(len(omega_range)):
        z = np.exp(1j * omega_range[omega])

        # 零点からの距離の積
        for q in zeros:
            H_mag[omega] *= np.abs(z - q)
            H_phase[omega] += np.angle(z - q)

        # 極からの距離の積
        for p in poles:
            H_mag[omega] /= np.abs(z - p)
            H_phase[omega] -= np.angle(z - p)

    return omega_range, H_mag, H_phase


# 例: 2次バターワースLPFの幾何学的周波数応答
b, a = butter(2, 0.2)   # デジタル周波数で 0.2π (正規化)
zeros, poles, gain = tf2zpk(b, a)

omega, H_mag, H_phase = freq_response_from_zpk(zeros, poles, gain)

# scipy.signal.freqz による結果と比較
w, h = freqz(b, a, worN=1024)

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

axes[0].plot(omega / np.pi, 20 * np.log10(H_mag + 1e-12), "b--",
             linewidth=2, label="Geometric (from z-plane)")
axes[0].plot(w / np.pi, 20 * np.log10(np.abs(h) + 1e-12), "r",
             linewidth=1, alpha=0.7, label="scipy.signal.freqz")
axes[0].set_xlabel("Normalized Frequency [× π rad/sample]")
axes[0].set_ylabel("Magnitude [dB]")
axes[0].set_title("Frequency Response: Geometric vs freqz")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(-60, 5)

axes[1].plot(omega / np.pi, np.degrees(np.unwrap(H_phase)), "b--",
             linewidth=2, label="Geometric")
axes[1].plot(w / np.pi, np.degrees(np.unwrap(np.angle(h))), "r",
             linewidth=1, alpha=0.7, label="scipy.signal.freqz")
axes[1].set_xlabel("Normalized Frequency [× π rad/sample]")
axes[1].set_ylabel("Phase [degrees]")
axes[1].set_title("Phase Response")
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

from scipy.signal import residuez

# H(z) = (1 + 0.5z^-1) / (1 - 0.9z^-1 + 0.2z^-2)
b = np.array([1.0, 0.5])
a = np.array([1.0, -0.9, 0.2])

# residuez は X(z) = sum(r_k / (1 - p_k z^-1)) + k の形で返す
r, p, k = residuez(b, a)

print("部分分数展開: H(z) = sum(r_k / (1 - p_k * z^-1)) + k(z)")
for i, (ri, pi) in enumerate(zip(r, p)):
    print(f"  r[{i}] = {ri:.4f},  p[{i}] = {pi:.4f}  (|p| = {abs(pi):.4f})")
print(f"  直達項 k = {k}")

# 逆Z変換（インパルス応答の解析計算）
N = 20
n = np.arange(N)
h_analytical = np.zeros(N, dtype=complex)
for ri, pi in zip(r, p):
    h_analytical += ri * pi**n
h_analytical = h_analytical.real

# scipy.signal.lfilter による数値計算と比較
from scipy.signal import lfilter, unit_impulse
impulse = unit_impulse(N)
h_numerical = lfilter(b, a, impulse)

fig, ax = plt.subplots(figsize=(10, 4))
ax.stem(n, h_numerical, linefmt="b-", markerfmt="bo", basefmt="k-", label="lfilter (numerical)")
ax.stem(n, h_analytical, linefmt="r--", markerfmt="rx", basefmt="k-", label="Partial fraction (analytical)")
ax.set_xlabel("Sample n")
ax.set_ylabel("Amplitude")
ax.set_title("Impulse Response: Partial Fraction vs Direct Computation")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print(f"\n最大誤差: {np.max(np.abs(h_numerical - h_analytical)):.2e}")

from scipy.signal import tf2sos, sos2tf, zpk2sos

# バターワース4次フィルタを各形式で表現
b, a = butter(4, 0.3)

# BA → ZPK
zeros, poles, gain = tf2zpk(b, a)

# ZPK → BA（逆変換）
b_r, a_r = zpk2tf(zeros, poles, gain)

# BA → SOS
sos = tf2sos(b, a)

# SOS → BA
b_sos, a_sos = sos2tf(sos)

print("=== 4次バターワースフィルタの表現形式比較 ===")
print(f"\nBA形式 (b): {b}")
print(f"BA形式 (a): {a}")
print(f"\n極 (poles): {poles}")
print(f"零点 (zeros): {zeros}")
print(f"ゲイン: {gain:.6f}")
print(f"\nSOS形式:\n{sos}")

# 全極の絶対値（安定性確認）
print(f"\n全極の絶対値（すべて < 1 ならば安定）:")
for i, p in enumerate(poles):
    print(f"  p[{i}] = {p:.4f},  |p| = {abs(p):.4f} {'✓' if abs(p) < 1 else '✗'}")