如何在C++中使用FFTW库计算3D数组的FFT？

Question

如何在C++中使用FFTW库计算3D数组的FFT？

3

我在MATLAB中有一个三维数组，如下所示：

val(:,:,1) =

   1.0000 + 1.0000i   2.0000 + 2.0000i
   3.0000 + 3.0000i   4.0000 + 4.0000i


val(:,:,2) =

   5.0000 + 5.0000i   6.0000 + 6.0000i
   7.0000 + 7.0000i   8.0000 + 8.0000i

我在MATLAB中进行一维FFT，如下所示：

clear
close all
%3D ARRAY Init
A3D(:,:,1)=[1+1i,2+2i;3+3i,4+4i];
A3D(:,:,2)=[5+5i,6+6i;7+7i,8+8i];
n=size(A3D,1);
res=fft(A3D,n,1)

FFT的结果是：

val(:,:,1) =

   4.0000 + 4.0000i   6.0000 + 6.0000i
  -2.0000 - 2.0000i  -2.0000 - 2.0000i


val(:,:,2) =

  12.0000 +12.0000i  14.0000 +14.0000i
  -2.0000 - 2.0000i  -2.0000 - 2.0000i

我用C++写了一个程序，像下面这样转换输入到行主序格式，并发送到fftw_plan_dft_1d进行FFT，但我的结果与MATLAB的结果不同。有人能帮我解决这个问题吗？

#include <cstdio>

#include <fftw3.h>

int main()
{
    fftw_complex *in, *out;
    fftw_plan p;

    int N = 8;

    in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
    for (int i=0;i<8;i++){
        in[i][0]=i+1;
        in[i][1]=i+1;
    }


    for(int i = 0; i < N; i++)
{
    printf("%f\t%f\n", in[i][0],in[i][1]);
}

    p = fftw_plan_dft_1d(8, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(p); /* repeat as needed */

    for(int i = 0; i < N; i++)
{
    printf("%f\t%f\n", out[i][0],out[i][1]);
}


    fftw_destroy_plan(p);
    fftw_free(in); fftw_free(out);

    return 0;
}

- Alexanov

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ted Lyngmo · Answer 1

用户应注意，FFTW计算未归一化的DFT，其指数的符号由fftw_create_plan的dir参数给出。因此，计算正向变换后跟反向变换（或反之亦然）会导致原始数组按n缩放。

因此，您应该进行归一化。以下是一个示例，其中包括归一化和执行反向变换以获取原始值，您将注意到这将起作用（但不进行归一化将无法起作用）。您期望的值并非由Matlab生成，因此我无法解释该部分。

#include <fftw3.h>

#include <iomanip>
#include <ios>
#include <iostream>

std::ostream& operator<<(std::ostream& os, const fftw_complex& c) {
    return os << std::setw(10) << c[0] << " + " << std::setw(10) << c[1] << 'i';
}

void display(const char* title, const fftw_complex* arr, size_t N) {
    std::cout << title << ":\n";
    for(size_t i = 0; i < N; i++) {
        std::cout << arr[i] << '\n';
    }
}

void normalize(fftw_complex* arr, size_t N) {
    for(size_t bin = 0; bin < N; ++bin) {
        arr[bin][0] /= static_cast<double>(N);
        arr[bin][1] /= static_cast<double>(N);
    }
}

int main() {
    fftw_complex *in, *out;
    fftw_plan p, p2;
    std::cout << std::setprecision(6) << std::fixed;

    size_t N = 8;

    in = fftw_alloc_complex(N);
    out = fftw_alloc_complex(N);

    for(size_t i = 1; i <= N; i++) {
        in[i - 1][0] = static_cast<double>(i);
        in[i - 1][1] = static_cast<double>(i);
    }

    display("in", in, N);

    p = fftw_plan_dft_1d(static_cast<int>(N), in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    p2 = fftw_plan_dft_1d(static_cast<int>(N), out, in, FFTW_BACKWARD, FFTW_ESTIMATE);

    fftw_execute(p);

    display("out", out, N);

    normalize(out, N);
    display("normalized", out, N);

    fftw_execute(p2);
    display("backward", in, N);

    fftw_destroy_plan(p2);
    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
}

输出：

in:
  1.000000 +   1.000000i
  2.000000 +   2.000000i
  3.000000 +   3.000000i
  4.000000 +   4.000000i
  5.000000 +   5.000000i
  6.000000 +   6.000000i
  7.000000 +   7.000000i
  8.000000 +   8.000000i
out:
 36.000000 +  36.000000i
-13.656854 +   5.656854i
 -8.000000 +   0.000000i
 -5.656854 +  -2.343146i
 -4.000000 +  -4.000000i
 -2.343146 +  -5.656854i
  0.000000 +  -8.000000i
  5.656854 + -13.656854i
normalized:
  4.500000 +   4.500000i
 -1.707107 +   0.707107i
 -1.000000 +   0.000000i
 -0.707107 +  -0.292893i
 -0.500000 +  -0.500000i
 -0.292893 +  -0.707107i
  0.000000 +  -1.000000i
  0.707107 +  -1.707107i
backward:
  1.000000 +   1.000000i
  2.000000 +   2.000000i
  3.000000 +   3.000000i
  4.000000 +   4.000000i
  5.000000 +   5.000000i
  6.000000 +   6.000000i
  7.000000 +   7.000000i
  8.000000 +   8.000000i