Roman Rdgz Roman Rdgz - 7 days ago 5
C++ Question

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?

Answer

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.