[질문] c 로 fft코드를 작성해

아래는 C로 작성된 FFT 코드입니다.

fft.h 파일:

c
#ifndef FFT_H
#define FFT_H

int fft(long N, double XR[], double XI[]);
char *cplxStr(double re, double im);

#endif

fft.c 파일:

c
#include <math.h>
#include <stdio.h>
#include "fft.h"

#define PI 3.14159265358979
int fft(long N, double XR[], double XI[]) {
    if ((N != 0) && ((N & (N - 1)) != 0))
        return -1;

    long N2 = N >> 1;
    double WR[N2], WI[N2];
    double T = 2 * PI / N;

    for (long i = 0; i < N2; ++i) {
        WR[i] = cos(T * i);
        WI[i] = -sin(T * i);
    }

    int log2N = (int)log2(N);
    double tmp;

    for (long n = 1; n < N - 1; ++n) {
        long m = 0;

        for (int i = 0; i < log2N; ++i) {
            m |= ((n >> i) & 1) << (log2N - i - 1);
        }

        if (n < m) {
            tmp = XR[n];
            XR[n] = XR[m];
            XR[m] = tmp;
        }
    }

    for (int loop = 0; loop < log2N; ++loop) {
        long regionSize = 1 << (loop + 1);
        long kJump = 1 << (log2N - loop - 1);
        long half = regionSize >> 1;

        for (register long i = 0; i < N; i += regionSize) {
            long blockEnd = i + half - 1;
            long k = 0;
            double TR, TI;

            for (register long j = i; j <= blockEnd; ++j) {
                TR = WR[k] * XR[j + half] - WI[k] * XI[j + half];
                TI = WI[k] * XR[j + half] + WR[k] * XI[j + half];
                XR[j + half] = XR[j] - TR;
                XI[j + half] = XI[j] - TI;
                XR[j] = XR[j] + TR;
                XI[j] = XI[j] + TI;
                k += kJump;
            }
        }
    }

    return 0;
}

char *cplxStr(double re, double im) {
    static char s[50];

    if (re == 0 && im == 0) {
        sprintf(s, "0");
    } else if (re == 0) {
        sprintf(s, "%.12fi", im);
    } else if (im == 0) {
        sprintf(s, "%.12f", re);
    } else if (im > 0) {
        sprintf(s, "%.12f + %.12fi", re, im);
    } else {
        sprintf(s, "%.12f - %.12fi", re, -im);
    }

    return s;
}

test_fft.c 파일:

c
#include <stdio.h>
#include <time.h>
#include <math.h>
#include "fft.h"

int test_simple(void) {
    const long N = 8;
    double XR[8] = {0, 0, 0, 0, 1, 1, 1, 1};
    double XI[8];

    printf("** test_simple_fft **\n");
    int result = fft(N, XR, XI);

    for (long i = 0; i < N; ++i) {
        printf("%s, ", cplxStr(XR[i], XI[i]));
    }

    printf("\n");

    return 0;
}

#define CNT 65536
#define PI 3.14159265358979
int test64k(void) {
    long N = CNT;
    static double XR[CNT];
    static double XI[CNT];

    printf("** test_measure_time **\n");

    for (long i = 0; i < N; ++i) {
        XR[i] = cos(2 * PI * 0.004 * i);
    }

    double beforeT, diffT;
    int result = -1;

    beforeT = clock();
    result = fft(N, XR, XI);
    diffT = clock() - beforeT;

    printf("%.0fms\n", diffT);

    if (result < 0) {
        printf("FFT calculation Error\n");
    }

    for (long i = 0; i < 8; ++i) {
        printf("%s ,", cplxStr(XR[i], XI[i]));
    }

    printf("\n");

    for (long i = N - 8; i < N; ++i) {
        printf("%s ,", cplxStr(XR[i], XI[i]));
    }

    return 0;
}

int main() {
    test_simple();
    test64k();

    return 0;
}

위의 코드는 FFT를 구현한 C 코드입니다. fft.h 파일에는 FFT 함수와 cplxStr 함수의 선언이 포함되어 있고, fft.c 파일에는 실제 FFT 알고리즘이 구현되어 있습니다. test_fft.c 파일은 FFT 함수를 테스트하는 예제 코드입니다.

참고 URL: https://infograph.tistory.com/358