你自己实现一个FFT的C函数

我自己实现一个快速傅立叶变换（FFT）的C函数，通常会采用Cooley-Tukey算法，这是一种分治策略，将大数组分解成较小的部分然后递归处理。以下是一个简单的C语言版本的单通道一维离散傅立叶变换（DFT）的实现，适用于复数输入：

#include <complex.h>
#include <stdio.h>

#define N 1024 // 变换的长度

// 带符号长整型用于存储复数分量，以支持较大的数据范围
typedef long double Complex;

void fft(Complex *in, Complex *out) {
    if (N <= 1) {
        out = in;
        return;
    }

    Complex t;
    Complex *even = in + 0, *odd = in + N / 2;

    // 分治法：递归地计算偶数和奇数索引部分的FFT
    fft(even, out + 0);
    fft(odd, out + N / 2);

    for (int k = 0; k < N / 2; ++k) {
        t = cexp(-2 * M_PI * I * k / N) * odd[k]; // 计算乘积
        out[k] = even[k] + t; // 将结果相加到偶部
        out[k + N / 2] = even[k] - t; // 将结果相减到奇部
    }
}

// 为了使用，需要定义常数I为虚数单位（i）
const double I = sqrt(-1.0);

int main() {
    Complex input[N], output[N];

    // 初始化输入数组
    for (int i = 0; i < N; ++i) {
        input[i] = (Complex){sin(i * 2 * M_PI / N), cos(i * 2 * M_PI / N)};
    }

    fft(input, output);

    printf("Input array:\n");
    for (int i = 0; i < N; ++i) {
        printf("%g+%gi\n", creal(input[i]), cimag(input[i]));
    }

    printf("\nOutput array (after FFT):\n");
    for (int i = 0; i < N; ++i) {
        printf("%g+%gi\n", creal(output[i]), cimag(output[i]));
    }

    return 0;
}