(C++) FFT 알고리즘 구현

왜 직접 짜보나

FFT는 라이브러리를 쓰면 된다. FFTW는 빠르고 검증돼 있다. 그런데 두 가지 이유로 직접 짜볼 일이 생긴다.

하나는 라이선스다. FFTW는 GPL이라 상용 프로그램에 넣으려면 별도 라이선스를 사야 한다. 다른 하나는 배포다. 보드에 올릴 프로그램에 외부 의존성을 하나라도 줄이고 싶을 때가 있다.

무엇보다 안에서 뭘 하는지 한 번은 알아야 결과를 해석할 수 있다. 빈 번호를 주파수로 바꾸거나, 창을 걸었을 때 진폭이 왜 변하는지 같은 건 블랙박스로 두면 계속 걸린다.

DFT를 반으로 쪼개기

DFT의 정의는 이렇다.

\[X(k) = \sum_{n=0}^{N-1} x(n)\, e^{-j 2\pi k n / N}, \quad k = 0, 1, \dots, N-1\]

• $X(k)$는 $k$번째 복소 계수
• $x(n)$은 $n$번째 복소 입력 샘플
• $N$은 샘플 개수

정의대로 계산하면 $k$마다 $N$번 곱하니 $O(N^2)$다. $N = 512$면 26만 번이다.

Cooley-Tukey는 이 합을 짝수 인덱스와 홀수 인덱스로 나눈다. 각각이 크기 $N/2$의 DFT가 되고, 두 결과를 회전 인자로 합치면 원래 결과가 나온다. 이걸 재귀적으로 반복하면 $O(N \log N)$이 된다. 512면 4600번 정도다.

$N$이 2의 거듭제곱이라는 가정이 여기 들어간다. 계속 반으로 나눠야 하기 때문이다.

재귀 구현

#include <iostream>
#include <vector>
#include <cmath>
#include <complex>

typedef complex<double> Complex;
const double PI = 3.14159265358979323846;

vector<Complex> fft(vector<Complex>& a) {
    int n = a.size();
    if (n == 1)
        return a;
    
    vector<Complex> a_even(n / 2), a_odd(n / 2);
    for (int i = 0; i < n / 2; i++) {
        a_even[i] = a[i * 2];
        a_odd[i] = a[i * 2 + 1];
    }
    vector<Complex> y_even = fft(a_even);
    vector<Complex> y_odd = fft(a_odd);
    vector<Complex> y(n);
    for (int k = 0; k < n / 2; k++) {
        Complex t = polar(1.0, -2.0 * PI * k / n) * y_odd[k];
        y[k] = y_even[k] + t;
        y[k + n / 2] = y_even[k] - t;
    }
    return y;
}

크기가 1이면 그대로 돌려주고(그 자체가 DFT다), 아니면 짝수/홀수로 나눠 각각 재귀 호출하고, 회전 인자를 곱해 합친다.

마지막 루프의 두 줄이 버터플라이다. $k$번째와 $k + N/2$번째를 한 번에 만든다. $e^{-j2\pi(k+N/2)/N} = -e^{-j2\pi k/N}$이라 부호만 바뀌기 때문에 곱셈을 한 번만 하면 된다.

IFFT는 켤레를 이용한다

vector<Complex> ifft(vector<Complex>& a) {
    int n = a.size();
    for (int i = 0; i < n; i++) {
        a[i] = conj(a[i]);
    }
    vector<Complex> y = fft(a);
    for (int i = 0; i < n; i++) {
        y[i] = conj(y[i]) / double(n);
    }
    return y;
}

역변환을 따로 짜지 않아도 된다. 입력을 켤레로 바꾸고 FFT를 돌린 뒤 결과를 다시 켤레로 바꾸고 $N$으로 나누면 IFFT가 된다. 정변환 코드 하나만 유지하면 되니 관리가 편하다.

여기 함정이 하나 있다. 입력 a가 파괴된다. 첫 루프가 제자리에서 켤레를 취하기 때문에, 호출한 쪽의 데이터가 바뀐 채로 남는다.

vector<Complex> spec = fft(x);
vector<Complex> back = ifft(spec);
// 이 시점에 spec 의 내용은 켤레로 바뀌어 있다

spec을 나중에 다시 쓰려고 하면 값이 달라져 있다. 이런 건 한참 뒤에야 발견된다. 인자를 값으로 받거나 const 참조 + 내부 복사로 바꿔야 한다.

fft 쪽도 vector<Complex>&로 받는데 정작 수정은 안 한다. 비-const 참조라 상수나 임시 객체를 못 넘긴다. const vector<Complex>&가 맞다.

확인

int main() {
    vector<Complex> a = {1, 2, 3, 4};
    vector<Complex> y = fft(a);
    for (int i = 0; i < y.size(); i++)
        cout << y[i] << " ";
    cout << endl;
    vector<Complex> x = ifft(y);
    for (int i = 0; i < x.size(); i++)
        cout << x[i] << " ";
    cout << endl;
    return 0;
}

FFT를 하고 다시 IFFT를 하면 원래 값이 나와야 한다. 이 왕복 검사가 FFT 구현이 맞는지 보는 제일 빠른 방법이다.

다만 완전히 같지는 않다. 부동소수점 연산이라 $10^{-15}$ 수준의 오차가 남는다. 자동 검사를 만들 때는 정확한 일치가 아니라 허용 오차로 비교해야 한다.

한 가지 더 확인하면 좋은 게 있다. 알려진 주파수의 정현파를 넣어서 그 빈에만 값이 서는지 보는 것이다. 왕복 검사는 정변환과 역변환이 서로 반대이기만 하면 통과하기 때문에, 회전 인자의 부호가 통째로 뒤집혀 있어도 지나간다. 그러면 나중에 도플러의 접근/이탈 방향이 반대로 나온다.

// 빈 3에 정현파 하나
const int N = 8;
vector<Complex> x(N);
for (int n = 0; n < N; ++n)
    x[n] = std::polar(1.0, 2.0 * PI * 3 * n / N);
// fft(x) 의 결과는 인덱스 3 에만 크기 N, 나머지는 0 이어야 한다

재귀 구현의 대가

동작은 하는데 실전에 넣기에는 몇 가지가 걸린다.

호출마다 벡터를 새로 잡는다. 한 번의 fft 호출에서 a_even, a_odd, y 세 개를 할당하고, 재귀 호출이 반환한 y_even, y_odd까지 받는다. $N = 512$면 재귀 노드가 1023개라 할당이 수천 번 일어난다. 프레임마다 FFT를 도는 코드에서는 이 비용이 계산보다 클 수 있다.

회전 인자를 매번 계산한다. polar(1.0, -2.0 * PI * k / n)이 안쪽 루프에서 매번 불린다. sin과 cos은 싼 연산이 아니다. 값은 $N$과 $k$로 정해지니 미리 표로 만들어두면 된다.

$N$이 2의 거듭제곱이 아니면 조용히 틀린다. 예를 들어 $N = 6$이면 n/2 = 3이라 짝수 3개, 홀수 3개로 나뉘고, 그 3짜리가 또 n/2 = 1로 나뉘면서 샘플 하나가 사라진다. 크래시도 안 나고 에러도 없이 다른 값이 나온다.

검사를 넣어야 한다.

if (n & (n - 1)) throw std::invalid_argument("FFT size must be a power of two");

n & (n-1)이 0이면 2의 거듭제곱이다. 비트가 하나만 서 있는 수에서 1을 빼면 그 아래 비트가 전부 뒤집히기 때문이다.

길이가 안 맞으면 0으로 채워(제로 패딩) 다음 2의 거듭제곱까지 늘린다. 이때 스펙트럼의 빈 간격이 촘촘해지지만 분해능이 좋아지는 건 아니다. 보간이 세밀해지는 것뿐이다.

실전에서는 반복 버전

위 문제들을 없애려면 재귀를 펼쳐서 반복문으로 바꾼다. 입력을 비트 역순으로 재배치한 다음, 크기 2부터 $N$까지 두 배씩 늘려가며 제자리에서 버터플라이를 돌린다.

• 추가 메모리가 없다 (제자리 연산)
• 회전 인자를 미리 표로 만들어 재사용한다
• 함수 호출 오버헤드가 없다

같은 크기의 FFT를 반복해서 돌리는 상황이면 차이가 크다. 실제로 쓰는 형태의 코드는 CW 수신 체인 부록에 Fft 클래스로 정리해뒀다. 생성자에서 비트 역순 인덱스와 회전 인자를 한 번 계산해두고 프레임마다 재사용한다.

회전 인자를 표로 두는 게 성능뿐 아니라 정확도 문제이기도 하다. 교과서 예제처럼 안쪽 루프에서 w *= wl로 회전 인자를 곱해 나가면, $N$번 곱하는 동안 위상 오차가 누적된다. 채널 간 위상차를 0.1° 수준에서 읽어야 하는 작업에서는 이게 그대로 각도 오차가 된다.

다른 FFT 알고리즘

Cooley-Tukey가 가장 널리 쓰이지만 유일한 건 아니다.

• Radix-2: 여기서 구현한 것. $N$이 2의 거듭제곱일 때만 동작한다
• Mixed-Radix: 여러 기수를 섞어 쓴다. 2의 거듭제곱이 아닌 크기에 쓸 수 있다
• Prime Factor Algorithm: $N$이 서로소인 인수들의 곱일 때 효율적이다
• Rader: $N$이 소수일 때 쓴다. 순환 합성곱으로 변환한다
• Bluestein: 크기에 제약이 없다. 처프 z-변환을 이용하는데, 대신 상수 인자가 크다
• Winograd: 작은 크기의 DFT를 곱셈 수 기준으로 최적화한 것. 큰 변환을 쪼갠 뒤 조각에 쓴다

FFTW가 빠른 이유 중 하나가 이것들을 크기에 따라 골라 쓰고, 실행 시간에 실제로 재보고 계획을 세우기 때문이다. 직접 짠 radix-2와 성능 차이가 몇 배 나는 게 이상한 일이 아니다.

실무에서는 대부분 $N$을 2의 거듭제곱으로 잡아버리는 것으로 해결된다. 프레임 길이는 내가 정하는 값이니 512나 1024로 두면 radix-2 하나로 끝난다.

실수 입력이면 절반만 계산해도 된다

측정 데이터가 실수라면 스펙트럼이 켤레 대칭이다.

\[X(N-k) = \overline{X(k)}\]

절반은 계산하지 않아도 나머지 절반에서 나온다. 실수 전용 FFT는 이걸 이용해 계산량과 메모리를 절반으로 줄인다.

다만 IQ 신호를 다루면 입력이 처음부터 복소수라 이 최적화가 적용되지 않는다. 그리고 복소 FFT에서는 상위 절반 빈이 음의 주파수를 뜻해서, 이걸 잘못 읽으면 멀어지는 표적이 초고속으로 다가오는 표적으로 뒤집힌다. 그 부분은 CW 수신 체인 2절에 적어뒀다.

정리하면

• 재귀 구현은 구조를 이해하기 좋지만 호출마다 벡터를 할당하고 삼각함수를 다시 계산한다
• IFFT를 제자리 켤레로 만들면 입력이 파괴된다. 호출한 쪽 데이터가 바뀐다
• $N$이 2의 거듭제곱이 아니면 에러 없이 틀린 값이 나온다. n & (n-1)로 검사한다
• 왕복 검사만으로는 회전 인자의 부호 실수를 못 잡는다. 알려진 주파수의 정현파도 같이 확인한다
• 같은 크기를 반복해서 돌리면 반복 버전 + 회전 인자 테이블이 성능과 정확도 양쪽에서 낫다

FFTW로 실제 측정 데이터를 처리한 이야기는 iFFT & FFT 알고리즘 구현에 있다.