为什么 FFTW 输入非零,输出 returns 全零(C++)
why is FFTW input non-zero, output returns all zeros (C++)
好的,所以我已经在以下问题上工作了一段时间,但似乎无法找到答案。简而言之,我正在对雷达信号进行层析成像,但似乎无法让 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
};
如有任何帮助,我们将不胜感激。
编辑:代码应该更完整一些。
混淆可能源于行
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];
}
为了快速而粗暴地思考这个问题,您先用 fftw_plan_dft_1d
创建 fftw_plan
,然后将有效输入数组的内容修改为fftw_plan_dft_1d
正如 SleuthEye 提到的那样,然后 然后 调用 fftw_execute
.
好的,所以我已经在以下问题上工作了一段时间,但似乎无法找到答案。简而言之,我正在对雷达信号进行层析成像,但似乎无法让 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
};
如有任何帮助,我们将不胜感激。
编辑:代码应该更完整一些。
混淆可能源于行
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];
}
为了快速而粗暴地思考这个问题,您先用 fftw_plan_dft_1d
创建 fftw_plan
,然后将有效输入数组的内容修改为fftw_plan_dft_1d
正如 SleuthEye 提到的那样,然后 然后 调用 fftw_execute
.