FFT용 윈도우 함수 구현 예제 소스 코드
본 글은 고속 푸리에 변환(Fast Fourier Transform,FFT)에 윈도우 함수(window
function)를 사용하는 이유와 구현된 소스 코드를 싣고 있다.
스펙트럼 누설(Spectral leakage)
신호처리에서 신호를 주파수 도메인에서 분석하기 위해 고속 푸리에 변환을 많이
사용한다. 이 FFT(Fast Fourier Transform) 변환할 때 종종 window를 사용하여
신호를 한번 가공하는 것을 볼 수 있다. 이 window를 왜 사용하는 것인가?
그 이유는 스펙트럼 누설 현상이 발생하기 때문이다. 스펙트럼 누설(Spectral
leakage)이라는 것은 신호를 스펙트럼 분석했을 때, 원래의 신호에는 포함되어 있지
않은 주파수 성분이 관측되는 현상을 말한다
그런데, 재밌는 사실은 스펙트럼 누설 현상은 이 window 때문에 발생한다.
아래 그림과 같이 만약 전체 신호를 FFT를 통해 주파수 분석을 하게 되면 스펙트럼
누설 현상이 나타나지 않는다. 하지만, 전체가 아닌 일정 부분을 샘플링(다르게
말하면, rectangular window을 사용) 하여 FFT를 수행하면 녹색 화살표의 주파수
같은 스펙트럼 누설이 나타난다.
실제 우리가 분석해야 할 신호들은 전체를 한꺼번에 FFT를 취할 수 없기 때문에,
어쩔 수 없이 일부분을 잘라 주파수 분석을 해야 한다. 그럼 이 스펙트럼 누설을
어떻게 해야 하나?
다행이도 스펙트럼 누설을 줄이기 위한 window function들이 존재한다.
window function
다양한 window function들 중 몇 가지를 구현하여 스펙트럼 누설이 사라지는지에
대해 테스트했다.
테스트에서 fft는 kiss fft를 사용하였으며, fft 수행 코드는 아래와 같다.
#include "kiss_fftr.h"
/*
sample buffer length is N
freq_spectrum buffer length is N/2
*/
int do_fft(int N, double *sample, double *freq_spectrum)
{
kiss_fft_cpx *ffr_outdata = NULL;
kiss_fftr_cfg fft_ctx = NULL;
if (!N || !sample ||
!freq_spectrum)
return -1;
fft_ctx = kiss_fftr_alloc(N, 0,
NULL, NULL);
if (fft_ctx == NULL)
return -1;
ffr_outdata =
(kiss_fft_cpx*)malloc(N*sizeof(kiss_fft_cpx));
if (ffr_outdata == NULL)
{
free(ffr_outdata);
return -1;
}
memset(ffr_outdata, 0, N *
sizeof(kiss_fft_cpx));
// fft
kiss_fftr(fft_ctx, sample,
ffr_outdata);
for (int i = 0; i < N / 2;i++)
{
freq_spectrum[i] +=
sqrt(pow(ffr_outdata[i].r, 2.0) + pow(ffr_outdata[i].i, 2.0));
}
safe_free(ffr_outdata);
kiss_fftr_free(fft_ctx);
return 0;
}
테스트에 사용된 신호는 sin파를 사용하였다.
double *get_signal(int size, double sample_rate)
{
double *signal = NULL;
double angle_tab =
((sample_rate/2)*PI) / (sample_rate);
signal =
(double*)malloc(sizeof(double)*size);
if (!signal)
return NULL;
for (int i = 0; i < size; i++)
{
signal[i] = sin(i*angle_tab);
}
return signal;
}
신호에 window적용은 아래 코드와 같이 각 샘플들을 곱해주면 된다.
for (int i = 0; i < signal_size; i++)
{
sample[i] = signal[i] * window[i];
}
아래 결과 그래프들의 좌측은 각 윈도우와 윈도우를 적용한 신호를 보여주고 있고,
우측은 FFT결과를 보여준다. FFT결과의 파란색 선은 해당 윈도우를 적용하지 않은
경우이고, 주황색은 윈도우를 적용한 결과이다.
전반적으로 스펙트럼 누설이 사라져 보이며, 주파수의 amplitude가 윈도우를
사용하지 않은 경우 대비 작이진 것을 볼 수 있다.
1. hamming window
double *alloc_win_hamming(int N)
{
int half, i, idx;
double *w = NULL;
w =
(double*)malloc(N*sizeof(double));
if (!w)
return 0;
memset(w, 0, N*sizeof(double));
if (N % 2 == 0)
// even
{
half = N / 2;
for (i = 0; i<half; i++)
w[i] = 0.54 - 0.46 * cos(2 * PI*(i)
/ (N - 1));
idx = half - 1;
for (i = half; i<N; i++) {
w[i] = w[idx];
idx--;
}
}
else
// odd
{
half = (N + 1) / 2;
for (i = 0; i<half; i++)
w[i] = 0.54 - 0.46 * cos(2 * PI*(i)
/ (N - 1));
idx = half - 2;
for (i = half; i<N; i++) {
w[i] = w[idx];
idx--;
}
}
return w;
}
2. hann window
double *alloc_win_hann(int N)
{
int half = 0;
int i = 0;
int idx = 0;
double *w = NULL;
w =
(double*)malloc(N*sizeof(double));
if (!w)
return 0;
memset(w, 0, N*sizeof(double));
if (N % 2 == 0)
// even
{
half = N / 2;
for (i = 0; i<half; i++)
w[i] = 0.5*(1.0 - cos(2 * PI*(i) /
(N - 1)));
idx = half - 1;
for (i = half; i<N; i++) {
w[i] = w[idx];
idx--;
}
}
else
// odd
{
half = (N + 1) / 2;
for (i = 0; i<half; i++)
w[i] = 0.5*(1.0 - cos(2 * PI*(i) /
(N - 1)));
idx = half - 2;
for (i = half; i<N; i++) {
w[i] = w[idx];
idx--;
}
}
return w;
}
3. blackman window
double *alloc_win_balckman(int N)
{
int i = 0;
double *w = NULL;
double alpha = 0.16;
w =
(double*)malloc(N*sizeof(double));
if (!w)
return 0;
memset(w, 0, N*sizeof(double));
while (i < N)
{
w[i] = (1 - alpha) / 2;
w[i] -= 0.5*cos((2 * PI*i) / (N -
1));
w[i] += (alpha / 2)*cos((4 * PI*i) /
(N - 1));
i++;
}
return w;
}
4. flattop window
double *alloc_win_flattop(int N)
{
int i = 0;
double *w = NULL;
w =
(double*)malloc(N*sizeof(double));
if (!w)
return 0;
memset(w, 0, N*sizeof(double));
while (i < N)
{
double a = 2.0*PI / (double)(N - 1);
w[i] = (1) - (1.93)*cos(a*i) +
(1.29)*cos(2 * a*i) - (0.388)*cos(3 * a*i) + (0.028)*cos(4 * a*i);
i++;
}
return w;
}
5. gaussian window
double *alloc_win_gaussian(int N)
{
int i = 0;
double *w = NULL;
double sigma = 0.5;// 0.1*N;
w =
(double*)malloc(N*sizeof(double));
if (!w)
return 0;
memset(w, 0, N*sizeof(double));
for (i = 0; i<N; i++)
{
double s = (i - (N - 1) / 2) /
(sigma*(N - 1) / 2);
s = s*s;
w[i] = exp(-s / 2);
}
return w;
}
댓글
댓글 쓰기