Использование FFTW только для преобразования в частотную область и обратно дает результат, отличный от источника

Question

В моей тестовой папке я использую FFTW только для преобразования изображения в частотную область и обратно в пространственный домен. Но то, что я получаю, - результат, отличный от того, с чего я начал.

Чтобы получить некоторые вопросы из пути:

• Причина Я создаю план FFTW на каждом вызове fft() и 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);, результат будет правильным (так, результат будет неправильным только в первый раз).

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 пытается много ДПФ алгоритмов выбрать самый быстрый один. Проблема в том, что входной массив in.data перезаписан. Следовательно, после fftwf_execute(plan);, result.data не является ДПФ in.data, как это было до создания плана. Согласно документации FFTW на planner flags:

Важно: планировщик перезаписывает массив ввода при не планируют, если сохраненный план (см мудрости) доступен для этой проблемы, так что вы должны инициализировать входные данные после создания плана , Исключением являются только флаги FFTW_ESTIMATE и FFTW_WISDOM_ONLY, как указано ниже.

документация обеспечивает решение: с помощью флага FFTW_ESTIMATE гарантирует, что входные выходные массивы/не перезаписываются во время планирования.

В данном конкретном случае не ожидается, что использование FFTW_ESTIMATE вместо FFTW_MEASURE вызовет значительное увеличение времени вычисления. В самом деле, как много DFT вычисляются во время создания плана, создание fftw_plan с использованием FFTW_MEASURE будет намного медленнее, чем создание с использованием FFTW_ESTIMATE. Тем не менее, если необходимо выполнить много ДПФ того же размера, план должен быть создан один раз для всех с использованием FFTW_MEASURE и сохранен. При необходимости можно применять new-array execution functions.

Я думаю, чем вы уже знаете о r2c и c2r трансформирует, направленные на сокращение дискового пространства и времени вычислений почти 2.