1

I'm trying to implement the Matlab fft2() function in C using the FFTW3 library.

However, I've got different results.

Considering the next matrix: 

Z=[ 
    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765     
    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695     
    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765     
    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695     
    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765     
    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695     
    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765     
    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695     
    .... 
]

And using the following code: 

Z-> Double* 

fftw_complex* fft2; 
fft2 = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*Samples*(Lines)); 

fftw_plan p1; 

p1 = fftw_plan_dft_r2c_2d(Lines,Samples, Z, fft2, FFTW_ESTIMATE); 
fftw_execute(p1);

The results with Matlab: 

fft2= [ 
 5534,25859596829 + 0,00000000000000i     186,747610745237 - 529,515274347496i
 42,6452471730436 - 321,074636721419i    -21,4495750160608 - 190,407528614266i
-50,3875107145668 - 50,5480303619799i     30,1151029075525 + 378,240946095017i
-196,295569635431 + 228,972218925794i     35,6434356803659 - 5,46216875816971i
 36,2702126322693 - 38,5502177293316i     18,5093049539101 - 33,4608602804025i
     .... 
 ]

The results with my C code: 

5534.260423 + 0.000000 i           186.731496 + -529.495788 i 
  42.655319 + -321.068356 i        -21.425010 + -190.382717 i 
 -50.277195 + -50.384210 i          29.909846 + 377.823957 i 
 -195.767224 + 228.693862 i         35.241375 + -5.315382 i 
 36.134134 + -38.527643 i           18.406395 + -33.467351 i 
    .... 
]

What am I doing wrong?

Shai
  • 111,146
  • 38
  • 238
  • 371
user3789508
  • 19
  • 1
  • 3
    You have wrong expectations. Computer calculations are in general not exact. – undur_gongor Jun 30 '14 at 08:33
  • 2
    @JoachimPileborg No. The MATLAB engine is written in Fortran and/or C using good old fashioned IEEE754 floating point types. – David Heffernan Jun 30 '14 at 08:43
  • @DavidHeffernan MATLAB uses double, not float – mch Jun 30 '14 at 08:47
  • 3
    @mch I did not say otherwise. Or are you suggesting that the IEEE754 double precision type is not a floating point type? – David Heffernan Jun 30 '14 at 08:49
  • @undur_gongor You're correct, but OP's discrepancy is large enough to warrant concern IMO. If the same types and similar algorithms are used, then the results should only be off on the order of machine precision. It's possible the version of FFTW you're using is single precision (I'm pretty sure the MATLAB version is double precision). But from my experiences, even though MATLAB claims to used FFTW, I think they have altered the code somehow (maybe optimized the codelets better). I noted a marked different in speed [here](http://stackoverflow.com/questions/15301426/fftw-vs-matlab-fft). – JustinBlaber Jun 30 '14 at 16:08

2 Answers2

1

You aren't doing anything wrong in your C implementation but you can't be sure that you are comparing like with like.

There are varying reasons why you could achieve different results between FFTW in C and fft2() in MATLAB.

  • MATLAB uses its own specifically optimised version of FFTW, compiled with multithreading support and SSE/AVX vectorised instructions. You can find details of it with version('-fftw') in MATLAB.
  • MATLAB uses it owns abstraction (libmwmfl_fft) on top of FFTW to do away with the planner routines and expose simple functions like fft() and fft2(). You can't be sure that it will choose the same planner routine that you have in your question.
  • FFTW uses heuristics to determine the optimal algorithm for computing the FFT for the data size and type you've specified. These heuristics could vary from run to run on the same data and the same machine especially with FFTW_ESTIMATE as it uses a less rigourous testing procedure that FFTW_MEASURE, FFTW_PATIENT or FFTW_EXHAUSTIVE.
  • Floating point operations are not necessarily associative and could result in different values being computed depending on the algorithm FFTW decides to use to compute the FFT. Depending on the algorithm chosen, floating point round-off error could be propagating through some of the computations resulting in lesser accuracy.

I've looked at the relative error for the data you've posted in your question and for some values it is 1e^-5 whilst the mean relative error is 0.2341.

On the whole though given the nature of floating point calculations who is to say which computes the correct values? FFTW in C or MATLAB?

IKavanagh
  • 6,089
  • 11
  • 42
  • 47
0

Before comparing Matlab's fft2 vs FFTW, you should verify reproducible results between execution given the same input, for both Matlab and FFTW. Because the default FFTW will not produce the same result every time for the same input. I know this is disturbing, (and will probably induce a lot of negative comments on this answer.) The differences are small.

I ran into this issue myself with FFTW. I use MD5 checksums of the output to verify no difference. I was able to reproduce equal checksums (same output) by changing the default Planner Flags for FFTW. This issue is covered in Q3.8 of the FFTW FAQ.

Once you have deterministic FFTW working, do the same for Matlab (since it just wraps FFTW, and offers options to modify the FFTW plan).

Then compare FFTW and Matlab results.

Tyson Hilmer
  • 741
  • 7
  • 25