我正在尝试对一些真实数据执行二维傅里叶变换。 我正在使用 FFTW 库来做到这一点,因为它比犰狳的库快得多。
给定一个简单的 (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 ( : "在许多实际应用中,输入数据 [i] 是纯实数,在这种情况下,DFT 输出满足"埃尔米特"冗余:out[i] 是 out[n-i] 的共轭。可以利用这些情况,在速度和内存使用方面实现大约两倍的改进。
因此,矩阵 CCC 需要某种操作来检索 Hermet 冗余,但我在数学上太菜鸟了,无法理解该操作是什么。 谁能帮我?
此外,犰狳以 col 主要格式存储数据,以行主要格式存储 FFTW,根据文档,只要您将行/列维度以相反的顺序传递给计划函数,这应该无关紧要?
谢谢你的关注。
附上我的代码:
#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;
}
`
你有一个真正的函数 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与犰狳行少校不同时,你也应该在你的代码中反映它的后果:
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
逆变换到实数是完整的。
为了补充Kaveh的答案(接受的答案(,为了完全恢复Hermetian冗余,有必要按照他的答案中所述执行DFT,然后选择一个忽略零频率的子矩阵。做一个从左到右的翻转,一个向上向下的翻转,然后对所得矩阵进行复共轭。希望这对其他人有所帮助。这是代码
#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;
}