FFT 알고리즘
FFT를 설명하기 위하여,
1) 푸리에 급수와 변환에 관한 공식들을 정리하고,
2) 이산 시간 푸리에 변환에 관하여 설명하고,
3) 고속 푸리에 변환에 관하여 설명하고자 한다.
참고사이트 :
| 푸리에 급수
주기 함수 f(t)는 주기가 T 일 때, 다음과 같은 삼각급수로 표현될 수 있다.
n이 증가함에 따라, 푸리에급수는 본 함수 f(x)에 가까워진다.
이해를 돕기 위해 다음 주기함수를 푸리에급수로 표현한 예를 살펴보자.
본래의 함수를 -180부터 180까지 그려보면 다음과 같다.
이것의 푸리에급수의 계수는 암산해보면 다음과 같다.
즉, 푸리에급수는 다음과 같은 함수가 계속적으로 더해져서 만들어진다.
만약 이러한 주기 신호가 1초 걸렸다면,
n=1 에 더해지는 함수는 1 Hz를 나타내고,
n=3 에 더해지는 함수는 3 Hz를 나타낸다.
| 복소 푸리에급수
푸리에급수는 오일러(Euler) 공식에 의해 지수함수 형태로 표현될 수 있다.
| 푸리에 변환
| 이산 푸리에 변환 (Descrete Fourier tranform : DFT)
함수 f(x)를 주기가 1인 주기함수라고 하고,
일정하게 등간격으로 0 부터 N-1까지의 구간에서 f(x)값을 N개 측정(추출)했다고 하자.
잘못된 신호(aliasing)를 피하기 위해서는 n 은 N/2보다 훨씬 작은 값을 사용해야 한다.
| 고속 푸리에 변환 (Fast Fourier Transform:FFT)
FFT 소스:
struct complex { double r,i; }; -
int FFT(complex f[], int N) { int i,j,k,it,xp,xp2,j1,j2,r; double w,wr,wi,dr1,dr2,di1,di2,tr,ti,arg; -
/* 데이터 수는 2 이상 */ if(N < 2) return -1; r = log((double)N)/log(2.0); /* n=2^r */ -
/* N이 2의 제곱수인가 검사 */ j = 1; for(i=0; i < r; i++) j *= 2; if(fabs(n-j) > 1.0E-6) return -2; -
/* FFT */ xp2 = N; for(it=0; it < r; it++) { xp = xp2; xp2 /= 2; /* 절반 */ w = PI / xp2; for(k=0; k < xp2; k++) { arg = k * w; wr = cos(arg); wi = -sin(arg); i = k - xp; for(j=xp; j <= N; j+=xp) { j1 = j + i; /* j1 = k,...,N-k+xp */ j2 = j1 + xp2; /* j1에서 xp2만틈 떨어진 */ dr1 = f[j1].r; dr2 = f[j2].r; di1 = f[j1].i; di2 = f[j2].i; tr = dr1 - dr2; ti = di1 - di2; f[j1].r = dr1 + dr2; f[j1].i = di1 + di2; /* (a+bi)(c+di) = (ac-bd)+(bc+ad)i */ f[j2].r = tr*wr - ti*wi; /* (ac-bd) */ f[j2].i = ti*wr + tr*wi; /* (bc+ad) */ } } } return 0; } -
|
|
댓글
댓글 쓰기