Transformée de Fourier 2D avec Eigen et FFTW

Question

J'essaye de faire une Transformée de Fourier 2d à valeur réelle avec FFTW. Mes données sont stockées dans une matrice propre de taille dynamique. Voici la classe wrapper je l'ai écrit:

FFT2D.h:

#include <Eigen> 

class FFT2D { 
public: 

    enum FFT_TYPE {FORWARD=0, REVERSE=1}; 
    FFT2D(EMatrix &input, EMatrix &output, FFT_TYPE type_ = FORWARD); 
    ~FFT2D(); 

    void execute(); 

private: 
    EMatrix& input; 
    EMatrix& output; 
    fftw_plan plan; 
    FFT_TYPE type; 
}; 

FFT2D.cpp:

#include "FFT2D.h" 
#include <fftw3.h> 
#include "Defs.h" 


FFT2D::FFT2D(EMatrix &input_, EMatrix &output_, FFT_TYPE type_) 
     : type(type_), input(input_), output(output_) { 

    if (type == FORWARD) 
     plan = fftw_plan_dft_2d((int) input.rows(), (int) input.cols(), 
       (fftw_complex *) &input(0), (fftw_complex *) &output(0), 
       FFTW_FORWARD, FFTW_ESTIMATE); 
    else 
     // placeholder for ifft-2d code, unwritten 
} 


FFT2D::~FFT2D() { 
    fftw_destroy_plan(plan); 
} 

void FFT2D::execute() { 
    fftw_execute(plan); // seg-fault here 
} 

Et une définition pour EMatrix:

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> EMatrix; 

Le problème est, est-ce que je reçois une erreur seg dans FFT2D::execute(). Je sais que je mets quelque chose en erreur dans le constructeur, et j'ai essayé plusieurs façons différentes, mais je n'arrive pas à progresser dans ce domaine.

choses que j'ai essayé comprennent: le changement EMatrix typedef à Eigen::ColMajor, en passant (fftw_complex *) input.data()-fftw_plan_dft_2d, en utilisant différents plans FFTW (fftw_plan_dft_r2c_2d). Mon C++ est (clairement) rouillé, mais à la fin de la journée ce dont j'ai besoin est de faire un FT 2D sur une matrice propre de valeurs réelles 2D de doubles. Merci d'avance.

Answer 1

Le problème majeur ici est qu'il n'existe pas de "transformée de Fourier à valeur réelle". Il est juste une transformée de Fourier de quelque chose avec une partie imaginaire zéro, mais les zéros doivent encore être là, comme vous pouvez le voir fftw_complex Définition:

typedef double fftw_complex[2];

Cela est logique que la sortie peut (et probablement) avoir une partie imaginaire non nulle.
_{La sortie aura cependant des propriétés symétriques, c'est-à-dire dans le cas d'une transformation 1D, il s'agira d'une fonction paire. Par conséquent, le cast (fftw_complex *) &input(0) ne fonctionne pas vraiment - FFTW attend deux fois plus de valeurs double que vous lui passez.}

La solution consiste à imbriquer les données brutes de votre matrice avec des zéros, et il existe un certain nombre de façons de le faire. Quelques exemples:

• Vous pouvez copier la matrice entière dans un nouveau tableau avant de la transmettre à FFTW, en ajoutant les zéros en cours.
• Vous pouvez réserver l'espace pour les zéros dans la matrice elle-même - de cette façon, vous serez en mesure d'éviter la copie, mais il faudra probablement beaucoup de refactoring :)
• La meilleure façon que je peux penser est d'utiliser std::complex<double> en tant que scalaire. Cela nuira quelque peu à votre avis de «FFT à valeur réelle», mais encore une fois, il n'y a pratiquement rien de tel en premier lieu. Au lieu de cela, vous serez en mesure de garder toutes vos opérations en valeur réelle telles qu'elles sont, et la mise en page de std::complex s'adaptera parfaitement à fftw_complex.

Il pourrait y avoir d'autres choses à considérer ici, comme ordre de stockage (fonctionne FFTW sur les tableaux dans l'ordre des lignes grande, donc les matrices Eigen doivent être conformes) et la validité de l'accès linéaire aux données de la matrice Eigen (semble OK moi).