Armadillo 2D FFTW 实数到复数 DFT Hermetian 对称性

Question

我正在尝试对一些真实数据执行二维傅里叶变换。我正在使用 FFTW 库来执行此操作，因为它比 Armadillo 的库快得多。

给定一个简单的 (4x4) 起始矩阵： AAA:

    0        0        0        0
    0   1.0000   2.0000   3.0000
    0   2.0000   4.0000   6.0000
    0   3.0000   6.0000   9.0000

`

如果我在犰狳中使用内置的 FFT，输出如下所示：

BBB： (+3.600e+01,+0.000e+00) (-1.200e+01,+1.200e+01) (-1.200e+01,+0.000e+00) (-1.200e+01,-1.200e+01) (-1.200e+01,+1.200e+01) (+0.000e+00,-8.000e+00) (+4.000e+00,-4.000e+00) (+8.000e+00,+0.000e+00) (-1.200e+01,+0.000e+00) (+4.000e+00,-4.000e+00) (+4.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00) (-1.200e+01,-1.200e+01) (+8.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00) (+0.000e+00,+8.000e+00)

但是如果我使用 FFTW，我会得到：

CCC： (+3.600e+01,+0.000e+00) (+0.000e+00,-8.000e+00) (+4.000e+00,+0.000e+00) (0,0) (-1.200e+01,+1.200e+01) (+4.000e+00,-4.000e+00) (-1.200e+01,-1.200e+01) (0,0) (-1.200e+01,+0.000e+00) (-1.200e+01,+0.000e+00) (+8.000e+00,+0.000e+00) (0,0) (-1.200e+01,+1.200e+01) (+4.000e+00,-4.000e+00) (+4.000e+00,+4.000e+00) (0,0

对矩阵 BBB 和 CCC 执行各自的 IFFT 准确给出起始矩阵 AAA。

根据文档：（http://www.fftw.org/fftw3_doc/One_002dDimensional-DFTs-of-Real-Data.html#One_002dDimensional-DFTs-of-Real-Data）： “在许多实际应用中，输入数据 in[i] 是纯实数，在这种情况下 DFT 输出满足“Hermitian”冗余：out[i] 是 out[n-i] 的共轭。可以利用这些情况在速度和内存使用方面实现大约两倍的改进。”

因此，矩阵 CCC 需要某种运算来检索 Hermetian 冗余，但我在数学上太菜鸟了，无法理解该运算是什么。谁能帮我解决这个问题？

此外，Armadillo 以 col 主要格式存储数据，以行主要格式存储 FFTW，根据文档，只要您以相反的顺序将 row/col 维度传递给计划功能？

感谢观看。

附上我的代码：

#include <iostream>
#include <fftw3.h>
#include "armadillo"

using namespace arma;
using namespace std;


int main(int argc, char** argv)
{

    mat AAA=zeros(4,4);
    mat IBB=zeros(4,4);
    cx_mat BBB(4,4);



    for (int xx=0;xx<=3;xx++){
        for ( int yy=0;yy<=3;yy++){

    AAA(xx,yy)= xx*yy;
        }
    }


    cx_mat CCC (4,4);
    cx_mat CCCC(4,4);
    mat ICC =zeros(4,4);



     fftw_plan plan=fftw_plan_dft_r2c_2d(4, 4,(double(*))&AAA(0,0), (double(*)[2])&CCC(0,0), FFTW_ESTIMATE);
     fftw_plan plan2=fftw_plan_dft_c2r_2d(4, 4,(double(*)[2])&CCCC(0,0), (double(*))&ICC(0,0), FFTW_ESTIMATE);

     //Perform Armadillo FFT (Correct output)
     BBB=fft2(AAA);
     //Perform armadillo IFFT
     IBB=real(ifft2(BBB));




    //Perform FFTW- FFT
    fftw_execute(plan);
    //Allocate fourier array to another array as imput array is destroyed
    CCCC=CCC;

    //Perform FFTW- IFFT on newly allocated array
    fftw_execute(plan2);
    //Must re-normalise the array by the number of elements
    ICC=ICC/(4*4);

    //myst rescale by the number of elements in the array

BBB.print("BBB:");
CCC.print("CCC:");  

IBB.print("IBB:");
ICC.print("ICC:");




    return 0;
}
`

Answer 1

你有一个真正的函数A:

0        0        0        0
0   1.0000   2.0000   3.0000
0   2.0000   4.0000   6.0000
0   3.0000   6.0000   9.0000

实函数的傅立叶变换是厄密函数。这意味着频谱的实部是偶数，X(iw) = X(-iw)，而频谱的虚部是奇数，X(iw)=-X(-iw)。换句话说

imag(BBB) == -imag(BBB') && real(BBB) == real(BBB')

但在这种情况下，我知道仅从上对角线上的系数，例如我能够重建 BBB 以进行逆变换。

fftw_plan_dft_r2c_2d 也解释了为什么 CCC 需要 nd/2+1 x nd（1 列填充）来存储输出。所以你可以像这样声明 CCC 和 CCCC：

cx_mat CCC (4,3);
cx_mat CCCC(4,3);

但是正如您自己提到 FFTW 与 Armadillo row-major 不同，您还应该在您的代码中反映其后果：

cx_mat CCC (3,4);
cx_mat CCCC(3,4);

突然间你的结果看起来完全不同了：

BBB:
(+3.600e+01,+0.000e+00) (-1.200e+01,+1.200e+01) (-1.200e+01,+0.000e+00) (-1.200e+01,-1.200e+01)
(-1.200e+01,+1.200e+01) (+0.000e+00,-8.000e+00) (+4.000e+00,-4.000e+00) (+8.000e+00,+0.000e+00)
(-1.200e+01,+0.000e+00) (+4.000e+00,-4.000e+00) (+4.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00)
(-1.200e+01,-1.200e+01) (+8.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00) (+0.000e+00,+8.000e+00)
CCC:
(+3.600e+01,+0.000e+00) (-1.200e+01,+1.200e+01) (-1.200e+01,+0.000e+00) (-1.200e+01,-1.200e+01)
(-1.200e+01,+1.200e+01) (+0.000e+00,-8.000e+00) (+4.000e+00,-4.000e+00) (+8.000e+00,+0.000e+00)
(-1.200e+01,+0.000e+00) (+4.000e+00,-4.000e+00) (+4.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00)
IBB:
    0        0        0        0
    0   1.0000   2.0000   3.0000
    0   2.0000   4.0000   6.0000
    0   3.0000   6.0000   9.0000
ICC:
    0        0        0        0
    0   1.0000   2.0000   3.0000
    0   2.0000   4.0000   6.0000
    0   3.0000   6.0000   9.0000

根据 CCC 中的内容，您可以重建剩余的系数，并且您的逆变换是正确的。有什么不明白的可以私信我

重构可以这样进行，例如：

(+3.6e+01,+0.0e+00) (-1.2e+01,+1.2e+01) (-1.2e+01,+0.0e+00 (-1.2e+01,-1.2e+01)
(-1.2e+01,+1.2e+01) (+0.0e+00,-8.0e+00) (+4.0e+00,-4.0e+00 (+8.0e+00,+0.0e+00)
(-1.2e+01,+0.0e+00) (+4.0e+00,-4.0e+00) (+4.0e+00,+0.0e+00 (+4.0e+00,+4.0e+00)
conj(CCC(1,2))      conj(CCC(4,2))      conj(CCC(4,3))     conj(CCC(2,2))

并且CCCC完成了对实数的逆变换。

Answer 2

只是添加到 Kaveh 的答案（已接受的答案）中，为了完全恢复 Hermetian 冗余，有必要按照他的答案中所述执行 DFT，然后 select 忽略零频率的子矩阵。从左到右翻转，上下翻转，然后对结果矩阵进行复共轭。希望这对其他人有帮助。这是代码

#include <iostream>
#include <fftw3.h>
#include "armadillo"

using namespace arma;
using namespace std;


int main(int argc, char** argv)
{

    mat AAA=zeros(6,6);
    mat IBB=zeros(6,6);
    cx_mat ccgin(6,6);
    cx_mat ccgout(6,6);
    cx_mat BBB(6,6);



    for (int xx=0;xx<=5;xx++){
        for ( int yy=0;yy<=5;yy++){

    AAA(xx,yy)= xx*(xx+yy);
        }
    }


    cx_mat CCC (4,6);
    cx_mat CCCC(4,6);
    mat ICC =zeros(6,6);

    cx_mat con(3,5);



     fftw_plan plan=fftw_plan_dft_r2c_2d(6, 6,(double(*))&AAA(0,0), (double(*)[2])&CCC(0,0), FFTW_ESTIMATE);
     fftw_plan plan2=fftw_plan_dft_c2r_2d(6, 6,(double(*)[2])&CCCC(0,0), (double(*))&ICC(0,0), FFTW_ESTIMATE);

     //Perform Armadillo FFT (Correct output)
     BBB=fft2(AAA);

     //Perform armadillo IFFT
     IBB=real(ifft2(BBB));




    //Perform FFTW- FFT
    fftw_execute(plan);
    //Allocate fourier array to another array as imput array is destroyed
    CCCC=CCC;

    //Perform FFTW- IFFT on newly allocated array
    fftw_execute(plan2);
    //Must re-normalise the array by the number of elements
    ICC=ICC/(6*6);

    //must rescale by the number of elements in the array

BBB.print("BBB:");
CCC.print("CCC:");  

IBB.print("IBB:");
ICC.print("ICC:");


//Recover Hermetian redundancy
con=fliplr(flipud(conj(CCC(span(1,3),span(1,5)))));
con.print("fliplr(flipud(conj(CCC(span(1,3),span(1,5)))));:");


    return 0;
}

Armadillo 2D FFTW 实数到复数 DFT Hermetian 对称性

Armadillo 2D FFTW real to complex DFT Hermetian symmetry

2d

fftw

armadillo

symmetry