오디오 필터(19)-주파수 샘플 방법으로 FIR 필터 구하기

---

주파수 샘플링

주파수 샘플링 방법으로 FIR 필터를 구하는 방법은 참으로 기막힌 기법입니다. 시간축에서 계속 변화하는 아날로그 신호를 샘플링으로 디지털 신호로 변화하는 방법이 있었습니다.  이 방법을 주파수 축에서 적용하는 해서 필터를 구하는 방법이 주파수 샘플링 방법입니다.

이  글에서는 주파수 샘플링 방법으로 원하는 특성의 필터를 쉽게 구하는 방법에 대해서 알아보겠습니다. 사실 이 방법이 실제로 필터를 구하기 더 쉽습니다. 

• 먼저 주파수 축에서 구현될 필터 모양을 결정해야 합니다. 수식으로 해도 되고 실제 원하는 필터의 모양에서 필터의 크기에 관한 정보를 수집해야 합니다.
• 그리고 주파수 축에서 각주파수로 [0, 2$\pi$] 사이에 필터의 모양에 따라 N개의 샘플을 등간격으로 채취합니다. 이를 H(k)로 정의해서 k =0...(N-1)의 순서대로 필터의 높이를 기지고 있습니다. 
• 이때 $\pi$를 기준으로 좌우대칭에 되도록 필터의 모양을 복사합니다. 그러니 실제 구현 필터는 [0, $\pi$] 안에 있어야 합니다. 그래야 나이퀴스트 조건을 만족하겠죠.
• 샘플링한 수량에 따라 아래의 수식에 의해 필터 계수를 구하면 끝입니다. 

$$h(n)=\left\{\begin{matrix}
\frac{1}{N}\left [ H(0)+2\sum_{k=1}^{\frac{N-1}{2}}Re(H(k)e^{\frac{j2\pi kn}{N}}) \right ], N=odd\\ 
\frac{1}{N}\left [ H(0)+2\sum_{k=1}^{\frac{N}{2}-1}Re(H(k)e^{\frac{j2\pi kn}{N}}) \right ], N=even
\end{matrix}\right.$$

---

필터 구하기

아래 그림은 주파수 축에서 주파수 샘플링 방법으로 저역 통과 필터 구하기 위해 그린 그림입니다. 저번 글에서 인과성에 대해서 설명한 적 있습니다. 이를 고려하면 필터는 [0, 2$\pi$]에서 형성되어야 하고 $\pi$를 기준으로 저역 부분이 대칭 복사되어  $2\pi$ 부분에 나타납니다. 즉 $-\omega$부분 까지 전부 표시합니다. 

주파수 축에서 N개의 샘플링을 취하면 이때 샘플링 위치의 각주파수는 $\omega = \frac{2\pi}{N}k$ 가 되며, 이를 수집해서 H(k)를 만듭니다. 원하는 필터를 좀 더 정밀하게 반영하고 싶다면 샘플링 숫자 N을 늘려야겠죠.

위의 에서 H(k)는 다음과 같이 구해집니다. 여기서 N=11입니다. 

$$H(k)=\left\{\begin{matrix}
1; k=0,1,2,8,9,10\\ 
0; k=3,4,5,6,7
\end{matrix}\right.$$

H(k)를 구했으니 위의 식을 이용해서 h(n)을 구하면 됩니다. 너무 쉽죠. 다만 $H_d(\omega)$를 어떻게 구할지가 문제인데, 이상적인 필터 형태를 사용해도 되고, 기존의 IIR 필터에서 구현한 필터의 주파수 응답을 사용해도 됩니다.

주파수 샘플링 방법은 저역 통과 필터 등의 필터 특성에 구해를 받지 않고 원하는 형태로 만들 수 있는 장점도 있습니다. 다만 세밀한 필터 변화 특성과 차단 추파수 특성을 반영하려면 많은 주파수 축 샘플링을 해야 합니다. 

---

저역통과 필터 구현

다음의 GNU Octave의 스크립트는 주파수 샘플링 방법으로 저역 통과 필터를 구현하는 내용입니다.

clear all;
pkg load signal;

N=61;
if(mod(N,2)==0)
  N =N+1;
endif
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

## FIR LPF
function [h,w]=LPF(fs,fc,N)
  H = zeros(N,1);
  h = zeros(N,1);
  w = zeros(N,1);
  dw = 2*pi/(N-1);
  w(1,1)=0;
  for k=2:N
    w(k,1)=w(k-1,1)+dw; 
  endfor
  M=(N-1)/2;
  wc = 2*pi*fc/fs;
  for k=0:M
    if(w(k+1,1)<=wc)
      H(k+1,1) = 1;
      kk = (2*M)+k;
      H(kk+1,1) = 1;
    endif
  endfor
  for n =0:N-1
    tmp = 0;
    for k=1:(N-1)/2
      tmp = tmp+(H(k+1,1)*cos(2*pi*k*(M-n)/N));
    endfor
    h(n+1)=(H(0+1,1) + (2*tmp))/N;
  endfor
endfunction

## FIR 주파수 응답 
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

[h,w]=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;

아래 그림은 주파수 샘플 수 N=5와 N=61의 결과입니다. N=61 이상이 되어야 $f_c$에 접근하는 필터 구현이 됩니다.

아래 그림은 N=201일때의 필터 모양으로 필터의 경사부는 매우 급격하지만 통과 대역에 리플이 많이 있습니다. 

---

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

---

위의 내용을 참조용으로만 사용해주세요. 무단 도용이나 무단 복제는 불허합니다.

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