就地二维快速傅里叶变换的问题

Question

我正在使用 fftw 库实现一个 2D fft 示例，以了解如何使用它，并且我遇到了就地 r2c 案例的问题。

据我了解，我正在按照应该完成的方式进行操作：i）分配正确大小的内存（

N_x * 2 * (N_y/2 + 1)

，我从here获得）和ii）创建计划将输入作为 double* 传递，并将其作为输出传递，将其转换为 (fftw_complex*) 类型。然而，结果与我实现的异地 r2c 和就地 c2c 都不匹配。

运行这三种情况并进行比较的代码如下：

#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

int main(){
    // define parameters
    int N_x = 20;
    int N_y = 20;
    int N_y_complex = N_y/2 + 1;

    double a = 2;
    double b = 2;

    double range_x = 10;
    double range_y = 10;

    double dx = range_x/(N_x);
    double dy = range_y/(N_y);

    ////////////////////////////////////////// OUT-OF-PLACE R2C TRANSFORM //////////////////////////////////////////

    // declare functions
    double *in = (double*)fftw_malloc(N_x * N_y * sizeof(double));
    fftw_complex *out = (fftw_complex*)fftw_malloc(N_x * N_y_complex * sizeof(fftw_complex));

    for(int i = 0; i < N_x; i++){
        double x = i*dx;
        for(int j = 0; j < N_y; j++){
            double y = j*dy;
            
            if(fabs(x) <= a/2 && fabs(y) <= b/2){
                in[i*N_y + j] = 1;
            }
            else if(fabs(range_x - x) <= a/2 && fabs(range_y - y) <= b/2){
                in[i*N_y + j] = 1;
            }
            else if(fabs(x) <= a/2 && fabs(range_y - y) <= b/2){
                in[i*N_y + j] = 1;
            }
            else if(fabs(range_x - x) <= a/2 && fabs(y) <= b/2){
                in[i*N_y + j] = 1;
            }
            else{
                in[i*N_y + j] = 0;
            }
        }
    }

    // Check if the function is even
    bool is_even = true;
    for (int i = 0; i < N_x / 2; i++) {
        for (int j = 0; j < N_y; j++) {
            int i_neg = (N_x - i) % N_x;
            int j_neg = (N_y - j) % N_y;

            if (in[i * N_y + j] != in[i_neg * N_y + j_neg]) {
                cout << "Symmetry error: f(" << i * dx << "," << j * dy << ") = " << in[i * N_y + j] 
                     << " but f(" << i_neg * dx << "," << j_neg * dy << ") = " << in[i_neg * N_y + j_neg] << endl;
                is_even = false;
                cout << i << ", " << j << " || " << i_neg << ", " << j_neg  <<  endl;
            }
        }
    }

    cout << "Checking f(x,y) parity: " << endl;
    if (is_even) {
        cout << "   The function is even." << endl;
    } else {
        cout << "   The function is not even." << endl;
    }

    fftw_plan p = fftw_plan_dft_r2c_2d(N_x,N_y,in,out,FFTW_ESTIMATE);
    fftw_execute(p);

    ////////////////////////////////////////// IN-PLACE C2C FORWARD TRANSFORM //////////////////////////////////////////

    fftw_complex *in2 = (fftw_complex*)fftw_malloc(N_x * N_y * sizeof(fftw_complex));

    for(int i = 0; i < N_x; i++){
        for(int j = 0; j < N_y; j++){
            in2[i*N_y + j][0] = in[i*N_y + j];
            in2[i*N_y + j][1] = 0;
        }
    }

    fftw_plan p_inPlace_c2c = fftw_plan_dft_2d(N_x, N_y, in2, in2, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(p_inPlace_c2c);

    bool success = false;
    double tol = 1e-10;

    for(int i = 0; i < N_x; i++){
        for(int j = 0; j < N_y/2; j++){

            if(abs(out[i * N_y_complex + j][0] - in2[i * N_y_complex + j][0]) > tol && abs(out[i * N_y_complex + j][1] - in2[i * N_y_complex + j][1]) > tol){
                success = false;
                goto skiploop;
            }
            else{
                success = true;
            }
        }
    }

skiploop:

    cout << "Checking if in-place (forward C2C) and out-of-place results are matching: " << endl;
    if(success == true){
        cout << "   Matching results for in-place c2c and out-of-place" << endl;
    }
    else{
        cout << "   Something went wrong" << endl;
    }

    ////////////////////////////////////////// IN-PLACE R2C TRANSFORM //////////////////////////////////////////

    // Calculate padded size: N_x * 2 * (N_y/2 + 1) = N_x * 2 * N_y_complex
    int padded_size = N_x * 2 * N_y_complex;

    // Allocate enough memory to hold complex data
    double *in3 = (double*)fftw_malloc(padded_size * sizeof(double));

    // Assign values for the real function to be transformed
    for(int i = 0; i < N_x; i++){
        for(int j = 0; j < N_y; j++){
            in3[i*N_y + j] = in[i*N_y + j];
        }
    }

    // Create plan casting in3 to fftw_complex type
    fftw_plan p_inPlace_r2c = fftw_plan_dft_r2c_2d(N_x, N_y, in3, (fftw_complex*)in3, FFTW_ESTIMATE);

    // Execute plan
    fftw_execute(p_inPlace_r2c);

    // printing first real value for debugging:
    cout << "out-of-place r2c: " << out[0][0] << " || in-place c2c: " << in2[0][0] << " || in-place r2c: " << ((fftw_complex*)in3)[0][0] << endl;

    // Check final results compared to out-of-place transform
    for(int i = 0; i < N_x; i++){
        for(int j = 0; j < N_y/2; j++){

            if(abs(out[i * N_y_complex + j][0] - ((fftw_complex*)in3)[i * N_y_complex + j][0]) > tol && abs(out[i * N_y_complex + j][1] - ((fftw_complex*)in3)[i * N_y_complex + j][1]) > tol){
                success = false;
                goto skiploop2;
            }
            else{
                success = true;
            }
        }
    }

skiploop2:

    cout << "Checking if in-place (R2C) and out-of-place results are matching: " << endl;
    if(success == true){
        cout << "   Matching results for in-place r2c and out-of-place" << endl;
    }
    else{
        cout << "   Something went wrong" << endl;
    }

    fftw_free(in);
    fftw_free(out);
    fftw_free(in2);
    fftw_free(in3);
    fftw_destroy_plan(p);
    fftw_destroy_plan(p_inPlace_c2c);
    fftw_destroy_plan(p_inPlace_r2c);

}

我为获取所有三种情况显示的第一个实际值而输入的打印语句：

out-of-place r2c: 25 || in-place c2c: 25 || in-place r2c: -2.58976e+284

在这种情况下我错过了什么？我将其编译为

g++ -I/path/to/include -o fft2d fft2d.cpp -L/path/to/lib -lfftw3

，因为我正在使用的机器在自定义目录中有 fftw 库。

Answer 1

感谢@TimRoberts 为我澄清了这一点！我原以为额外的填充仅在用新的复数值覆盖数组时才重要，所以我只是在

in3

的末尾留下了所有额外的空间。

因此，在将 2D 数组表示为 1D 行主数组时，我需要修改这部分代码，以确保所有额外的填充都放置在与每行末尾相对应的正确位置：

    // Assign values for the real function to be transformed
    int extra_row_padding;
    for(int i = 0; i < N_x; i++){
        for(int j = 0; j < N_y; j++){
            if(N_y % 2 == 0){
                extra_row_padding = 2;
            }
            else{
                extra_row_padding = 1;
            }
            in3[i*(N_y + extra_row_padding) + j] = in[i * N_y + j];
        }
    }

这是完整的工作示例：

#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

int main(){
    // define parameters
    int N_x = 20;
    int N_y = 20;
    int N_y_complex = N_y/2 + 1;

    double a = 2;
    double b = 2;

    double range_x = 10;
    double range_y = 10;

    double dx = range_x/(N_x);
    double dy = range_y/(N_y);

    ////////////////////////////////////////// OUT-OF-PLACE R2C TRANSFORM //////////////////////////////////////////

    // declare functions
    double *in = (double*)fftw_malloc(N_x * N_y * sizeof(double));
    fftw_complex *out = (fftw_complex*)fftw_malloc(N_x * N_y_complex * sizeof(fftw_complex));

    for(int i = 0; i < N_x; i++){
        double x = i*dx;
        for(int j = 0; j < N_y; j++){
            double y = j*dy;
            
            if(fabs(x) <= a/2 && fabs(y) <= b/2){
                in[i*N_y + j] = 1;
            }
            else if(fabs(range_x - x) <= a/2 && fabs(range_y - y) <= b/2){
                in[i*N_y + j] = 1;
            }
            else if(fabs(x) <= a/2 && fabs(range_y - y) <= b/2){
                in[i*N_y + j] = 1;
            }
            else if(fabs(range_x - x) <= a/2 && fabs(y) <= b/2){
                in[i*N_y + j] = 1;
            }
            else{
                in[i*N_y + j] = 0;
            }
        }
    }

    // Check if the function is even
    bool is_even = true;
    for (int i = 0; i < N_x / 2; i++) {
        for (int j = 0; j < N_y; j++) {
            int i_neg = (N_x - i) % N_x;
            int j_neg = (N_y - j) % N_y;

            if (in[i * N_y + j] != in[i_neg * N_y + j_neg]) {
                cout << "Symmetry error: f(" << i * dx << "," << j * dy << ") = " << in[i * N_y + j] 
                     << " but f(" << i_neg * dx << "," << j_neg * dy << ") = " << in[i_neg * N_y + j_neg] << endl;
                is_even = false;
                cout << i << ", " << j << " || " << i_neg << ", " << j_neg  <<  endl;
            }
        }
    }

    cout << "Checking f(x,y) parity: " << endl;
    if (is_even) {
        cout << "   The function is even." << endl;
    } else {
        cout << "   The function is not even." << endl;
    }

    fftw_plan p = fftw_plan_dft_r2c_2d(N_x,N_y,in,out,FFTW_ESTIMATE);
    fftw_execute(p);

    ////////////////////////////////////////// IN-PLACE C2C FORWARD TRANSFORM //////////////////////////////////////////

    fftw_complex *in2 = (fftw_complex*)fftw_malloc(N_x * N_y * sizeof(fftw_complex));

    for(int i = 0; i < N_x; i++){
        for(int j = 0; j < N_y; j++){
            in2[i*N_y + j][0] = in[i*N_y + j];
            in2[i*N_y + j][1] = 0;
        }
    }

    fftw_plan p_inPlace_c2c = fftw_plan_dft_2d(N_x, N_y, in2, in2, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(p_inPlace_c2c);

    bool success = false;
    double tol = 1e-10;

    for(int i = 0; i < N_x; i++){
        for(int j = 0; j < N_y/2; j++){

            if(abs(out[i * N_y_complex + j][0] - in2[i * N_y_complex + j][0]) > tol && abs(out[i * N_y_complex + j][1] - in2[i * N_y_complex + j][1]) > tol){
                success = false;
                goto skiploop;
            }
            else{
                success = true;
            }
        }
    }

skiploop:

    cout << "Checking if in-place (forward C2C) and out-of-place results are matching: " << endl;
    if(success == true){
        cout << "   Matching results for in-place c2c and out-of-place" << endl;
    }
    else{
        cout << "   Something went wrong" << endl;
    }

    ////////////////////////////////////////// IN-PLACE R2C TRANSFORM //////////////////////////////////////////

    // Calculate padded size: N_x * 2 * (N_y/2 + 1) = N_x * 2 * N_y_complex
    int padded_size = N_x * 2 * N_y_complex;

    // Allocate enough memory to hold complex data
    double *in3 = (double*)fftw_malloc(padded_size * sizeof(double));

    // Assign values for the real function to be transformed
    int extra_row_padding;
    for(int i = 0; i < N_x; i++){
        for(int j = 0; j < N_y; j++){
            if(N_y % 2 == 0){
                extra_row_padding = 2;
            }
            else{
                extra_row_padding = 1;
            }
            in3[i*(N_y + extra_row_padding) + j] = in[i * N_y + j];
        }
    }

    // Create plan casting in3 to fftw_complex type
    fftw_plan p_inPlace_r2c = fftw_plan_dft_r2c_2d(N_x, N_y, in3, (fftw_complex*)in3, FFTW_ESTIMATE);

    // Execute plan
    fftw_execute(p_inPlace_r2c);

    // printing first real value for debugging:
    cout << "out-of-place r2c: " << out[0][0] << " || in-place c2c: " << in2[0][0] << " || in-place r2c: " << ((fftw_complex*)in3)[0][0] << endl;

    // Check final results compared to out-of-place transform
    for(int i = 0; i < N_x; i++){
        for(int j = 0; j < N_y/2; j++){

            if(abs(out[i * N_y_complex + j][0] - ((fftw_complex*)in3)[i * N_y_complex + j][0]) > tol && abs(out[i * N_y_complex + j][1] - ((fftw_complex*)in3)[i * N_y_complex + j][1]) > tol){
                success = false;
                goto skiploop2;
            }
            else{
                success = true;
            }
        }
    }

skiploop2:

    cout << "Checking if in-place (R2C) and out-of-place results are matching: " << endl;
    if(success == true){
        cout << "   Matching results for in-place r2c and out-of-place" << endl;
    }
    else{
        cout << "   Something went wrong" << endl;
    }

    fftw_free(in);
    fftw_free(out);
    fftw_free(in2);
    fftw_free(in3);
    fftw_destroy_plan(p);
    fftw_destroy_plan(p_inPlace_c2c);
    fftw_destroy_plan(p_inPlace_r2c);

}