C语言如何使用FFT

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

typedef struct {
    double real;
    double imag;
} complex;

// 复数加法
complex add(complex a, complex b) {
    complex result;
    result.real = a.real + b.real;
    result.imag = a.imag + b.imag;
    return result;
}

// 复数减法
complex subtract(complex a, complex b) {
    complex result;
    result.real = a.real - b.real;
    result.imag = a.imag - b.imag;
    return result;
}

// 复数乘法
complex multiply(complex a, complex b) {
    complex result;
    result.real = a.real * b.real - a.imag * b.imag;
    result.imag = a.real * b.imag + a.imag * b.real;
    return result;
}

// 检查是否为2的幂
int is_power_of_two(int n) {
    return (n > 0) && ((n & (n - 1)) == 0);
}

// 位反转函数
unsigned int reverse_bits(unsigned int x, int log2n) {
    unsigned int n = 0;
    for (int i = 0; i < log2n; i++) {
        n <<= 1;
        n |= (x & 1);
        x >>= 1;
    }
    return n;
}

// 位反转排列
void bit_reverse(complex *x, int n) {
    int log2n = log2(n);
    for (int i = 0; i < n; i++) {
        int j = reverse_bits(i, log2n);
        if (i < j) {
            complex temp = x[i];
            x[i] = x[j];
            x[j] = temp;
        }
    }
}

// FFT函数
void fft(complex *x, int n) {
    if (!is_power_of_two(n)) {
        fprintf(stderr, "错误：输入长度必须是2的幂。\n");
        return;
    }

    bit_reverse(x, n);

    int log2n = log2(n);
    for (int s = 1; s <= log2n; s++) {
        int m = 1 << s; // 当前块大小：2^s
        int half_m = m >> 1;
        // 基础旋转因子：W_m = e^(-2πi/m)
        double theta = -2 * M_PI / m;
        complex w_m = {cos(theta), sin(theta)};

        for (int k = 0; k < n; k += m) { // 处理每个块
            complex w = {1.0, 0.0}; // W_m^0 = 1
            for (int j = 0; j < half_m; j++) {
                complex t = multiply(w, x[k + j + half_m]);
                complex u = x[k + j];
                x[k + j] = add(u, t);
                x[k + j + half_m] = subtract(u, t);
                w = multiply(w, w_m); // 更新旋转因子：W_m^(j+1)
            }
        }
    }
}

// 示例：计算4点FFT
int main() {
    int n = 4;
    complex x[] = {{1, 0}, {2, 0}, {3, 0}, {4, 0}};

    fft(x, n);

    printf("FFT结果：\n");
    for (int i = 0; i < n; i++) {
        printf("x[%d] = %.2f + %.2fi\n", i, x[i].real, x[i].imag);
    }

    return 0;
}