地震数据处理二维快速傅里叶变换的c语言代码

以下是一份基于C语言的地震数据处理二维快速傅里叶变换代码，供参考：

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

#define PI 3.14159265358979323846

int main()
{
    int Nx = 256; // X维数据点数
    int Ny = 256; // Y维数据点数
    
    // 生成数据
    double **data = (double **)malloc(Nx * sizeof(double *));
    for (int i = 0; i < Nx; i++) {
        data[i] = (double *)malloc(Ny * sizeof(double));
        for (int j = 0; j < Ny; j++) {
            data[i][j] = sin(i * j * PI / 128.0);
        }
    }
    
    // 二维快速傅里叶变换
    double complex **fft2 = (double complex **)malloc(Nx * sizeof(double complex *));
    for (int i = 0; i < Nx; i++) {
        fft2[i] = (double complex *)malloc(Ny * sizeof(double complex));
        for (int j = 0; j < Ny; j++) {
            fft2[i][j] = 0.0;
        }
    }
    double complex *fft1 = (double complex *)malloc(Ny * sizeof(double complex));
    for (int i = 0; i < Nx; i++) {
        for (int j = 0; j < Ny; j++) {
            fft1[j] = data[i][j];
        }
        fft(fft1, Ny);
        for (int j = 0; j < Ny; j++) {
            fft2[i][j] = fft1[j];
        }
    }
    for (int j = 0; j < Ny; j++) {
        double complex *fft1 = (double complex *)malloc(Nx * sizeof(double complex));
        for (int i = 0; i < Nx; i++) {
            fft1[i] = fft2[i][j];
        }
        fft(fft1, Nx);
        for (int i = 0; i < Nx; i++) {
            fft2[i][j] = fft1[i];
        }
    }
    
    // 输出结果
    for (int i = 0; i < Nx; i++) {
        for (int j = 0; j < Ny; j++) {
            printf("%f + %fi ", creal(fft2[i][j]), cimag(fft2[i][j]));
        }
        printf("\n");
    }
    
    // 释放内存
    for (int i = 0; i < Nx; i++) {
        free(data[i]);
        free(fft2[i]);
    }
    free(data);
    free(fft2);
    free(fft1);
    
    return 0;
}

void fft(double complex *x, int N)
{
    if (N <= 1) {
        return;
    }
    double complex *even = (double complex *)malloc(N / 2 * sizeof(double complex));
    double complex *odd = (double complex *)malloc(N / 2 * sizeof(double complex));
    for (int i = 0; i < N / 2; i++) {
        even[i] = x[2 * i];
        odd[i] = x[2 * i + 1];
    }
    fft(even, N / 2);
    fft(odd, N / 2);
    for (int k = 0; k < N / 2; k++) {
        double complex t = cexp(-2 * PI * I * k / N) * odd[k];
        x[k] = even[k] + t;
        x[k + N / 2] = even[k] - t;
    }
    free(even);
    free(odd);
}

需要注意的是，这份代码没有实现数据的读取和输出，需要自行添加。另外，这份代码中的快速傅里叶变换实现是递归形式的Cooley-Tukey算法，可以对任意长度的一维数据进行快速傅里叶变换。