2D fourier transform with Eigen and FFTW

Question

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.

score 0 · Answer 1 · answered Mar 23 '17 at 07:42

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).

Thanks for the informative post. I guess my language was sloppy, but what I meant was the input to the FFT was "real-valued" ie, has all zeros for the imaginary part like you said. My followup question is in regards to your third bullet point: What do you mean by "use `std::complex` as a scalar"? Can I just drop this into the `EMatrix` typedef I'm using (as above)? — halp_me, Mar 23 '17 at 15:52
@halp_me yeah, you can just replace `double` with `std::complex` in the typedef. However if you do have interleaving zeroes already, it's perhaps simplier to just keep it this way. The reason I thought you don't have imaginary part is that you pass `input.rows()` and `input.cols()` to FFTW - I am almost sure you want to divide one of them by 2 for this to work — Ap31, Mar 23 '17 at 16:33

2D fourier transform with Eigen and FFTW

1 Answers1