FFT : FFTW Matlab FFT2 mystery

Question

I inherited an old fortran code with fft subroutine and I am unable to trace the source of that program. The only thing I know is that there is a call to ff2prp() and the call to fft2() to perform 2D forward and backward DFT. In order to know what the code is doing I take DFT of a 4x4 2D array (matrix) and the results are very different from Matlab and FFTW results.

Question: Can someone tell what the code is doing by looking at the output. The input and output are both real arrays

Input array

 0.20000     0.30000     1.00000     1.20000
 0.00000    12.00000     5.00000     1.30000
 0.30000     0.30000     1.00000     1.40000
 0.00000     0.00000     0.00000     1.50000

After forward FFT with fft2() fortran routine

 0.16875    -0.01875    -0.05000     0.05625
 0.00000    12.00000     5.00000     1.30000
 0.30000     0.30000     1.00000     1.40000
 0.00000     0.00000     0.00000     1.50000

Matlab output performing DCT: dct2(input)

    6.3750   -0.8429   -3.4250   -2.4922
    2.4620    0.6181   -2.6356   -0.9887
   -4.2750   -0.9798    4.2250    2.2730
   -4.8352   -1.2387    5.0695    3.4819

Output from C++ code using FFTW library. DCT from FFTW

(6.3750, 0.00)  (-0.8429, 0.00) (-3.4250, 0.00) (-2.4922, 0.00) 
(2.4620, 0.00)  (0.6181, 0.00)  (-2.6356, 0.00) (-0.9887, 0.00) 
(-4.2750, 0.00) (-0.9798, 0.00) (4.2250, 0.00)  (2.2730, 0.00)  
(-4.8352, 0.00) (-1.2387, 0.00) (5.0695, 0.00)  (3.4819, 0.00)

Forward FFT with Matlab - fft2(input)

  25.5000 + 0.0000i  -6.5000 - 7.2000i -10.5000 + 0.0000i  -6.5000 + 7.2000i
  -0.3000 -16.8000i -12.3000 + 4.8000i   0.1000 + 6.8000i  12.1000 + 5.2000i
 -14.1000 + 0.0000i   3.5000 +11.2000i   9.1000 + 0.0000i   3.5000 -11.2000i
  -0.3000 +16.8000i  12.1000 - 5.2000i   0.1000 - 6.8000i -12.3000 - 4.8000i

Forward FFT with FFTW

(25.50, 0.00)   (-6.50, -7.20)  (-10.50, 0.00)  (-6.50, 7.20)   
(-0.30, -16.80) (-12.30, 4.80)  (0.10, 6.80)    (12.10, 5.20)   
(-14.10, 0.00)  (3.50, 11.20)   (9.10, 0.00)    (3.50, -11.20)  
(-0.30, 16.80)  (12.10, -5.20)  (0.10, -6.80)   (-12.30, -4.80)

As you cane see the output of Matlab and FFTW agrees with each other but not with the output of the fortran code. I would like to use FFTW but the results are different due to FFT. I can't figure out what FFT the fortran program is doing. Can anyone tell by looking at the output.

Answer 1

As far as I can tell, fft2 seems to have computed the 1D FFT of the first row (leaving all 3 others unchanged), with the result being scaled by 1/16 and packed in r0, r2, r1, i1 format.

In other words the output can be constructed in Matlab using:

input = [0.2 0.3 1 1.2;0 12 5 1.3;0.3 0.3 1 1.4;0 0 0 1.5];
N = size(input,2);
A = fft(input(1,:))/16;
B = reshape([real(A);imag(A)],1,2*N);
B(2) = B(N+1);
output = [B(1:N);A(2:size(input,1),:)];

If you have some reasons to believe that fft2 should be computing 2D FFTs, then there might some problem in the way you pass your data to this routine which causes incorrect results. Also, additional test cases (or how you call ff2prp) for different sizes inputs may provide more insight about the choice of scaling factor (eg. is it 1/N^2 or 1/4N, or something else).