用C语言设计FFT

实现FFT算法的C语言代码相对较长，这里只提供一个简单的例子，仅供参考。

#include <stdio.h>
#include <math.h>

#define PI 3.1415926535897932384626

int reverse(int N, int n) {
    int j, p = 0;
    for (j = 1; j <= log2(N); j++) {
        if (n & (1 << (int)(log2(N) - j)))
            p |= 1 << (j - 1);
    }
    return p;
}

void FFT(double *x_real, double *x_imag, int N) {
    int i, j, k;
    int m, p, q, r;
    double arg, cosarg, sinarg;
    double *p_real, *p_imag, *q_real, *q_imag, *t_real, *t_imag;

    p_real = (double *)malloc(N * sizeof(double));
    p_imag = (double *)malloc(N * sizeof(double));
    q_real = (double *)malloc(N * sizeof(double));
    q_imag = (double *)malloc(N * sizeof(double));
    t_real = (double *)malloc(N * sizeof(double));
    t_imag = (double *)malloc(N * sizeof(double));

    for (i = 0; i < N; i++) {
        p_real[i] = x_real[reverse(N, i)];
        p_imag[i] = x_imag[reverse(N, i)];
    }

    for (m = 1; m <= log2(N); m++) {
        for (k = 0; k < pow(2, m - 1); k++) {
            for (q = k; q < N; q += pow(2, m)) {
                j = q + pow(2, m - 1);
                arg = -2 * PI * k / pow(2, m);
                cosarg = cos(arg);
                sinarg = sin(arg);
                t_real[q] = p_real[q] + cosarg * p_real[j] - sinarg * p_imag[j];
                t_imag[q] = p_imag[q] + cosarg * p_imag[j] + sinarg * p_real[j];
                t_real[j] = p_real[q] - cosarg * p_real[j] + sinarg * p_imag[j];
                t_imag[j] = p_imag[q] - cosarg * p_imag[j] - sinarg * p_real[j];
            }
            for (q = 0; q < N; q++) {
                p_real[q] = t_real[q];
                p_imag[q] = t_imag[q];
            }
        }
    }

    for (i = 0; i < N; i++) {
        x_real[i] = p_real[i];
        x_imag[i] = p_imag[i];
    }

    free(p_real);
    free(p_imag);
    free(q_real);
    free(q_imag);
    free(t_real);
    free(t_imag);
}

int main() {
    int N = 8;
    double x_real[8] = {1, 2, 3, 4, 5, 6, 7, 8};
    double x_imag[8] = {0};

    FFT(x_real, x_imag, N);

    for (int i = 0; i < N; i++) {
        printf("X[%d] = %f + %fi\n", i, x_real[i], x_imag[i]);
    }

    return 0;
}

这段代码实现的是一个基于蝴蝶算法的FFT。需要注意的是，这里的 FFT 函数输入的实部和虚部数组都要有 N 个元素，而且 N 必须是 2 的幂次方。在代码中，reverse 函数是用来翻转位序的，以便后续做蝴蝶算法的时候能够按位位置正确地计算。在 FFT 函数中，先将输入数组按位翻转，然后进行蝴蝶算法的迭代计算，最后将计算结果写回输入数组中。