소스 10-7은 고속 이산 푸리에 변환 알고리즘을 C/C++로 구현한 FFT1d 함수이다. 이 함수의 인자는 DFT1d 함수와 동일하다. 즉, 입력 데이터의 실수부와 허수부가 각각 re와 im으로 전달되고, 데이터의 개수는 N에 저장된다. 네 번째 인자 dir의 값이 +1이면 순방향 FFT를 수행하고, -1이면 역방향 FFT를 수행한다. 이산 푸리에 변환된 결과는 다시 re와 im에 저장된다. FFT1d 함수의 선언은 앞서 10.2절에서 소스 10-3에 나타내었다. 참고로 입력 데이터 순서 바꾸기 코드는 전통적인 비트 반전 방법 대신 좀 더 효율적이고 널리 사용되고 있는 코드를 사용하였다.
#include <algorithm> void FFT1d(double* re, double* im, int N, int dir) { register int i, j, k; //------------------------------------------------------------------------- // 입력 데이터 순서 바꾸기 //------------------------------------------------------------------------- int n2 = N >> 1; int nb = 0; while (N != (1 << nb)) nb++; for (i = 0, j = 0; i < N - 1; i++) { if (i < j) { std::swap(re[i], re[j]); std::swap(im[i], im[j]); } k = n2; while (k <= j) { j -= k; k >>= 1; } j += k; } //------------------------------------------------------------------------- // 버터플라이(Butterfly) 알고리즘 //------------------------------------------------------------------------- int i1, l, l1, l2; double c1, c2, t1, t2, u1, u2, z; c1 = -1.0; c2 = 0.0; l2 = 1; for (l = 0; l < nb; l++) { l1 = l2; l2 <<= 1; u1 = 1.0; u2 = 0.0; for (j = 0; j < l1; j++) { for (i = j; i < N; i += l2) { i1 = i + l1; t1 = u1 * re[i1] - u2 * im[i1]; t2 = u1 * im[i1] + u2 * re[i1]; re[i1] = re[i] - t1; im[i1] = im[i] - t2; re[i] += t1; im[i] += t2; } z = u1 * c1 - u2 * c2; u2 = u1 * c2 + u2 * c1; u1 = z; } c2 = sqrt((1.0 - c1) / 2.0); if (dir == 1) // Forward c2 = -c2; c1 = sqrt((1.0 + c1) / 2.0); } if (dir == -1) // IDFT { for (i = 0; i < N; i++) { re[i] /= static_cast<double>(N); im[i] /= static_cast<double>(N); } } }