Another Matlab fft2 != FFTW

Hello,

I have something strange...

There is a matrix of 5x5, padded with zeroes to 9x9:
Columns 1 through 7

   -0.1449   -0.0888    0.0095   -0.0816    0.1081         0         0
    0.0013   -0.1255   -0.0517   -0.0423    0.0090         0         0
   -0.0156    0.2098    0.2978   -0.1366   -0.1449         0         0
    0.0103    0.0943    0.0068    0.0311   -0.0550         0         0
    0.0711    0.0695   -0.0577    0.0572   -0.0314         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0

  Columns 8 through 9

         0         0
         0         0
         0         0
         0         0
         0         0
         0         0
         0         0
         0         0
         0         0

In matlab, an fft2 yields the following:
 Columns 1 through 3

   0.0000             0.2733 - 0.1159i  -0.2439 - 0.4495i
  -0.4673 - 0.1859i  -0.6651 - 0.2544i  -0.6534 + 0.4118i
  -0.3924 + 0.2796i  -0.4806 + 0.5024i   0.3644 + 0.6977i
  -0.1653 + 0.2693i  -0.0195 + 0.3936i   0.5265 - 0.1896i
   0.1352 + 0.2382i   0.2658 + 0.0587i   0.0051 - 0.4679i
   0.1352 - 0.2382i  -0.1037 - 0.5388i  -0.5974 - 0.2733i
  -0.1653 - 0.2693i  -0.7101 - 0.2598i  -0.5318 + 0.5770i
  -0.3924 - 0.2796i  -0.7927 + 0.2416i   0.2980 + 0.5578i
  -0.4673 + 0.1859i  -0.2157 + 0.7054i   0.4214 - 0.1171i

  Columns 4 through 6

  -0.3748 + 0.1383i  -0.0044 + 0.1137i  -0.0044 - 0.1137i
  -0.2228 + 0.3993i   0.1978 + 0.4462i   0.2247 - 0.4734i
   0.3786 - 0.0615i   0.2339 - 0.1726i  -0.3105 - 0.3324i
  -0.3444 - 0.5698i  -0.2941 + 0.1250i   0.0202 - 0.0751i
  -0.5753 - 0.0160i  -0.0069 + 0.2812i  -0.0116 - 0.2485i
  -0.4691 + 0.2022i  -0.0116 + 0.2485i  -0.0069 - 0.2812i
   0.0517 + 0.3219i   0.0202 + 0.0751i  -0.2941 - 0.1250i
  -0.0253 - 0.4140i  -0.3105 + 0.3324i   0.2339 + 0.1726i
  -0.5867 - 0.0759i   0.2247 + 0.4734i   0.1978 - 0.4462i

  Columns 7 through 9

  -0.3748 - 0.1383i  -0.2439 + 0.4495i   0.2733 + 0.1159i
  -0.5867 + 0.0759i   0.4214 + 0.1171i  -0.2157 - 0.7054i
  -0.0253 + 0.4140i   0.2980 - 0.5578i  -0.7927 - 0.2416i
   0.0517 - 0.3219i  -0.5318 - 0.5770i  -0.7101 + 0.2598i
  -0.4691 - 0.2022i  -0.5974 + 0.2733i  -0.1037 + 0.5388i
  -0.5753 + 0.0160i   0.0051 + 0.4679i   0.2658 - 0.0587i
  -0.3444 + 0.5698i   0.5265 + 0.1896i  -0.0195 - 0.3936i
   0.3786 + 0.0615i   0.3644 - 0.6977i  -0.4806 - 0.5024i
  -0.2228 - 0.3993i  -0.6534 - 0.4118i  -0.6651 + 0.2544i

Now, I have following program:
#include <complex>
#include <iostream>
#include <fftw3.h>
#include <algorithm>

double data [81] = {
-0.144892,  -0.0888301,  0.00954482,  -0.0816278, 0.108072,
0,  0,  0,  0,
0.00130495,  -0.125482,  -0.0517049,  -0.0422596,  0.00901076,
0,  0,  0,  0,
-0.0155715,  0.209789,  0.297849,  -0.136606,  -0.144892,
0,  0,  0,  0,
0.010323,  0.0943392,  0.0067982,  0.0311058,  -0.0550466,
0,  0,  0,  0,
0.0710997,  0.069528,  -0.0577017,  0.0572293,  -0.0313799,
0,  0,  0,  0,
0,  0,  0,  0,  0,  0,  0,  0,  0,
0,  0,  0,  0,  0,  0,  0,  0,  0,
0,  0,  0,  0,  0,  0,  0,  0,  0,
0,  0,  0,  0,  0,  0,  0,  0,  0};


int main()
{
	double in [81];
	std::complex<double> out2d [81];

	fftw_plan p2d = fftw_plan_dft_r2c_2d(9, 9, in,
				(fftw_complex*)out2d,
					   FFTW_ESTIMATE);

	std::copy(data, data+81, in);

	for (int x = 0; x < 9; ++x)
	{
		for (int y = 0; y < 9; ++y)
			std::cout << in[9*x+y] << "   ";
		std::cout << std::endl;
	}
	std::cout << std::endl;

	fftw_execute(p2d);
	fftw_destroy_plan(p2d);

	for (int x = 0; x < 9; ++x)
	{
		for (int y = 0; y < 9; ++y)
			std::cout << out2d[9*x+y] << "   ";
		std::cout << std::endl;
	}
	std::cout << std::endl;

	return 0;
}

And compute the 2D FFT with it. Then the computed values are:
(-3.7e-07,0)   (0.273315,-0.115934)   (-0.243932,-0.449487)
(-0.374841,0.138284)   (-0.00435327,0.113728)   (-0.467347,-0.18594)
(-0.665096,-0.254414)   (-0.653422,0.411802)   (-0.222791,0.399276)
(0.197808,0.446171)   (-0.392351,0.279649)   (-0.480622,0.502412)
(0.364381,0.697738)   (0.37858,-0.0614555)   (0.233949,-0.172564)
(-0.16532,0.269268)   (-0.0194797,0.393635)   (0.526491,-0.189582)
(-0.344401,-0.569772)   (-0.29414,0.124964)   (0.135219,0.238206)
(0.265819,0.0586845)   (0.00508875,-0.46789)   (-0.575287,-0.0159722)
(-0.00689842,0.281154)   (0.135219,-0.238206)   (-0.103669,-0.538752)
(-0.597382,-0.273305)   (-0.469052,0.202163)   (-0.0115958,0.248459)
(-0.16532,-0.269268)   (-0.710068,-0.259756)   (-0.531826,0.577016)
(0.0516525,0.32187)   (0.0202206,0.0750899)   (-0.392351,-0.279649)
(-0.792725,0.241556)   (0.298028,0.557757)   (-0.0253494,-0.414048)
(-0.310451,0.33237)   (-0.467347,0.18594)   (-0.215681,0.705421)
(0.421411,-0.117124)   (-0.58673,-0.075928)   (0.224719,0.473378)
(0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)
(0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)
(0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)
(0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)

... which contains a lot of zeroes that are not in matlab's result.

Anyone an idea of what is going on here ?

Thanks in advance,

Fokko Beekhof

... the answer to my own question is that in FFTW's real-to-complex
transforms and vice versa, the number of elements is not equal as in my
code.

And this is written in FFTW's documentation.

Fokko Beekhof

Fokko Beekhof wrote:
> Hello,
> 
> I have something strange...
> 
> There is a matrix of 5x5, padded with zeroes to 9x9:
> Columns 1 through 7
> 
>    -0.1449   -0.0888    0.0095   -0.0816    0.1081         0         0
>     0.0013   -0.1255   -0.0517   -0.0423    0.0090         0         0
>    -0.0156    0.2098    0.2978   -0.1366   -0.1449         0         0
>     0.0103    0.0943    0.0068    0.0311   -0.0550         0         0
>     0.0711    0.0695   -0.0577    0.0572   -0.0314         0         0
>          0         0         0         0         0         0         0
>          0         0         0         0         0         0         0
>          0         0         0         0         0         0         0
>          0         0         0         0         0         0         0
> 
>   Columns 8 through 9
> 
>          0         0
>          0         0
>          0         0
>          0         0
>          0         0
>          0         0
>          0         0
>          0         0
>          0         0
> 
> In matlab, an fft2 yields the following:
>  Columns 1 through 3
> 
>    0.0000             0.2733 - 0.1159i  -0.2439 - 0.4495i
>   -0.4673 - 0.1859i  -0.6651 - 0.2544i  -0.6534 + 0.4118i
>   -0.3924 + 0.2796i  -0.4806 + 0.5024i   0.3644 + 0.6977i
>   -0.1653 + 0.2693i  -0.0195 + 0.3936i   0.5265 - 0.1896i
>    0.1352 + 0.2382i   0.2658 + 0.0587i   0.0051 - 0.4679i
>    0.1352 - 0.2382i  -0.1037 - 0.5388i  -0.5974 - 0.2733i
>   -0.1653 - 0.2693i  -0.7101 - 0.2598i  -0.5318 + 0.5770i
>   -0.3924 - 0.2796i  -0.7927 + 0.2416i   0.2980 + 0.5578i
>   -0.4673 + 0.1859i  -0.2157 + 0.7054i   0.4214 - 0.1171i
> 
>   Columns 4 through 6
> 
>   -0.3748 + 0.1383i  -0.0044 + 0.1137i  -0.0044 - 0.1137i
>   -0.2228 + 0.3993i   0.1978 + 0.4462i   0.2247 - 0.4734i
>    0.3786 - 0.0615i   0.2339 - 0.1726i  -0.3105 - 0.3324i
>   -0.3444 - 0.5698i  -0.2941 + 0.1250i   0.0202 - 0.0751i
>   -0.5753 - 0.0160i  -0.0069 + 0.2812i  -0.0116 - 0.2485i
>   -0.4691 + 0.2022i  -0.0116 + 0.2485i  -0.0069 - 0.2812i
>    0.0517 + 0.3219i   0.0202 + 0.0751i  -0.2941 - 0.1250i
>   -0.0253 - 0.4140i  -0.3105 + 0.3324i   0.2339 + 0.1726i
>   -0.5867 - 0.0759i   0.2247 + 0.4734i   0.1978 - 0.4462i
> 
>   Columns 7 through 9
> 
>   -0.3748 - 0.1383i  -0.2439 + 0.4495i   0.2733 + 0.1159i
>   -0.5867 + 0.0759i   0.4214 + 0.1171i  -0.2157 - 0.7054i
>   -0.0253 + 0.4140i   0.2980 - 0.5578i  -0.7927 - 0.2416i
>    0.0517 - 0.3219i  -0.5318 - 0.5770i  -0.7101 + 0.2598i
>   -0.4691 - 0.2022i  -0.5974 + 0.2733i  -0.1037 + 0.5388i
>   -0.5753 + 0.0160i   0.0051 + 0.4679i   0.2658 - 0.0587i
>   -0.3444 + 0.5698i   0.5265 + 0.1896i  -0.0195 - 0.3936i
>    0.3786 + 0.0615i   0.3644 - 0.6977i  -0.4806 - 0.5024i
>   -0.2228 - 0.3993i  -0.6534 - 0.4118i  -0.6651 + 0.2544i
> 
> Now, I have following program:
> #include <complex>
> #include <iostream>
> #include <fftw3.h>
> #include <algorithm>
> 
> double data [81] = {
> -0.144892,  -0.0888301,  0.00954482,  -0.0816278, 0.108072,
> 0,  0,  0,  0,
> 0.00130495,  -0.125482,  -0.0517049,  -0.0422596,  0.00901076,
> 0,  0,  0,  0,
> -0.0155715,  0.209789,  0.297849,  -0.136606,  -0.144892,
> 0,  0,  0,  0,
> 0.010323,  0.0943392,  0.0067982,  0.0311058,  -0.0550466,
> 0,  0,  0,  0,
> 0.0710997,  0.069528,  -0.0577017,  0.0572293,  -0.0313799,
> 0,  0,  0,  0,
> 0,  0,  0,  0,  0,  0,  0,  0,  0,
> 0,  0,  0,  0,  0,  0,  0,  0,  0,
> 0,  0,  0,  0,  0,  0,  0,  0,  0,
> 0,  0,  0,  0,  0,  0,  0,  0,  0};
> 
> 
> int main()
> {
> 	double in [81];
> 	std::complex<double> out2d [81];
> 
> 	fftw_plan p2d = fftw_plan_dft_r2c_2d(9, 9, in,
> 				(fftw_complex*)out2d,
> 					   FFTW_ESTIMATE);
> 
> 	std::copy(data, data+81, in);
> 
> 	for (int x = 0; x < 9; ++x)
> 	{
> 		for (int y = 0; y < 9; ++y)
> 			std::cout << in[9*x+y] << "   ";
> 		std::cout << std::endl;
> 	}
> 	std::cout << std::endl;
> 
> 	fftw_execute(p2d);
> 	fftw_destroy_plan(p2d);
> 
> 	for (int x = 0; x < 9; ++x)
> 	{
> 		for (int y = 0; y < 9; ++y)
> 			std::cout << out2d[9*x+y] << "   ";
> 		std::cout << std::endl;
> 	}
> 	std::cout << std::endl;
> 
> 	return 0;
> }
> 
> And compute the 2D FFT with it. Then the computed values are:
> (-3.7e-07,0)   (0.273315,-0.115934)   (-0.243932,-0.449487)
> (-0.374841,0.138284)   (-0.00435327,0.113728)   (-0.467347,-0.18594)
> (-0.665096,-0.254414)   (-0.653422,0.411802)   (-0.222791,0.399276)
> (0.197808,0.446171)   (-0.392351,0.279649)   (-0.480622,0.502412)
> (0.364381,0.697738)   (0.37858,-0.0614555)   (0.233949,-0.172564)
> (-0.16532,0.269268)   (-0.0194797,0.393635)   (0.526491,-0.189582)
> (-0.344401,-0.569772)   (-0.29414,0.124964)   (0.135219,0.238206)
> (0.265819,0.0586845)   (0.00508875,-0.46789)   (-0.575287,-0.0159722)
> (-0.00689842,0.281154)   (0.135219,-0.238206)   (-0.103669,-0.538752)
> (-0.597382,-0.273305)   (-0.469052,0.202163)   (-0.0115958,0.248459)
> (-0.16532,-0.269268)   (-0.710068,-0.259756)   (-0.531826,0.577016)
> (0.0516525,0.32187)   (0.0202206,0.0750899)   (-0.392351,-0.279649)
> (-0.792725,0.241556)   (0.298028,0.557757)   (-0.0253494,-0.414048)
> (-0.310451,0.33237)   (-0.467347,0.18594)   (-0.215681,0.705421)
> (0.421411,-0.117124)   (-0.58673,-0.075928)   (0.224719,0.473378)
> (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)
> (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)
> (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)
> (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)   (0,0)
> 
> ... which contains a lot of zeroes that are not in matlab's result.
> 
> Anyone an idea of what is going on here ?
> 
> Thanks in advance,
> 
> Fokko Beekhof
>