IIR ν•„ν„° 섀계

Previous FIR ν•„ν„° 섀계(FIR Filter Design)
πŸ“‹ Topic list
Next 닀쀑λ₯  μ‹ ν˜Έ 처리(Multirate Signal Processing)

IIR ν•„ν„° 섀계

ν•™μŠ΅ λͺ©ν‘œ

  • 고전적인 μ•„λ‚ λ‘œκ·Έ μ›ν˜• ν•„ν„°(λ²„ν„°μ›ŒμŠ€(Butterworth), 체비쇼프(Chebyshev), νƒ€μ›ν˜•(Elliptic))의 νŠΉμ„± 이해
  • μ•„λ‚ λ‘œκ·Έ-디지털 λ³€ν™˜μ„ μœ„ν•œ 쌍일차 λ³€ν™˜(bilinear transform)κ³Ό μž„νŽ„μŠ€ λΆˆλ³€ 방법(impulse invariance) λ§ˆμŠ€ν„°
  • 주파수 μ˜μ—­ κ·œκ²©μœΌλ‘œλΆ€ν„° ν•„ν„° 차수 κ²°μ • 방법 ν•™μŠ΅
  • Python의 scipy.signal λͺ¨λ“ˆμ„ μ‚¬μš©ν•œ IIR 디지털 ν•„ν„° 섀계
  • κ·Ή-영점 뢄석(pole-zero analysis)을 ν†΅ν•œ ν•„ν„° μ•ˆμ •μ„± 검증
  • ν•„ν„° μœ ν˜• 비ꡐ 및 νŠΉμ • μ‘μš©μ— μ ν•©ν•œ 섀계 선택

λͺ©μ°¨

  1. IIR ν•„ν„° 섀계 μ†Œκ°œ
  2. μ•„λ‚ λ‘œκ·Έ μ›ν˜• ν•„ν„°
  3. λ²„ν„°μ›ŒμŠ€ ν•„ν„°
  4. 체비쇼프 Type I ν•„ν„°
  5. 체비쇼프 Type II ν•„ν„°
  6. νƒ€μ›ν˜•(μ½”μ–΄) ν•„ν„°
  7. μ•„λ‚ λ‘œκ·Έ-디지털 λ³€ν™˜
  8. 쌍일차 λ³€ν™˜
  9. μž„νŽ„μŠ€ λΆˆλ³€ 방법
  10. μ™„μ „ν•œ IIR 섀계 절차
  11. μ•ˆμ •μ„± 뢄석
  12. ν•„ν„° μœ ν˜• 비ꡐ
  13. Python κ΅¬ν˜„
  14. μ—°μŠ΅ 문제

1. IIR ν•„ν„° 섀계 μ†Œκ°œ

1.1 IIR ν•„ν„° ꡬ쑰

IIR(λ¬΄ν•œ μž„νŽ„μŠ€ 응닡(Infinite Impulse Response)) ν•„ν„°λŠ” 순방ν–₯(feedforward)κ³Ό ν”Όλ“œλ°±(feedback) 경둜λ₯Ό λͺ¨λ‘ κ°–μŠ΅λ‹ˆλ‹€. 전달 ν•¨μˆ˜λŠ” λ‹€ν•­μ‹μ˜ λΉ„μœ¨λ‘œ λ‚˜νƒ€λƒ…λ‹ˆλ‹€:

$$H(z) = \frac{B(z)}{A(z)} = \frac{\sum_{k=0}^{M} b_k z^{-k}}{1 + \sum_{k=1}^{N} a_k z^{-k}}$$

μ°¨λΆ„ 방정식(difference equation)은 λ‹€μŒκ³Ό κ°™μŠ΅λ‹ˆλ‹€:

$$y[n] = \sum_{k=0}^{M} b_k x[n-k] - \sum_{k=1}^{N} a_k y[n-k]$$

1.2 섀계 접근법

FIR ν•„ν„°(이산 μ˜μ—­μ—μ„œ 직접 섀계)와 달리, IIR ν•„ν„°λŠ” 일반적으둜 λ‹€μŒκ³Ό 같이 μ„€κ³„ν•©λ‹ˆλ‹€:

  1. 잘 μ•Œλ €μ§„ νŠΉμ„±μ„ κ°–λŠ” μ•„λ‚ λ‘œκ·Έ μ›ν˜• $H_a(s)$μ—μ„œ μ‹œμž‘
  2. μ•„λ‚ λ‘œκ·Έ ν•„ν„°λ₯Ό 디지털 ν•„ν„° $H(z)$둜 λ³€ν™˜
β”Œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”
β”‚                 IIR Filter Design Pipeline                       β”‚
β”‚                                                                  β”‚
β”‚  Digital Specs     Analog Specs     Analog Filter    Digital     β”‚
β”‚  β”Œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”     β”Œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”     β”Œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”    β”Œβ”€β”€β”€β”€β”€β”€β”€β”€β”  β”‚
β”‚  β”‚ Ο‰p, Ο‰s   β”‚ ──▢ β”‚ Ξ©p, Ξ©s   β”‚ ──▢ β”‚ Ha(s)    β”‚ ──▢│ H(z)   β”‚  β”‚
β”‚  β”‚ δ₁, Ξ΄β‚‚   β”‚     β”‚ δ₁, Ξ΄β‚‚   β”‚     β”‚          β”‚    β”‚        β”‚  β”‚
β”‚  β””β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”˜     β””β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”˜     β””β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”˜    β””β”€β”€β”€β”€β”€β”€β”€β”€β”˜  β”‚
β”‚                                                                  β”‚
β”‚  Step 1:           Step 2:          Step 3:         Step 4:     β”‚
β”‚  Specify digital   Pre-warp to      Design analog   Apply       β”‚
β”‚  requirements      analog specs     prototype       BLT/IIM    β”‚
β””β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”˜

1.3 μ•„λ‚ λ‘œκ·Έ μ›ν˜•μ—μ„œ μ‹œμž‘ν•˜λŠ” 이유

  • μˆ˜μ‹­ λ…„μ˜ 이둠: λ²„ν„°μ›ŒμŠ€(Butterworth), 체비쇼프(Chebyshev), νƒ€μ›ν˜• ν•„ν„° μ„€κ³„λŠ” 잘 μ •λ¦½λ˜μ–΄ 있음
  • λ‹«νžŒ ν˜•νƒœμ˜ ν•΄: 극점(pole)/영점(zero) μœ„μΉ˜λ₯Ό ν•΄μ„μ μœΌλ‘œ 계산 κ°€λŠ₯
  • 졜적 νŠΉμ„±: 각 μœ ν˜•μ€ νŠΉμ • μ˜λ―Έμ—μ„œ 졜적 (평탄도, λ“±λ¦¬ν”Œ, μ΅œμ†Œ 차수)
  • 주파수 λ³€ν™˜: μ €μ—­ 톡과(lowpass) μ›ν˜•μ„ κ³ μ—­ 톡과(highpass), λŒ€μ—­ 톡과(bandpass), λŒ€μ—­ 저지(bandstop)둜 λ³€ν™˜ κ°€λŠ₯

2. μ•„λ‚ λ‘œκ·Έ μ›ν˜• ν•„ν„°

2.1 μ•„λ‚ λ‘œκ·Έ ν•„ν„° 규격

μ•„λ‚ λ‘œκ·Έ μ €μ—­ 톡과 ν•„ν„° κ·œκ²©μ€ λ‹€μŒμœΌλ‘œ κ΅¬μ„±λ©λ‹ˆλ‹€:

  • ν†΅κ³ΌλŒ€μ—­ 엣지 주파수(passband edge frequency) $\Omega_p$
  • μ €μ§€λŒ€μ—­ 엣지 주파수(stopband edge frequency) $\Omega_s$
  • ν†΅κ³ΌλŒ€μ—­ λ¦¬ν”Œ(passband ripple) $R_p$ (dB) λ˜λŠ” ν—ˆμš©μ˜€μ°¨ $\epsilon$
  • μ €μ§€λŒ€μ—­ 감쇠(stopband attenuation) $A_s$ (dB)

λ¦¬ν”Œ νŒŒλΌλ―Έν„° $\epsilon$κ³Ό ν†΅κ³ΌλŒ€μ—­ λ¦¬ν”Œμ˜ 관계:

$$R_p = 10\log_{10}(1 + \epsilon^2) \quad \text{dB}$$

2.2 μ •κ·œν™”λœ μ €μ—­ 톡과 μ›ν˜•

λͺ¨λ“  고전적 μ„€κ³„λŠ” ν†΅κ³ΌλŒ€μ—­ 엣지가 $\Omega_p = 1$ rad/s인 μ •κ·œν™”λœ μ €μ—­ 톡과 μ›ν˜•(normalized lowpass prototype)μ—μ„œ μ‹œμž‘ν•©λ‹ˆλ‹€. 이후 μ›ν˜•μ„ 주파수 μŠ€μΌ€μΌλ§ν•˜κ³  μ›ν•˜λŠ” ν•„ν„° μœ ν˜•μœΌλ‘œ λ³€ν™˜ν•©λ‹ˆλ‹€.

선택도 인수(selectivity factor):

$$k = \frac{\Omega_p}{\Omega_s} < 1$$

$k$κ°€ μž‘μ„μˆ˜λ‘ μ „μ΄λŒ€μ—­(transition band)이 ν†΅κ³ΌλŒ€μ—­μ— λΉ„ν•΄ λ„“μ–΄μ Έ 섀계가 μ‰¬μ›Œμ§‘λ‹ˆλ‹€.

3. λ²„ν„°μ›ŒμŠ€ ν•„ν„°

3.1 크기 응닡

λ²„ν„°μ›ŒμŠ€(Butterworth) ν•„ν„°λŠ” ν†΅κ³ΌλŒ€μ—­μ—μ„œ μ΅œλŒ€ 평탄 크기 응닡(maximally flat magnitude response)을 μ œκ³΅ν•©λ‹ˆλ‹€. 제곱 크기 응닡은:

$$|H_a(j\Omega)|^2 = \frac{1}{1 + \left(\Omega / \Omega_c\right)^{2N}}$$

μ—¬κΈ°μ„œ $N$은 ν•„ν„° 차수, $\Omega_c$λŠ” $-3$ dB 차단 μ£ΌνŒŒμˆ˜μž…λ‹ˆλ‹€.

νŠΉμ„±: - $|H_a(j\Omega)|^2$의 $2N-1$μ°¨κΉŒμ§€ λͺ¨λ“  미뢄값이 $\Omega = 0$μ—μ„œ 0 - 단쑰 κ°μ†Œ(monotonically decreasing) 크기 - 차단 μ£ΌνŒŒμˆ˜μ—μ„œ $|H_a(j\Omega_c)|^2 = 1/2$ ($-3$ dB) - λ‘€μ˜€ν”„(rolloff) 속도: μ €μ§€λŒ€μ—­μ—μ„œ $-20N$ dB/decade

3.2 극점 μœ„μΉ˜

λ²„ν„°μ›ŒμŠ€ ν•„ν„°μ˜ 극점은 $s$-ν‰λ©΄μ—μ„œ λ°˜μ§€λ¦„ $\Omega_c$인 원 μœ„μ— λ“±κ°„κ²©μœΌλ‘œ μœ„μΉ˜ν•©λ‹ˆλ‹€:

$$s_k = \Omega_c \exp\left[j\frac{\pi}{2N}(2k + N - 1)\right], \quad k = 0, 1, \ldots, 2N-1$$

μ•ˆμ •μ (인과적) ν•„ν„°λ₯Ό μœ„ν•΄ μ’Œλ°˜ν‰λ©΄(left-half-plane) 극점($\text{Re}(s_k) < 0$)만 μ„ νƒν•©λ‹ˆλ‹€.

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

def butterworth_poles(N, Omega_c=1.0):
    """Compute Butterworth filter poles in the s-plane."""
    poles = []
    for k in range(2 * N):
        s_k = Omega_c * np.exp(1j * np.pi * (2 * k + N - 1) / (2 * N))
        if np.real(s_k) < 0:  # Left-half plane only
            poles.append(s_k)
    return np.array(poles)

# Visualize pole locations for different orders
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

for ax, N in zip(axes, [2, 4, 8]):
    poles = butterworth_poles(N)
    theta = np.linspace(0, 2 * np.pi, 200)

    # Unit circle
    ax.plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.3)
    # Poles
    ax.plot(np.real(poles), np.imag(poles), 'rx', markersize=12, markeredgewidth=2)
    ax.axhline(0, color='k', linewidth=0.5)
    ax.axvline(0, color='k', linewidth=0.5)
    ax.set_xlabel('Real')
    ax.set_ylabel('Imaginary')
    ax.set_title(f'Butterworth N={N} Poles')
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-1.5, 0.5)
    ax.set_ylim(-1.5, 1.5)

plt.tight_layout()
plt.savefig('butterworth_poles.png', dpi=150)
plt.close()

3.3 차수 κ²°μ •

규격 $(\Omega_p, \Omega_s, R_p, A_s)$κ°€ 주어지면 μ΅œμ†Œ λ²„ν„°μ›ŒμŠ€ μ°¨μˆ˜λŠ”:

$$N \geq \frac{\log\left(\frac{10^{A_s/10} - 1}{10^{R_p/10} - 1}\right)}{2\log(\Omega_s / \Omega_p)}$$

def butterworth_order(Omega_p, Omega_s, Rp_dB, As_dB):
    """Compute minimum Butterworth filter order."""
    numerator = np.log10((10**(As_dB / 10) - 1) / (10**(Rp_dB / 10) - 1))
    denominator = 2 * np.log10(Omega_s / Omega_p)
    N = int(np.ceil(numerator / denominator))
    return N

# Example
N = butterworth_order(1.0, 2.0, 1.0, 40)
print(f"Minimum Butterworth order: N = {N}")
# Compare with scipy
N_scipy, Wn = signal.buttord(1.0, 2.0, 1.0, 40, analog=True)
print(f"scipy.signal.buttord: N = {N_scipy}, Wn = {Wn:.4f}")

3.4 크기 응닡 μ‹œκ°ν™”

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Different orders
orders = [1, 2, 4, 8, 16]
for N in orders:
    b, a = signal.butter(N, 1.0, analog=True)
    w, H = signal.freqs(b, a, worN=np.logspace(-1, 1, 1000))
    H_dB = 20 * np.log10(np.abs(H))

    axes[0].plot(w, H_dB, linewidth=1.5, label=f'N={N}')
    axes[1].plot(w, np.abs(H), linewidth=1.5, label=f'N={N}')

axes[0].set_xlabel('Frequency (rad/s)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('Butterworth: Magnitude (dB)')
axes[0].set_xlim(0.1, 10)
axes[0].set_ylim(-80, 5)
axes[0].axhline(-3, color='gray', linestyle=':', alpha=0.5, label='-3 dB')
axes[0].axvline(1.0, color='gray', linestyle=':', alpha=0.5)
axes[0].set_xscale('log')
axes[0].legend()
axes[0].grid(True, which='both', alpha=0.3)

axes[1].set_xlabel('Frequency (rad/s)')
axes[1].set_ylabel('|H(jΞ©)|')
axes[1].set_title('Butterworth: Magnitude (Linear)')
axes[1].set_xlim(0, 3)
axes[1].axhline(1/np.sqrt(2), color='gray', linestyle=':', alpha=0.5, label='-3 dB')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('butterworth_responses.png', dpi=150)
plt.close()

4. 체비쇼프 Type I ν•„ν„°

4.1 크기 응닡

체비쇼프(Chebyshev) Type I ν•„ν„°λŠ” λ“±λ¦¬ν”Œ ν†΅κ³ΌλŒ€μ—­(equiripple passband)κ³Ό 단쑰 μ €μ§€λŒ€μ—­μ„ κ°–μŠ΅λ‹ˆλ‹€:

$$|H_a(j\Omega)|^2 = \frac{1}{1 + \epsilon^2 T_N^2(\Omega / \Omega_p)}$$

μ—¬κΈ°μ„œ $T_N(\cdot)$λŠ” $N$μ°¨ 제1μ’… 체비쇼프 닀항식(Chebyshev polynomial of the first kind)μž…λ‹ˆλ‹€:

$$T_N(x) = \begin{cases} \cos(N \cos^{-1}(x)), & |x| \leq 1 \\ \cosh(N \cosh^{-1}(x)), & |x| > 1 \end{cases}$$

νŠΉμ„±: - λ¦¬ν”Œ $\epsilon$의 λ“±λ¦¬ν”Œ ν†΅κ³ΌλŒ€μ—­ - ν†΅κ³ΌλŒ€μ—­ λ¦¬ν”Œ: $R_p = 10\log_{10}(1 + \epsilon^2)$ dB - 동일 μ°¨μˆ˜μ—μ„œ λ²„ν„°μ›ŒμŠ€λ³΄λ‹€ κΈ‰κ²©ν•œ 전이 - 단쑰 κ°μ†Œ μ €μ§€λŒ€μ—­

4.2 체비쇼프 닀항식

처음 λͺ‡ 개의 체비쇼프 닀항식:

$$T_0(x) = 1, \quad T_1(x) = x, \quad T_2(x) = 2x^2 - 1$$

$$T_3(x) = 4x^3 - 3x, \quad T_4(x) = 8x^4 - 8x^2 + 1$$

점화식(recurrence relation): $T_{N+1}(x) = 2x \cdot T_N(x) - T_{N-1}(x)$

4.3 차수 κ²°μ •

$$N \geq \frac{\cosh^{-1}\left(\sqrt{\frac{10^{A_s/10} - 1}{10^{R_p/10} - 1}}\right)}{\cosh^{-1}(\Omega_s / \Omega_p)}$$

4.4 극점 μœ„μΉ˜

체비쇼프 Type I ν•„ν„°μ˜ 극점은 $s$-ν‰λ©΄μ—μ„œ 타원 μœ„μ— μœ„μΉ˜ν•©λ‹ˆλ‹€:

$$s_k = \sigma_k + j\omega_k$$

μ—¬κΈ°μ„œ:

$$\sigma_k = -\sinh\left(\frac{1}{N}\sinh^{-1}\left(\frac{1}{\epsilon}\right)\right) \sin\left(\frac{(2k-1)\pi}{2N}\right)$$

$$\omega_k = \cosh\left(\frac{1}{N}\sinh^{-1}\left(\frac{1}{\epsilon}\right)\right) \cos\left(\frac{(2k-1)\pi}{2N}\right)$$

$k = 1, 2, \ldots, N$에 λŒ€ν•΄.

def chebyshev1_analysis(N, Rp_dB=1.0):
    """Analyze Chebyshev Type I filter."""
    epsilon = np.sqrt(10**(Rp_dB / 10) - 1)

    # Pole locations
    poles = []
    for k in range(1, N + 1):
        theta_k = (2 * k - 1) * np.pi / (2 * N)
        sigma_k = -np.sinh(np.arcsinh(1 / epsilon) / N) * np.sin(theta_k)
        omega_k = np.cosh(np.arcsinh(1 / epsilon) / N) * np.cos(theta_k)
        poles.append(sigma_k + 1j * omega_k)

    return np.array(poles), epsilon

# Compare Butterworth and Chebyshev
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for N in [2, 4, 6, 8]:
    # Chebyshev Type I (1 dB ripple)
    b, a = signal.cheby1(N, 1.0, 1.0, analog=True)
    w, H = signal.freqs(b, a, worN=np.logspace(-1, 1, 1000))
    axes[0].plot(w, 20 * np.log10(np.abs(H)), linewidth=1.5, label=f'N={N}')

axes[0].set_xlabel('Frequency (rad/s)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('Chebyshev Type I (Rp = 1 dB)')
axes[0].set_xlim(0.1, 10)
axes[0].set_ylim(-80, 5)
axes[0].set_xscale('log')
axes[0].legend()
axes[0].grid(True, which='both', alpha=0.3)

# Pole-zero plot for N=4
poles_butter = butterworth_poles(4)
poles_cheby, _ = chebyshev1_analysis(4, 1.0)

theta = np.linspace(0, 2 * np.pi, 200)
axes[1].plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.3, label='Unit circle')
axes[1].plot(np.real(poles_butter), np.imag(poles_butter), 'bx',
             markersize=12, markeredgewidth=2, label='Butterworth')
axes[1].plot(np.real(poles_cheby), np.imag(poles_cheby), 'r+',
             markersize=12, markeredgewidth=2, label='Chebyshev I')
axes[1].set_xlabel('Real')
axes[1].set_ylabel('Imaginary')
axes[1].set_title('Pole Locations: Butterworth vs Chebyshev I (N=4)')
axes[1].set_aspect('equal')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(-1.5, 0.5)
axes[1].set_ylim(-1.5, 1.5)

plt.tight_layout()
plt.savefig('chebyshev1_analysis.png', dpi=150)
plt.close()

5. 체비쇼프 Type II ν•„ν„°

5.1 크기 응닡

체비쇼프 Type II(μ—­ 체비쇼프(inverse Chebyshev)) ν•„ν„°λŠ” ν‰νƒ„ν•œ ν†΅κ³ΌλŒ€μ—­κ³Ό λ“±λ¦¬ν”Œ μ €μ§€λŒ€μ—­μ„ κ°–μŠ΅λ‹ˆλ‹€:

$$|H_a(j\Omega)|^2 = \frac{1}{1 + \frac{1}{\epsilon^2 T_N^2(\Omega_s / \Omega)}}$$

νŠΉμ„±: - 단쑰 κ°μ†Œ(평탄) ν†΅κ³ΌλŒ€μ—­ - μ €μ§€λŒ€μ—­μ—μ„œ λ“±λ¦¬ν”Œ λ™μž‘ - 극점과 영점 λͺ¨λ‘ 쑴재 (Type I은 극점만) - 영점이 $j\Omega$ 좕에 μœ„μΉ˜ν•˜μ—¬ μ €μ§€λŒ€μ—­μ— λ…ΈμΉ˜(null) 생성

5.2 Type I과의 비ꡐ

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

N = 5

# Chebyshev Type I (1 dB ripple)
b1, a1 = signal.cheby1(N, 1.0, 1.0, analog=True)
w1, H1 = signal.freqs(b1, a1, worN=np.logspace(-1, 1, 2000))

# Chebyshev Type II (40 dB stopband attenuation)
b2, a2 = signal.cheby2(N, 40, 1.0, analog=True)
w2, H2 = signal.freqs(b2, a2, worN=np.logspace(-1, 1, 2000))

# Butterworth (same order)
b_bw, a_bw = signal.butter(N, 1.0, analog=True)
w_bw, H_bw = signal.freqs(b_bw, a_bw, worN=np.logspace(-1, 1, 2000))

# dB plot
for w, H, name, style in [(w_bw, H_bw, 'Butterworth', 'b-'),
                            (w1, H1, 'Chebyshev I', 'r-'),
                            (w2, H2, 'Chebyshev II', 'g-')]:
    axes[0].plot(w, 20 * np.log10(np.abs(H) + 1e-15), style,
                 linewidth=1.5, label=name)

axes[0].set_xlabel('Frequency (rad/s)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title(f'Filter Comparison (N={N})')
axes[0].set_xlim(0.1, 10)
axes[0].set_ylim(-80, 5)
axes[0].set_xscale('log')
axes[0].legend()
axes[0].grid(True, which='both', alpha=0.3)

# Passband detail (linear scale)
for w, H, name, style in [(w_bw, H_bw, 'Butterworth', 'b-'),
                            (w1, H1, 'Chebyshev I', 'r-'),
                            (w2, H2, 'Chebyshev II', 'g-')]:
    axes[1].plot(w, np.abs(H), style, linewidth=1.5, label=name)

axes[1].set_xlabel('Frequency (rad/s)')
axes[1].set_ylabel('|H(jΞ©)|')
axes[1].set_title('Passband Detail')
axes[1].set_xlim(0, 2)
axes[1].set_ylim(0, 1.2)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('chebyshev2_comparison.png', dpi=150)
plt.close()

6. νƒ€μ›ν˜•(μ½”μ–΄) ν•„ν„°

6.1 크기 응닡

νƒ€μ›ν˜•(Elliptic, Cauer) ν•„ν„°λŠ” ν†΅κ³ΌλŒ€μ—­κ³Ό μ €μ§€λŒ€μ—­ λͺ¨λ‘μ— λ¦¬ν”Œμ„ λΆ„μ‚°μ‹œμΌœ 주어진 κ·œκ²©μ— λŒ€ν•΄ μ΅œμ†Œ 차수λ₯Ό λ‹¬μ„±ν•©λ‹ˆλ‹€:

$$|H_a(j\Omega)|^2 = \frac{1}{1 + \epsilon^2 R_N^2(\Omega / \Omega_p)}$$

μ—¬κΈ°μ„œ $R_N(\cdot)$은 유리 체비쇼프(μ•Όμ½”λΉ„ 타원(Jacobian elliptic)) ν•¨μˆ˜μž…λ‹ˆλ‹€.

νŠΉμ„±: - ν†΅κ³ΌλŒ€μ—­κ³Ό μ €μ§€λŒ€μ—­ λͺ¨λ‘ λ“±λ¦¬ν”Œ - 주어진 μ°¨μˆ˜μ— λŒ€ν•΄ κ°€μž₯ κΈ‰κ²©ν•œ 전이 - 차수 λŒ€ 규격 λ©΄μ—μ„œ κ°€μž₯ 효율적 - 극점과 영점 λͺ¨λ‘ 쑴재

6.2 차수 μš°μœ„

λ™μΌν•œ κ·œκ²©μ— λŒ€ν•΄ ν•„μš”ν•œ ν•„ν„° μ°¨μˆ˜λŠ”:

$$N_\text{Elliptic} \leq N_\text{Chebyshev} \leq N_\text{Butterworth}$$

def compare_filter_orders(Rp_dB, As_dB, wp, ws):
    """Compare required orders for different filter types."""
    # Butterworth
    N_butter, Wn_butter = signal.buttord(wp, ws, Rp_dB, As_dB, analog=True)

    # Chebyshev Type I
    N_cheby1, Wn_cheby1 = signal.cheb1ord(wp, ws, Rp_dB, As_dB, analog=True)

    # Chebyshev Type II
    N_cheby2, Wn_cheby2 = signal.cheb2ord(wp, ws, Rp_dB, As_dB, analog=True)

    # Elliptic
    N_ellip, Wn_ellip = signal.ellipord(wp, ws, Rp_dB, As_dB, analog=True)

    print(f"Specifications: Rp={Rp_dB} dB, As={As_dB} dB, wp={wp}, ws={ws}")
    print(f"{'Filter Type':<20} {'Order N':<10} {'Natural Freq':>12}")
    print("-" * 45)
    print(f"{'Butterworth':<20} {N_butter:<10} {Wn_butter:>12.4f}")
    print(f"{'Chebyshev I':<20} {N_cheby1:<10} {Wn_cheby1:>12.4f}")
    print(f"{'Chebyshev II':<20} {N_cheby2:<10} {Wn_cheby2:>12.4f}")
    print(f"{'Elliptic':<20} {N_ellip:<10} {Wn_ellip:>12.4f}")

    return N_butter, N_cheby1, N_cheby2, N_ellip

orders = compare_filter_orders(1.0, 60, 1.0, 1.5)

6.3 νƒ€μ›ν˜• ν•„ν„° μ‹œκ°ν™”

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

# Same specifications, different filter types
Rp = 1.0  # dB
As = 60   # dB
wp = 2 * np.pi * 1000  # rad/s
ws = 2 * np.pi * 1300  # rad/s

filters = [
    ('Butterworth', signal.buttord, signal.butter),
    ('Chebyshev I', signal.cheb1ord, signal.cheby1),
    ('Chebyshev II', signal.cheb2ord, signal.cheby2),
    ('Elliptic', signal.ellipord, signal.ellip),
]

for ax, (name, ord_func, design_func) in zip(axes.flat, filters):
    N, Wn = ord_func(wp, ws, Rp, As, analog=True)

    if 'Chebyshev I' in name:
        b, a = design_func(N, Rp, Wn, analog=True)
    elif 'Chebyshev II' in name:
        b, a = design_func(N, As, Wn, analog=True)
    elif 'Elliptic' in name:
        b, a = design_func(N, Rp, As, Wn, analog=True)
    else:
        b, a = design_func(N, Wn, analog=True)

    w, H = signal.freqs(b, a, worN=np.linspace(0, 3 * ws, 5000))
    H_dB = 20 * np.log10(np.abs(H) + 1e-15)

    ax.plot(w / (2 * np.pi), H_dB, 'b-', linewidth=1.5)
    ax.axhline(-Rp, color='g', linestyle='--', alpha=0.5, label=f'-{Rp} dB')
    ax.axhline(-As, color='r', linestyle='--', alpha=0.5, label=f'-{As} dB')
    ax.axvline(wp / (2 * np.pi), color='g', linestyle=':', alpha=0.5)
    ax.axvline(ws / (2 * np.pi), color='r', linestyle=':', alpha=0.5)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Magnitude (dB)')
    ax.set_title(f'{name} (N={N})')
    ax.set_ylim(-80, 5)
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

plt.suptitle(f'Filter Type Comparison (Rp={Rp} dB, As={As} dB)', fontsize=13)
plt.tight_layout()
plt.savefig('filter_type_comparison.png', dpi=150)
plt.close()

7. μ•„λ‚ λ‘œκ·Έ-디지털 λ³€ν™˜

7.1 방법 κ°œμš”

μ•„λ‚ λ‘œκ·Έ ν•„ν„° $H_a(s)$λ₯Ό 디지털 ν•„ν„° $H(z)$둜 λ³€ν™˜ν•˜λŠ” μ„Έ 가지 μ£Όμš” 방법:

β”Œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”
β”‚            Analog-to-Digital Conversion Methods                  β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”¬β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Method           β”‚ Mapping s β†’ z                                β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Bilinear         β”‚ s = (2/T)(z-1)/(z+1)                        β”‚
β”‚ Transform        β”‚ - No aliasing                                β”‚
β”‚                  β”‚ - Frequency warping (correctable)            β”‚
β”‚                  β”‚ - Most widely used                           β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Impulse          β”‚ h[n] = T Β· h_a(nT)                           β”‚
β”‚ Invariance       β”‚ - Preserves impulse response shape           β”‚
β”‚                  β”‚ - Aliasing in frequency domain               β”‚
β”‚                  β”‚ - Only for bandlimited filters (LP, BP)      β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Matched          β”‚ Map poles: s_k β†’ z_k = e^(s_k T)             β”‚
β”‚ Z-Transform      β”‚ Map zeros: same mapping                      β”‚
β”‚                  β”‚ - Simple but no formal optimality            β”‚
β””β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”΄β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”˜

8. 쌍일차 λ³€ν™˜

8.1 맀핑

쌍일차 λ³€ν™˜(bilinear transform, BLT)은 $s$-평면 전체λ₯Ό $z$-ν‰λ©΄μœΌλ‘œ λ‹€μŒκ³Ό 같이 λ§€ν•‘ν•©λ‹ˆλ‹€:

$$s = \frac{2}{T} \cdot \frac{z - 1}{z + 1}$$

λ˜λŠ” λ™λ“±ν•˜κ²Œ:

$$z = \frac{1 + (T/2)s}{1 - (T/2)s}$$

μ—¬κΈ°μ„œ $T$λŠ” μƒ˜ν”Œλ§ μ£ΌκΈ°μž…λ‹ˆλ‹€.

8.2 μ£Όμš” νŠΉμ„±

μ•ˆμ •μ„± 보쑴: $s$-ν‰λ©΄μ˜ μ’Œλ°˜ν‰λ©΄μ΄ $z$-ν‰λ©΄μ˜ λ‹¨μœ„μ› λ‚΄λΆ€λ‘œ λ§€ν•‘λ©λ‹ˆλ‹€. 이λ₯Ό 톡해 μ•ˆμ •μ μΈ μ•„λ‚ λ‘œκ·Έ ν•„ν„°λŠ” 항상 μ•ˆμ •μ μΈ 디지털 ν•„ν„°λ₯Ό μƒμ„±ν•©λ‹ˆλ‹€.

주파수 맀핑: ν—ˆμˆ˜μΆ•($s = j\Omega$)은 λΉ„μ„ ν˜• 주파수 μ›Œν•‘(frequency warping)으둜 λ‹¨μœ„μ›($z = e^{j\omega}$)으둜 λ§€ν•‘λ©λ‹ˆλ‹€:

$$\omega = 2\arctan\left(\frac{\Omega T}{2}\right) \quad \Leftrightarrow \quad \Omega = \frac{2}{T}\tan\left(\frac{\omega}{2}\right)$$

8.3 주파수 μ›Œν•‘

Bilinear Transform Frequency Warping:

Ξ© (analog)     Ο‰ (digital)
   ↑              ↑
   β”‚     β•±       π│─────────────────╱
   β”‚   β•±          β”‚              β•±
   β”‚  β•±           β”‚            β•±
   β”‚ β•±            β”‚         β•±      Nonlinear
   β”‚β•±             β”‚      β•±        compression
 0 β”œβ”€β”€β”€β†’ Ο‰     0 β”œβ”€β”€β”€β•±β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β†’ Ξ©
   0    Ο€         0

- Low frequencies: Ξ© β‰ˆ Ο‰ (nearly linear)
- High frequencies: compressed into [0, Ο€]
- Ξ© = ∞ maps to Ο‰ = Ο€

# Visualize frequency warping
T = 1.0  # Sampling period
omega = np.linspace(0, np.pi - 0.01, 1000)
Omega = (2 / T) * np.tan(omega / 2)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].plot(Omega, omega / np.pi, 'b-', linewidth=2, label='BLT mapping')
axes[0].plot(Omega, Omega * T / np.pi, 'r--', linewidth=1.5, label='Linear (ideal)')
axes[0].set_xlabel('Analog Frequency Ξ© (rad/s)')
axes[0].set_ylabel('Digital Frequency Ο‰/Ο€')
axes[0].set_title('Bilinear Transform: Frequency Warping')
axes[0].set_xlim(0, 20)
axes[0].set_ylim(0, 1.1)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Warping effect at different frequencies
freq_analog = np.array([0.1, 0.5, 1.0, 2.0, 5.0, 10.0])
freq_digital = 2 * np.arctan(freq_analog * T / 2)

axes[1].bar(range(len(freq_analog)),
            freq_digital / (freq_analog * T) * 100 - 100,
            tick_label=[f'{f:.1f}' for f in freq_analog])
axes[1].set_xlabel('Analog Frequency Ξ© (rad/s)')
axes[1].set_ylabel('Warping Error (%)')
axes[1].set_title('Frequency Warping Error')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('bilinear_warping.png', dpi=150)
plt.close()

8.4 사전 μ›Œν•‘(Pre-warping)

μž„κ³„ μ£ΌνŒŒμˆ˜κ°€ μ˜¬λ°”λ₯΄κ²Œ λ§€ν•‘λ˜λ„λ‘ ν•˜κΈ° μœ„ν•΄ μ•„λ‚ λ‘œκ·Έ κ·œκ²©μ„ 사전 μ›Œν•‘(pre-warp)ν•©λ‹ˆλ‹€:

$$\Omega_p' = \frac{2}{T}\tan\left(\frac{\omega_p}{2}\right), \quad \Omega_s' = \frac{2}{T}\tan\left(\frac{\omega_s}{2}\right)$$

그런 λ‹€μŒ 사전 μ›Œν•‘λœ 주파수둜 μ•„λ‚ λ‘œκ·Έ ν•„ν„°λ₯Ό μ„€κ³„ν•˜κ³  쌍일차 λ³€ν™˜μ„ μ μš©ν•©λ‹ˆλ‹€.

8.5 단계별 BLT 섀계

def iir_design_blt(wp_digital, ws_digital, Rp_dB, As_dB, fs, ftype='butter'):
    """
    Complete IIR design using bilinear transform with pre-warping.

    Parameters:
        wp_digital: digital passband edge (Hz)
        ws_digital: digital stopband edge (Hz)
        Rp_dB: passband ripple (dB)
        As_dB: stopband attenuation (dB)
        fs: sampling frequency (Hz)
        ftype: 'butter', 'cheby1', 'cheby2', or 'ellip'

    Returns:
        b, a: digital filter coefficients
        N: filter order
    """
    T = 1 / fs

    # Step 1: Pre-warp digital frequencies to analog
    Omega_p = 2 * fs * np.tan(np.pi * wp_digital / fs)
    Omega_s = 2 * fs * np.tan(np.pi * ws_digital / fs)

    print(f"Pre-warped frequencies: Ξ©p = {Omega_p:.2f}, Ξ©s = {Omega_s:.2f} rad/s")

    # Step 2: Determine order
    if ftype == 'butter':
        N, Wn = signal.buttord(Omega_p, Omega_s, Rp_dB, As_dB, analog=True)
        # Step 3: Design analog prototype
        ba, aa = signal.butter(N, Wn, analog=True)
    elif ftype == 'cheby1':
        N, Wn = signal.cheb1ord(Omega_p, Omega_s, Rp_dB, As_dB, analog=True)
        ba, aa = signal.cheby1(N, Rp_dB, Wn, analog=True)
    elif ftype == 'cheby2':
        N, Wn = signal.cheb2ord(Omega_p, Omega_s, Rp_dB, As_dB, analog=True)
        ba, aa = signal.cheby2(N, As_dB, Wn, analog=True)
    elif ftype == 'ellip':
        N, Wn = signal.ellipord(Omega_p, Omega_s, Rp_dB, As_dB, analog=True)
        ba, aa = signal.ellip(N, Rp_dB, As_dB, Wn, analog=True)

    print(f"Analog filter order: N = {N}")

    # Step 4: Bilinear transform
    b, a = signal.bilinear(ba, aa, fs)

    return b, a, N

# Example
fs = 8000
wp = 1000  # Hz
ws = 1500  # Hz
Rp = 1.0   # dB
As = 60    # dB

b, a, N = iir_design_blt(wp, ws, Rp, As, fs, ftype='ellip')
print(f"\nDigital filter coefficients:")
print(f"b = {b}")
print(f"a = {a}")

# Verify
w, H = signal.freqz(b, a, worN=4096, fs=fs)
H_dB = 20 * np.log10(np.abs(H) + 1e-15)

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

# Magnitude
axes[0, 0].plot(w, H_dB, 'b-', linewidth=1.5)
axes[0, 0].axhline(-Rp, color='g', linestyle='--', alpha=0.5, label=f'-{Rp} dB')
axes[0, 0].axhline(-As, color='r', linestyle='--', alpha=0.5, label=f'-{As} dB')
axes[0, 0].axvline(wp, color='g', linestyle=':', alpha=0.5)
axes[0, 0].axvline(ws, color='r', linestyle=':', alpha=0.5)
axes[0, 0].set_xlabel('Frequency (Hz)')
axes[0, 0].set_ylabel('Magnitude (dB)')
axes[0, 0].set_title(f'Elliptic IIR (N={N}): Magnitude')
axes[0, 0].set_ylim(-80, 5)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Phase
axes[0, 1].plot(w, np.unwrap(np.angle(H)) * 180 / np.pi, 'r-', linewidth=1.5)
axes[0, 1].set_xlabel('Frequency (Hz)')
axes[0, 1].set_ylabel('Phase (degrees)')
axes[0, 1].set_title('Phase Response')
axes[0, 1].grid(True, alpha=0.3)

# Group delay
w_gd, gd = signal.group_delay((b, a), w=4096, fs=fs)
axes[1, 0].plot(w_gd, gd, 'g-', linewidth=1.5)
axes[1, 0].set_xlabel('Frequency (Hz)')
axes[1, 0].set_ylabel('Group Delay (samples)')
axes[1, 0].set_title('Group Delay (Non-constant!)')
axes[1, 0].grid(True, alpha=0.3)

# Pole-zero plot
z_zeros, z_poles, k = signal.tf2zpk(b, a)
theta = np.linspace(0, 2 * np.pi, 200)
axes[1, 1].plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.3)
axes[1, 1].plot(np.real(z_zeros), np.imag(z_zeros), 'bo', markersize=10, label='Zeros')
axes[1, 1].plot(np.real(z_poles), np.imag(z_poles), 'rx', markersize=10,
                markeredgewidth=2, label='Poles')
axes[1, 1].set_xlabel('Real')
axes[1, 1].set_ylabel('Imaginary')
axes[1, 1].set_title('Pole-Zero Plot')
axes[1, 1].set_aspect('equal')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.suptitle('IIR Filter Design via Bilinear Transform', fontsize=13)
plt.tight_layout()
plt.savefig('iir_blt_design.png', dpi=150)
plt.close()

9. μž„νŽ„μŠ€ λΆˆλ³€ 방법

9.1 원리

μž„νŽ„μŠ€ λΆˆλ³€(impulse invariance) 방법은 디지털 ν•„ν„°μ˜ μž„νŽ„μŠ€ 응닡을 μ•„λ‚ λ‘œκ·Έ μž„νŽ„μŠ€ μ‘λ‹΅μ˜ μƒ˜ν”Œκ³Ό λ™μΌν•˜κ²Œ μ„€μ •ν•©λ‹ˆλ‹€:

$$h[n] = T \cdot h_a(nT)$$

μ—¬κΈ°μ„œ $h_a(t)$λŠ” μ•„λ‚ λ‘œκ·Έ μž„νŽ„μŠ€ 응닡, $T$λŠ” μƒ˜ν”Œλ§ μ£ΌκΈ°μž…λ‹ˆλ‹€.

9.2 맀핑

μ•„λ‚ λ‘œκ·Έ ν•„ν„°κ°€ λΆ€λΆ„ λΆ„μˆ˜ μ „κ°œ(partial fraction expansion)λ₯Ό κ°€μ§ˆ λ•Œ:

$$H_a(s) = \sum_{k=1}^{N} \frac{A_k}{s - s_k}$$

디지털 ν•„ν„°λŠ”:

$$H(z) = T \sum_{k=1}^{N} \frac{A_k}{1 - e^{s_k T} z^{-1}}$$

각 μ•„λ‚ λ‘œκ·Έ 극점 $s_k$λŠ” 디지털 극점 $z_k = e^{s_k T}$둜 λ§€ν•‘λ©λ‹ˆλ‹€.

9.3 에일리어싱 문제

디지털 ν•„ν„°μ˜ 주파수 응닡은 μ•„λ‚ λ‘œκ·Έ μ‘λ‹΅μ˜ μ΄λ™λœ λ³΅μ‚¬λ³Έλ“€μ˜ ν•©μž…λ‹ˆλ‹€:

$$H(e^{j\omega}) = \frac{1}{T} \sum_{k=-\infty}^{\infty} H_a\left(j\frac{\omega - 2\pi k}{T}\right)$$

이둜 인해 $H_a(j\Omega)$κ°€ λŒ€μ—­ μ œν•œ(bandlimited)λ˜μ§€ μ•ŠμœΌλ©΄ 에일리어싱(aliasing)이 λ°œμƒν•˜λ―€λ‘œ, μž„νŽ„μŠ€ λΆˆλ³€ 방법은 κ³ μ—­ 톡과 및 λŒ€μ—­ 저지 필터에 μ ν•©ν•˜μ§€ μ•ŠμŠ΅λ‹ˆλ‹€.

9.4 κ΅¬ν˜„ 예제

def impulse_invariance(ba, aa, fs):
    """
    Apply impulse invariance method to convert analog filter to digital.

    Parameters:
        ba, aa: analog filter coefficients (transfer function form)
        fs: sampling frequency

    Returns:
        bd, ad: digital filter coefficients
    """
    T = 1 / fs

    # Convert to zeros, poles, gain form
    z_a, p_a, k_a = signal.tf2zpk(ba, aa)

    # Partial fraction expansion
    residues, poles, direct = signal.residue(ba, aa)

    # Map poles: z_k = exp(s_k * T)
    z_poles = np.exp(poles * T)

    # Construct digital transfer function
    # H(z) = T * sum(A_k / (1 - e^(s_k*T) * z^-1))
    bd = np.array([0.0])
    ad = np.array([1.0])

    for A_k, z_k in zip(residues, z_poles):
        # Each term: T * A_k / (1 - z_k * z^-1)
        bd_k = np.array([T * A_k])
        ad_k = np.array([1, -z_k])

        # Combine fractions
        bd_new = np.convolve(bd, ad_k) + np.convolve(bd_k, ad)
        ad_new = np.convolve(ad, ad_k)
        bd = bd_new
        ad = ad_new

    # Ensure real coefficients
    bd = np.real(bd)
    ad = np.real(ad)

    return bd, ad

# Compare BLT and impulse invariance for a Butterworth lowpass
fs = 8000
N = 4
fc = 1000  # Cutoff frequency (Hz)

# Analog prototype
Omega_c = 2 * np.pi * fc
ba, aa = signal.butter(N, Omega_c, analog=True)

# Method 1: Bilinear transform
bd_blt, ad_blt = signal.bilinear(ba, aa, fs)

# Method 2: Impulse invariance
bd_ii, ad_ii = impulse_invariance(ba, aa, fs)

# Compare frequency responses
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

w_blt, H_blt = signal.freqz(bd_blt, ad_blt, worN=4096, fs=fs)
w_ii, H_ii = signal.freqz(bd_ii, ad_ii, worN=4096, fs=fs)

# Also plot the analog response for reference
w_a, H_a = signal.freqs(ba, aa, worN=np.linspace(0, 2 * np.pi * fs / 2, 4096))

axes[0].plot(w_a / (2 * np.pi), 20 * np.log10(np.abs(H_a) + 1e-15),
             'k--', linewidth=1.5, label='Analog', alpha=0.5)
axes[0].plot(w_blt, 20 * np.log10(np.abs(H_blt) + 1e-15),
             'b-', linewidth=1.5, label='Bilinear Transform')
axes[0].plot(w_ii, 20 * np.log10(np.abs(H_ii) + 1e-15),
             'r-', linewidth=1.5, label='Impulse Invariance')
axes[0].set_xlabel('Frequency (Hz)')
axes[0].set_ylabel('Magnitude (dB)')
axes[0].set_title('BLT vs Impulse Invariance')
axes[0].set_ylim(-80, 5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Near Nyquist detail
axes[1].plot(w_blt, 20 * np.log10(np.abs(H_blt) + 1e-15),
             'b-', linewidth=1.5, label='BLT')
axes[1].plot(w_ii, 20 * np.log10(np.abs(H_ii) + 1e-15),
             'r-', linewidth=1.5, label='Impulse Invariance')
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_ylabel('Magnitude (dB)')
axes[1].set_title('Detail Near Nyquist (aliasing visible)')
axes[1].set_xlim(2000, 4000)
axes[1].set_ylim(-60, -20)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('impulse_invariance.png', dpi=150)
plt.close()

10. μ™„μ „ν•œ IIR 섀계 절차

10.1 단계별 절차

β”Œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”
β”‚                  IIR Filter Design Steps                         β”‚
β”‚                                                                  β”‚
β”‚  1. Specify digital requirements:                                β”‚
β”‚     - Passband/stopband edges (Ο‰p, Ο‰s or fp, fs)               β”‚
β”‚     - Passband ripple Rp (dB)                                    β”‚
β”‚     - Stopband attenuation As (dB)                               β”‚
β”‚                                                                  β”‚
β”‚  2. Choose filter type:                                          β”‚
β”‚     - Butterworth: maximally flat                                β”‚
β”‚     - Chebyshev I: equiripple passband                           β”‚
β”‚     - Chebyshev II: equiripple stopband                          β”‚
β”‚     - Elliptic: minimum order                                    β”‚
β”‚                                                                  β”‚
β”‚  3. Choose conversion method (usually BLT)                       β”‚
β”‚                                                                  β”‚
β”‚  4. Pre-warp critical frequencies                                β”‚
β”‚                                                                  β”‚
β”‚  5. Determine minimum analog filter order                        β”‚
β”‚                                                                  β”‚
β”‚  6. Design analog prototype                                      β”‚
β”‚                                                                  β”‚
β”‚  7. Apply analog-to-digital transformation                       β”‚
β”‚                                                                  β”‚
β”‚  8. Verify specifications                                        β”‚
β”‚     - Check passband ripple                                      β”‚
β”‚     - Check stopband attenuation                                 β”‚
β”‚     - Check stability (all poles inside unit circle)             β”‚
β”‚                                                                  β”‚
β”‚  9. Implement (Direct Form II, SOS cascade, etc.)                β”‚
β””β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”˜

10.2 scipy.signal 직접 섀계 ν•¨μˆ˜ μ‚¬μš©

SciPyλŠ” 사전 μ›Œν•‘μ„ λ‚΄λΆ€μ μœΌλ‘œ μ²˜λ¦¬ν•˜λŠ” κ³ μˆ˜μ€€ ν•¨μˆ˜λ₯Ό μ œκ³΅ν•©λ‹ˆλ‹€:

def design_all_types(fs, f_pass, f_stop, Rp, As):
    """Design IIR lowpass filter using all four types."""
    results = {}

    # Direct digital design (scipy handles pre-warping internally)
    # Butterworth
    N_b, Wn_b = signal.buttord(f_pass, f_stop, Rp, As, fs=fs)
    b_b, a_b = signal.butter(N_b, Wn_b, fs=fs)
    results['Butterworth'] = (b_b, a_b, N_b)

    # Chebyshev Type I
    N_c1, Wn_c1 = signal.cheb1ord(f_pass, f_stop, Rp, As, fs=fs)
    b_c1, a_c1 = signal.cheby1(N_c1, Rp, Wn_c1, fs=fs)
    results['Chebyshev I'] = (b_c1, a_c1, N_c1)

    # Chebyshev Type II
    N_c2, Wn_c2 = signal.cheb2ord(f_pass, f_stop, Rp, As, fs=fs)
    b_c2, a_c2 = signal.cheby2(N_c2, As, Wn_c2, fs=fs)
    results['Chebyshev II'] = (b_c2, a_c2, N_c2)

    # Elliptic
    N_e, Wn_e = signal.ellipord(f_pass, f_stop, Rp, As, fs=fs)
    b_e, a_e = signal.ellip(N_e, Rp, As, Wn_e, fs=fs)
    results['Elliptic'] = (b_e, a_e, N_e)

    return results

# Design and compare
fs = 16000
f_pass = 2000
f_stop = 2500
Rp = 0.5
As = 60

results = design_all_types(fs, f_pass, f_stop, Rp, As)

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

for ax, (name, (b, a, N)) in zip(axes.flat, results.items()):
    w, H = signal.freqz(b, a, worN=4096, fs=fs)
    H_dB = 20 * np.log10(np.abs(H) + 1e-15)

    ax.plot(w, H_dB, 'b-', linewidth=1.5)
    ax.axhline(-Rp, color='g', linestyle='--', alpha=0.5, label=f'Rp = -{Rp} dB')
    ax.axhline(-As, color='r', linestyle='--', alpha=0.5, label=f'As = -{As} dB')
    ax.axvline(f_pass, color='g', linestyle=':', alpha=0.5)
    ax.axvline(f_stop, color='r', linestyle=':', alpha=0.5)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Magnitude (dB)')
    ax.set_title(f'{name} (N={N})')
    ax.set_ylim(-80, 5)
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

plt.suptitle(f'IIR Filter Comparison (fp={f_pass} Hz, fs_stop={f_stop} Hz)', fontsize=13)
plt.tight_layout()
plt.savefig('iir_all_types.png', dpi=150)
plt.close()

# Print summary
print(f"\n{'Filter Type':<15} {'Order':<8} {'Actual Rp (dB)':<16} {'Actual As (dB)':<16}")
print("-" * 55)
for name, (b, a, N) in results.items():
    w, H = signal.freqz(b, a, worN=4096, fs=fs)
    H_mag = np.abs(H)
    f = w

    # Passband ripple
    pass_idx = f <= f_pass
    Rp_actual = -20 * np.log10(np.min(H_mag[pass_idx]))

    # Stopband attenuation
    stop_idx = f >= f_stop
    As_actual = -20 * np.log10(np.max(H_mag[stop_idx]))

    print(f"{name:<15} {N:<8} {Rp_actual:<16.4f} {As_actual:<16.4f}")

11. μ•ˆμ •μ„± 뢄석

11.1 μ•ˆμ •μ„± κΈ°μ€€

디지털 IIR ν•„ν„°λŠ” $H(z)$의 λͺ¨λ“  극점이 λ‹¨μœ„μ› 내뢀에 μ—„κ²©νžˆ μœ„μΉ˜ν•  λ•Œλ§Œ μ•ˆμ •(stable)ν•©λ‹ˆλ‹€:

$$|z_k| < 1, \quad \forall k$$

11.2 수치적 μ•ˆμ •μ„±μ„ μœ„ν•œ 2μ°¨ μ„Ήμ…˜(SOS)

직접 ν˜•μ‹(direct form)으둜 κ΅¬ν˜„λœ κ³ μ°¨ IIR ν•„ν„°λŠ” κ³„μˆ˜ μ–‘μžν™”(coefficient quantization) 및 수치 정밀도 λ¬Έμ œκ°€ λ°œμƒν•  수 μžˆμŠ΅λ‹ˆλ‹€. 해결책은 2μ°¨ μ„Ήμ…˜(second-order sections, biquads)의 쒅속 μ—°κ²°(cascade)둜 κ΅¬ν˜„ν•˜λŠ” κ²ƒμž…λ‹ˆλ‹€:

$$H(z) = \prod_{i=1}^{L} \frac{b_{0i} + b_{1i}z^{-1} + b_{2i}z^{-2}}{1 + a_{1i}z^{-1} + a_{2i}z^{-2}}$$

def stability_analysis(b, a, title=""):
    """Analyze IIR filter stability."""
    # Poles and zeros
    zeros, poles, gain = signal.tf2zpk(b, a)

    # Check stability
    pole_mags = np.abs(poles)
    is_stable = np.all(pole_mags < 1.0)
    max_pole_mag = np.max(pole_mags) if len(poles) > 0 else 0

    print(f"Filter: {title}")
    print(f"  Number of poles: {len(poles)}")
    print(f"  Number of zeros: {len(zeros)}")
    print(f"  Maximum pole magnitude: {max_pole_mag:.6f}")
    print(f"  Stable: {is_stable}")
    print(f"  Stability margin: {1.0 - max_pole_mag:.6f}")

    # Pole-zero plot
    fig, ax = plt.subplots(figsize=(8, 8))
    theta = np.linspace(0, 2 * np.pi, 200)
    ax.plot(np.cos(theta), np.sin(theta), 'k-', alpha=0.3, label='Unit circle')
    ax.plot(np.real(zeros), np.imag(zeros), 'bo', markersize=10,
            label=f'Zeros ({len(zeros)})')
    ax.plot(np.real(poles), np.imag(poles), 'rx', markersize=10,
            markeredgewidth=2, label=f'Poles ({len(poles)})')

    # Color poles by stability
    for p in poles:
        color = 'green' if np.abs(p) < 1 else 'red'
        circle = plt.Circle((np.real(p), np.imag(p)), 0.03,
                            fill=True, color=color, alpha=0.3)
        ax.add_patch(circle)

    ax.set_xlabel('Real')
    ax.set_ylabel('Imaginary')
    ax.set_title(f'Pole-Zero Plot: {title}\n(Stable: {is_stable}, '
                 f'Max |pole| = {max_pole_mag:.4f})')
    ax.set_aspect('equal')
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-1.5, 1.5)
    ax.set_ylim(-1.5, 1.5)

    plt.tight_layout()
    return fig, is_stable

# Design a high-order filter and check stability
N = 12
b_high, a_high = signal.butter(N, 0.3)  # 12th-order Butterworth
fig, stable = stability_analysis(b_high, a_high, f"Butterworth N={N}")
plt.savefig('stability_analysis.png', dpi=150)
plt.close()

# Convert to SOS form for better numerical stability
sos = signal.tf2sos(b_high, a_high)
print(f"\nSOS representation: {sos.shape[0]} second-order sections")
print(f"SOS coefficients:\n{sos}")

11.3 수치적 비ꡐ: tf vs SOS

def compare_tf_vs_sos(order=20):
    """Compare direct form vs SOS implementation for high-order filter."""
    # Design high-order filter
    b, a = signal.butter(order, 0.1)
    sos = signal.butter(order, 0.1, output='sos')

    # Test signal
    np.random.seed(42)
    x = np.random.randn(1000)

    # Filter using tf form
    y_tf = signal.lfilter(b, a, x)

    # Filter using SOS form
    y_sos = signal.sosfilt(sos, x)

    # Compare
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))

    # Frequency response comparison
    w_tf, H_tf = signal.freqz(b, a, worN=4096)
    w_sos, H_sos = signal.sosfreqz(sos, worN=4096)

    axes[0, 0].plot(w_tf / np.pi, 20 * np.log10(np.abs(H_tf) + 1e-15),
                     'b-', linewidth=1.5, label='tf form')
    axes[0, 0].plot(w_sos / np.pi, 20 * np.log10(np.abs(H_sos) + 1e-15),
                     'r--', linewidth=1.5, label='SOS form')
    axes[0, 0].set_xlabel('Normalized Frequency (Γ—Ο€)')
    axes[0, 0].set_ylabel('Magnitude (dB)')
    axes[0, 0].set_title('Frequency Response')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)

    # Output comparison
    axes[0, 1].plot(y_tf[:200], 'b-', alpha=0.7, label='tf form')
    axes[0, 1].plot(y_sos[:200], 'r--', alpha=0.7, label='SOS form')
    axes[0, 1].set_xlabel('Sample')
    axes[0, 1].set_ylabel('Output')
    axes[0, 1].set_title('Output Signal')
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3)

    # Difference
    diff = np.abs(y_tf - y_sos)
    axes[1, 0].semilogy(diff, 'k-')
    axes[1, 0].set_xlabel('Sample')
    axes[1, 0].set_ylabel('|y_tf - y_sos|')
    axes[1, 0].set_title(f'Numerical Difference (max: {np.max(diff):.2e})')
    axes[1, 0].grid(True, alpha=0.3)

    # Pole magnitudes
    zeros_tf, poles_tf, _ = signal.tf2zpk(b, a)
    axes[1, 1].plot(np.abs(poles_tf), 'rx', markersize=10, markeredgewidth=2)
    axes[1, 1].axhline(1.0, color='k', linestyle='--', alpha=0.5, label='Unit circle')
    axes[1, 1].set_xlabel('Pole Index')
    axes[1, 1].set_ylabel('|pole|')
    axes[1, 1].set_title(f'Pole Magnitudes (N={order})')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)

    plt.suptitle(f'Direct Form vs SOS (Order={order})', fontsize=13)
    plt.tight_layout()
    plt.savefig('tf_vs_sos.png', dpi=150)
    plt.close()

compare_tf_vs_sos(order=20)

12. ν•„ν„° μœ ν˜• 비ꡐ

12.1 μš”μ•½ ν‘œ

β”Œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”
β”‚                     IIR Filter Type Comparison                          β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”¬β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”¬β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”¬β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”¬β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚              β”‚ Butterworth β”‚ Chebyshev Iβ”‚ Chebyshev IIβ”‚ Elliptic      β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Passband     β”‚ Maximally   β”‚ Equiripple β”‚ Monotonic  β”‚ Equiripple    β”‚
β”‚              β”‚ flat        β”‚            β”‚            β”‚               β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Stopband     β”‚ Monotonic   β”‚ Monotonic  β”‚ Equiripple β”‚ Equiripple    β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Transition   β”‚ Widest      β”‚ Medium     β”‚ Medium     β”‚ Sharpest      β”‚
β”‚ band         β”‚             β”‚            β”‚            β”‚               β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Order for    β”‚ Highest     β”‚ Medium     β”‚ Medium     β”‚ Lowest        β”‚
β”‚ given specs  β”‚             β”‚            β”‚            β”‚               β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Group delay  β”‚ Most        β”‚ Less       β”‚ Better thanβ”‚ Most          β”‚
β”‚              β”‚ uniform     β”‚ uniform    β”‚ Type I     β”‚ non-uniform   β”‚
β”œβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”Όβ”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€
β”‚ Use case     β”‚ General,    β”‚ Sharp      β”‚ Flat       β”‚ Minimum order,β”‚
β”‚              β”‚ smooth      β”‚ cutoff,    β”‚ passband,  β”‚ tight specs   β”‚
β”‚              β”‚ response    β”‚ ok ripple  β”‚ ok stopbandβ”‚               β”‚
β””β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”΄β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”΄β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”΄β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”΄β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”€β”˜

12.2 μ’…ν•© μ‹œκ°μ  비ꡐ

def comprehensive_comparison():
    """Complete comparison of all four IIR filter types."""
    fs = 16000
    fp = 2000
    fstop = 3000
    Rp = 1.0
    As = 50

    filters = {}

    # Butterworth
    N, Wn = signal.buttord(fp, fstop, Rp, As, fs=fs)
    sos = signal.butter(N, Wn, fs=fs, output='sos')
    filters['Butterworth'] = (sos, N)

    # Chebyshev I
    N, Wn = signal.cheb1ord(fp, fstop, Rp, As, fs=fs)
    sos = signal.cheby1(N, Rp, Wn, fs=fs, output='sos')
    filters['Chebyshev I'] = (sos, N)

    # Chebyshev II
    N, Wn = signal.cheb2ord(fp, fstop, Rp, As, fs=fs)
    sos = signal.cheby2(N, As, Wn, fs=fs, output='sos')
    filters['Chebyshev II'] = (sos, N)

    # Elliptic
    N, Wn = signal.ellipord(fp, fstop, Rp, As, fs=fs)
    sos = signal.ellip(N, Rp, As, Wn, fs=fs, output='sos')
    filters['Elliptic'] = (sos, N)

    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    colors = {'Butterworth': 'blue', 'Chebyshev I': 'red',
              'Chebyshev II': 'green', 'Elliptic': 'purple'}

    # Magnitude response overlay
    ax = axes[0, 0]
    for name, (sos, N) in filters.items():
        w, H = signal.sosfreqz(sos, worN=4096, fs=fs)
        ax.plot(w, 20 * np.log10(np.abs(H) + 1e-15),
                color=colors[name], linewidth=1.5, label=f'{name} (N={N})')
    ax.axhline(-Rp, color='gray', linestyle='--', alpha=0.5)
    ax.axhline(-As, color='gray', linestyle='--', alpha=0.5)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Magnitude (dB)')
    ax.set_title('Magnitude Response')
    ax.set_ylim(-70, 5)
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

    # Passband detail
    ax = axes[0, 1]
    for name, (sos, N) in filters.items():
        w, H = signal.sosfreqz(sos, worN=4096, fs=fs)
        ax.plot(w, 20 * np.log10(np.abs(H) + 1e-15),
                color=colors[name], linewidth=1.5, label=f'{name} (N={N})')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Magnitude (dB)')
    ax.set_title('Passband Detail')
    ax.set_xlim(0, fp * 1.2)
    ax.set_ylim(-2, 0.5)
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

    # Group delay
    ax = axes[1, 0]
    for name, (sos, N) in filters.items():
        b, a = signal.sos2tf(sos)
        w_gd, gd = signal.group_delay((b, a), w=4096, fs=fs)
        ax.plot(w_gd, gd, color=colors[name], linewidth=1.5, label=f'{name}')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Group Delay (samples)')
    ax.set_title('Group Delay')
    ax.set_xlim(0, fs / 2)
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

    # Step response
    ax = axes[1, 1]
    step = np.ones(200)
    for name, (sos, N) in filters.items():
        y = signal.sosfilt(sos, step)
        ax.plot(y, color=colors[name], linewidth=1.5, label=f'{name}')
    ax.set_xlabel('Sample')
    ax.set_ylabel('Amplitude')
    ax.set_title('Step Response')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

    plt.suptitle('Comprehensive IIR Filter Comparison', fontsize=14)
    plt.tight_layout()
    plt.savefig('comprehensive_comparison.png', dpi=150)
    plt.close()

comprehensive_comparison()

13. Python κ΅¬ν˜„

13.1 μ™„μ „ν•œ IIR 섀계 μ›Œν¬ν”Œλ‘œμš°

def complete_iir_workflow(fs, filter_type, band_type, freqs, Rp, As):
    """
    Complete IIR filter design and analysis workflow.

    Parameters:
        fs: sampling frequency (Hz)
        filter_type: 'butter', 'cheby1', 'cheby2', 'ellip'
        band_type: 'lowpass', 'highpass', 'bandpass', 'bandstop'
        freqs: (f_pass, f_stop) or ((fp1, fp2), (fs1, fs2)) for bandpass/stop
        Rp: passband ripple (dB)
        As: stopband attenuation (dB)

    Returns:
        sos: second-order sections
        info: design information
    """
    f_pass, f_stop = freqs

    # Order determination
    ord_funcs = {
        'butter': signal.buttord,
        'cheby1': signal.cheb1ord,
        'cheby2': signal.cheb2ord,
        'ellip': signal.ellipord,
    }
    N, Wn = ord_funcs[filter_type](f_pass, f_stop, Rp, As, fs=fs)

    # Design
    design_funcs = {
        'butter': lambda: signal.butter(N, Wn, btype=band_type, fs=fs, output='sos'),
        'cheby1': lambda: signal.cheby1(N, Rp, Wn, btype=band_type, fs=fs, output='sos'),
        'cheby2': lambda: signal.cheby2(N, As, Wn, btype=band_type, fs=fs, output='sos'),
        'ellip': lambda: signal.ellip(N, Rp, As, Wn, btype=band_type, fs=fs, output='sos'),
    }
    sos = design_funcs[filter_type]()

    info = {
        'filter_type': filter_type,
        'band_type': band_type,
        'order': N,
        'Wn': Wn,
        'Rp': Rp,
        'As': As,
    }

    return sos, info

# Example: Bandpass filter
fs = 44100
sos_bp, info_bp = complete_iir_workflow(
    fs=fs,
    filter_type='ellip',
    band_type='bandpass',
    freqs=([800, 3000], [500, 3500]),
    Rp=0.5,
    As=50
)

print(f"Bandpass filter: {info_bp}")

# Apply to signal
t = np.arange(0, 0.5, 1/fs)
# Speech-like signal: mix of tones + noise
x = (0.5 * np.sin(2*np.pi*200*t) +    # Low freq
     np.sin(2*np.pi*1000*t) +           # Mid freq (passband)
     0.8 * np.sin(2*np.pi*2000*t) +     # Mid freq (passband)
     0.3 * np.sin(2*np.pi*5000*t) +     # High freq
     0.2 * np.random.randn(len(t)))      # Noise

y = signal.sosfilt(sos_bp, x)

fig, axes = plt.subplots(3, 1, figsize=(14, 10))

# Time domain
axes[0].plot(t[:2000] * 1000, x[:2000], 'b-', alpha=0.7, label='Input')
axes[0].plot(t[:2000] * 1000, y[:2000], 'r-', linewidth=1.5, label='Filtered')
axes[0].set_xlabel('Time (ms)')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Time Domain: Bandpass Filtering')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Spectrum
N_fft = len(x)
freqs_fft = np.fft.rfftfreq(N_fft, 1/fs)
X = np.abs(np.fft.rfft(x)) / N_fft
Y = np.abs(np.fft.rfft(y)) / N_fft

axes[1].plot(freqs_fft, 20*np.log10(X + 1e-15), 'b-', alpha=0.7, label='Input')
axes[1].plot(freqs_fft, 20*np.log10(Y + 1e-15), 'r-', linewidth=1.5, label='Filtered')
axes[1].axvspan(800, 3000, alpha=0.1, color='green', label='Passband')
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_ylabel('Magnitude (dB)')
axes[1].set_title('Frequency Domain')
axes[1].set_xlim(0, 8000)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# Filter frequency response
w, H = signal.sosfreqz(sos_bp, worN=4096, fs=fs)
axes[2].plot(w, 20*np.log10(np.abs(H) + 1e-15), 'b-', linewidth=1.5)
axes[2].axvspan(800, 3000, alpha=0.1, color='green', label='Passband')
axes[2].set_xlabel('Frequency (Hz)')
axes[2].set_ylabel('Magnitude (dB)')
axes[2].set_title(f'Filter Response (Elliptic N={info_bp["order"]})')
axes[2].set_xlim(0, 8000)
axes[2].set_ylim(-70, 5)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('iir_bandpass_demo.png', dpi=150)
plt.close()

13.2 SOSλ₯Ό μ‚¬μš©ν•œ μ‹€μ‹œκ°„ IIR 필터링

class RealtimeIIRFilter:
    """Real-time IIR filter using second-order sections."""

    def __init__(self, sos):
        """
        Initialize with SOS coefficients.

        Parameters:
            sos: second-order sections array (L x 6)
        """
        self.sos = np.array(sos, dtype=np.float64)
        self.n_sections = self.sos.shape[0]
        # State variables: 2 delay elements per section
        self.state = np.zeros((self.n_sections, 2))

    def process_sample(self, x):
        """Process a single input sample."""
        for i in range(self.n_sections):
            b0, b1, b2, a0, a1, a2 = self.sos[i]

            # Direct Form II Transposed
            y = b0 * x + self.state[i, 0]
            self.state[i, 0] = b1 * x - a1 * y + self.state[i, 1]
            self.state[i, 1] = b2 * x - a2 * y

            x = y  # Output of this section is input to next

        return y

    def process_block(self, x_block):
        """Process a block of samples."""
        y = np.zeros_like(x_block)
        for n in range(len(x_block)):
            y[n] = self.process_sample(x_block[n])
        return y

    def reset(self):
        """Reset filter state."""
        self.state = np.zeros((self.n_sections, 2))

# Demonstration
sos = signal.butter(6, 1000, fs=8000, output='sos')
filt = RealtimeIIRFilter(sos)

# Process in blocks (simulating real-time)
np.random.seed(42)
x = np.random.randn(1000)
block_size = 64

y_realtime = np.zeros_like(x)
for i in range(0, len(x), block_size):
    block = x[i:i + block_size]
    y_realtime[i:i + len(block)] = filt.process_block(block)

# Compare with scipy batch processing
y_batch = signal.sosfilt(sos, x)

error = np.max(np.abs(y_realtime - y_batch))
print(f"Maximum error between real-time and batch: {error:.2e}")

14. μ—°μŠ΅ 문제

μ—°μŠ΅ 1: λ²„ν„°μ›ŒμŠ€ ν•„ν„° 섀계

λ‹€μŒ 규격의 디지털 λ²„ν„°μ›ŒμŠ€ μ €μ—­ 톡과 ν•„ν„°λ₯Ό μ„€κ³„ν•˜μ„Έμš”: - μƒ˜ν”Œλ§ 주파수: $f_s = 10000$ Hz - ν†΅κ³ΌλŒ€μ—­ 엣지: $f_p = 1500$ Hz, λ¦¬ν”Œ $\leq 0.5$ dB - μ €μ§€λŒ€μ—­ 엣지: $f_s = 2000$ Hz, 감쇠 $\geq 40$ dB

(a) λ²„ν„°μ›ŒμŠ€ 차수 곡식을 μ‚¬μš©ν•˜μ—¬ ν•„μš”ν•œ ν•„ν„° 차수λ₯Ό κ³„μ‚°ν•˜μ„Έμš”.

(b) 디지털 주파수λ₯Ό μ•„λ‚ λ‘œκ·Έ 주파수둜 사전 μ›Œν•‘ν•˜μ„Έμš”.

(c) μ•„λ‚ λ‘œκ·Έ λ²„ν„°μ›ŒμŠ€ μ›ν˜•μ„ μ„€κ³„ν•˜κ³  쌍일차 λ³€ν™˜μ„ μ μš©ν•˜μ„Έμš”.

(d) 섀계가 κ·œκ²©μ„ μΆ©μ‘±ν•˜λŠ”μ§€ κ²€μ¦ν•˜μ„Έμš”. 크기 응닡, μœ„μƒ 응닡, κ΅° 지연(group delay), κ·Ή-영점 λ‹€μ΄μ–΄κ·Έλž¨μ„ κ·Έλ¦¬μ„Έμš”.

(e) μˆ˜λ™ 섀계λ₯Ό scipy.signal.butter와 λΉ„κ΅ν•˜μ„Έμš”. κ²°κ³Όκ°€ λ™μΌν•œκ°€μš”?

μ—°μŠ΅ 2: 체비쇼프 ν•„ν„° 비ꡐ

규격: $f_s = 8000$ Hz, $f_p = 1000$ Hz, $f_\text{stop} = 1200$ Hz, $R_p = 1$ dB, $A_s = 50$ dB에 λŒ€ν•΄:

(a) 체비쇼프 Type Iκ³Ό Type II ν•„ν„°λ₯Ό λͺ¨λ‘ μ„€κ³„ν•˜μ„Έμš”.

(b) 단일 κ·Έλž˜ν”„μ—μ„œ 주파수 응닡을 λΉ„κ΅ν•˜μ„Έμš”.

(c) ν†΅κ³ΌλŒ€μ—­μ—μ„œ μ–΄λŠ 것이 κ΅° 지연 평탄도가 더 μ’‹μ€κ°€μš”? μ •λŸ‰μ μœΌλ‘œ λ³΄μ—¬μ£Όμ„Έμš”.

(d) 두 ν•„ν„°λ₯Ό μ²˜ν”„(chirp) μ‹ ν˜Έ(100 Hzμ—μ„œ 4000 Hz μŠ€μœ•)에 μ μš©ν•˜κ³  좜λ ₯을 λΉ„κ΅ν•˜μ„Έμš”.

μ—°μŠ΅ 3: μ˜€λ””μ˜€μš© νƒ€μ›ν˜• ν•„ν„°

μ˜€λ””μ˜€ μ•ˆν‹°μ—μΌλ¦¬μ–΄μ‹±(anti-aliasing)을 μœ„ν•œ νƒ€μ›ν˜• μ €μ—­ 톡과 ν•„ν„°λ₯Ό μ„€κ³„ν•˜μ„Έμš”: - μž…λ ₯ μƒ˜ν”Œλ§ 속도: 96 kHz (48 kHz둜 λ‹€μš΄μƒ˜ν”Œλ§ μ˜ˆμ •) - ν†΅κ³ΌλŒ€μ—­: 0~20 kHz, λ¦¬ν”Œ $\leq 0.01$ dB - μ €μ§€λŒ€μ—­: 24 kHz~48 kHz, 감쇠 $\geq 96$ dB (16λΉ„νŠΈ 정밀도)

(a) μ΅œμ†Œ ν•„ν„° 차수λ₯Ό κ²°μ •ν•˜μ„Έμš”.

(b) SOS ν˜•μ‹μœΌλ‘œ κ΅¬ν˜„ν•˜μ„Έμš”.

(c) μ„ ν˜• 및 dB μŠ€μΌ€μΌ λͺ¨λ‘μ—μ„œ 크기 응닡을 κ·Έλ¦¬μ„Έμš”.

(d) κ΅° 지연을 κ³„μ‚°ν•˜κ³  κ·Έλ¦¬μ„Έμš”. μ˜€λ””μ˜€μ—μ„œ κ΅° 지연 변동이 ν—ˆμš© κ°€λŠ₯ν•œκ°€μš”?

μ—°μŠ΅ 4: μˆ˜μž‘μ—… 쌍일차 λ³€ν™˜

1μ°¨ μ•„λ‚ λ‘œκ·Έ μ €μ—­ 톡과 ν•„ν„°κ°€ μ£Όμ–΄μ‘Œμ„ λ•Œ:

$$H_a(s) = \frac{\Omega_c}{s + \Omega_c}$$

(a) 쌍일차 λ³€ν™˜ $s = \frac{2}{T}\frac{z-1}{z+1}$을 μ μš©ν•˜μ—¬ $H(z)$λ₯Ό μœ λ„ν•˜μ„Έμš”.

(b) $\Omega_c = 2\pi \times 1000$ rad/s, $f_s = 8000$ Hz에 λŒ€ν•΄ 디지털 ν•„ν„° κ³„μˆ˜ $b_0, b_1, a_0, a_1$을 κ³„μ‚°ν•˜μ„Έμš”.

(c) $H(z)$의 주파수 응닡과 μ›Œν•‘λœ μ•„λ‚ λ‘œκ·Έ 응닡을 λΉ„κ΅ν•˜μ—¬ κ²€μ¦ν•˜μ„Έμš”.

(d) λ™μΌν•œ 차단 주파수의 2μ°¨ λ²„ν„°μ›ŒμŠ€ 필터에 λŒ€ν•΄ λ°˜λ³΅ν•˜μ„Έμš”.

μ—°μŠ΅ 5: μ•ˆμ •μ„± 쑰사

(a) 16μ°¨ 체비쇼프 Type I ν•„ν„°($R_p = 3$ dB, $\omega_c = 0.1\pi$)λ₯Ό tf와 sos ν˜•μ‹μœΌλ‘œ μ„€κ³„ν•˜μ„Έμš”.

(b) 두 ν‘œν˜„μ— λŒ€ν•΄ 극점 크기λ₯Ό κ³„μ‚°ν•˜κ³  ν‘œμ‹œν•˜μ„Έμš”.

(c) 두 ν‘œν˜„μ„ μ‚¬μš©ν•˜μ—¬ 백색 작음 μ‹ ν˜Έλ₯Ό ν•„ν„°λ§ν•˜μ„Έμš”. tf ν˜•μ‹ 좜λ ₯이 μœ νš¨ν•œκ°€μš”?

(d) 체비쇼프 Type I ν•„ν„°μ—μ„œ tf ν˜•μ‹μ΄ 일반적으둜 μ–΄λŠ μ°¨μˆ˜μ—μ„œ λ¬΄λ„ˆμ§€κΈ° μ‹œμž‘ν•˜λ‚˜μš”? 차수 4, 8, 12, 16, 20, 24둜 μ‹€ν—˜ν•΄λ³΄μ„Έμš”.

μ—°μŠ΅ 6: IIR vs FIR 비ꡐ

규격: $f_s = 16000$ Hz, $f_p = 3000$ Hz, $f_\text{stop} = 3500$ Hz, $A_s = 60$ dB에 λŒ€ν•΄:

(a) IIR(νƒ€μ›ν˜•)κ³Ό FIR(Parks-McClellan) ν•„ν„° λͺ¨λ‘ μ„€κ³„ν•˜μ„Έμš”.

(b) λ‹€μŒμ„ λΉ„κ΅ν•˜μ„Έμš”: ν•„ν„° 차수, μƒ˜ν”Œλ‹Ή 계산 λΉ„μš©, κ΅° 지연, 크기 응닡.

(c) μŒμ„±κ³Ό μœ μ‚¬ν•œ μ‹ ν˜Έ(100, 200, ..., 5000 Hz의 고쑰파 ν•©μ‚°)λ₯Ό ν•„ν„°λ§ν•˜μ„Έμš”. μ‹œκ°„ μ˜μ—­ 좜λ ₯을 μ‹œκ°μ μœΌλ‘œ λΉ„κ΅ν•˜κ³  SNR κ°œμ„ μ„ κ³„μ‚°ν•˜μ„Έμš”.

(d) IIR을 FIR보닀 선택해야 ν•˜λŠ” μ‘μš© μ‹œλ‚˜λ¦¬μ˜€λŠ” 무엇인지, κ·Έ λ°˜λŒ€λŠ” 무엇인지 μ„€λͺ…ν•˜μ„Έμš”?

μ—°μŠ΅ 7: λŒ€μ—­ 저지(λ…ΈμΉ˜) ν•„ν„°

심전도(ECG) μ‹ ν˜Έμ—μ„œ 60 Hz 전원선 κ°„μ„­(powerline interference)을 μ œκ±°ν•˜κΈ° μœ„ν•œ 쒁은 λŒ€μ—­ 저지 ν•„ν„°λ₯Ό μ„€κ³„ν•˜μ„Έμš”: - μƒ˜ν”Œλ§ 주파수: 500 Hz - 제거 λŒ€μƒ: 59-61 Hz - ν†΅κ³ΌλŒ€μ—­: 0-55 Hz 및 65-250 Hz - ν†΅κ³ΌλŒ€μ—­ λ¦¬ν”Œ: $\leq 0.1$ dB - μ €μ§€λŒ€μ—­ 감쇠: $\geq 30$ dB

(a) νƒ€μ›ν˜• μ›ν˜•μ„ μ‚¬μš©ν•˜μ—¬ μ„€κ³„ν•˜μ„Έμš”. μ–΄λ–€ μ°¨μˆ˜κ°€ ν•„μš”ν•œκ°€μš”?

(b) 60 Hz 주변을 ν™•λŒ€ν•œ 뷰둜 크기 응닡을 κ·Έλ¦¬μ„Έμš”.

(c) 60 Hz μ˜€μ—Όμ΄ μžˆλŠ” ν•©μ„± ECG μ‹ ν˜Έλ₯Ό μƒμ„±ν•˜κ³  ν•„ν„°μ˜ 효과λ₯Ό μ‹œμ—°ν•˜μ„Έμš”.

(d) λ²„ν„°μ›ŒμŠ€ λŒ€ νƒ€μ›ν˜• κ΅¬ν˜„μ˜ κ΅° 지연을 κ³„μ‚°ν•˜κ³  λΉ„κ΅ν•˜μ„Έμš”.

μ°Έκ³ λ¬Έν—Œ

  1. Oppenheim, A. V., & Schafer, R. W. (2010). Discrete-Time Signal Processing (3rd ed.). Pearson. Chapter 7.
  2. Proakis, J. G., & Manolakis, D. G. (2007). Digital Signal Processing (4th ed.). Pearson. Chapter 11.
  3. Parks, T. W., & Burrus, C. S. (1987). Digital Filter Design. Wiley.
  4. Antoniou, A. (2006). Digital Signal Processing: Signals, Systems, and Filters. McGraw-Hill.
  5. SciPy Documentation -- Filter design functions: https://docs.scipy.org/doc/scipy/reference/signal.html
  6. Smith, S. W. (1997). The Scientist and Engineer's Guide to Digital Signal Processing. Available free at http://www.dspguide.com/

탐색

Previous FIR ν•„ν„° 섀계(FIR Filter Design)
πŸ“‹ Topic list
Next 닀쀑λ₯  μ‹ ν˜Έ 처리(Multirate Signal Processing)
← → to navigate between lessons