写一个实现二维离散傅里叶变换的C语言代码

以下为实现二维离散傅里叶变换的C语言代码：

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

#define PI 3.14159265358979323846

/* 复数结构体 */
typedef struct {
    double real;    // 实部
    double imag;    // 虚部
} complex;

/* 二维离散傅里叶变换 */
void dft2d(complex **input, complex **output, int width, int height) {
    int i, j, u, v;
    double cos_val, sin_val;
    complex sum, temp;

    // 对每一个频率进行计算
    for (u = 0; u < height; u++) {
        for (v = 0; v < width; v++) {
            sum.real = 0.0;
            sum.imag = 0.0;
            for (i = 0; i < height; i++) {
                for (j = 0; j < width; j++) {
                    // 计算旋转因子
                    cos_val = cos(2 * PI * (u * i / (double)height + v * j / (double)width));
                    sin_val = -sin(2 * PI * (u * i / (double)height + v * j / (double)width));
                    // 计算加权和
                    temp.real = input[i][j].real * cos_val - input[i][j].imag * sin_val;
                    temp.imag = input[i][j].real * sin_val + input[i][j].imag * cos_val;
                    sum.real += temp.real;
                    sum.imag += temp.imag;
                }
            }
            output[u][v].real = sum.real;
            output[u][v].imag = sum.imag;
        }
    }
}

int main() {
    int width = 4;
    int height = 4;
    int i, j;
    complex **input, **output;

    // 分配内存空间
    input = (complex **)malloc(sizeof(complex *) * height);
    output = (complex **)malloc(sizeof(complex *) * height);
    for (i = 0; i < height; i++) {
        input[i] = (complex *)malloc(sizeof(complex) * width);
        output[i] = (complex *)malloc(sizeof(complex) * width);
    }

    // 初始化输入矩阵
    input[0][0].real = 1; input[0][0].imag = 1;
    input[0][1].real = 2; input[0][1].imag = 2;
    input[0][2].real = 3; input[0][2].imag = 3;
    input[0][3].real = 4; input[0][3].imag = 4;
    input[1][0].real = 5; input[1][0].imag = 5;
    input[1][1].real = 6; input[1][1].imag = 6;
    input[1][2].real = 7; input[1][2].imag = 7;
    input[1][3].real = 8; input[1][3].imag = 8;
    input[2][0].real = 9; input[2][0].imag = 9;
    input[2][1].real = 10; input[2][1].imag = 10;
    input[2][2].real = 11; input[2][2].imag = 11;
    input[2][3].real = 12; input[2][3].imag = 12;
    input[3][0].real = 13; input[3][0].imag = 13;
    input[3][1].real = 14; input[3][1].imag = 14;
    input[3][2].real = 15; input[3][2].imag = 15;
    input[3][3].real = 16; input[3][3].imag = 16;

    // 进行二维离散傅里叶变换
    dft2d(input, output, width, height);

    // 输出结果
    for (i = 0; i < height; i++) {
        for (j = 0; j < width; j++) {
            printf("(%f, %f) ", output[i][j].real, output[i][j].imag);
        }
        printf("\n");
    }

    // 释放内存空间
    for (i = 0; i < height; i++) {
        free(input[i]);
        free(output[i]);
    }
    free(input);
    free(output);

    return 0;
}

该代码中，我们定义了一个复数结构体，用于存储实部和虚部。在主函数中，我们定义了一个 $4 \times 4$ 的输入矩阵，并对其进行二维离散傅里叶变换。输出结果为一个 $4 \times 4$ 的矩阵，其中每个元素为一个复数，表示对应频率的幅度和相位。

在 dft2d 函数中，我们使用了四重循环，对每个频率进行计算。内部两重循环用于遍历输入矩阵，计算加权和；外部两重循环用于遍历频率，计算旋转因子和加权和。旋转因子的计算使用了欧拉公式：$$ e^{ix} = \cos x + i\sin x $$ 因此，我们可以用 cos 和 sin 函数计算出旋转因子的实部和虚部。最后，将加权和存储到输出矩阵中。

需要注意的是，在本代码中，输入矩阵和输出矩阵都是二维复数数组，因此在分配内存空间时需要进行嵌套的 malloc 调用。同时，计算旋转因子时需要将整数除法结果强制转换为 double 类型，否则可能会出现精度误差。