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;
}



댓글

댓글 쓰기

이 블로그의 인기 게시물

간단한 cfar 알고리즘에 대해

이미지
Constant False Alarm Rate (CFAR)은 기본적으로 테스트하고자 하는 위치의 cell과 주변 cell의 관계를 보고 테스트 cell이 Target 인지 아닌지를 구분하는 알고리즘이다. 주변 cell과의 관계를 어떤 방식으로 비교하는가에 따라 여러가지 알고리즘이 존재한다. CA-CFAR, GO-CFAR, SO-CFAR, OS-CFAR, VI-CFAR, OSVI-CFAR등이 있다. 이 많은 알고리즘 중 간단하고 기본적인 CA-CFAR과 OS-CFAR에 대해 알아보자. CA-CFAR (Cell Average CFAR)   CA-CFAR의 원리는 위 그림과 같다. 좌우의 Reference Cell의 평균에 Scale Factor T를 곱하고, ADT (Average Decision Threshold)를 합하여 CUT가 타겟인지 아닌지를 구분한다. ADT는 신호의 white noise를 고려하기 위해 더해준다. Scale Factor T와 ADT는 아래와 같은 공식으로 사용할 수 있다. 사용하려는 시스템에 따라 공식은 달라질 수도 있다.   Pfa는 False Alarm 발생 확률을 나타내고, M과N는각 좌우의 Reference Cell의 개수를 나타낸다. OS-CFAR (Order Static CFAR)   OS-CFAR의 원리는 위 그림과 같다. 좌우의 reference cell중 K번째 큰 값의 cell을 Reference Cell의 대표 값으로 정해 CUT가 타겟인지 아닌지를 구분한다. 보통 전제 reference cell의 크기 순에서 (3*N)/4 번째 cell이 성능이 좋다고 알려져 있다. CFAR 사용 예 아래 사진은 위 두 cfar을 사용한 예다. 파란색은 신호이며, 붉은 색은 char을 사용해 만든 threshold이다. 이 threshold보다 큰 값의 신호가 있는 위치에 실제 Target이 있다고 볼 수 있
자세한 내용 보기

아두이노(arduino) 심박센서 (heart rate sensor) 심박수 측정 example code

이미지
하드 정리 중 예전에 만든 아두이노에서 심박센서를 사용하여 심박수 측정하는 소스 코드가 있어, 이를 블로그 글로 남긴다. 인터넷 쇼핑몰에서 쉽게 구할 수 있는 저렴함 센서가 사용 되었으며, 연결은 아래 그림과 같이 연결 되었다.      심박센서의 하트 부위에 손가락을 올려 놓고, 아날로그 A0값을 읽으면 심박의 움직임에 따라 A0값이 변하는 동작이 간단한 센서이다. 아래 그림의 왼쪽은 일정 시간 간격 마다 A0값을 읽어 시리얼 플로터로 본 심박에 따른 데이터이고, 오른쪽은 평균 필터를 사용해 노이즈를 제거한 후 시리얼 플로터로 본 데이터이다.   심박수는 수집된 데이터에서 threshold를 넘는 peak들의 시간 간격을 측정하여 계산한다.  peak위치를 찾는 방법은 여러 가지 방법이 있으나, 소스에 구현된 방법은 신호의 기울기가 +에서 0또는 -로 바뀌는 위치를 peak으로 선정하도록 하였다. 심장이 t시간동안 N회 뛰는 경우 심박수는 heart rate bpm = (N/t)*60 을 표현 할 수 있다.   threshold는 데이터를 보고 임의로 정했으며, 측정 데이터에 따라 자동으로 정하고 싶을 때는 cfar 같은 방법을 사용하면 된다.  소스는 아래와 같다. #define heartbit_adc_threshold      511 #define meanfilter_size             5 #define heartbit_tick_buffer_size   10 #define measure_buffer_size       3 float meanfilter[meanfilter_size] = {0,};
자세한 내용 보기

쉽게 설명한 파티클 필터(particle filter) 동작 원리와 예제

이미지
파티클 필터(particle filter)는 칼만 필터(kalman filter) 와 마찬가지로 노이즈가 있는 환경에서 측정된 데이터를 필터를 사용해 실제 위치를 추정하는 도구다. 파티클 필터(particle filter)는 보통 가우시안 분포가 아닌 측정 데이터를 다루기 위해 사용된다고 한다. 물론 가우시안 분포의 데이터에서 사용하지 말라는 건 아니다.  본 포스트에서는 파티클 필터(particle filter)의 어려운 수학적인 내용은 제하고, 쉽게 예를 들며 필터의 동작원리에 대해 알아보았다. 파티클 필터(particle filter)에 대해 검색해 보면 아래와 같은 그림을 많이 보게 된다. 아래 그림은 파티클 필터의 estimation cycle을 도식화한 것이다. <’Real-Time Tracking of Multiple Moving Objects Using Particle Filters and Probabilistic Data Association’ original scientific paper 中> 위 그림의 검은색 원들은 particle을 의미한다. 이 particle에는 보통 위치 데이터와 weight가 포함된다. 검은색 원이 큰 것이 있고 작은 것이 있는 이유는 weight가 크고 작음을 의마한다. typedef struct _particle_t {    int x;    int y;    float weight; }particle_t; 노이즈가 있는 환경에서 레이저 센서나 레이더등을 이용하여 물체의 위치를 측정할 때를 가정해 순서대로 이 필터가 동작하는 과정에 대해 알아 보겠다. (순서는 위 그림과 조금 다르다) 1. 초기 상태 측정 범위안에 랜덤 혹은 일정한 간경으로 particle을 뿌려 놓는다. 이 떄 particle의 weight는 모두 0으로 한다.   for(i
자세한 내용 보기

mkfs.fat Device or resource busy 에러 해결법

이미지
 리눅스에서 SD 카드를 포맷할 때 아래와 같이 "Device or resource busy"에러를 출력하며 포맷이 되지 않는 경우가 있다. mkfs.fat: unable to open /dev/mmcblk0p1: Device or resource busy 해결법을 찾아보았다. 이 문제의 원인은  1. 파티션 dev path 가 맞지 않거나,  2. umount가 완전히 이루어 지지 않았거나,  3. mbr파티션에 문제가 있어 발생한다고 한다.  파티션에 문제가 있는 경우 fdisk를 사용하면 문제를 해결할 수 있다. umount에 문제가 있는 경우, umount할 때 -l 옵션이 사용되는 경우 umount 명령이 끝나도 커널 내부에서 리소스가 계속 사용될 가능성이 있다고 한다.  # umount -l /mnt/mmc umount2("/mnt/mmc",MNT_DETACH); 이런 경우 일정 시간 후 다시 시도하면 에러 없이 포맷이 되는 경우가 있다. 또, umount 가 되지 않을 때 "lsof" , "fuser" 명령어를 사용하여 어떤 프로세스가 mmc를 사용하고 있는지도 확인 가능하다. * lsof에서 사용 내역이 없다고 나와도, fuser에서 나올수 있다.  lsof, fuser 명령에서도 mmc를 사용하는 프로세스가 없는 경우에도 "Device or resource busy"에러가 발생할 경우, umount에 -lf 옵션을 두어 강제로 언마운트 후 SD 카드를 포맷하면 "Device or resource busy"에러가 나타나지 않을 확률이 높아진다. # umount -lf /mnt/mmc umount2("/mnt/mmc",MNT_DETACH | MNT_FORCE); * fat32에서는 "-f","
자세한 내용 보기

windows에서 간단하게 크롬캐스트(Chromecast)를 통해 윈도우 화면 미러링 방법

이미지
 Windows 10에서는 미라캐스트를 사용한 무선 디스플레이연결을 지원한다. 하지만, 만약 크롬캐스트(Chromecast)를 사용해 윈도우 화면을 미러링 하고자 할 때, 구글 크롬 브라우저를 사용하면 쉽고 편하게 할 수 있다. https://www.google.co.kr/intl/ko/chrome/ VLC에서도 비디오 화면을 미러링 할 수 있다고 하는데, 잘 동작하지 않고 윈도우 화면 전체를 미러링하지 않아 크롬 브라우저를 사용하고 있다. 크롬 캐스트 윈도우 화면 전송 방법 0. 크롬 캐스트 디바이스와 windows pc는 같은 네트웍에 연결되어 있어야 한다. (같은 공유기에 연결되어 있어야 한다.) 1. 크롬 브라우저의 우측 상단의 메뉴 > 전송을 클릭하면, 현재 사용가능한 크롬캐스트 디바이스 항목을 볼 수 있다. 2. 편의에 따라 탭 전송 아이콘에서 우측 버튼을 클릭하여 항상 표시되게 설정할 수도 있다. 3.  윈도우 화면 전체를 전송하고자 할 때는 탭 전송의 소스에서 '데스크톱 전송'을 선택하고, 바로 위의 크롬캐스트 디바이스를 클릭하면 된다. ** 소스는 미러링 중에는 변경되지 않는다. 먼저 소스를 선택하고 미러링해야 한다.  4. 화면 전송을 끝내고자 할 때는, 탭 전송에서 화면 전송중인 크롬캐스트 디바이스를 클릭하거나, 크롬 브라우저를 종료하면 된다.
자세한 내용 보기