Python으로 배우는 DSP/Part 3. 디지털 필터 설계와 구현

[17] FIR 필터 설계: 윈도우 기법을 이용한 설계와 노이즈 제거 실전

multimedia 2026. 5. 30. 10:00
반응형

들어가며

지난 16편에서는 Z-변환을 통해 시스템의 전달함수 $H(z)$를 구하고, 복소평면 위에서 극점(Pole)과 영점(Zero)의 배치가 필터의 주파수 응답과 안정성에 어떤 영향을 미치는지 시각적으로 확인했습니다. 이론적으로 시스템을 분석하는 렌즈를 장착했다면, 이제는 그 렌즈를 활용해 우리가 원하는 특성의 필터를 직접 깎아내고 만들어볼 차례입니다.

디지털 필터는 크게 두 가지로 나뉩니다. 출력이 오직 현재와 과거의 입력값에만 의존하는 FIR(Finite Impulse Response) 필터와, 과거의 출력값까지 피드백으로 사용하는 IIR(Infinite Impulse Response) 필터입니다. 이번 17편에서는 FIR 필터에 집중합니다.

 FIR 필터는 피드백(Feedback) 구조가 없어 극점이 항상 원점에 위치하므로 구조적으로 100% 안정적이며, 신호의 파형을 왜곡시키지 않는 '선형 위상(Linear Phase)' 특성을 가집니다. 오디오, 영상, 의료 신호처리 등 위상 왜곡이 허용되지 않는 분야에서 FIR 필터가 압도적으로 선호되는 이유가 바로 여기에 있습니다.

이번 글에서는 FIR 필터의 선형 위상 구조를 결정하는 Type I ~ IV 분류 체계를 먼저 이해하고, 윈도우 기법(Window Method)을 이용한 설계 원리를 단계별로 익힙니다. 이어서 Hamming, Hann, Kaiser 등 다양한 윈도우 함수를 적용하여 저역통과(LPF), 고역통과(HPF), 대역통과(BPF) 필터를 직접 설계하고, 실제 잡음이 섞인 신호에 필터를 적용하여 노이즈 제거 성능을 검증하는 실전 실습까지 함께 진행하겠습니다.

 

1. FIR 필터란 무엇인가? 

FIR 필터의 출력 y[n]은 다음과 같이 표현됩니다.

$$y[n] = \sum_{k=0}^{N-1} h[k] \cdot x[n-k]$$

여기서 h[k]는 필터 계수(Filter Coefficients), N은 필터의 탭 수(Tap 수, 또는 차수+1)입니다. 이 식은 16편에서 배운 컨볼루션 연산 그 자체입니다. 즉, FIR 필터는 임펄스 응답 h[k]와 입력 신호 x[n]의 컨볼루션으로 정의됩니다.

FIR 필터의 핵심 특성을 정리하면 다음과 같습니다.

특성 내용
안정성 극점이 항상 원점에만 위치하므로 무조건 안정
위상 특성 선형 위상(Linear Phase) 보장 가능
계수 대칭성 대칭 또는 반대칭 구조로 선형 위상 구현
계산량 IIR 대비 동일 성능에 더 많은 탭 수 필요

 

2. 선형 위상을 위한 FIR 필터의 4가지 유형 (Type I ~ IV)  

FIR 필터가 사랑받는 가장 큰 이유는 신호의 위상을 왜곡하지 않는 선형 위상(Linear Phase)을 설계할 수 있기 때문입니다. 이를 위해서는 필터의 임펄스 응답 $h[n]$이 정중앙을 기준으로 대칭(Symmetric)이거나 반대칭(Anti-symmetric)이어야 합니다.

여기에 필터의 탭(Tap) 수, 즉 길이 $N$이 홀수인지 짝수인지에 따라 FIR 필터는 다음과 같이 4가지 유형으로 나뉩니다. 이 분류는 단순한 학문적 구분이 아니라, 어떤 타입이 어떤 필터(LPF, HPF, BPF 등)를 설계할 수 있는지를 직접 제한하기 때문에, 필터를 설계할 때 이 유형을 제대로 알지 못하면 원하는 주파수 특성을 얻을 수 없으므로 매우 중요합니다.

2.1. 분류기준

타입 대칭성 탭 수 N 차수 M=N-1
Type I 대칭: h[n] = h[N-1-n] 홀수 짝수
Type II 대칭: h[n] = h[N-1-n] 짝수 홀수
Type III 반대칭: h[n] = -h[N-1-n] 홀수 짝수
Type IV 반대칭: h[n] = -h[N-1-n] 짝수 홀수

 

Type I (대칭, 홀수 탭)

가장 범용적인 구조입니다. 대칭 중심이 정수 인덱스에 위치하며, $\omega=0$ 및 $\omega=\pi$에서 응답값에 제약이 없습니다. 저역통과, 고역통과, 대역통과, 대역저지 등 모든 필터 유형에 사용 가능합니다. 윈도우 기법 설계 시 특별한 이유가 없다면 Type I을 기본으로 선택합니다.

$$H(e^{j\omega}) = e^{-j\omega M/2} \cdot A(\omega)$$

여기서 $A(\omega)$는 실수 함수(Zero-phase Response)이고, $e^{-j\omega M/2}$가 선형 위상 성분입니다.

Type II (대칭, 짝수 탭)

대칭 구조이므로 선형 위상은 보장되지만, 대칭 중심이 두 샘플 사이의 중간점에 위치합니다. 이로 인해 $\omega=\pi$(Nyquist 주파수)에서 항상 $H=0$이 됩니다. 따라서 고역통과(HPF)나 대역차단(BSF) 필터에는 사용할 수 없습니다.

⚠️ firwin으로 짝수 탭의 HPF를 설계하면 경고 또는 오류가 발생하는 이유가 바로 Type II의 이 제약 때문입니다.

Type III (반대칭, 홀수 탭)

$h[n] = -h[N-1-n]$ 구조이며, 중심 계수 $h[M/2]$는 반드시 0입니다. 직류 성분($\omega=0$)과 Nyquist 주파수($\omega=\pi$) 모두에서 $H=0$이 됩니다. 저역통과와 고역통과 필터에는 사용할 수 없으며, 주로 힐버트 변환기(Hilbert Transformer)나 미분기(Differentiator) 설계에 활용됩니다. 위상은 선형이지만 90° 위상 이동이 추가됩니다. 

Type IV (반대칭, 짝수 탭)

Type III의 짝수 길이 버전입니다. $\omega=0$에서 항상 $H=0$이므로 저역통과와 대역저지 필터에는 사용할 수 없습니다. $\omega=\pi$에서는 제약이 없어 고역통과 응용이 가능하며, Type III와 마찬가지로 힐버트 변환기와 미분기 설계에 활용됩니다.

4가지 타입 요약 및 적용 가능 필터

타입 LPF HPF BPF BSF 힐버트/미분기 주요 용도
Type I 가능 가능 가능 가능 — 범용 (가장 널리 사용)
Type II 가능 불가능 가능 불가능 — 저역통과, 대역통과
Type III 불가능 불가능 가능 불가능 가능 힐버트 변환기, 미분기
Type IV 불가능 가능 가능 불가능 가능 힐버트 변환기, 미분기

2.2. Python 실습 1: 다양한 FIR 필터 Type 

Python의 scipy.signal 라이브러리는 firwin 함수를 통해 복잡한 수식 없이도 윈도우 기법 기반의 FIR 필터를 설계할 수 있게 해 줍니다. 타입별 임펄스 응답 및 주파수 응답을 아래 코드를 통해 시각화하였습니다. 

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

plt.rcParams['font.family'] = 'Malgun Gothic'
plt.rcParams['axes.unicode_minus'] = False

fs = 8000
cutoff = 1500 / (fs / 2)

# 각 타입의 필터 계수 생성
# Type I: 대칭, 홀수 탭 (N=51, LPF)
h1 = signal.firwin(51, cutoff, window='hamming')

# Type II: 대칭, 짝수 탭 (N=50, LPF — HPF는 불가)
h2 = signal.firwin(50, cutoff, window='hamming')

# Type III: 반대칭, 홀수 탭 — 힐버트 변환기
h3 = signal.firwin(51, [0.05, 0.95], window='hamming', pass_zero=False)
# 반대칭 강제 적용
M = len(h3)
mid = M // 2
h3_anti = h3.copy()
for i in range(mid):
    h3_anti[i] = h3[i]
    h3_anti[M - 1 - i] = -h3[i]
h3_anti[mid] = 0.0

# Type IV: 반대칭, 짝수 탭 — HPF 계열
h4 = signal.firwin(50, [0.05, 0.95], window='hamming', pass_zero=False)
M4 = len(h4)
mid4 = M4 // 2
h4_anti = h4.copy()
for i in range(mid4):
    h4_anti[i] = h4[i]
    h4_anti[M4 - 1 - i] = -h4[i]

filters = [
    ('Type I\n(대칭, 홀수, LPF)', h1, 'steelblue'),
    ('Type II\n(대칭, 짝수, LPF)', h2, 'tomato'),
    ('Type III\n(반대칭, 홀수, BPF계열)', h3_anti, 'seagreen'),
    ('Type IV\n(반대칭, 짝수, BPF계열)', h4_anti, 'darkorchid'),
]

fig, axes = plt.subplots(4, 2, figsize=(13, 14))
fig.suptitle('FIR 필터 Type I ~ IV: 임펄스 응답 및 주파수 응답', fontsize=13, fontweight='bold')

for i, (title, h, color) in enumerate(filters):
    # 임펄스 응답 (계수 구조)
    axes[i, 0].stem(h, linefmt=color, markerfmt=f'o', basefmt='gray')
    axes[i, 0].set_title(f'{title} — 임펄스 응답 h[n]')
    axes[i, 0].set_xlabel('샘플 인덱스')
    axes[i, 0].set_ylabel('계수값')
    axes[i, 0].grid(True, alpha=0.3)

    # 주파수 응답
    w, H = signal.freqz(h, fs=fs, worN=4096)
    axes[i, 1].plot(w, 20 * np.log10(np.abs(H) + 1e-10), color=color)
    axes[i, 1].set_title(f'{title} — 주파수 응답')
    axes[i, 1].set_xlabel('주파수 (Hz)')
    axes[i, 1].set_ylabel('크기 (dB)')
    axes[i, 1].set_xlim(0, fs / 2)
    axes[i, 1].set_ylim(-80, 5)
    axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('fir_type_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

Type I과 Type II의 임펄스 응답이 중심을 기준으로 좌우 대칭인 반면, Type III와 Type IV는 반대칭(중심에 대해 부호가 반전)임을 시각적으로 확인할 수 있습니다. 또한 Type II의 주파수 응답은 $\omega=\pi$ ($4000\text{ Hz}$)에서 반드시 0이 됨을 그래프에서 확인하십시오. 

다음 그림은 위 예제의 타입별 Pole-Zero Plot입니다. FIR 필터는 피드백이 없으므로 극점(Pole)은 항상 원점($z=0$)에만 위치합니다. 차수가 $N-1$이면 원점에 $N-1$중 극점이 놓이며, 이것이 FIR 필터가 무조건 안정한 이유입니다. 필터의 특성은 오직 영점(Zero)의 배치로 결정됩니다.

이 그래프에서 확인할 수 있는 타입별 핵심포인트는 다음과 같습니다.

Type I (좌상) — 극점 $\times 50$이 원점에, 영점은 단위원 내외에 켤레 대칭으로 고르게 분포합니다. $z=\pm1$ 어디에도 구조적 영점이 없으므로 모든 필터 유형에 적용 가능합니다.

Type II (우상) — $z=-1$ ($\omega=\pi$)에 구조적 영점이 반드시 존재합니다(주황색 주석). 계수의 짝수 대칭 구조에서 수학적으로 필연적으로 도출되는 것으로, 이 때문에 HPF·BSF 설계가 불가능합니다.

Type III (좌하) — $z=+1$ ($\omega=0$)에 구조적 영점이 존재합니다(파란색 주석). 반대칭 홀수 구조이므로 $\omega=\pi$에도 영점이 위치하여 $H(0)=H(\pi)=0$이 됩니다.

Type IV (우하) — $z=+1$ ($\omega=0$)에 구조적 영점이 존재합니다. Type III와 달리 짝수 탭이므로 $\omega=\pi$의 제약은 없지만, LPF·BSF에는 사용할 수 없습니다.

실무에서의 선택 원칙은 다음과 같습니다.

  • 일반 필터(LPF, HPF, BPF, BSF) 설계: Type I을 기본으로 사용하고, 탭 수가 홀수인지 항상 확인합니다.
  • 짝수 탭을 써야 할 상황: LPF 또는 BPF라면 Type II도 무방하나, HPF·BSF에는 절대 사용하지 않습니다.
  • 힐버트 변환기, 90° 위상 이동기: Type III 또는 Type IV를 사용합니다.
  • 미분기(Differentiator): Type III(홀수 탭) 또는 Type IV(짝수 탭)를 활용합니다.

 

3. 윈도우 함수를 이용한 필터 

3.1. 이상적인 필터의 한계와 깁스 현상

우리가 꿈꾸는 '이상적인 저역통과 필터(Ideal LPF)'는 특정 차단 주파수($\omega_c$) 이하의 신호는 완벽하게 통과시키고, 그 이상의 주파수는 칼로 자른 듯이 완벽하게 차단하는 벽돌(Brick-wall) 형태의 주파수 응답을 가집니다.

이 이상적인 주파수 응답을 역 이산 시간 푸리에 변환(IDTFT)하여 시간 영역의 임펄스 응답 $h_{ideal}[n]$으로 변환하면, 양쪽으로 끝없이 진동하며 이어지는 $\text{sinc}$ 함수 형태가 도출됩니다.

$$h_{ideal}[n] = \frac{\sin(\omega_c n)}{\pi n}$$

문제는 컴퓨터나 DSP 칩에서 이 필터를 연산하려면 무한한 길이의 메모리와 시간이 필요하다는 점입니다. 따라서 실제 구현을 위해서는 임펄스 응답의 길이를 유한하게(Finite) 잘라내야(Truncation) 합니다.

하지만 무한한 신호를 중간에서 갑자기 싹둑 잘라버리면, 주파수 영역에서는 통과 대역과 차단 대역의 경계 부근에서 파도처럼 출렁이는 깁스 현상(Gibbs Phenomenon)이 발생합니다. 이는 필터의 성능을 심각하게 저하시키는 원인이 됩니다.

3.2. 윈도우 기법(Windowing Method)의 원리

갑작스러운 절단으로 인한 충격을 완화하기 위해 고안된 방법이 바로 윈도우 기법입니다. 잘라낸 임펄스 응답에 양끝이 부드럽게 0으로 수렴하는 '윈도우 함수 $w[n]$'을 곱해주는 방식입니다.

$$h[n] = h_{ideal}[n] \cdot w[n]$$

11편에서 다루었던 다양한 윈도우 함수들이 여기서 맹활약합니다.

  • Rectangular (사각 윈도우): 윈도우를 적용하지 않고 그냥 자른 것과 같습니다. 전이 대역은 좁지만 리플(Ripple)이 매우 큽니다.
  • Hann & Hamming (한, 해밍 윈도우): 코사인 함수 기반으로 양끝을 부드럽게 깎아줍니다. 리플을 효과적으로 억제하여 실무에서 가장 기본적으로 널리 쓰입니다.
  • Kaiser (카이저 윈도우): $\beta$ 파라미터를 조절하여 전이 대역의 폭과 차단 대역의 감쇠량 사이의 트레이드오프를 설계자가 직접 미세 조정할 수 있는 강력한 윈도우입니다.

3.3. 주요 윈도우 함수 비교

대표적인 윈도우 함수들과 그 특성을 비교합니다.

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

plt.rcParams['font.family'] = 'Malgun Gothic'
plt.rcParams['axes.unicode_minus'] = False

N = 51  # 윈도우 길이 (홀수: Type I 보장)

windows = {
    'Rectangular': np.ones(N),
    'Hann':        signal.windows.hann(N),
    'Hamming':     signal.windows.hamming(N),
    'Blackman':    signal.windows.blackman(N),
    'Kaiser (β=8)': signal.windows.kaiser(N, beta=8),
}

fig, axes = plt.subplots(len(windows), 2, figsize=(12, 14))
fig.suptitle('윈도우 함수 비교: 시간 영역 vs 주파수 영역', fontsize=14, fontweight='bold')

for i, (name, w) in enumerate(windows.items()):
    # 시간 영역
    axes[i, 0].plot(w, color='steelblue')
    axes[i, 0].set_title(f'{name} — 시간 영역')
    axes[i, 0].set_xlabel('샘플')
    axes[i, 0].set_ylabel('진폭')
    axes[i, 0].set_ylim(-0.1, 1.1)
    axes[i, 0].grid(True, alpha=0.3)

    # 주파수 영역 (dB)
    W      = np.fft.fft(w, 2048)
    freq   = np.fft.fftfreq(2048)
    mag_db = 20 * np.log10(np.abs(W) / np.abs(W).max() + 1e-10)
    axes[i, 1].plot(freq[:1024], mag_db[:1024], color='tomato')
    axes[i, 1].set_title(f'{name} — 주파수 응답')
    axes[i, 1].set_xlabel('정규화 주파수')
    axes[i, 1].set_ylabel('크기 (dB)')
    axes[i, 1].set_ylim(-120, 5)
    axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('window_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

주요 윈도우 함수의 성능 지표는 아래와 같습니다.

윈도우 저지대역 감쇄 (dB) 전이대역 폭 특징
Rectangular ~21 가장 좁음 리플 심각, 실사용 거의 없음
Hann ~44 중간 범용 목적에 적합
Hamming ~53 중간 오디오 처리에 널리 사용
Blackman ~74 넓음 고감쇄 요구 시 사용
Kaiser 가변 ($\beta$로 제어) 가변 가장 유연, 설계 자유도 최고

💡 Kaiser 윈도우의 $\beta$ 파라미터: $\beta$ 값이 클수록 저지대역 감쇠는 증가하지만 전이대역이 넓어집니다. $\beta=0$이면 Rectangular, $\beta \approx 5.4$이면 Hamming과 유사한 특성을 보입니다.

 

4. FIR 필터 설계 

4.1. 필터 설계 단계

윈도우 기법으로 FIR 필터를 설계하는 절차는 다음과 같습니다.

① 사양(Specification) 정의

  • 필터 유형: LPF / HPF / BPF / BSF
  • 차단 주파수 ($f_c$), 전이대역 폭 ($\Delta f$)
  • 요구 저지대역 감쇠 (dB)

② FIR 타입 및 윈도우 함수 선택

설계하려는 필터 유형에 맞는 FIR 타입을 먼저 결정합니다. 일반 필터라면 Type I(홀수 탭)을 기본으로 선택하고, 요구 감쇠량을 기준으로 적절한 윈도우를 선택합니다.

③ 필터 차수 N 결정

Kaiser 윈도우를 사용할 경우 차수 추정 공식이 있습니다.

$$N \approx \frac{A - 7.95}{2.285 \cdot \Delta\omega}$$

여기서 $A$는 요구 감쇠량(dB), $\Delta\omega = 2\pi\Delta f / f_s$입니다.

④ 이상적 임펄스 응답 계산 후 윈도우 적용

⑤ 주파수 응답 검증 및 반복 수정

4.2. 저역통과 FIR 필터 설계

scipy.signal.firwin을 이용하면 단 몇 줄로 FIR 필터를 설계할 수 있습니다. 기본적으로 firwin은 홀수 탭을 사용할 때 Type I 필터를 생성합니다.

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

plt.rcParams['font.family'] = 'Malgun Gothic'
plt.rcParams['axes.unicode_minus'] = False

fs     = 8000   # 샘플링 주파수 (Hz)
fc     = 1000   # 차단 주파수 (Hz)
N      = 101    # 필터 탭 수 (홀수 → Type I 보장)

cutoff_norm = fc / (fs / 2)

window_list = ['rectangular', 'hann', 'hamming', 'blackman', ('kaiser', 8.0)]
labels      = ['Rectangular', 'Hann', 'Hamming', 'Blackman', 'Kaiser (β=8)']
colors      = ['gray', 'steelblue', 'tomato', 'seagreen', 'darkorchid']

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))
fig.suptitle(f'저역통과 FIR 필터 비교 (fc={fc}Hz, N={N}, Type I)', fontsize=13, fontweight='bold')

for win, label, color in zip(window_list, labels, colors):
    h = signal.firwin(N, cutoff_norm, window=win)
    w, H = signal.freqz(h, fs=fs, worN=8000)
    ax1.plot(w, 20 * np.log10(np.abs(H) + 1e-10), label=label, color=color)
    ax2.plot(w, np.abs(H), label=label, color=color)

for ax in (ax1, ax2):
    ax.axvline(fc, color='black', linestyle='--', linewidth=1.2, label=f'fc={fc}Hz')
    ax.set_xlabel('주파수 (Hz)')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0, fs / 2)

ax1.set_title('주파수 응답 (dB)')
ax1.set_ylabel('크기 (dB)')
ax1.set_ylim(-120, 5)

ax2.set_title('주파수 응답 (선형)')
ax2.set_ylabel('크기')

plt.tight_layout()
plt.savefig('lpf_window_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

출력된 그래프에서 Rectangular 윈도우의 깁스 현상으로 인한 리플이 뚜렷하게 관찰되는 반면, Blackman 및 Kaiser 윈도우에서는 저지대역이 매우 깔끔하게 감쇠됨을 확인할 수 있습니다.

4.3. 고역통과 및 대역통과 필터 설계 

HPF와 BPF 설계 시에도 탭 수(홀수/짝수)에 따른 타입 제약을 반드시 고려해야 합니다. 

fs = 8000

# --- 고역통과 필터 (HPF, Type I) ---
# HPF는 반드시 홀수 탭(Type I)을 사용해야 합니다.
# 짝수 탭(Type II)은 ω=π에서 H=0이므로 HPF 설계가 불가합니다.
h_hpf = signal.firwin(101, 500 / (fs/2), window='hamming', pass_zero=False)

# --- 대역통과 필터 (BPF, Type I) ---
h_bpf = signal.firwin(101, [800 / (fs/2), 2000 / (fs/2)],
                      window='hamming', pass_zero=False)

# --- 대역저지 필터 (BSF, Type I) ---
# BSF도 Type II(짝수 탭)를 사용할 수 없으므로 홀수 탭을 사용합니다.
h_bsf = signal.firwin(101, [800 / (fs/2), 2000 / (fs/2)],
                      window='hamming', pass_zero=True)

# ⚠️ Type II로 HPF 시도 시 경고 발생 예시 (주석 참고)
# h_hpf_wrong = signal.firwin(100, 500/(fs/2), window='hamming', pass_zero=False)
# → ValueError: Type II filter cannot have a highpass or bandstop response.

filters = {
    '고역통과 (HPF, fc=500Hz, Type I)': h_hpf,
    '대역통과 (BPF, 800~2000Hz, Type I)': h_bpf,
    '대역저지 (BSF, 800~2000Hz, Type I)': h_bsf,
}

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
fig.suptitle('FIR 필터 유형별 주파수 응답 (모두 Type I)', fontsize=13, fontweight='bold')

for ax, (title, h) in zip(axes, filters.items()):
    w, H = signal.freqz(h, fs=fs, worN=8000)
    ax.plot(w, 20 * np.log10(np.abs(H) + 1e-10), color='steelblue')
    ax.set_title(title, fontsize=10)
    ax.set_xlabel('주파수 (Hz)')
    ax.set_ylabel('크기 (dB)')
    ax.set_xlim(0, fs / 2)
    ax.set_ylim(-80, 5)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('filter_types.png', dpi=150, bbox_inches='tight')
plt.show()

💡 firwin 내부의 타입 자동 감지: scipy.signal.firwin은 내부적으로 탭 수와 필터 유형을 조합하여 Type II로는 구현 불가한 HPF·BSF가 요청될 경우 명시적인 오류를 발생시킵니다. 이는 잘못된 필터가 조용히 생성되는 상황을 방지하는 안전장치입니다.

 

5. 실전 실습: 노이즈가 섞인 신호 복원하기

이제 실제 상황을 가정한 노이즈 제거 실습을 진행합니다. 순수한 1kHz 정현파에 광대역 백색 잡음과 특정 주파수(3.5kHz)의 고조파 간섭이 섞인 신호를 생성하고, Type I FIR 저역통과 필터로 복원합니다.

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

plt.rcParams['font.family'] = 'Malgun Gothic'
plt.rcParams['axes.unicode_minus'] = False

np.random.seed(42)

# ── 신호 생성 ──────────────────────────────────────────
fs       = 8000
duration = 0.5
t        = np.linspace(0, duration, int(fs * duration), endpoint=False)

f_clean  = 1000   # 원본 신호 주파수 (Hz)
f_noise  = 3500   # 고조파 간섭 주파수 (Hz)
snr_db   = 10     # 목표 SNR (dB)

clean    = np.sin(2 * np.pi * f_clean * t)
harmonic = 0.5 * np.sin(2 * np.pi * f_noise * t)
noise_std = np.sqrt(np.mean(clean**2) / (10 ** (snr_db / 10)))
wgn      = np.random.normal(0, noise_std, len(t))

noisy = clean + harmonic + wgn

# ── FIR 저역통과 필터 설계 (Type I: 홀수 탭) ──────────
N      = 101                  # 홀수 → Type I 보장
cutoff = 1800 / (fs / 2)
h_lpf  = signal.firwin(N, cutoff, window='hamming')

# 필터링 (위상 왜곡 방지: filtfilt 사용)
filtered = signal.filtfilt(h_lpf, 1.0, noisy)

# ── SNR 계산 ───────────────────────────────────────────
def calc_snr(clean_sig, noisy_sig):
    noise = noisy_sig - clean_sig
    return 10 * np.log10(np.mean(clean_sig**2) / (np.mean(noise**2) + 1e-10))

snr_before = calc_snr(clean, noisy)
snr_after  = calc_snr(clean, filtered)
print(f"필터링 전 SNR : {snr_before:.2f} dB")
print(f"필터링 후 SNR : {snr_after:.2f} dB")
print(f"SNR 개선량    : {snr_after - snr_before:.2f} dB")

# ── 시각화 ─────────────────────────────────────────────
fig, axes = plt.subplots(3, 2, figsize=(14, 10))
fig.suptitle('FIR Type I 필터를 이용한 노이즈 제거 실전', fontsize=13, fontweight='bold')

plot_range = slice(0, 400)

signals = [
    ('원본 신호 (clean)', clean,    'steelblue'),
    ('잡음 신호 (noisy)', noisy,    'tomato'),
    ('복원 신호 (filtered)', filtered, 'seagreen'),
]

for i, (label, sig, color) in enumerate(signals):
    axes[i, 0].plot(t[plot_range] * 1000, sig[plot_range], color=color, linewidth=0.8)
    axes[i, 0].set_title(f'{label} — 시간 영역')
    axes[i, 0].set_xlabel('시간 (ms)')
    axes[i, 0].set_ylabel('진폭')
    axes[i, 0].grid(True, alpha=0.3)

    freqs    = np.fft.rfftfreq(len(sig), 1/fs)
    spectrum = 20 * np.log10(np.abs(np.fft.rfft(sig)) / len(sig) + 1e-10)
    axes[i, 1].plot(freqs, spectrum, color=color, linewidth=0.8)
    axes[i, 1].set_title(f'{label} — 주파수 스펙트럼')
    axes[i, 1].set_xlabel('주파수 (Hz)')
    axes[i, 1].set_ylabel('크기 (dB)')
    axes[i, 1].set_xlim(0, fs / 2)
    axes[i, 1].set_ylim(-80, 0)
    axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('noise_removal_result.png', dpi=150, bbox_inches='tight')
plt.show()
필터링 전 SNR : 4.57 dB
필터링 후 SNR : 13.71 dB
SNR 개선량    : 9.14 dB

필터링 후 스펙트럼에서 3.5kHz 고조파 성분이 완전히 제거되고 백색 잡음도 크게 감쇠되어, 필터링 전 4.57 dB에서 필터링 후 13.71 dB로, SNR이 약 9.14 dB 향상되었음을 확인할 수 있습니다. 

 

6. filtfilt vs lfilter: 위상 왜곡 문제

위 코드에서 signal.filtfilt를 사용한 이유가 있습니다.

함수 방식 위상 왜곡 지연 실시간 처리
lfilter 단방향 (인과적) 있음 반 차수 지연 가능
filtfilt 양방향 (비인과적) 없음 0 불과 (오프라인 전용)

FIR 필터는 선형 위상 특성을 설계 단계에서 보장하지만, lfilter를 사용하면 필터 차수의 절반에 해당하는 그룹 지연(Group Delay)이 발생합니다. 오프라인 분석이라면 filtfilt로 이 지연을 완전히 제거할 수 있습니다. 반면 실시간 처리에서는 반드시 lfilter를 사용하고 지연을 별도로 보상해야 합니다.

h = signal.firwin(101, 0.3, window='hamming')
demo = np.sin(2 * np.pi * 0.05 * np.arange(500)) + 0.5 * np.random.randn(500)

out_lfilter  = signal.lfilter(h, 1.0, demo)
out_filtfilt = signal.filtfilt(h, 1.0, demo)

plt.figure(figsize=(10, 4))
plt.plot(demo,        alpha=0.4, label='입력 신호',          color='gray')
plt.plot(out_lfilter,  label='lfilter (그룹 지연 있음)', color='tomato')
plt.plot(out_filtfilt, label='filtfilt (지연 없음)',     color='steelblue')
plt.title('lfilter vs filtfilt: 그룹 지연 비교')
plt.xlabel('샘플')
plt.ylabel('진폭')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('filtfilt_vs_lfilter.png', dpi=150, bbox_inches='tight')
plt.show()

 

7. Kaiser 윈도우를 이용한 정밀 설계

Kaiser 윈도우는 $\beta$ 파라미터 하나로 감쇠량과 전이대역 폭을 동시에 제어할 수 있어 가장 실용적인 윈도우입니다. scipy.signal.kaiserord 함수는 공학적 사양을 입력하면 최적의 $N$과 $\beta$를 자동으로 계산해 줍니다. 

from scipy.signal import kaiserord, firwin

fs      = 8000
ripple  = 60.0   # 요구 저지대역 감쇠 (dB)
width   = 200.0  # 전이대역 폭 (Hz)
cutoff  = 1000.0 # 차단 주파수 (Hz)

N_opt, beta_opt = kaiserord(ripple, width / (fs / 2))

# N을 홀수로 보정 → Type I 보장
if N_opt % 2 == 0:
    N_opt += 1

h_kaiser = firwin(N_opt, cutoff / (fs / 2), window=('kaiser', beta_opt))

print(f"자동 계산된 탭 수 (N) : {N_opt}")
print(f"자동 계산된 β         : {beta_opt:.4f}")

w, H = signal.freqz(h_kaiser, fs=fs, worN=8000)
plt.figure(figsize=(8, 4))
plt.plot(w, 20 * np.log10(np.abs(H) + 1e-10), color='darkorchid')
plt.axvline(cutoff, color='red',  linestyle='--', label=f'fc={cutoff}Hz')
plt.axhline(-ripple, color='gray', linestyle=':', label=f'-{ripple}dB 감쇠선')
plt.title(f'Kaiser 윈도우 LPF, Type I (N={N_opt}, β={beta_opt:.2f})')
plt.xlabel('주파수 (Hz)')
plt.ylabel('크기 (dB)')
plt.xlim(0, fs / 2)
plt.ylim(-90, 5)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('kaiser_optimized.png', dpi=150, bbox_inches='tight')
plt.show()
자동 계산된 탭 수 (N) : 147
자동 계산된 β         : 5.6533

kaiserord는 "저지대역 감쇠 $60\text{ dB}$ 이상, 전이대역 $200\text{ Hz}$ 이내"라는 사양으로부터 $N=147$, $\beta=5.6533$을 자동으로 도출하였습니다. 결과로 얻어진 $N=146$이 짝수였으므로 $1$을 더해 홀수($147$)로 보정하여 Type I 구조를 확실히 보장하였습니다.

그래프에서 확인할 수 있는 주요 특징은 다음과 같습니다.

첫째, $1000\text{ Hz}$ 차단 주파수를 기준으로 전이대역이 약 $200\text{ Hz}$ 폭 안에서 급격하게 감소하고 있습니다. 둘째, 저지대역의 첫 번째 리플 피크가 약 $-63\text{ dB}$로 요구 사양($-60\text{ dB}$)을 충족하며, 이후 리플은 $-75 \sim -90\text{ dB}$ 수준으로 더욱 깊이 감쇠됩니다. 셋째, 저지대역 전체에 걸쳐 $-60\text{ dB}$ 기준선 아래를 유지하고 있어 설계 사양이 완전히 만족되었음을 확인할 수 있습니다.

$147$개의 탭은 결코 적지 않습니다. 이처럼 FIR 필터는 엄격한 전이대역 사양을 만족할수록 탭 수가 급격히 증가하며, 이것이 다음 편에서 다룰 IIR 필터의 존재 이유가 됩니다. FIR 필터는 저지대역 감쇠를 높이려면 전이대역이 넓어지고, 전이대역을 좁히려면 탭 수(계산량)가 늘어납니다. 이 세 가지 요소 — 감쇠량, 전이대역 폭, 탭 수 — 는 항상 상호 절충 관계에 있습니다. Kaiser 윈도우의 kaiserord가 유용한 이유는 이 트레이드오프를 수학적으로 정량화하여, 사용자가 명시한 사양을 가장 효율적으로 만족하는 최소 자원의 필터를 제안해 주기 때문입니다.

 

[Insight] 

우리가 일상에서 가장 쉽게 볼 수 있는 FIR필터는 주가 차트의 5일선, 20일선, 60일선 등입니다. 일반적으로 FIR 필터 설계는 제약 사항을 인지하고 타협점을 찾아가는 과정입니다.

단순히 탭 수(길이)를 짝수로 잘못 설정하는 것만으로도 고역통과 필터(HPF)가 망가져버리는 Type II의 수학적 한계가 존재합니다. 또한, 주파수를 더 날카롭게 자르고 싶거나 차단 대역의 노이즈를 더 완벽하게 누르고 싶어 탭 수를 늘리다 보면, 연산량과 메모리 사용량이 극심해져 실시간 임베디드 기기에서는 치명적인 병목 현상을 일으키기도 합니다.

무조건적으로 안정적이고 위상을 왜곡하지 않는다는 완벽한 장점 이면에는, 이처럼 필터의 유형(Type I ~ IV)을 섬세하게 고려해야 하는 수학적 디테일과 '비싼 연산 비용'이라는 현실적인 제약이 숨어 있습니다. 

또한, 매우 좁은 전이대역이나 급격한 감쇠 사양을 만족하기 위한 수백 개의 탭의 사용은 연산 비용을 폭발적으로 증가시키기 때문에 다음 편에서 다룰 IIR 필터의 존재 이유가 생겨납니다. 

 

다음 글 예고 

오늘 우리는 선형 위상을 보장하는 FIR 필터의 4가지 유형을 이해하고, 윈도우 기법을 활용하여 실전 노이즈 제거까지 수행해 보았습니다.

FIR 필터가 "안정성과 선형 위상"을 보장하는 대신 많은 탭 수를 요구한다면, 연산 자원이 넉넉하지 않은 환경에서, 아주 적은 수의 계수만으로도 예리하게 주파수를 깎아내야 한다면 어떻게 해야 할까요? 다음 시간에는 [18. IIR 필터 설계: Butterworth, Chebyshev 필터 및 안정성 분석] 편을 통해 피드백(Feedback)을 허용하여 극점을 Z-평면 곳곳에 활용하는 IIR 필터의 세계로 들어갑니다. 

 

📌 이전 연재 글 보기

 

 

반응형

'Python으로 배우는 DSP/Part 3. 디지털 필터 설계와 구현'의 다른글

  • 현재글[17] FIR 필터 설계: 윈도우 기법을 이용한 설계와 노이즈 제거 실전

관련글

티스토리툴바