I'm trying to perform a 2d Fourier transform of some real data. I'm using the FFTW library to do this as it's significantly faster than Armadillo's library.
Given a simple (4x4) start matrix: 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
`
If I use the built in FFT in armadillo, the output looks as follows:
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)
But if I use FFTW I get:
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
Performing the respective IFFT on both matrix BBB and CCC give exactly the starting matrix AAA.
According to the documentation: ( http://www.fftw.org/fftw3_doc/One_002dDimensional-DFTs-of-Real-Data.html#One_002dDimensional-DFTs-of-Real-Data ) : “In many practical applications, the input data in[i] are purely real numbers, in which case the DFT output satisfies the “Hermitian” redundancy: out[i] is the conjugate of out[n-i]. It is possible to take advantage of these circumstances in order to achieve roughly a factor of two improvement in both speed and memory usage.”
Therefore, matrix CCC requires some kind of operation to retrieve the Hermetian redundancy but I'm too mathematically noob to understand what that operation is. Could anyone help me with that?
Also, Armadillo stores data in a col major format and FFTW in a row major format, according to the doc this shouldn't matter as long as you pass the the row/col dimensions in reverse order to the plan function?
Thanks for looking.
Attached is my code:
#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;
}
`
You have a real function 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
Fourier transform of a real function is hermitian. Meaning that the real part of the spectrum is even, X(iw) = X(-iw)
, and the imaginary part of the spectrum is odd, X(iw)=-X(-iw)
. In other words
imag(BBB) == -imag(BBB') && real(BBB) == real(BBB')
But under these circumstances, I know that from only the coefficients on the upper diagonal for example I am able to reconstruct BBB for inverse transformation.
fftw_plan_dft_r2c_2d
also explains that this is why CCC
needs to be nd/2+1 x nd (1 column of padding) to store the output. So you are safe to declare CCC and CCCC like so:
cx_mat CCC (4,3);
cx_mat CCCC(4,3);
But as you mention yourself that FFTW is unlike Armadillo row-major, you should also reflect the consequences of that in your code:
cx_mat CCC (3,4);
cx_mat CCCC(3,4);
And suddenly your result looks completely different:
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
From what is in CCC you can reconstruct the remaining coefficients and your inverse transform is correct. Hit me if anything remains unclear.
The reconstruction can be done for example as follows:
(+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))
And CCCC
is complete for inverse transform to real.