基于FFTW实现二维复数数组行维度4096点FFT及与MATLAB对比

首先明确：你需要对每行独立执行1D的4096点FFT，直接用FFTW的1D FFT例程循环处理每行就好，完全不需要2D FFT（2D是针对整个二维数组的二维变换，和你的需求不匹配）。

下面是具体实现步骤、代码，以及和MATLAB对比的要点：

核心思路拆解 #

你的输入是1001行×144列的复数数组，但每行只有144个点，要做4096点FFT必须先补零到4096个元素。FFTW的1D规划可以复用，循环处理每行即可——这样既高效又节省内存（不需要一次性分配1001×4096的大数组）。

完整代码示例 #

#include <fftw3.h>
#include <complex.h>
#include <stdio.h>
#include <string.h>

#define ROWS 1001
#define COLS 144
#define FFT_SIZE 4096

int main() {
    // 初始化你的输入数组（按定义赋值）
    double complex reshaped_calib[ROWS][COLS];
    for (int i = 0; i < ROWS; i++) {
        for (int j = 0; j < COLS; j++) {
            reshaped_calib[i][j] = 1.0 + 0.1 * I;
        }
    }

    // 分配FFTW单行输入/输出缓冲区（复用处理所有行）
    fftw_complex *input_buf = fftw_malloc(sizeof(fftw_complex) * FFT_SIZE);
    fftw_complex *output_buf = fftw_malloc(sizeof(fftw_complex) * FFT_SIZE);

    // 创建1D正向FFT规划：FFTW_MEASURE会生成最优执行方案（首次运行慢，后续复用效率高）
    fftw_plan fft_plan = fftw_plan_dft_1d(FFT_SIZE, input_buf, output_buf, FFTW_FORWARD, FFTW_MEASURE);

    // 分配结果数组：1001行 × 4096列（如果你的需求确实是4096×144，那应该是对每列做FFT，需调整循环逻辑）
    fftw_complex *fft_result = fftw_malloc(sizeof(fftw_complex) * ROWS * FFT_SIZE);

    // 循环处理每一行
    for (int row = 0; row < ROWS; row++) {
        // 1. 复制当前行到输入缓冲区，剩余位置补零
        memcpy(input_buf, reshaped_calib[row], sizeof(double complex) * COLS);
        for (int col = COLS; col < FFT_SIZE; col++) {
            input_buf[col][0] = 0.0; // 实部补零
            input_buf[col][1] = 0.0; // 虚部补零
        }

        // 2. 执行FFT
        fftw_execute(fft_plan);

        // 3. 将结果存入结果数组的对应行
        memcpy(fft_result + row * FFT_SIZE, output_buf, sizeof(fftw_complex) * FFT_SIZE);
    }

    // --- 可选：将结果写入二进制文件，用于和MATLAB对比 ---
    FILE *out_file = fopen("c_fft_result.bin", "wb");
    fwrite(fft_result, sizeof(double), 2 * ROWS * FFT_SIZE, out_file); // 每个复数占2个double
    fclose(out_file);

    // 清理FFTW资源
    fftw_destroy_plan(fft_plan);
    fftw_free(input_buf);
    fftw_free(output_buf);
    fftw_free(fft_result);

    return 0;
}

关键细节说明 #

1. 为什么不用2D FFT？
   2D FFT是对整个二维数组先做行FFT再做列FFT（或反过来），这不是你要的——你只需要单独对每行做1D变换，循环调用1D例程更直接，也避免了不必要的计算。

2. 内存与类型兼容

   • FFTW的fftw_complex本质是double[2]（实部+虚部），和C标准的double complex二进制完全兼容，所以可以用memcpy直接复制，不用手动拆分实虚部。
   • 用单行缓冲区复用处理所有行，比一次性分配1001×4096的输入数组更节省内存。

3. 和MATLAB结果对比的步骤
   在MATLAB中执行以下代码，生成对比基准并验证C代码的结果：

   % 生成和C代码一致的输入数组
   reshaped_calib = 1.0 + 0.1*1i * ones(1001, 144);
   % 对每行做4096点FFT（和C代码逻辑一致）
   matlab_result = fft(reshaped_calib, 4096, 2);
   
   % 读取C代码输出的二进制文件
   fid = fopen('c_fft_result.bin', 'r');
   c_raw = fread(fid, [FFT_SIZE*2, ROWS], 'double');
   fclose(fid);
   % 转换为复数数组（C是行优先，所以转置后和MATLAB对齐）
   c_result = c_raw(1:FFT_SIZE,:) + 1i*c_raw(FFT_SIZE+1:end,:);
   c_result = c_result';
   
   % 计算最大误差（应该在浮点精度范围内，~1e-15级别）
   max_error = max(max(abs(c_result - matlab_result)))
   

4. 关于输出格式的疑问
   你提到最终输出是4096×144，但按“对每行做4096点FFT”的逻辑，输出应该是1001×4096（1001行FFT结果，每行4096个频率点）。如果确实需要4096×144，那应该是对每列做4096点FFT——只需要把代码中的循环从行改为列，调整数组维度即可。

内容的提问来源于stack exchange，提问作者Makau