CODEWITHSUNDEEP
c++fftw
Roman Rdgz
Roman Rdgz

Reputation: 13254

FFTW output vector size is wrong after executing the FFT plan

I am having trouble executing an FFT plan with fftw library. I am declaring fftIn and fftOut vectors with size 2048, and defining a fftwf_plan_dft_1d with them.

But when I execute the plan, suddenly the fftOut vector gets all wrong. In my code, I can see while debugging that its size goes from 2048 to 0. When I wrote the small executable example below, I just get bad sizes such as 17293858104588369909.

Of course, when I try to access to the first item of the vector, a SIGSEGV happens.

My code:

#include <complex>
#include <iostream>
#include <fftw3.h>
#include <vector>
using namespace std;

typedef std::vector<std::complex<float>> fft_vector;

int main() {
    const unsigned int VECTOR_SIZE = 2048;
    fft_vector* fftIn = new fft_vector(VECTOR_SIZE);
    fft_vector* fftOut = new fft_vector(VECTOR_SIZE);

    fftwf_plan fftPlan = fftwf_plan_dft_1d(VECTOR_SIZE, reinterpret_cast<fftwf_complex*>(fftIn), reinterpret_cast<fftwf_complex*>(fftOut), FFTW_FORWARD, FFTW_ESTIMATE);

    fftwf_execute(fftPlan);

    std::cout << fftOut->size() << std::endl;
    std::cout << fftOut->at(0).real() << fftOut->at(0).imag() << std::endl;

    return 0;
}

Of course, I know fftIn vector is empty in this example, but output is broken when it is not anyway. In this case, the SIGSEGV happens in the second cout as described before.

My full code has threads (but FFT happens all within the same thread, so race conditions should not apply), and that was one of the reasons for trying to isolate the code in this small example, just in case, but it seems to have something wrong anyway.

Any idea?

Upvotes: 0

Views: 934

Answers (2)

Andrew
Andrew

Reputation: 518

The main issue is that you are passing a vector to it, which won't work; you need to pass the vector contents: &((*fftIn)[0]) and &((*fftOut)[0]) or something comparable. As it is, you are telling fftw to stomp on the vector object's metadata (including the length, which explains why it was sometimes 0 and sometimes gibberish). It was literally writing to the start of the vector structure because that's what the pointer points to. Also it was using fftIn's metadata as part of the input to the fft, which is also not what you want.

You might consider using fftw_complex and fftw_malloc instead, which will ensure that your data is stored aligned the way fftw needs: http://www.fftw.org/doc/SIMD-alignment-and-fftw_005fmalloc.html You can put fftOut in a vector afterwards if you really need it. This will ensure that you get the advantages of SIMD, if available, and avoid any compiler- or platform-specific behavior with complex types and memory allocation. Using fftw's type and allocator means your code will always work optimally and as you expect, regardless of which platform you run it on and which compiler you use.

Here is an example of your transform from fftw's documentation (http://www.fftw.org/doc/Complex-One_002dDimensional-DFTs.html):

#include <fftw3.h>
...
{
    fftw_complex *in, *out;
    fftw_plan p;
    ...
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    ...
    fftw_execute(p); /* repeat as needed */
    ...
    fftw_destroy_plan(p);
    fftw_free(in); fftw_free(out);
}

Try to match this example and your code should do what you expect.

EDIT: If you must use C++ complex and vector, this should work:

const unsigned int VECTOR_SIZE = 2048;
fft_vector* fftIn = new fft_vector(VECTOR_SIZE);
fft_vector* fftOut = new fft_vector(VECTOR_SIZE);

fftwf_plan fftPlan = fftwf_plan_dft_1d(VECTOR_SIZE, 
    reinterpret_cast<fftwf_complex*>(&(*fftIn)[0]),
    reinterpret_cast<fftwf_complex*>(&(*fftOut)[0]),
    FFTW_FORWARD, FFTW_ESTIMATE);

fftwf_execute(fftPlan);

std::cout << fftOut->size() << std::endl;
std::cout << fftOut->at(0).real() << fftOut->at(0).imag() << std::endl;

Notice that you have to dereference the vector pointer in order to get the [] operator to work; getting index 0 gives you the first element of the actual vector data and so the address of that element is the address (pointer) you need to give to fftw_plan.

Upvotes: 3

Jean-Fran&#231;ois Fabre
Jean-Fran&#231;ois Fabre

Reputation: 140196

you cannot reinterpret_cast a vector to a pointer on complex elements: you're accessing the vector inner variables (like the size), not necessarily the vector raw data, which probably explains why the object is trashed in output, and all the bad things that happen next.

On the other hand, casting the address of the vector data works.

Just do (no need to new the vectors too, use standard declaration):

fft_vector fftIn(VECTOR_SIZE);
fft_vector fftOut(VECTOR_SIZE);
// load your data in fftIn
fftwf_plan_dft_1d(VECTOR_SIZE, static_cast<const fftwf_complex*>(&fftIn[0]), static_cast<fftwf_complex*>(&fftOut[0]), FFTW_FORWARD, FFTW_ESTIMATE);

Upvotes: 1

Related Questions