具有复数的 C++ FFTW（快速傅里叶变换）

Question

我正在尝试在 C++ 中对复数使用 FFT。问题是当使用 fftw_complex 数字时，我无法生成内积或具有共同语法的复杂向量的总和。

这里是一个简化的简单程序：

#include <complex.h>
#include <fftw3.h>


int main(void){

int n=1024;
    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);
//  out = 1*in;   /!\
//  out = in+0;  / ! \

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

    return 0;
}

注释的两行不起作用：

i) 第一个在编译时给出：

test.cpp:13:10: error: invalid operands of types ‘int’ and ‘fftw_complex* {aka __complex__ double*}’ to binary ‘operator*’

ii) 第二个在执行期间给出：

*** glibc detected *** ./prog: corrupted double-linked list: 0x09350030 ***

根据 fftw.org [link] 的文档，应该接受用于操作 fftw_complex 的标准语法。为什么不使用 fftw_complex 的向量？

例如，我想计算：

r = i*gamma_p*w/w0*m*exp(-l*z);

其中i是虚数； gamma_p、w0、z为实数； w、m 和 l 是复数向量。

我试过了：

    fftw_complex* i_NUMBER = (fftw_complex*) fftw_malloc(1*sizeof(fftw_complex));
memset(i_NUMBER, 0, n*sizeof(fftw_complex));
    i_NUMBER[0][0]=0;
    i_NUMBER[0][1]=1;   //#real(i_NUMBER)=0; imag(i_NUMBER)=1; i.e. i_NUMBER=i;
fftw_complex* r = (fftw_complex*) fftw_malloc(n*sizeof(fftw_complex));

for (int i=0; i<n; i++){
    r[i] = i_NUMBER[0]*gamma_p*w[i]/w0*m[i]*pow(e_NUMBER, -l[i]*z);
}

但是编译器给出了关于操作数和操作类型的错误。

gnlse.h:58:22: error: invalid operands of types ‘fftw_complex {aka double [2]}’ and ‘double’ to binary ‘operator*’
gnlse.h:58:61: error: wrong type argument to unary minus

欢迎任何帮助，谢谢！

Answer 1

您必须编写显式循环才能使用向量。

行out = in+0;会产生两个指向同一位置的指针。当你释放它两次时，你会得到损坏的堆

更新

我想我知道可能是什么问题。默认 FFTW 模式是 C89。使用复数的自动计算几乎需要 C99 模式或使用 <complex> header 的 C++。我猜你选择了一些例子，修改了它并想知道为什么它不起作用。

在 C89 模式下 fftw_complex 只是一个双精度 [2]。因此，C89 模式下的任何表达式都需要对实部和虚部进行大量手动处理，这正是我所建议的。不能混入表达式 coerced-to-pointer 复杂和双精度、函数并期望自动处理。

在 C99 模式下，您可以获得原生复合体，这使得在表达式中使用它们成为 C89 的真正享受。

C++ 本机复合体应该也能正常工作。

有关详细说明，请参阅 http://www.fftw.org/doc/Complex-numbers.html。

我猜您现在实际上处于 C89 模式。您必须使用 C99 编译示例（通常在 linux /usr/bin/c99 上），或者移动到 C++ 并沿线使用 C++ 原生复杂

#include <complex>

...

std::complex<double>* in  = new std::complex<double>[n];
std::complex<double>* out = new std::complex<double>[n];
...
p = fftw_plan_dft_1d(n,
                     reinterpret_cast<fftw_complex*>(in),
                     reinterpret_cast<fftw_complex*>(out),
                     FFTW_FORWARD, FFTW_ESTIMATE);
...

具有复数的 C++ FFTW（快速傅里叶变换）

C++ FFTW (fast fourier transforms) with complex numbers

c++

fft

fftw