CODEWITHSUNDEEP
c++ffteigenfftw
halp_me
halp_me

Reputation: 46

2D fourier transform with Eigen and FFTW

I'm trying to do a real-valued 2d Fourier Transform with FFTW. My data is stored in a dynamically sized Eigen Matrix. Here's the wrapper class I wrote:

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
}

And a definition for EMatrix:

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

The problem is, is I'm getting a seg fault in FFT2D::execute(). I know I'm setting something up wrong in the constructor, and I've tried a number of different ways, but I can't seem to make any headway on this.

Things I've tried include: changing EMatrix typedef to Eigen::ColMajor, passing (fftw_complex *) input.data() to fftw_plan_dft_2d, using different fftw plans (fftw_plan_dft_r2c_2d).

My C++ is (clearly) rusty, but at the end of the day what I need is to do a 2D FT on a real-valued 2D Eigen Matrix of doubles. Thanks in advance.

Upvotes: 2

Views: 2585

Answers (2)

Steve
Steve

Reputation: 217

If your input is real, you can use:

fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1, double *in, fftw_complex *out, unsigned flags);

(which avoids you converting the reals to complexes).

Upvotes: 0

Ap31
Ap31

Reputation: 3324

The major problem here is that there is no such thing as "real-valued Fourier transform". It's just a Fourier transform of something with zero imaginary part, but the zeroes still have to be there, as you can see from fftw_complex definition:

typedef double fftw_complex[2];

This makes sense as the output can (and probably will) have non-zero imaginary part.
The output will have some symmetric properties though, i.e. in case of 1D transform it would be an even function.

As a result the (fftw_complex *) &input(0) cast doesn't really work - FFTW expects twice as many double values as you pass to it.

The solution is to interleave your matrix raw data with zeroes, and there is a number of ways to do that. Few examples:

  • You could copy the whole matrix into a new array before passing it to FFTW, adding the zeroes in process.
  • You could reserve the space for zeroes in the matrix itself - this way you'll be able to avoid copying, but it will probably require a lot of refactoring:)
  • The best way I can think of is to use std::complex<double> as a scalar. This will somewhat hurt your notice of "real-valued FFT", but again there is hardly such thing in the first place. Instead you'll be able to keep all your real-value operations as they are, and the layout of std::complex will fit fftw_complex perfectly.

There could be some other things to consider here, like storage order (FFTW operates on arrays in row-major order, so Eigen matrices should comply) and the validity of linear access to Eigen matrix data (seems OK to me).

Upvotes: 0

Related Questions