Почему вход FFTW не равен нулю, вывод возвращает все нули (C++)

Question

Хорошо, поэтому я некоторое время работаю над следующей проблемой и не могу найти ответ. Короче говоря, я реализую томографическую визуализацию для радиолокационного сигнала и не могу заставить FFTW давать мне результат, отличный от нулей. Я реализую через C++ и используя библиотеку FFTW3.

В настоящее время вход представляет собой смоделированный точечный рассеиватель в начале сцены, который имеет фазовый отклик всех реальных (см. Переменную rxphase). Из-за используемой пропускной способности и т. Д. IFFT отфильтрованного принятого сигнала должен быть довольно большим (я еще не масштабировал его), но я получаю нули для каждого импульса.

выдержка из моего кода:

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 
}; 

любая помощь будет принята с благодарностью.

EDIT: код должен быть немного более полным.

Answer 1

Путаница, вероятно, происходит с линией

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

В то время как это обновляет локальную переменную in точкам на область памяти, где хранятся ваши данные (буфер хранения в filteredData вектора), он не обновить указатель, хранящийся внутри, на plan. Область памяти, которая FFTW знает и использует в качестве входного буфера, таким образом, до сих пор один вы первоначально наделенным

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

и указаны в качестве аргумента

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

Обратите внимание, что эта область памяти остается инициализируется (который может оказаться нулями) на протяжении всей вашей программы.

Вы должны вместо этого написать непосредственно в in входного буфера с:

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(); 
} 

Излишне говорить, что в этом случае in = reinterpret_cast<fftw_complex*>(&filteredData[0]) линии также должны быть удалены.

В качестве альтернативы, если вы используете компилятор, где fftw_complex и std::complex<double> являются двоичными совместимы (см FFTW documentation of Complex numbers type), вы можете предпочесть:

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]; 
} 

Answer 2

Для быстрого и грязного способа мышления об этом, вы создаете fftw_plan с fftw_plan_dft_1dсначала, измените содержимое эффективного входного массива на fftw_plan_dft_1d, как указано в описании SleuthEye, и , затем наберите fftw_execute.