仅使用FFTW转换到频域并返回，结果与信号源完全不同

在我的测试用例中，我正在使用FFTW只是将图像转换到频域并返回到空间域。但是我最终得到的结果与我刚开始的结果截然不同。

要解决一些问题：

• 我为每次调用fft()和创建FFTW计划ifft()而不是缓存它们的原因是为了简化。那不应该是问题
• 为了简单起见，在将实值数组传递给fftw之前将其转换为复数数组的原因。我发现这种方法比使用真正复杂的DFT函数更容易出错。

我的代码（包括实用程序功能，宏和类）如下。它使用libcinder（其旧版本）来处理矢量类和图形，但是即使您不了解libcinder，它也将是不言而喻的。请注意代码长度，我已经尽力使它最小化，但是某些实用程序（例如Array2D类）占用了相当大的空间。

main.cpp：

#include "precompiled.h"
typedef std::complex<float> Complexf;
using namespace ci;

void createConsole()
{
    AllocConsole();
    std::fstream* fs = new std::fstream("CONOUT$");
    std::cout.rdbuf(fs->rdbuf());
}

#define forxy(image) \
    for(Vec2i p(0, 0); p.x < image.w; p.x++) \
        for(p.y = 0; p.y < image.h; p.y++)
template<class T>
class ArrayDeleter
{
public:
    ArrayDeleter(T* arrayPtr)
    {
        refcountPtr = new int(0);
        (*refcountPtr)++;

        this->arrayPtr = arrayPtr;
    }

    ArrayDeleter(ArrayDeleter const& other)
    {
        arrayPtr = other.arrayPtr;
        refcountPtr = other.refcountPtr;
        (*refcountPtr)++;
    }

    ArrayDeleter const& operator=(ArrayDeleter const& other)
    {
        reduceRefcount();

        arrayPtr = other.arrayPtr;
        refcountPtr = other.refcountPtr;
        (*refcountPtr)++;

        return *this;
    }

    ~ArrayDeleter()
    {
        reduceRefcount();
    }

private:
    void reduceRefcount()
    {
        (*refcountPtr)--;
        if(*refcountPtr == 0)
        {
            delete refcountPtr;
            fftwf_free(arrayPtr);
        }
    }

    int* refcountPtr;
    T* arrayPtr;
};

template<class T>
struct Array2D
{
    T* data;
    int area;
    int w, h;
    ci::Vec2i Size() const { return ci::Vec2i(w, h); }
    ArrayDeleter<T> deleter;

    Array2D(Vec2i s) : deleter(Init(s.x, s.y)) { }
    Array2D() : deleter(Init(0, 0)) { }

    T* begin() { return data; }
    T* end() { return data+w*h; }

    T& operator()(Vec2i const& v) { return data[v.y*w+v.x]; }

private:
    T* Init(int w, int h) {
        data = (T*)fftwf_malloc(w * h * sizeof(T));
        area = w * h;
        this->w = w;
        this->h = h;
        return data;
    }
};

Array2D<Complexf> fft(Array2D<float> in)
{
    Array2D<Complexf> in_complex(in.Size());
    forxy(in)
    {
        in_complex(p) = Complexf(in(p));
    }

    Array2D<Complexf> result(in.Size());

    auto plan = fftwf_plan_dft_2d(in.h, in.w,
        (fftwf_complex*)in_complex.data, (fftwf_complex*)result.data, FFTW_FORWARD, FFTW_MEASURE);
    fftwf_execute(plan);
    auto mul = 1.0f / sqrt((float)result.area);
    forxy(result)
    {
        result(p) *= mul;
    }
    return result;
}

Array2D<float> ifft(Array2D<Complexf> in)
{
    Array2D<Complexf> result(in.Size());
    auto plan = fftwf_plan_dft_2d(in.h, in.w,
        (fftwf_complex*)in.data, (fftwf_complex*)result.data, FFTW_BACKWARD, FFTW_MEASURE);
    fftwf_execute(plan);

    Array2D<float> out_real(in.Size());
    forxy(in)
    {
        out_real(p) = result(p).real();
    }
    auto mul = 1.0f / sqrt((float)out_real.area);
    forxy(out_real)
    {
        out_real(p) *= mul;
    }
    return out_real;
}

gl::Texture uploadTex(Array2D<float> a)
{
    gl::Texture tex(a.w, a.h);
    tex.bind();
    glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, a.w, a.h, GL_LUMINANCE, GL_FLOAT, a.data);
    return tex;
}

struct SApp : ci::app::AppBasic {
    gl::Texture texSource;
    gl::Texture texbackAndForthed;
    void setup()
    {
        createConsole();

        Array2D<float> source(Vec2i(400, 400));
        forxy(source) {
            float dist = p.distance(source.Size() / 2);
            if(dist < 100)
                source(p) = 1.0f;
            else
                source(p) = 0.0f;
        }
        printSum("source", source);
        texSource = uploadTex(source);
        setWindowSize(source.w, source.h);
        auto backAndForthed = backAndForth(source);
        printSum("backAndForthed", backAndForthed);
        //backAndForthed = backAndForth(loaded);
        texbackAndForthed = uploadTex(backAndForthed);
    }
    void printSum(std::string description, Array2D<float> arr) {
        float sum = std::accumulate(arr.begin(), arr.end(), 0.0f);
        std::cout << "array '" << description << "' has sum " << sum << std::endl;
    }
    void draw()
    {
        gl::clear(Color(0, 0, 0));
        gl::draw(texSource, Rectf(0,0, getWindowWidth() / 2, getWindowHeight() /2));
        gl::draw(texbackAndForthed, Rectf(getWindowWidth() / 2, getWindowWidth() / 2, getWindowWidth(), getWindowHeight()));
    }
    Array2D<float> backAndForth(Array2D<float> in) {
        auto inFd = fft(in);
        auto inResult = ifft(inFd);
        return inResult;
    }
};

CINDER_APP_BASIC(SApp, ci::app::RendererGl)

precompiled.h：

#include <complex>
#include <cinder/app/AppBasic.h>
#include <cinder/gl/Texture.h>
#include <cinder/gl/gl.h>
#include <cinder/Vector.h>
#include <fftw3.h>
#include <numeric>

控制台输出：

array 'source' has sum 31397
array 'backAndForthed' has sum -0.110077

图形输出：

如您所见，右下角的圆圈有点暗和渐变。

注意：如果取消注释第二backAndForthed = backAndForth(loaded);行，则结果是正确的（因此，仅在第一次时结果是错误的）。