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

以下是一个简单的实现二维傅里叶变换的C语言代码，使用了FFT算法。

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

#define PI 3.14159265358979323846

void fft2d(complex double *data, int nx, int ny, int sign)
{
    int i, j, k, l, m, n;
    
    complex double *tmp = malloc(nx * ny * sizeof(complex double));
    
    for (i = 0; i < nx; i++) {
        for (j = 0; j < ny; j++) {
            tmp[i*ny+j] = 0;
            for (k = 0; k < nx; k++) {
                for (l = 0; l < ny; l++) {
                    double angle = 2.0 * PI * (double)(i*k)/nx + 2.0 * PI * (double)(j*l)/ny;
                    complex double c = cos(angle) - sign * I * sin(angle);
                    tmp[i*ny+j] += data[k*ny+l] * c;
                }
            }
        }
    }
    
    for (i = 0; i < nx; i++) {
        for (j = 0; j < ny; j++) {
            data[i*ny+j] = tmp[i*ny+j];
        }
    }
    
    free(tmp);
}

int main()
{
    int nx = 4, ny = 4;
    
    complex double data[] = {
        1, 2, 3, 4,
        5, 6, 7, 8,
        9, 10, 11, 12,
        13, 14, 15, 16
    };
    
    fft2d(data, nx, ny, 1);
    
    for (int i = 0; i < nx; i++) {
        for (int j = 0; j < ny; j++) {
            printf("%f+%fi ", creal(data[i*ny+j]), cimag(data[i*ny+j]));
        }
        printf("\n");
    }
    
    return 0;
}

这个代码假设输入数据是一个大小为`nx` $\times$ `ny`的复数矩阵，输出也是一个相同大小的矩阵。变换的方向由`sign`决定，`sign=1`表示正向变换，`sign=-1`表示逆向变换。在这个代码中，我们使用了暴力的FFT算法，它的时间复杂度是$O(N^4)$，其中$N=$`nx` $\times$ `ny`。对于大型数据，可以使用更高效的算法，例如快速傅里叶变换（FFT）。