오디오 필터(18) - 윈도우 방법을 이용한 FIR 필터 구하기

Audio Processing

오디오 필터(18) - 윈도우 방법을 이용한 FIR 필터 구하기

gigasound 2021. 10. 29. 16:49

윈도우 방법을 이용한 FIR 필터 구하기

지난 글에서 이상적인 필터의 모양을 기반으로 필터 $h(n)$을 구하는 방법을 알아봤습니다. 그런데 $h(n)$은 시간이 -t이고 주파수가 $-\omega $인 조건까지 사용해야 하기 때문에 공학적으로 사용할 수 없습니다. 이를 비 인과성(non-causalty)라고 합니다. 비 인과성은 현재, 과거, 미래의 모든 신호가 필요합니다. -t와 $-\omega $는 미래의 정보입니다.

그런데 우리는 미래의 시간을 만들 수 없습니다. 그러니 인과성(causalty)만 고려해야 합니다. 즉 현재와 과거의 신호만으로 공학적으로 의미 있는 신호처리를 해야 합니다.

이 글에서는 지난 글의 비인과성 조건의 필터를 인과성 조건으로 변경해서 저역 통과 필터를 만들어 보겠습니다.

Sinc 함수의 모양

지난 글에서 저역통과 필터는 주파수 영역에서 구형파 형태입니다. 이를 시간영역으로 역 푸리에 변환하면 Sinc()가 됩니다. Sinc()는 아래 그림과 같이 나타나며 $-\infty \leq t \leq \infty$ 에서 신호가 존재합니다. 즉 미래의 시간에서 과거의 시간 영역까지의 구간이 있어야 Sinc()를 표시할 수 있습니다.

그렇다면 Sinc()의 모든 부분을 과거로 옮기고(shift) 일부만 사용하도록 유도하면 되겠네요. 먼저 위의 그림을 2초 정도 이동해 보겠습니다. 그러면 아래 그림과 같이 Sinc()의 모양이 우측으로 이동하면서 모양은 유지합니다.

시간축에서 이동은 푸리에 변환으로 보면 다음과 같습니다.

$$x(t-c) \leftrightarrow X(e^{j\omega}) e^{-jc}$$

즉 h(n)으로 표시되는 신호를 h(n-M)으로 M 만큼 시간축 이동을 하면 $H(e^{j\omega}) e^{-jc}$로 수정해서 연산하면 됩니다. 그러면 저번 글에서 역 푸리에 변환은 다음과 같이 수정이 되며 저역 주파수의 h(n)도 수정됩니다. 여기서 $H_d(n)$ 은 $H_d(n-M)$로 적어야 정상이지만, 신호처리를 위한 함수 표현으로 $H_d(n)$을 사용합니다. 그러면 n=0으로 시작하고 M 만큼 시간 이동된 신호를 표시한 것으로 간주됩니다.

$$H_d(n)=\frac{1}{2\pi}\int_{\pi}^{-\pi}H_d(e^{j\omega})e^{j\omega(n-M)}d\omega$$

위의 시간 이동 방법을 보면 Sinc()의 제일 높은 신호를 기준으로 좌우가 대칭입니다. 그러므로 이 좌우 대칭(symmetric) 특징을 그대로 살려 주는 게 좋겠습니다. 사실 symmetric을 유지해야 필터의 위상 특성이 선형(linear phase)이 되도록 유도할 수 있습니다.

만약 우리가 N개의 샘플로 FIR 필터를 구현한다면 $M=\frac{N-1}{2}$을 만족하도록 유도합니다. 그러면 구간 [0, N-1]의 중심이 M이 되고, M 만큼 이동한 Sinc는 N구간에서 M을 기준으로 좌우 대칭을 유지하게 됩니다. 여기서 중요한 것은 N은 항상 정수이며 대칭성을 확보하기 위해서는 반드시 홀수여야만 합니다.

FIR 저역 통과 필터

그러면 다음과 같이 시간 이동된 저역 통과 필터를 구할 수 있습니다

$$h_{d,LPF} (n)=\frac{sin(\omega_c (n-M))}{\pi (n-M)} $$

이때 n=M이면 로비탈 정리로 값을 구해야 하며 $h(M)= \omega_c/(\pi) = 2f_c/f_s$로 구해집니다.

이제 남은 것은 그래도 여전히 $-\infty \leq t \leq \infty$인 조건의 Sinc()를 유한 조건의 영역으로 제한하는 것입니다. 이때 유용하게 사용되는 것이 윈도우 함수입니다.

FIR 필터에 윈도우함수 적용

FFT를 다루면서 윈도우 함수(Window function)를 다룬 적이 있습니다. 푸리에 변환이 무한한 시간의 신호를 주파수 축으로 변환하기 때문에 현실적이지 않기 때문에 유한한 구역만 적당히 잘라내서 푸리에 변환해서 사용할 수 있도록 윈도우를 사용 한다고 설명했습니다. 그런데 윈도우 함수는 윈도우를 형성하는 구간의 양 끝의 신호를 0으로 하고 나머지 구간 또한 모두 0을 처리해서 마치 윈도우 구간만 신호가 있는 것처럼 속이는 방법입니다. 이 방법을 Sinc()에도 그대로 사용하겠습니다.

구해진 $h_{d,LPF} (n)$에 원도우 함수 $W(n)$ 곱해서 무한 조건의 n을 유한 조건 $0 \leq n \leq W(N-1)$로 변경합니다.

$$h_{LPF} (n)=h_{d,LPF} (n)W(n)$$

즉 윈도우 구간 밖의 n에 대해서는 신호가 모두 0으로 만들어 미래 신호와 너무 많은 과거 신호 조건을 제한하여 사용할 수 있도록 합니다.

M은 필터의 기울기(slope)와 관련됩니다. 윈도우의 종류는 FIR 필터의 저지 대역이 형성되는 신호의 크기와 관련됩니다.

FIR 저역 통과 필터 구하기 예제

$f_s =44100Hz$, $f_c = 1000Hz$의 LPF, N=5를 구현해 보겠습니다. 그러면 $M=(5-1)/2 = 2$가 되며, $h_d(M)=H(2)=2*f_c/f_s=0.04535$입니다. 벌써 필터 계수 하나를 구했습니다.

필터의 시간축 특성이 M번째를 기준으로 좌우 대칭이기 때문에 h_d(0)=h_d(4), h_d(1)=h_d(3)입니다. 그러니 N=5인 조건에서는 h_(0), h_(1)만 더 구하면 됩니다. 그리고 여기에 해밍 윈도우를 적용하면 다음과 같이 필터 계수들을 구할 수 있습니다.

전달 함수는 다음과 같이 구해집니다.

$$H(z)=h(0)+h(1)z^{-1}+h(2)z^{-2}+h(3)z^{-3}+h(4)z^{-4}$$

이때 z와 f의 관계를 이용해서 필터의 진폭 응답을 구할 수 있습니다.

$$H(e^{j\omega})=H(z)\mid_{z=e^{j\omega}} =\sum_{n=0}^{N-1}h(n)e^{-j\omega}$$

$$\left | H(e^{j\omega}) \right |=h(M)+\sum_{n=0}^{M-1}h(n)cos(\omega (n-M))$$

아래의 GNU Octave 스크립은 윈도우 방식의 FIR 저역 통과 필터를 구하는 내용입니다.

clear all;
pkg load signal;

N=51;
fs=44100;
fc=1000;
vect_size = 1000;

function freq= Make_Lin_Freq(f_min, f_max, vect_size)
  freq = zeros(vect_size,i);
  df = (f_max-f_min)/(vect_size-1);
  tmp = f_min;
  freq(1,1)=f_min;
  for k=2:vect_size
    freq(k,1)=freq(k-1,1)+df; 
  endfor
endfunction

function W= Hanning_window(N)
  W=zeros(N,1);
  for i =0:N-1
    W((i+1),1)=0.5*(1-cos(2*pi*i/(N-1)));
  endfor
endfunction

function h = LPF(fs,fc,N)
  h = zeros(N,1);
  wc = 2*pi*fc/fs;
  M= (N-1)/2;
  W = Hanning_window(N);
  for n = 0: M
    if (n!=M)
      h(n+1,1) = sin(wc*(n-M))/pi/(n-M) *  W(n+1,1);
      nn = (2*M)-n;
      h(nn+1,1) = h(n+1,1);
    else
      h(n+1,1) = wc/pi * W(n+1,1);
    endif
  endfor
endfunction

function mag=Calc_Mag(fs,h,N,freq)
  M = (N-1)/2;
  mag = zeros(size(h,1),1);
  c = zeros(M,1);
  for n =0: (size(freq,1)-1)
    w = 2*pi*freq(n+1,1)/fs;
    tmp =0;
    for k=0:M-1
      tmp += h(k+1)*cos(w*(k-M));
    endfor
    mag(n+1,1) =20*log10(abs( h(M+1)+(2*tmp))); 
  endfor
endfunction

if(mod(N,2)==0)
  N =N+1;
endif

h = LPF(fs,fc,N);
freq= Make_Lin_Freq(20, 20000, vect_size);
mag=Calc_Mag(fs,h,N,freq);
hold on;
semilogx(freq,mag);
ylim([20, 20000]);
ylim([-80, 10]);
grid on;

아래의 그림은 $f_c=1000Hz$에서 N에 따라서 생성된 필터의 모양입니다. N이 작으면 원하는 필터 모양을 갖추지 못할 뿐만 아니라 $f_c$를 구현하지도 못합니다. 더군다나 $f_c$가 낮으면 큰 N을 사용해야 합니다. N=51이면 IIR 필터에 비해 매우 많은 연산을 사용해야 합니다.

Hanning 윈도우를 사용해서 최소 저지 대역은 -44dB 근처에서 발견됩니다.

FIR 고역 통과 필터

고역통과 필터는 저역통과 필터에 -1을 곱해서 구하면 됩니다. 이는 IDF 과정의 적분에서 적분 구간이 저역통과 필터와 고역 통과 필터가 다르기 때문에 그렇습니다.

나머지 필터 계수를 구하는 방법은 저역통과 필터와 같습니다. 다만 $h(M)= (\pi-\omega_c)/\pi = (1-2f_c)/f_s$가 됩니다.

$$h_{d,HPF} (n)=\frac{-sin(\omega_c (n-M))}{\pi (n-M)} $$

FIR 대역 통과 필터

대역 통과 필터도 다음과 같이 구하면 됩니다. 여기서 다만 $h(M)= (\pi-(\omega_{higher}-\omega_{lower}))/\pi$가 됩니다.

$$h_{d,BPF} (n)=\frac{sin(\omega_{higher} (n-M))}{\pi (n-M)}-\frac{sin(\omega_{lower} (n-M))}{\pi (n-M)}$$

광고좀 꾹 눌러주시면 고맙겠습니다.

기타 문의 사항은 gigasound@naver.com에 남겨 주시면 고맙겠습니다.

저작자표시 비영리 변경금지