使用 C 进行离散傅立叶变换及其逆变换

离散傅立叶变换（Discrete Fourier Transform，DFT）是从一个离散时间序列（或离散信号）中提取频域信息的一种方法。它可以为信号分解成一系列不同频率的正弦函数和余弦函数。同样地，离散傅立叶逆变换（Inverse Discrete Fourier Transform，IDFT）可以将这些不同的频域成分合成回原始的时间序列。

在 C 语言中，我们可以使用 Fast Fourier Transform（FFT）算法来进行离散傅立叶变换和逆变换。FFT 是一种高效的算法，能够在 $O(n\log n)$ 的时间复杂度内完成计算。我们需要用到的库是 math.h 和 complex.h。

离散傅立叶变换

下面是一个示例程序，用于对离散时间序列进行 DFT：

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

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327
#endif

void DFT(int N, double complex x[N], double complex X[N]) {
    for (int k = 0; k < N; k++) {
        X[k] = 0;
        for (int n = 0; n < N; n++) {
            double theta = 2 * M_PI * k * n / N;
            X[k] += x[n] * cexp(-I * theta);
        }
    }
}

在这个程序中，我们定义了一个 DFT 函数，参数为 N、输入序列 x 和输出序列 X。在函数内部，我们首先循环遍历输出序列中的每一个点，然后对输入序列进行求和。最后，我们使用复数指数函数 $e^{-i \theta}$ 对求和结果进行调整，其中 $\theta$ 是傅立叶变换中的旋转因子。

接下来是一个示例程序，用于测试上述的离散傅立叶变换函数：

int main() {
    int N = 8;
    double complex x[] = {1, 2, 3, 4, 4, 3, 2, 1};
    double complex X[N];

    // Perform DFT
    DFT(N, x, X);

    // Print the result
    for (int k = 0; k < N; k++) {
        printf("X[%d] = %lf + %lfi\n", k, creal(X[k]), cimag(X[k]));
    }

    return 0;
}

在上述程序中，我们使用了一个长度为 8 的测试序列，将其作为输入传递给 DFT 函数，并将得到的傅立叶变换结果输出到控制台上。在输出结果中，我们得到了每个频率成分的幅值和相位。

离散傅立叶逆变换

下面是一个示例程序，用于对离散频率序列进行 IDFT：

void IDFT(int N, double complex X[N], double complex x[N]) {
    for (int n = 0; n < N; n++) {
        x[n] = 0;
        for (int k = 0; k < N; k++) {
            double theta = 2 * M_PI * k * n / N;
            x[n] += X[k] * cexp(I * theta);
        }
        x[n] /= N;
    }
}

在这个程序中，我们定义了一个 IDFT 函数，参数为 N、输入频率序列 X 和输出序列 x。在函数内部，我们首先循环遍历输出序列中的每一个点，然后对输入频率序列进行求和。最后，我们使用复数指数函数 $e^{i \theta}$ 对求和结果进行调整，并将结果除以 $N$，得到逆变换。

接下来是一个示例程序，用于测试上述的离散傅立叶逆变换函数：

int main() {
    int N = 8;
    double complex X[] = {4, -2+2*I, 0, -2-I, 0, -2+I, 0, -2-I};
    double complex x[N];

    // Perform IDFT
    IDFT(N, X, x);

    // Print the result
    for (int n = 0; n < N; n++) {
        printf("x[%d] = %lf + %lfi\n", n, creal(x[n]), cimag(x[n]));
    }

    return 0;
}

在上述程序中，我们使用了一个长度为 8 的测试频率序列，将其作为输入传递给 IDFT 函数，并将得到的逆变换结果输出到控制台上。在输出结果中，我们得到了原始信号的值。

总结

在本文中，我们介绍了如何在 C 语言中使用 FFT 算法进行离散傅立叶变换和逆变换。我们编写了两个示例程序，用于对离散时间序列和离散频率序列进行处理，并在控制台上输出结果。这些示例程序可以供程序员参考，帮助他们更好地理解离散傅立叶变换的原理和实现。