FFT와 윈도우

---

오디오 신호에서 주파수 정보를 보기

오디오 신호의 특성을 시간축에서 보다 주파수 축에서 알아보는 것이 더욱 편리합니다. 즉 주파수별로 신호의 크기와 위상의 특성으로 오디오가 어떤 정보를 가지고 있는지 알아봅니다. 이를 위해서 다른 글에서 약간 설명한 푸리에 변환을 사용합니다.  

또 현실적인 내용으로 푸리에 변환을 적용하도록 도와주는 각종 윈도우에 대해서도 알아 보겠습니다. 이 윈도우는 디지털 필터를 만들때에 중요한 역할을 하기 때문에 특성과 종류를 알아두는 것이 좋습니다. 

---

푸리에 변환과 주파수 해상도

푸리에 변환(Fourier transform)은 시간축의 신호를 주파수 축의 정보로 변환하는 방법으로 다음과 같이 정의됩니다. 이는 $f(t)$가 연속적인 신호 즉 오디오 신호에 적합한 변환 방식입니다. 

$$F(f)=\int_{-\infty  }^{\infty}e^{-j2\pi ft}f(t)dt$$

 

그런데 우리가 다루고자 하는 신호는 디지털 신호이기 때문에 이를 이산 신호를 위한 푸리에 변환 (discreate Fourier trasnform, DFT)으로 수식을 변경해서 사용합니다. 여기서 N은 프리에 변환을 위해 사용되는 오디오 샘플 수입니다. 즉 DFT는 아날로그 푸리에 변환과 다르게 유한한 샘플 신호에 대해서 변환을 수행합니다. 보다 현실적이죠. 그런데 DFT는 연산이 많이 필요합니다.

$$X[k]=\sum_{n=0}^{N-1} x[n] e^{\frac{-j2\pi kn}{N}}$$

여기서 문제가 발생합니다 푸리에 변화는 신호를 무한한 시간 동안을 대상으로 하기 때문에 유한한 시간의 샘플을 사용하면서 샘플의 끝 구간의 문제로 인해 오류 정보를 가지게 됩니다. 어차피 현실적으로 활용하려면 샘플의 수는 제한될 수밖에 없으니 오류를 해결할 방법을 만들면 됩니다.  

이를 해결하기 위해서 사용하는 방법이 윈도우(window) 함수입니다. 아래 글에서 다시 다루겠습니다. 

 

DFT를 빠르게 연산하도록 개선한 변환 방법이 고속 푸리에 변환 (fast Fourier transform, FFT)이라고 합니다. 그런데 빠른 연산을 위해 만든 FFT로 제한 조건이 있는데, FFT를 수행하는 $N=2^{m}$으로 제한됩니다. 여기서 m은 정수입니다.

N은 FFT 크기(FFT size)라고 하며 신호의 주파수를 얼마나 세밀하게 분석할 수 있는지 따지는 주파수 해상도(frequency resolution, FR) $FR=f_s / N$에 영향을 줍니다. 특히 낮은 주파수 부분에까지 좀 더 많은 정보를 얻고자 한다면 주파수 해상도를 높여야 하며(FR이 작은 숫자로 되게 함), 그러면 FFT의 N를 크게 설정해서 주파수 해상도를 조정해야 합니다. 그런데 N이 커지면 연산량이 증가하게 됩니다.  

DFT의 연산 복잡도는 $N^2$이며 FFT는 $NlogN$으로 알려져 있습니다. 

---

윈도우

윈도우(window)는 샘플 된 구간의 맨 앞쪽과 뒤쪽의 신호가 0으로 변하게 합니다. 그리고 나머지 구간을 모두 0으로 간주하도록 유도합니다. 그러면 푸리에 변환의 무한대 시간영역 조건을 샘플 된 신호 구간으로 좁힐 수 있습니다. 이려면 일단 유한구간에 대해서만 신경 쓰면 되기 때문에 푸리에 변환을 유용하게 사용할 수 있습니다. 물론 이 방법에도 다소의 문제가 있는데 이는 아래의 글에서 다시 다루겠습니다.

다음과 같이 DFT과정에 윈도우 함수 $w[n]$을 추가합니다. $w[n]$은 위에서 설명한 유한 구간의 조건을 만들어 줍니다. 

$$X[k]=\sum_{n=0}^{N-1} x[n]w[n] e^{\frac{-j2\pi kn}{N}}$$

시간축에서 보면 아래와 같이 DFT를 위한 구간의 시작과 끝나는 위치에 신호를 0으로 만들어 줍니다.

원도우 함수도 여러 종류가 있는데 오디오에서 많이 사용하는 원도우는 다음과 같습니다. 사용 용 도에 따라서 적합한 것을 사용자가 골라야 합니다. 

Blackman

사이드 로브가 가장 작음(아래 글 참조), 왜곡 측정에 주로 사용

$$w(n)=0.42-\left (0.5cos\left ( \frac{\pi k}{N} \right )  \right ) + \left (0.08cos\left ( \frac{4\pi k}{N} \right )  \right )$$

Hamming

음향 측정

$$w(n)=0.54-\left ( 0.46cos\left (\frac{2\pi k}{N}  \right ) \right )$$

Hanning

왜곡, 잡음 등의 음향기기 측정

$$w(n)=0.5\left (1-cos\left (\frac{2\pi k}{N}  \right )  \right )$$

 

다음은 대표적인 윈도우를 정현파에 적용한 내용입니다. 이 방법은 디지털 필터 중에 FIR 필터를 구할 때 중요한 역할을 합니다.

clear all

function time_vector = Make_Time_Vector(start_t, end_t, vect_size)
  time_vector = zeros(vect_size,1);
  dt = (end_t - start_t)/(vect_size-1);
  for i =1: vect_size
    time_vector(i,1)= start_t + ((i-1)*dt);
  endfor
endfunction

function sin_vector = Make_sine_Vector(start_t,end_t,freq,time_vector)
  N= size(time_vector,1);
  sin_vector = zeros(N,1);
  for i =0:(N-1)
    sin_vector(i+1,1)= sin(2*pi*freq*time_vector(i+1));
  endfor
endfunction

function W = Rectangular_Window(N);
  W = ones(N,1);
endfunction

function W = Triangular_Window(N)
  W = zeros(N,1);
  for i =0:(N-1)
    W(i+1,1)=1-abs((i-((N-1)/2))/(N/2));
  endfor
endfunction 
 
function W = Sine_Window(N)
  W = zeros(N,1);
  for i =0:(N-1)
    W(i+1,1)=sin(pi*i/(N-1));
  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 W = Hamming_Window(N)
  W = zeros(N,1);
  for i =0:(N-1)
    W(i+1,1)=0.54-(0.46*cos(2*pi*i/(N-1)));
  endfor
endfunction

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

function W = Blackman_Harris_Window(N)
  W = zeros(N,1);
  for i =0:(N-1)
    W(i+1,1)=0.35875-(0.48829*cos(2*pi*i/(N-1)))+\
             (0.14128*cos(4*pi*i/(N-1)))-(0.01168*cos(6*pi*i/(N-1)));
  endfor
endfunction

function W = Gaussian_Window(std_dev,N)
  if(std_dev>0.5)
    std_dev =0.5;
  endif
  W = zeros(N,1);
  for i =0:(N-1)
    m = -0.5*((i-((N-1)/2))/(std_dev*(N-1)/2))^2
    W(i+1,1)= exp(m);
  endfor
endfunction

start_t = 0;
end_t = 1;
vect_size =200;
N= 100;
W = Rectangular_Window(N);
#W = Triangular_Window(N);
#W = Sine_Window(N);
#W = Hanning_Window(N,vect_size);
#W = Hamming_Window(N);
#W = Blackman_Window(N);
#W = Blackman_Harris_Window(N);
#W = Gaussian_Window(0.1,N);
time_vector = Make_Time_Vector(start_t, end_t, vect_size);
sin_vector = Make_sine_Vector(0,1,3000,time_vector);
sin_w = sin_vector;
for i =1:vect_size
  if(i<100)
    sin_w(i,1) = sin_vector(i,1)* W(i,1);
  else
    sin_w(i,1)=0;
  endif
endfor

clf 
figure(1)
plot(W);
xlim([-10,N+10]);
ylim([-0.1,1.1]);
grid on;

figure(2);
plot(sin_vector,'g');
hold on;
plot(sin_w,'r');
hold off;
xlim([-10,vect_size+10]);
ylim([-1.1,1.1]);
grid on;

Rectanglure Window의 모양과 정현파에 적용한 결과입니다. 녹색은 정현파이고 적색은 정현파에 윈도우를 적용한 결과입니다. 이를 시간축에 본 내용입니다. 윈도우가 적용된 간의 시작과 끝의 신호는 0이며 나머지 구간도 0이 됩니다. 이를, 통해 유한한 시간축의 신호를 마치 무한한 시간축 신호인 것처럼 DFT과정에서 속여서 연산합니다.

그러먼 오류가 발생하는데 아래쪽을 보면 내용을 알 수 있습니다. 

Triangular Window의 입니다.

 

Sine Window 입니다.

Hanning Window 입니다.

Hamming Window 입니다. Hanning Window와 유사해 보이지만 윈도우의 시작과 끝부분이 쪼금 다릅니다.

Blackman Window 입니다.

Blackman-Harris Window 입니다.

Gaussian Window입니다. 표준편차 값(std_dev)에 따라 모양이 다릅니다. 

---

윈도우의 영향

1000Hz의 신호를 윈도우를 사용하기 않고 푸리에 변환을 하면 다음과 같이 나타납니다. 예상과 다르게 다른 주파수 성분이 같이 보입니다. 이를 로빙 오류(lobbing error)라고 합니다. 좌우에 로빙 오류가 펼쳐 있어 모양을 사이드 로브(side lob)라고 합니다. 사이드 로브에 의해 푸리에 변환의 결과에서 정확한 결과를 얻기 어렵습니다. 특히 여러 신호가 복합적으로 있는 경우는 더욱 어렵습니다. 

 

이를 개선하기 위해 Hamming 윈도우를 적용해 보겠습니다. 로빙 오류가 만이 개선되었습니다. 그렇다고 완전히 오류가 없는 것은 아닙니다. 

이번엔 Hanning 윈도우를 사용해 보겠습니다. 더욱 로빙 오류가 줄었지만 사실 완전시 사라지지는 않았습니다. 화면에만 안 보일 뿐입니다. 

---

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

---

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

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