Frequency Sampling FIR Filter Design

FIR Filter Design에는 여러가지가 있는데, 그 중에서 가장 직관적이고 간단한 방법이 Frequency Sampling 방식이다.

방법은 이렇다.

1. 원하는 형태의 주파수 도메인의 Magnitude Response $*$H(k)$*$를 만든다.

2. 만들어진 주파수 도메인 데이터를 Inverse DFT해서 시간축 필터 $*$h(n)$*$을 생성한다.

말했듯히 매우 간단하고 직관적이다.

간단하게 Low-Pass Filter를 하나 만들어보자.

Filter의 조건은 다음과 같다.

- Filter의 Order는 $*$N$*$는 30으로 하자. 그러면 Filter의 Tap수 $*$L$*$은 30 + 1 = 31이다.

- Cutoff 주파수 $*${ \omega  }_{ c }$*$는 0.25 라고 하자. 단위는 Normalized Frequency [frequency/sampling frequency]이다.

- Filter Coefficient는 실수이다.

- Stop-Band의 Gain은 0이고, Pass-Band의  Gain은 1으로 하자.

Step 1. 위 조건을 만족시키는 Magnitude Response를 생성한다.

위 조건을 만족시키는 $*$H(k)$*$는 아래 그래프와 같다. 말 그대로 주파수 축에다 이렇게 아래 그래프처럼 그려주면 된다.

우리가 일단 원하는것은 실수값을 가지는 FIR필터이므로 Magnitude Response는 0을 기준으로 대칭이 되어야 한다.

이 때에 주파수축 위치는 어떻게 잡아야 할까?

정답은 IDFT 공식에 있다.

$$h(n)\quad =\quad \cfrac { 1 }{ n } \sum _{ k=0 }^{ N-1 }{ H(k) } { e }^{ jn(\cfrac { 2\pi  }{ N } )k }$$

여기에서 $*$ \frac { 1 }{ N } k $*$가 Normalized Frequency를 나타내는 부분이다. 따라서, $*$ \frac { 1 }{ N } *\quad [0,\quad 1,\quad 2,\quad ...,\quad N-1] $*$ 에 해당하는 주파수에서 값을 가지도록 해야 한다. 

그런데, Sampling된 신호는 1 Normalized Frequency마다 반복 되므로, 주파수축을 $*$ \left\lfloor \cfrac { N }{ 2 } \right\rfloor $*$ 만큼 왼쪽으로 이동시키면 위와 같은 그래프를 얻을 수 있다. ($*$\left\lfloor \right\rfloor $*$는 내림을 의미한다.)

다음 step의 IDFT또한 주파수축과 시간축 모두 $*$ \left\lfloor \cfrac { N }{ 2 } \right\rfloor $*$만큼 이동된 공식을 사용할 것이다. 다음과 같다.

$$h(n)\quad =\quad \cfrac { 1 }{ n } \sum _{ k=0 }^{ N-1 }{ H(k) } { e }^{ jn(\cfrac { 2\pi  }{ N } )(k-\left\lfloor \cfrac { N }{ 2 }  \right\rfloor ) },\quad \quad \quad n=-\left\lfloor \cfrac { N }{ 2 }  \right\rfloor ,-\left\lfloor \cfrac { N }{ 2 }  \right\rfloor +1,...,-\left\lfloor \cfrac { N }{ 2 }  \right\rfloor +N-2,-\left\lfloor \cfrac { N }{ 2 }  \right\rfloor +N-1$$

Step 2. IDFT를 수행해서 시간축 Filter Coefficient $*$h(n)$*$을 생성한다.

위 IDFT공식에 의해서 생성된 $*$h(n)$*$를 그려보면 아래 그림과 같다.

자 그럼, 이 $*$h(n)$*$을 이용하면 Stop-Band에 있는 신호들이 $*$H(k)$*$에 나타난것 처럼 다 잘려질까? ㅋㅋㅋㅋ

맨 위의 그림처럼 주파수 0을 기준으로 대칭인 rectangular window를 ideal low-pass filter라고 한다. Rectangular window의 Fourier Transform은 Sinc Function ($*$sinc(x)=\frac { sin(\pi x) }{ \pi x } $*$) 인데, 이 Sinc Function은, 모든 x축 위치에서 값을 가지는 끝나지 않는 함수다. 

다시 말하면, ideal low-pass filter를 표현하기 위해서는 무한대의 tap수가 필요하다는 이야기다. 따라서, 고작 31tap으로는 우리가 처음에 조건으로 만들었던 저 필터의 조건을 만족하기는 쉽지 않을것 같다.

그렇다면 위 그림에 나타난 filter의 실제 성능은 어떻게 될까?

matlab의 freqz함수를 이용해서 저 filter의 magnitude/phase response를 구해보면 다음과 같다.

실제로는 stop-band attenuation이 대략 -30dB가 안되는것을 볼 수 있고, pass-band에서 큰 ripple이 있는것 (위 그림에서 대략 0.45 Normalized Frequency위치)또한 볼 수 있다.

Frequency Sampling FIR Filter Design의 단점이 여기서 나타난다. 디자인 단계에서는 결과를 쉽게 예측하기가 힘들다.

특히 filter의 성능에서 중요한 요소인 Stop-Band Attenuation, Ripple, Transition Band등을 마음대로 조절하기 어렵다. 

직관적이고, 구현이 쉬운 만큼, 세밀한 조정을 하기 어려운 것이다.

물론, 이러한 문제점을 해결하기 위해서 Magnitude Response를 생성할 때, 주어진 Tap수 보다 더 큰 Tap수에서 디자인하고 down-sampling을 하거나, Least-Square 방식으로 원하는 Magnitude Response를 가지도록 하거나, Window Function을 이용해서 Filter의 특성을 디자인 하는 방법 등, 매우 많은 방법들이 있다.

Matlab Code

	
% Frequency Sampling - Low-pass Filter Design
% n37jan@gmail.com

function h = FreqSamplingFIRLPF(order, cutoff)



N = order+1;
H = zeros([N,1]);

% Magnitude Response Design
halfN = floor(N/2);
freq = ( (-halfN : (N-1-halfN)))/N;
for idx = 0 : (N-1)
    if( abs(freq(idx+1)) < cutoff)
        H(idx+1) = 1;
    end
end



% Inverse DFT
timeOffset = 0;
if( rem(order, 2) ~= 0)
    timeOffset = 0.5;
end

h = zeros([N,1]);
for idxTime = 0 : (N-1)
    
    t = (idxTime - halfN + timeOffset);
    
    tempSum = 0;
    for idxFreq = 0 : (N-1)
        dw = H(idxFreq+1) * exp(1j*2*pi*freq(idxFreq+1)* t);
        tempSum = tempSum + dw;
    end
    
    h(idxTime+1) = tempSum;
end

% Making sure the outputs are real numbers.
h = real(h);


% Gain Normalize
m = sum(h(:));
h = h/m;