C 语言实现 DFT 和 FFT

fft.c

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

#define PI 3.14159265358979323846

// O(n^2)
void dft(double *input, double *real, double *imag, int n) {
    for (int k = 0; k < n; k++) {
        real[k] = 0;
        imag[k] = 0;
        for (int j = 0; j < n; j++) {
            double angle = 2 * PI * k * j / n;
            real[k] += input[j] * cos(angle);
            imag[k] -= input[j] * sin(angle);
        }
    }
}

// O(nlogn)
void fft(double *input, double *real, double *imag, int n) {
    if (n == 1) {
        real[0] = input[0];
        imag[0] = 0;
        return;
    }
    double even_real[n/2], even_imag[n/2];
    double odd_real[n/2], odd_imag[n/2];
    for (int i = 0; i < n/2; i++) {
        even_real[i] = input[2*i];
        even_imag[i] = 0;
        odd_real[i] = input[2*i+1];
        odd_imag[i] = 0;
    }
    fft(even_real, even_real, even_imag, n/2);
    fft(odd_real, odd_real, odd_imag, n/2);
    for (int k = 0; k < n/2; k++) {
        double angle = -2 * PI * k / n;
        double w_real = cos(angle);
        double w_imag = sin(angle);
        real[k] = even_real[k] + w_real * odd_real[k] - w_imag * odd_imag[k];
        imag[k] = even_imag[k] + w_real * odd_imag[k] + w_imag * odd_real[k];
        real[k+n/2] = even_real[k] - w_real * odd_real[k] + w_imag * odd_imag[k];
        imag[k+n/2] = even_imag[k] - w_real * odd_imag[k] - w_imag * odd_real[k];
    }
}

int main() {
    double input[] = {1, 2, 3, 4, 5, 6, 7, 8};
    int n = sizeof(input) / sizeof(double);
    double real[n], imag[n];
    fft(input, real, imag, n);
    for (int i = 0; i < n; i++) {
        printf("X[%d] = %f + %fi\n", i, real[i], imag[i]);
    }
    return 0;
}