The results of an Inverse Fourier Transform (IFT) using FFTW3 doesn't match the analytic form, or other IFT Methods, why?

Question

I am evaluating a few different methods of performing an inverse Fourier transform (IFT), numerically on 1D data that is both real and even. The function in question is f(x) = sech(A * x), the hyperbolic secant function. The analytic form of the IFT is f(x) = (pi/A) * sech(pi^2 x / A ), a scaled sech function.

I'm using CImg, as this is related to image processing task, but first I was teaching myself (numerically) about forward and inverse transforms. I've hit a snag that I can't seem to find an answer to, which is this -

I can't get FFTW3's output to match either the analytic form, the result of an inverse fourier transform from CImg, or even a hand rolled simple IFT snippet of code. For example, the following table shows the comparison of values that are the result of the analytic IFT, CImg's IFT (Which uses fftw3), FFTW3 directly, and my hand-written IFT. FFTW3 is close, but not quite right. Is it my fftw plan? Could it be my data, which is a 1D array of an odd number of pixels? Is it something else? Just trying to figure it out.

Because I know in advance that my data is real, and even, I'm trying to use the Real-to-Real interface in FFTW. The kind of FFTW that I'm attempting to use is the FFTW_REDFT00. The full comparison between the four methods is shown below, with the IFT's being log scaled for clarity (left axis) and the original function on the right axis.

The full listing is below, with sample output. It will produce the full table in the terminal output, separated by tabs. The FFTW version I have is 3.3.10. Using Debian 12, with g++ 12.2. Compile it with the following command

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;
}