C语言手写代码实现傅里叶变换

下面是一个简单的 C 语言程序，用于计算离散傅里叶变换（DFT）：

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

#define PI 3.14159265358979323846

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

complex complex_add(complex a, complex b) {
    complex c;
    c.real = a.real + b.real;
    c.imag = a.imag + b.imag;
    return c;
}

complex complex_sub(complex a, complex b) {
    complex c;
    c.real = a.real - b.real;
    c.imag = a.imag - b.imag;
    return c;
}

complex complex_mul(complex a, complex b) {
    complex c;
    c.real = a.real * b.real - a.imag * b.imag;
    c.imag = a.real * b.imag + a.imag * b.real;
    return c;
}

void dft(complex *x, int N) {
    complex X[N];
    for (int k = 0; k < N; k++) {
        X[k].real = 0.0;
        X[k].imag = 0.0;
        for (int n = 0; n < N; n++) {
            complex Wnk;
            Wnk.real = cos(2 * PI * n * k / N);
            Wnk.imag = -sin(2 * PI * n * k / N);
            X[k] = complex_add(X[k], complex_mul(x[n], Wnk));
        }
    }
    for (int k = 0; k < N; k++) {
        x[k] = X[k];
    }
}

int main() {
    complex x[] = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}};
    int N = sizeof(x) / sizeof(complex);
    dft(x, N);
    for (int n = 0; n < N; n++) {
        printf("%f + %fi\n", x[n].real, x[n].imag);
    }
    return 0;
}

这个程序使用了一个复数结构体 `complex`，定义了三个复数运算 `complex_add`、`complex_sub` 和 `complex_mul`，以及一个离散傅里叶变换函数 `dft`。

在 `dft` 函数中，我们先定义了一个临时的复数数组 `X`，用于存储变换后的结果。然后，我们对每个频率分量 `k`，计算出对应的旋转因子 `Wnk`，并使用 `X[k]` 累加所有时域信号 `x[n]` 乘上旋转因子的结果。最后，我们再把计算出来的结果赋值回原始的信号数组 `x` 中。

在 `main` 函数中，我们定义了一个长度为 4 的测试信号数组 `x`，并调用 `dft` 函数对它进行了变换。最后，我们输出了变换后的结果。