c++ 错误结果,当使用FFTW时,仅执行正向FFT和逆FFT或FFTW与ESTIMATE和MEASURE?

lh80um4z  于 2023-07-01  发布在  其他
关注(0)|答案(1)|浏览(130)

在使用FFTW时,我只对FFT进行了正向FFT,然后对array进行了逆FFT,但是当FFT的length具有特定值时,例如length==84,我得到了错误的结果
代码示例:

int main() {
  const int rows = 84;
  fftw_complex *in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * rows);
  memset(in, 0, sizeof(fftw_complex) * rows);
  in[0][0] = in[1][0] = in[2][0] = 1.0;
  fftw_plan plan_fft  = fftw_plan_dft_1d(rows, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_plan plan_ifft = fftw_plan_dft_1d(rows, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(plan_fft);
  fftw_execute(plan_ifft);
  for (int i = 0; i < rows; ++i) {
    std::cout << in[i][0] / rows << " ";
  }
  std::cout << std::endl;
  return 0;
}

如上面代码所示,wanted result为:
实际结果为:
当我将FFT的标志值从FFTW_ESTIMATE更改为FFTW_MEASURE时,结果正确。
为什么会发生这种情况?我已经困惑了很长一段时间,有人能帮助我吗?
只有当FFT length的值为例如84时,才会产生这样的误差,而其他值又是正常的,例如95

7d7tgy0s

7d7tgy0s1#

它使用浮点数学,所以它不太可能精确地给予01,而是一些非常接近的值。
你也可以尝试其他的planner flags

  • FFTW_ESTIMATE规定,不是使用不同算法的实际测量,而是使用简单的启发式来快速选择(可能是次优的)计划。有了这个标志,输入/输出阵列在计划期间不会被覆盖。
  • FFTW_MEASURE告诉FFTW通过实际计算几个FFT并测量它们的执行时间来找到优化的计划。根据您的计算机,这可能需要一些时间(通常为几秒钟)。FFTW_MEASURE是默认计划选项。
  • FFTW_PATIENT类似于FFTW_MEASURE,但考虑了更广泛的算法,并且通常会产生“更优化”的计划(特别是对于大型变换),但代价是规划时间延长了几倍(特别是对于大型变换)。
  • FFTW_EXHAUSTIVE类似于FFTW_PATIENT,但考虑了更广泛的算法,包括许多我们认为不太可能快速的算法,以产生最优的计划,但大大增加了规划时间。
  • FFTW_WISDOM_ONLY是一种特殊的计划模式,在这种模式下,只有在对给定问题有智慧时才创建计划,否则返回NULL计划。这可以与其他标志组合,例如FFTW_WISDOM_ONLY | FFTW_PATIENT仅当在FFTW_PATIENTFFTW_EXHAUSTIVE模式下创建的智慧可用时才创建计划。FFTW_WISDOM_ONLY标志用于需要检测智慧是否可用的用户;例如,如果智能不可用,则可能希望分配新的阵列用于规划,使得用户数据不被重写。

您还可以使用std::fixed来更容易地看到结果更接近您想要的结果。
此外,由于您使用的是C++,因此使用std::vector<std::complex<double>>作为容器,而不是fftw_malloc,并且仅在与FFTW库交互时转换为fftw_complex。这样就不必手动释放分配的内存-或者像当前程序那样发生内存泄漏。
重写为使用C++容器、不同的规划器标志和std::fixed

#include <fftw3.h>

#include <complex>
#include <iostream>
#include <vector>

// helpers to interact with C functions:
template<class T>
using two = T[2];

template<class T>
two<T>* to_c(std::complex<T>* in) {
    return reinterpret_cast<two<T>*>(in);
}

template<class T>
two<T>& to_c(std::complex<T>& in) {
    return reinterpret_cast<two<T>&>(in);
}
// ---------------------------------------

int main() {
    constexpr unsigned rows = 84;
    std::vector<std::complex<double>> in(rows);

    fftw_plan plan_fft = fftw_plan_dft_1d(
        rows, to_c(in.data()), to_c(in.data()), FFTW_FORWARD, FFTW_MEASURE);
    fftw_plan plan_ifft = fftw_plan_dft_1d(
        rows, to_c(in.data()), to_c(in.data()), FFTW_BACKWARD, FFTW_MEASURE);

    to_c(in[0])[0] = to_c(in[1])[0] = to_c(in[2])[0] = 1.0;
    // or, the C++ way:
    //in[0].real(1.0);
    //in[1].real(1.0);
    //in[2].real(1.0);

    fftw_execute(plan_fft);
    fftw_execute(plan_ifft);

    std::cout << std::fixed;

    for(unsigned i = 0; i < rows; ++i) {
        std::cout << in[i].real() / rows << '\n';
    }

    std::cout << '\n';
}

相关问题