使用 FFTW3 的逆傅里叶变换 （IFT） 的结果与解析形式或其他 IFT 方法不匹配，为什么？

问：

我正在评估几种不同的方法，以数字方式对真实和偶数的一维数据进行逆傅里叶变换 （IFT）。所讨论的函数是 ，双曲割线函数。IFT 的解析形式是 ，一个缩放的 sech 函数。f(x) = sech(A * x)f(x) = (pi/A) * sech(pi^2 x / A )

我正在使用 CImg，因为这与图像处理任务有关，但首先我（以数字方式）自学了正向和反向变换。我遇到了一个我似乎找不到答案的障碍，那就是这个——

我无法让 FFTW3 的输出与解析形式、CImg 的逆傅里叶变换结果相匹配，甚至无法与手卷的简单 IFT 代码片段相匹配。例如，下表显示了分析 IFT、CImg 的 IFT（使用 fftw3）、直接使用 FFTW3 和我手写的 IFT 的结果值的比较。FFTW3 很接近，但不太正确。这是我的fftw计划吗？会不会是我的数据，这是一个奇数像素的一维数组？是别的吗？只是想弄清楚。

因为我事先知道我的数据是真实的，甚至，我正在尝试在 FFTW 中使用实对实界面。我尝试使用的 FFTW 类型是 FFTW_REDFT00。四种方法之间的完整比较如下所示，为清晰起见，IFT 的对数缩放（左轴）和原始函数在右轴上。

完整列表如下，并带有示例输出。它将在终端输出中生成完整的表，由制表符分隔。我拥有的 FFTW 版本是 3.3.10。使用 Debian 12，g++ 12.2。使用以下命令编译它

g++ -std=c++17 main_so.cpp -o main2 -lX11 -pthread -lm -lfftw3_threads -lfftw3 -Wall -Wextra -Wconversion -pedantic-errors

#define cimg_use_fftw3

#include "CImg.h"
#include "fftw3.h"

#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <vector>

// alias for a image with double values
using img_t = cimg_library::CImg<double>;

// injecting this into namespace std for this example.
namespace std { double sech(double arg); }

// The different methods used to create the inverse fourier transforms
img_t create_sech_function(int N);
img_t create_analytic_IFT(int N);
img_t create_CImg_IFT(int N, const img_t& img);
img_t create_FFTW3_IFT(int N, const img_t& img);
img_t create_simple_numerical_IFT(int N, const img_t& img);

int main (int argc, char** argv)
{      
    // Get the data array dimension
    const int N = std::atoi( argv[1] );     
    const int NO2{ N/2 };

    // compute the various inverse fourier transforms
    img_t sech = create_sech_function(N);
    img_t analytic_IFT = create_analytic_IFT(N);
    img_t CImg_IFT = create_CImg_IFT(N, sech);
    img_t FFTW3_IFT = create_FFTW3_IFT(N, sech);
    img_t simple_IFT = create_simple_numerical_IFT(N, sech);

    // if yout want to display the images for ease of analysis
    bool display_images{ false };

    if(display_images) {
        sech.display("sech(a*x)");
        analytic_IFT.display("Analytic IFT");
        CImg_IFT.display("CImg IFT");
        FFTW3_IFT.display("FFTW3 IFT");
        simple_IFT.display("Simple 1D r to c IFT");
    }

    // print the data to the screen
    std::cout<< "Index\tsech(A*x)\tanalytic\tCImg\t\tSimple\t\tFFTW3\n";
    std::cout<< std::scientific << std::setprecision(7);
    for(int i = 0; i < N; ++i) {
        std::cout<< -NO2 + i << '\t'
                 << sech(i) << '\t'
                 << analytic_IFT(i) << '\t'
                 << CImg_IFT(i) << '\t'
                 << simple_IFT(i) << '\t'
                 << FFTW3_IFT(i) << '\n';
    }

    return EXIT_SUCCESS;
}

// cmath doesn't have a sech(x) function, so we'll define 
// our own.
namespace std {
    double sech(double arg) {
        return 1.0 / std::cosh(arg) ;
    }
}

img_t create_sech_function(int N) {

    img_t sech(N, 1, 1, 1, double{});  

    const double pi{ 4.0 * std::atan(1.0) };  
    const double freq{ 8.0 * pi / N };   
    const int NO2{ N/2 };
    
    // scale the sech function so that it 
    // always fits well within the bounds of our image
    for (int i = -NO2; i <= NO2; ++i) {
        sech(i+NO2) = std::sech(freq * i);      
    }

    return sech;
}

img_t create_analytic_IFT(int N) {

    img_t analytic(N, 1, 1, 1, double{});   

    const double pi{ 4.0 * std::atan(1.0) };  
    const double freq{ 8.0 * pi / N };   
    const double period{ 1.0 / freq};  
    const double IN{ 1.0/N };
    const double C0{ period * pi * IN };
    const double C1{ C0*pi };    
    const int NO2{ N/2 };

    // the analytic function is another sech function.
    // very similar to a Gaussian like function in that sense.
    for (int i = -NO2; i <= NO2; ++i) {
        analytic(i+NO2) = C0 * std::sech(C1 * i);
    }

    return analytic;
}

img_t create_CImg_IFT(int N, const img_t& img) {

    // Get the inverse transform through cimg.  
    // NOTE: because we've defined cimg_use_fftw3, cimg actually 
    // uses the real to complex interface, and creates essentially 
    // the same plan as the simple numerical version, which is 
    // why they match 1 to 1 in the data output.
    cimg_library::CImgList<double> CImg_IFT = img.get_FFT('x',true);
      
    img_t IFT_magnitude(N,1,1,1, double{});   
    for(int i = 0; i < N; ++i){        
         IFT_magnitude(i) = std::hypot( CImg_IFT(0)(i), CImg_IFT(1)(i) );
    }
   
    const int NO2{ N/2 };
    IFT_magnitude.shift(NO2,0,0,0,2);

    return IFT_magnitude;
}

img_t create_FFTW3_IFT(int N, const img_t& img) {

    fftw_plan P1;

    double* in = fftw_alloc_real(N);
    double* out = fftw_alloc_real(N);

    // using the real to real interface, just as before.
    P1 = fftw_plan_r2r_1d(N, in, out, FFTW_REDFT00, FFTW_MEASURE);   
   
    // populate the in array
    std::copy(img.begin(), img.end(), in);

    fftw_execute(P1);

    // Create a cimg container from the data, normalizing by N' = 2(N-1)
    img_t fftw3_IFT(out, N, 1, 1, 1);    
    fftw3_IFT /= 2.0*(N - 1);

    // deallocate the memory 
    fftw_destroy_plan(P1);
    fftw_free(in);
    fftw_free(out);
    fftw_cleanup();

    // create an empty container with the same dimensions.
    img_t fftw3_IFT_shift(fftw3_IFT, "xyzc", 0.0);

    // Get only the even frequencies.
    for(int i = 0; i < N; i += 2) {
        fftw3_IFT_shift(i/2) = fftw3_IFT(i);
    }

    // Shift the edges around to the center.
    const int NO2{ N/2 };
    fftw3_IFT_shift.shift(NO2,0,0,0,2);

    // reflect about the central value    
    const int NM1{ N-1 };
    for(int i = 0; i < NO2; ++i) {
        fftw3_IFT_shift(i) = fftw3_IFT_shift(NM1-i);
    }  

    fftw3_IFT_shift.abs();
    return fftw3_IFT_shift;
}

img_t create_simple_numerical_IFT(int N, const img_t& img) {

    // Create the data input and output arrays.  
    std::vector<double> data_in{img.begin(), img.end()};
    std::vector<std::complex<double>> data_out(data_in.size());

    const double IN{ 1.0/N };
    const int NO2{ N/2 };
    const double pi{ 4.0 * std::atan(1.0) }; 

    // Create our inverse Fourier transform function     
    auto F = [&](double k) {                 
         std::complex<double> val{};
         for(int i = 0; i < N; ++i) {
                 val += data_in.at(i) * std::exp(std::complex<double>{0.0,2*pi*k*i*IN});
         }
         return val;     
    };   

    // calculate the IFT here
    for(int i = 0; i < N; ++i) { 
        data_out.at(i) = F(i); 
    }

    // make a vector for the IFT magnitudes, and calculate it below.
    std::vector<double> data_mag(data_in.size());
    
    std::transform(data_out.cbegin(), data_out.cend(), data_mag.begin(),
         [&](const std::complex<double>& C) { return std::abs(C)*IN; });
        
    // Shift the data back to the center
    std::rotate(data_mag.begin(), data_mag.begin() + NO2 + 1, data_mag.end()); 
       
    // put the data into a cimg object
    img_t my_IFT(data_mag.data(),N,1,1,1);      

    return my_IFT;
}