why is FFTW input non-zero, output returns all zeros (C++)

Question

Ok, so I have been working at the following problem for a while now and can't seem to come to an answer. In short, I am implementing tomographic imaging for a radar signal, and can't seem to get FFTW to give me an output other than zeros. I am implementing via C++ and using the FFTW3 library.

The input for now is a simulated point scatterer at the origin of the scene, which has a phase response of all real ones (see rxphase variable). Because of bandwidth used etc, the IFFT of the filtered received signal should be fairly large (I haven't scaled it yet) but I am getting zeros for each pulse.

excerpt from my code:

void ifftshift(std::vector<phdata> &argin, std::size_t count){
  int k = 0;
  int c = (int)floor((float)count/2);

  if (count % 2 == 0)
  {
    for (k=0; k<c;k++)
    {
      complex<double> tmp = argin[k];
      argin[k] = argin[k+c];
      argin[k+c] = tmp;
    }
  }
  else
  {
    complex<double> tmp = argin[count-1];
    for (k=c-1; k>=0; k--)
    {
      argin[c+k+1] = argin[k];
      argin[k] = argin[c+k];
    }
    argin[c] = tmp;
  }
};

void main(){

  std::vector<complex<double> > filteredData;
  // Number of pulses across the aperture
  std::int pulses = 11;
  // Define the frequency vector
  std::vector<double> freq;
  std::double freqmin = 9 * pow(10,9);
  std::double freqmax = 11 * pow(10,9);
  std::vector<complex<double> > outData;
  std::vector<complex<double> > rxphase;

  for  (int i = 0; i<64; i++)
  {
    freq.push_back(freqmin + (((freqmax-freqmin)/64)*i));
  }

  // Create a (samples X pulses) array of complex doubles
  rxphase.assign(64, std::vector<complex<double> > (11, {1,0}) );

  fftw_complex *in, *out;
  fftw_plan plan;

  in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * image.Nfft );
  out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * image.Nfft );
  plan = fftw_plan_dft_1d(image.Nfft, in, out, FFTW_BACKWARD, FFTW_MEASURE);

  // iterate through pulses to back project across the aperture
  for ( int i=0; i<=pulses; i++ )
  {
    // Reset vectors for each pulse
    filteredData.clear();
    outData.clear();

    for ( int ii=0; ii<freq.size(); ii++ )
    {
      // Apply simple ramp filter
      filteredData.push_back(freq[ii]*rxphase[ii][i]);
    }

    // Shift the data to put the center frequency at DC position (zero index)
    ifftshift(filteredData, filteredData.size());

    for (int ii=0; ii<filteredData.size(); ii++)
    {
      std::cout << filteredData[ii] ;
    }
    std::cout << std::endl;
    // filteredData is what I expect it to be.

    // Recast the vector to the type needed for FFTW and compute FFT
    in = reinterpret_cast<fftw_complex*>(&filteredData[0]);

    for (int ii=0; ii<filteredData.size(); ii++)
    {
      std::cout << "(" << (in[ii])[0] << "," << (in[ii])[1] << ");";
    }
    std::cout << std::endl;
    // values for in are the same as for filteredData, as expected.

    fftw_execute(plan);

    for (int ii=0; ii<filteredData.size(); ii++)
    {
      std::cout << "(" << (out[ii])[0] << "," << (out[ii])[1] << ");";
    }
    std::cout << std::endl;
    // The values for out are all (0,0)
  }

  fftw_destroy_plan(plan);
  //fftw_free(in); error here
  //fftw_free(out); error here
};

any help would be greatly appreciated.

EDIT: code should be a little more complete.

Answer 1

For a quick and dirty way of thinking about this, you create fftw_plan with fftw_plan_dft_1d first, modify the contents of the effective input array to fftw_plan_dft_1d as SleuthEye mentions, and then call fftw_execute.

Answer 2

The confusion probably originates with the line

in = reinterpret_cast<fftw_complex*>(&filteredData[0]);

While this updates the local variable in to points to the area in memory where your data is stored (the storage buffer of the filteredData vector), it does not update the pointer kept internally by the plan. The area of memory that FFTW knows about and uses as input buffer is thus still the one you had originally allocated with

in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * image.Nfft );

and specified as argument to

fftw_plan_dft_1d(image.Nfft, in, out, FFTW_BACKWARD, FFTW_MEASURE);
//                           ^^

Note that this memory area remains uninitialized (which may happen to be zeros) throughout your program.

You should instead write directly to the in input buffer with:

for ( int ii=0; ii<freq.size(); ii++ )
{
    // Apply simple ramp filter
    std::complex<double> value = freq[ii]*rxphase[ii][i];
    in[ii][0] = value.real();
    in[ii][1] = value.imag();
}

Needless to say that in this case the in = reinterpret_cast<fftw_complex*>(&filteredData[0]) line should also be removed.

Alternatively if you are using a compiler where fftw_complex and std::complex<double> are binary compatible (see FFTW documentation of Complex numbers type), you may prefer:

std::complex<double>* filteredData = reinterpret_cast<std::complex<double>*>(in);
for ( int ii=0; ii<freq.size(); ii++ )
{
    // Apply simple ramp filter
    filteredData[ii] = freq[ii]*rxphase[ii][i];
}