푸리에 분석 - 고속 푸리에 변환 구현하기

이 페이지는 최병선 교수님의 Wavelet 해석을 바탕으로 작성되었습니다. 그러나 책과는 달리 해석적으로 엄밀한 내용을 추구하지는 않습니다. 이산 푸리에 변환의 구현은 여기를 참고 바랍니다.

FFT 구현하기

이번 페이지에서는 Cooley-Tukey 알고리즘을 C로 구현해보도록 하겠습니다. FFT 자체는 많은 수치해석 라이브러리에 구현되어 있습니다. 그렇지만 저는 공부하는 차원에서 직접 구현해본다면 더 깊이 있게 라이브러를 이해하고 쓸 수 있다고 생각하기 때문에, 이 작업이 그렇게 무의미하다고 생각하지는 않습니다(보통은 이걸 ‘바퀴의 재발명’이라고 부릅니다).

Cooley-Tukey 알고리즘은 데이터의 갯수가 2의 거듭제곱으로 떨어졌을 때 구현하기가 쉽습니다. 보통 데이터의 갯수는 크게 문제가 되지 않는데, 일반적으로 FFT를 사용하는 경우는 신호처리가 필요할 때이고, 신호처리가 필요하다면 FFT의 구현에 맞게 신호를 샘플링하면 되기 때문입니다.

실수와 복소수 정의하기

먼저, 실수를 정의하겠습니다. 그런데 사실상 float와 double 외에는 선택지가 없습니다. 여기서는 실수를 대표하는 자료형으로 double을 쓰도록 하겠습니다.

typedef double Real;

복소수 집합은 사실 2차원 실수 집합과 거의 똑같습니다(이런 걸 Isomorphism 이라고 합니다). 다만 x축은 실수, y축은 허수에 대응시켰을 뿐이지요. 그렇기 때문에 복소수는 실수 2개를 묶어두면 됩니다. 그러므로 다음과 같이 정의하는 게 합당할 것 같습니다.

typedef struct _complex
{
    Real re;
    Real im;
} Complex;

여기에서 re는 실수부, im은 허수부입니다.

복소수를 정의하였으니, 복소수의 연산을 위한 몇 가지 기본적인 함수들을 정의하겠습니다.

Real cAbsolute(Complex a);
Real cArgument(Complex a);
Complex cAdd(Complex a, Complex b);
Complex cSubtract(Complex a, Complex b);
Complex cMultiply(Complex a, Complex b);
Complex cDivide(Complex a, Complex b);
Complex cNegate(Complex a);
Complex cPolarToComplex(Real length, Real angle);

벡터 정의하기

FFT의 계산에는 벡터 연산이 들어갑니다. 그런데 C에는 벡터라는 건 없습니다(C++에는 vector라는 컨테이너가 있고 수학의 벡터처럼 쓸 수 있기는 합니다). 그래서 직접 정의해줘야 합니다.

벡터는 임의의 길이의 숫자의 묶음이니까 배열로 정의하면 적당할 것 같습니다. 그런데 벡터의 길이도 알면 좋을 것 같으니, 구조체로 배열과 배열의 길이를 같이 묶어주면 좋을 것 같습니다. 그런데 임의의 길이의 배열을 정적으로 선언할 수는 없습니다. 그렇기 때문에 배열의 포인터와 배열의 길이를 묶어줘야 합니다.

typedef struct _vec {
    unsigned long long length;
    Complex* arr;
} vec;
typedef vec* vecRef;

vecRef라는 것은 벡터의 레퍼런스를 의미합니다. 쉽게 말해 자바의 객체처럼 쓰려고 정의한 겁니다.

그 외에 vecRef를 사용하기 위한 몇 가지 함수들을 정의하겠습니다.

vecRef vecAlloc(unsigned long long length);
void   vecRelease(vecRef* arr);

void      vadd(vecRef a, vecRef b, vecRef* output);
void vsubtract(vecRef a, vecRef b, vecRef* output);
void vmultiply(vecRef a, vecRef b, vecRef* output);
void   vdivide(vecRef a, vecRef b, vecRef* output);

void vaddScalar(vecRef a, Complex c, vecRef* output);
void vsubScalar(vecRef a, Complex c, vecRef* output);
void vmulScalar(vecRef a, Complex c, vecRef* output);
void vdivScalar(vecRef a, Complex c, vecRef* output);

느린 푸리에 변환 구현하기

이제 준비는 끝났습니다. 고속 푸리에 변환을 구현하기에 앞서, 느린 버전의 이산 푸리에 변환을 먼저 구현해보도록 하겠습니다. 이산 푸리에 변환은 이전의 포스트에서 설명했듯이, 행렬을 곱하는 것이 전부입니다. \begin{align} F = \frac{1}{\sqrt{N}} \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_1 & \omega_1^2 & \cdots & \omega_1^{N-1} \\ 1 & \omega_2 & \omega_2^2 & \cdots & \omega_2^{N-1} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & \omega_{N-1} & \omega_{N-1}^2 & \cdots & \omega_{N-1}^{N-1} \\ \end{bmatrix} \end{align}

그런데 사실 구현의 관점에서 보면 저 행렬을 정의하고, 행렬의 곱셈 같은 연산을 굳이 정의해줄 필요가 없습니다. 왜냐하면 저 행렬 $F$ 를 만드는 규칙을 저희는 이미 알고 있거든요. 게다가 $N$ 이 매우 커진다면, 행렬 $F$ 를 보관하기 위해 할당해야 하는 메모리도 매우 커질 것입니다. 이미 규칙을 아는 상황에서는 낭비일 뿐입니다.

행렬의 곱셈을 명시적으로 구현하지 않은 DFT의 구현은 다음과 같습니다.

void dft(const vecRef input, vecRef *output)
{
    assert(input->length >= 2);
    Long N = input->length;
    
    // If output is null, we make it
    // If output is not null and if it is not power of 2, we reallocate it
    // Else, use output
    if (*output == NULL || (*output)->length != N)
    {
        vecRelease(output);
        *output = vecAlloc(N);
    }
    
    for (Long i = 0; i < N; i++)
    {
        Complex tmp = ComplexMake(0, 0);
        for (Long j = 0; j < N; j++)
        {
            Complex wj = cPolarToComplex(1.0, -(2.0*M_PI*i*j)/N);
            tmp = cAdd(tmp, cMultiply(wj, input->arr[j]));
        }
        (*output)->arr[i] = tmp;
    }
    vdivScalar(*output, ComplexMake(sqrt(N), 0.0), output);
}

고속 푸리에 변환 구현하기

이제 고속 푸리에 변환을 구현해보겠습니다. Cooley-Tukey 알고리즘은 다음의 수학적 사실에 기반합니다. \[ F_N = \frac{1}{\sqrt{N}} G_0^{N} G_1^{N} \cdots G_{q-1}^{N} Q_N \]

위 식을 분해해보면, 우리가 구현해야 할 부분은 두 개입니다: $Q_{N}$ 과 $G_{k}^{N}$ 입니다.

비트 반전 구현하기

$Q_{N}$ 은 벡터의 원소의 순서를 바꾸는 것으로 구현해야 합니다. 그러기 위해서는 각각의 원소의 인덱스를 이진수로 바꾸고, 그 이진수의 비트를 반전시켜야 합니다. 이진수를 반전시키는 알고리즘은 구현이 상당히 잘 되어 있습니다. 저는 스택오버플로우에서 답으로 올라온 구현 중 하나를 참고하였습니다.

// http://stackoverflow.com/a/34236981
Long reverse(const Long n, const Long k)
{
    Long r, i;
    for (r = 0, i = 0; i < k; ++i)
        r |= ((n >> i) & 1) << (k - i - 1);
    return r;
}

이 함수는 최대 몇 개의 비트까지 뒤집을 것인가를 선택할 수 있습니다. 저는 이 구현에서 정수형 자료형으로 64비트 unsigned long long을 썼는데, 만일 64비트를 다 뒤집게 되면 정말 끔찍한 일이 벌어지게 됩니다. 1을 뒤집으면 $2^{63}$ 이 나오게 되는데, 이 정도의 숫자라면 거의 100%의 확률로 배열의 범위를 넘어서 참조하게 됩니다(나머지 확률은 제가 진짜로 이 정도 크기의 배열을 메모리에 동적할당했을 경우인데, 이마저도 64비트 컴퓨터를 쓰는 제 노트북 환경에서는 택도 없습니다).

이제 우리는 $Q_{N}$ 을 구현할 수 있습니다.

// http://stackoverflow.com/a/34236981
Long reverse(const Long n, const Long k)
void permutate(const int q, const vecRef input, vecRef *output)
{
    for (Long i = 0; i < (1 << q); i++)
    {
        (*output)->arr[i] = input->arr[reverse(i, q)];
    }
}

복소수 연산 구현하기

이제 진짜 푸리에 변환을 구현할 차례입니다. 행렬 $G_{k}^{N}$ 를 구현하면 되는데, 일단 정의만 보면 \[ G_{k}^{N} = \underbrace{E_{N/2^{k}} \bigoplus E_{N/2^{k}} \bigoplus \cdots \bigoplus E_{N/2^{k}}}_{2^{k}\text{개}}\]

이 행렬이 어떻게 생겨먹었는지 감이 잘 와닿지 않습니다. 그래서 한 번 $G_2^8$ 을 써보았습니다. \[ \left[\begin{array}{rrrr|rrrr} \;\;\;1 & \;\;\;0 & \omega^{0} & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & \omega^{1} & 0 & 0 & 0 & 0 \\ 1 & 0 & -\omega^{0} & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & -\omega^{1} & 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 & \;\;\;1 & \;\;\;0 & \omega^{0} & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & \omega^{1} \\ 0 & 0 & 0 & 0 & 1 & 0 & -\omega^{0} & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & -\omega^{1} \end{array}\right]\]

이 행렬의 연산을 구현하는 방법 중 하나는 행렬을 직접 구현하는 것입니다. 행렬을 직접 구현하고, 곱하기도 구현하면 제일 편합니다. 그런데 이 방법은 그다지 좋은 방법이 아닌 게, 보면 알겠지만 0이 너무 많습니다. 그렇기 때문에 메모리가 낭비될 겁니다. 게다가, 이 행렬을 보면 알겠지만 반복되는 부분이 있습니다. 반복되는 부분만 적절히 구현해주면, 굳이 행렬을 정의하는 방식으로 구현하지 않아도 될 것 같습니다.

$G_2^8$ 은 $E_4$ 가 두 개 붙어있는 형태입니다. 그렇기 때문에 $E_4$ 만 구현해주면 됩니다. 이를 일반화하면, $G_k^N$ 을 구현하려면 $E_{N/2^{k}}$ 만 구현하면 됩니다.

감을 잡기 위해, $E_{N/2^{k}} \cdot v$ 를 한 번 써보겠습니다. \begin{align} \begin{bmatrix} 1 & 0 & \omega^{0} & 0 \\ 0 & 1 & 0 & \omega^{1} \\ 1 & 0 & -\omega^{0} & 0 \\ 0 & 1 & 0 & -\omega^{1} \end{bmatrix} \cdot \begin{bmatrix} v_0 \\ v_1 \\ v_2 \\ v_3 \\ \end{bmatrix} = \begin{bmatrix} v_0 + \omega^{0} v_2 \\ v_1 + \omega^{1} v_3 \\ v_0 - \omega^{0} v_2 \\ v_1 - \omega^{1} v_3 \\ \end{bmatrix} \end{align}

0번째 행과 2번째 행, 1번째 행과 3번째 행은 상당히 닮아 있습니다. 가운데에 $+$ 냐 $-$ 냐만 빼면 연산에 들어간 숫자들이 똑같습니다. 이 관찰을 일반화하면 다음과 같습니다: $E_{N/2^{k}}$ 를 구현할 땐 $0$ 부터 $(\frac{N}{2^{k+1}}-1)$ 행까지만 반복문을 돌면 된다.

이제 $G_k^N$ 을 구현하겠습니다. $E_{N/2^{k}}$ 를 하나의 블록이라고 생각하면, 이 블록은 총 $2^k$ 번 반복됩니다. 그러므로, $E_{N/2^{k}}$ 의 구현을 $2^k$ 번 반복하여 돌면 $G_k^N$ 의 구현이 끝납니다.

void computeFFT(const int q, const vecRef input, vecRef *output)
{
    const Long countBlock = 1 << q;
    const Long lenBlock = (*output)->length / countBlock;
    
    for (Long baseIdx = 0; baseIdx < lenBlock/2; baseIdx++)
    {
        const Complex w = cPolarToComplex(1.0, -(M_PI * baseIdx)/(lenBlock/2));
        for (Long block = 0; block < countBlock; block++)
        {
            const Long idx = baseIdx + block * lenBlock;
            const Complex a = input->arr[idx];
            const Complex b = input->arr[idx + lenBlock/2];
            const Complex wB = cMultiply(w, b);
            // Upper part
            (*output)->arr[idx] = cAdd(a, wB);
            // Lower part
            (*output)->arr[idx + lenBlock/2] = cSubtract(a, wB);
        }
    }
}

묶어주기

이제 핵심적인 부분은 모두 구현하였으므로, 이를 종합하여 FFT의 구현을 마무리해보겠습니다.

void fft(const vecRef input, vecRef *output)
{
    assert(input->length >= 2);
    Long N = input->length;
    if (checkPowerOfTwo(N) != 1)
    {
        N = 1 << (int)(round(log2(N))+1);
    }
    int rounds = (int)log2(N);
    
    // If output is null, we make it
    // If output is not null and if it is not power of 2, we reallocate it
    // Else, use output
    if (*output == NULL || (*output)->length != N)
    {
        vecRelease(output);
        *output = vecAlloc(N);
    }
    
    // 순서 섞기: 행렬 Q를 곱해줍니다
    permutate(rounds, input, output);
    // 계산을 몇 번 돌려야 하므로,
    // 임시 벡터를 하나 만듭니다
    vecRef tmp = vecAlloc(N);
    memcpy(tmp->arr, (*output)->arr, sizeof(Complex) * tmp->length);
    for (int q = rounds-1; q >= 0; q--)
    {
        vecSwap(output, &tmp);
        // 복소수 계산하기: 행렬 G를 곱해줍니다
        computeFFT(q, tmp, output);
    }
    vecRelease(&tmp);
    vdivScalar(*output, ComplexMake(sqrt(N), 0.0), output);
}

DFT vs FFT

이제 DFT와 FFT의 속도 차이를 비교해보겠습니다. 비교를 위한 코드는 다음과 같습니다.

typedef struct {
    clock_t fft;
    clock_t dft;
} elapsedtime;

elapsedtime compare(const unsigned long long q)
{
    vecRef input = vecAlloc(q);
    for (unsigned long long i = 0; i < q; i++)
    {
        Real f = 0.0;
        f += 1.0 * sin((2.0 * M_PI * i) / q);
        f += 4.0 * cos((4.0 * M_PI * i) / q);
        f -= 2.0 * cos((16.0 * M_PI * i) / q);
        input->arr[i] = ComplexMake(f, 0);
    }
    vecRef output = NULL;
    clock_t tStart = clock();
    fft(input, &output);
    clock_t tEnd = clock();
    
    vecRef output2 = NULL;
    clock_t tStart2 = clock();
    dft(input, &output2);
    clock_t tEnd2 = clock();
    
    vecRelease(&input);
    vecRelease(&output);
    vecRelease(&output2);
    
    elapsedtime e;
    e.dft = tEnd2 - tStart2;
    e.fft = tEnd - tStart;
    return e;
}

int main(int argc, const char * argv[])
{
    for(unsigned long long q = 1; q <= 10; q++)
    {
        clock_t tFFT = 0.0;
        clock_t tDFT = 0.0;
        for (int i = 0; i < 1000; i++)
        {
            elapsedtime t = compare(1 << q);
            tFFT += t.fft;
            tDFT += t.dft;
        }
        printf("N: %llu FFT: %f DFT: %f\n", q, tFFT/1000.0, tDFT/1000.0);
    }
    return 0;
}

테스트는 맥북프로 레티나(15인치, 2015년 Mid형, Intel i7 2.5GHz)에서 길이 $N = 2^1, 2^2, \cdots , 2^{10}$ 인 데이터를 대상으로 각각 1000번씩 수행하는 것으로 이루어졌습니다. 그 결과는 아래의 표에 있습니다. \begin{array}{c|cc} N & \text{FFT} & \text{DFT} \\ \hline 2^{1} & 0.841 & 0.627 \\ 2^{2} & 1.084 & 1.268 \\ 2^{3} & 1.695 & 3.828 \\ 2^{4} & 2.897 & 14.219 \\ 2^{5} & 6.058 & 57.979 \\ 2^{6} & 12.669 & 229.600 \\ 2^{7} & 29.593 & 933.523 \\ 2^{8} & 67.785 & 3451.708 \\ 2^{9} & 138.741 & 13810.502 \\ 2^{10} & 286.715 & 55147.116 \end{array}

역시 $O(N \log N)$ 의 위엄은 대단합니다. 데이터가 커지면 커질수록 FFT가 월등히 빠릅니다.

참고문헌

최병선, Wavelet해석, 세경사, 2001

Best Algorithm for Bit Reversal ( from MSB->LSB to LSB->MSB) in C