使用 FFTW 只是转换到频域并返回得到与源完全不同的结果

Question

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

解决一些问题：

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

我的代码（包含实用函数、宏和 classes）如下。它使用 libcinder（它的旧版本）来处理矢量 classes 和图形，但即使您不了解 libcinder，它也是不言自明的。请原谅代码长度，我已经努力让它变得非常小，但是像我的 Array2D class 这样的实用程序占用了相当多的 space.

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);，结果是正确的（所以，结果只是第一次错误）。

Answer 1

问题出在这里：

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 plan = fftwf_plan_dft_2d(in.h, in.w,
    (fftwf_complex*)in.data, (fftwf_complex*)result.data,
    FFTW_BACKWARD, FFTW_MEASURE);
fftwf_execute(plan);

使用标志 FFTW_MEASURE 意味着 FFTW 会尝试许多 DFT 算法以选择最快的一个。问题是输入数组 in.data 在途中被覆盖了。因此，在 fftwf_execute(plan); 之后，result.data 不是 in.data 的 DFT，因为它是在创建计划之前。根据 planner flags 上 FFTW 的文档：

Important: the planner overwrites the input array during planning unless a saved plan (see Wisdom) is available for that problem, so you should initialize your input data after creating the plan. The only exceptions to this are the FFTW_ESTIMATE and FFTW_WISDOM_ONLY flags, as mentioned below.

文档提供了一个解决方案：使用标志FFTW_ESTIMATE确保input/output数组在规划期间不被覆盖

在这种特殊情况下，使用 FFTW_ESTIMATE 而不是 FFTW_MEASURE 预计不会触发计算时间的大幅增加。实际上，由于在计划创建期间计算了很多 DFT，因此使用 FFTW_MEASURE 创建 fftw_plan 将比使用 FFTW_ESTIMATE 创建慢得多。然而，如果要执行许多相同大小的 DFT，则必须使用 FFTW_MEASURE 一次性创建计划并存储。如有需要，可申请new-array execution functions。

我猜你已经知道 r2c 和 c2r 转换，旨在将存储 space 和计算时间减少近 2 倍。

使用 FFTW 只是转换到频域并返回得到与源完全不同的结果

Using FFTW just to convert to frequency domain and back gives result quite different from the source

c++

signal-processing

fft

fftw