weird results from MKL FFT

MH_A_ · ‎07-24-2013

I am trying to perfrom 1D (in place) a complex to a complex fourier transform (forward and backward). In my code I added

Status = DftiCreateDescriptor(desc1,DFTI_DOUBLE,DFTI_COMPLEX,1, N)

Status = DftiCommitDescriptor(desc1)
Status = DftiComputeForward(desc1,coef)

               do j=1,n
                coefout(j)=coef(j)*(2.d0*pi*cmplx(0.d0,1.d0)*dble(j)/tlength)**2
               enddo

Status = DftiComputeBackward(desc1,coefout)
Status = DftiFreeDescriptor(desc1)

The program finish without any error but the results looks very weired.

The input data looks like

       0.145213596D-05   0.000000000D+00
      0.235852281D-05   0.000000000D+00
      0.379253834D-05   0.000000000D+00
       0.603777471D-05   0.000000000D+00
       0.951657967D-05   0.000000000D+00
      0.148505296D-04   0.000000000D+00
       0.229435209D-04   0.000000000D+00
       0.350941882D-04   0.000000000D+00
       0.531456135D-04   0.000000000D+00
      0.796813547D-04   0.000000000D+00

However after performing forward FFT I obtained

      -0.280573615D-21   0.154142831D-43
      -0.153507094D-21   0.154142831D-43
      -0.342774500D+06   0.980908925D-44
       0.226293969D+06 -0.980908925D-44
       0.466064128D+10   0.000000000D+00
       0.830659660D+20   0.000000000D+00
       0.279607085D+32   0.280259693D-44
      0.803789090D-34 -0.280259693D-44
      0.933711541D-28   0.000000000D+00
      0.741872405D-23   0.000000000D+00

While after backward FFT I got

       0.176396978D-37 -0.267003409D-40
      -0.786617726D-27 -0.303521247D-41
      -0.209912859D-17 -0.277344992D-40
      -0.898662859D+12 -0.461027195D-41
       0.113104136D+37 -0.283342549D-40
      -0.465754617D+22 -0.627921842D-41
      -0.118460876D+32 -0.286397380D-40
       0.209648241D+35 -0.788790906D-41
       0.274214524D-33 -0.286523497D-40
      0.218714569D+03 -0.941112049D-41

The results look very weired to me. do you have any idea how can I check if the fft code works probably?

Gennady_F_Intel · ‎07-24-2013

at the first glance all are correct. Could you check the a simple example of double-precision complex-to-complex in-place 1D ( "basic_dp_complex_dft_1d.c " ) from <mkl_root>\examples\dftc\source\ directory?

SergeyKostrov · ‎07-25-2013

Here are results of my verification ( with a modified basic_dp_complex_dft_1d.c test ) and I see that all results after backward DFT transform are multiplied by 10: ..\Tests>basic_dp_complex_dft_1d.exe Intel(R) Math Kernel Library Version 10.3.12 Product Build 20120831 for 32-bit applications Example basic_dp_complex_dft_1d Forward and backward double-precision complex in-place 1D FFTs Configuration parameters: DFTI_PRECISION = DFTI_DOUBLE DFTI_FORWARD_DOMAIN = DFTI_COMPLEX DFTI_DIMENSION = 1 DFTI_LENGTHS = {10} Create DFTI descriptor Commit DFTI descriptor Allocate input array Initialize input for forward transform Verify input data for forward FFT real=0.0000014521359600 imag=0.0000000000000000 real=0.0000023585228100 imag=0.0000000000000000 real=0.0000037925383400 imag=0.0000000000000000 real=0.0000060377747100 imag=0.0000000000000000 real=0.0000095165796700 imag=0.0000000000000000 real=0.0000148505296000 imag=0.0000000000000000 real=0.0000229435209000 imag=0.0000000000000000 real=0.0000350941882000 imag=0.0000000000000000 real=0.0000531456135000 imag=0.0000000000000000 real=0.0000796813547000 imag=0.0000000000000000 Compute forward transform Verify the result of forward FFT real=0.0002288727583900 imag=0.0000000000000000 real=0.0000315968694775 imag=0.0001279132334088 real=-0.0000276552848308 imag=0.0000726986815771 real=-0.0000415068627525 imag=0.0000402202220098 real=-0.0000460244304642 imag=0.0000182537898717 real=-0.0000471719816500 imag=0.0000000000000000 real=-0.0000460244304642 imag=-0.0000182537898717 real=-0.0000415068627525 imag=-0.0000402202220098 real=-0.0000276552848308 imag=-0.0000726986815771 real=0.0000315968694775 imag=-0.0001279132334088 Compute backward transform Verify the result of backward FFT real=0.0000145213596000 imag=0.0000000000000000 real=0.0000235852281000 imag=0.0000000000000000 real=0.0000379253834000 imag=0.0000000000000000 real=0.0000603777471000 imag=0.0000000000000000 real=0.0000951657967000 imag=0.0000000000000000 real=0.0001485052960000 imag=0.0000000000000000 real=0.0002294352090000 imag=0.0000000000000000 real=0.0003509418820000 imag=0.0000000000000000 real=0.0005314561350000 imag=0.0000000000000000 real=0.0007968135470000 imag=0.0000000000000000 Release the DFTI descriptor Free data array TEST PASSED ... For example, for a 1st element it looks like: Input .: real=0.0000014521359600 imag=0.0000000000000000 Output: real=0.0000145213596000 imag=0.0000000000000000 Is everything right?

MH_A_ · ‎07-25-2013

I think my problem is that in the mkl the indexes of the equation go from zero to n-1. However in my case the indexes go from -n/2+1 to n/2. Thus I need to perform some modification for the input and the output data before I use it.

SergeyKostrov · ‎07-25-2013

>>...I see that all results after backward DFT transform are multiplied by 10... I wonder if somebody from Intel could follow up on that isse?

Dmitry_B_Intel · ‎07-25-2013

Sergey Kostrov wrote:

>>...I see that all results after backward DFT transform are multiplied by 10...

I wonder if somebody from Intel could follow up on that isse?

Documented behaviour of Intel MKL FFT is that the scaling factor is 1, which means the forward+backward transforms result in multiplication by N, the size of the transform. If you need backward transform be the inverse of the forward, set the scale factor to 1/N for one of them, or 1/sqrt(N) for both. Other libraries follow this convention, e.g. FFTW uses factor 1, and you should incorporate the scale into surrounding computation instead of doing this while computing the FFT.

Thanks
Dima

Dmitry_B_Intel · ‎07-25-2013

MH A. wrote:

I think my problem is that in the mkl the indexes of the equation go from zero to n-1. However in my case the indexes go from -n/2+1 to n/2. Thus I need to perform some modification for the input and the output data before I use it.

If you posted a small self-contained test showing how do you set up the data arrays, compute the FFT, and use the results, it would greatly help in figuring out the reason why you don't get what you expect. It is also not a bad idea to look at the examples in mkl/examples/dftf (Fortran) directory and use them to verify if MKL really produces weird results for your case.

Thanks
Dima

SergeyKostrov · ‎07-26-2013

>>>>...Here are results of my verification ( with a modified basic_dp_complex_dft_1d.c test ) and I see that all results >>>>after backward DFT transform are multiplied by 10... >>... >>It is also not a bad idea to look at the examples in mkl/examples/dftf (Fortran) directory and use them to verify >>if MKL really produces weird results for your case... >>... I used an MKL sample basic_dp_complex_dft_1d.c for the verification 'MH A.' results. Thanks for the note regarding the scale factor and I finally fixed that issue.

SergeyKostrov · ‎07-26-2013

This is how output of the test case with 'MH A.' data looks like: ... Intel(R) Math Kernel Library Version 10.3.12 Product Build 20120831 for 32-bit applications Example basic_dp_complex_dft_1d Forward and backward double-precision complex in-place 1D DFTs Configuration parameters: DFTI_PRECISION = DFTI_DOUBLE DFTI_FORWARD_DOMAIN = DFTI_COMPLEX DFTI_DIMENSION = 1 DFTI_LENGTHS = {10} Create DFTI descriptor Commit DFTI descriptor Allocate input array Initialize input for forward transform Verify input data for forward DFT real=0.0000014521359600 imag=0.0000000000000000 real=0.0000023585228100 imag=0.0000000000000000 real=0.0000037925383400 imag=0.0000000000000000 real=0.0000060377747100 imag=0.0000000000000000 real=0.0000095165796700 imag=0.0000000000000000 real=0.0000148505296000 imag=0.0000000000000000 real=0.0000229435209000 imag=0.0000000000000000 real=0.0000350941882000 imag=0.0000000000000000 real=0.0000531456135000 imag=0.0000000000000000 real=0.0000796813547000 imag=0.0000000000000000 Compute forward transform Verify the result of forward DFT real=0.0002288727583900 imag=0.0000000000000000 real=0.0000315968694775 imag=0.0001279132334088 real=-0.0000276552848308 imag=0.0000726986815771 real=-0.0000415068627525 imag=0.0000402202220098 real=-0.0000460244304642 imag=0.0000182537898717 real=-0.0000471719816500 imag=0.0000000000000000 real=-0.0000460244304642 imag=-0.0000182537898717 real=-0.0000415068627525 imag=-0.0000402202220098 real=-0.0000276552848308 imag=-0.0000726986815771 real=0.0000315968694775 imag=-0.0001279132334088 Compute backward transform Verify the result of backward DFT real=0.0000014521359600 imag=0.0000000000000000 real=0.0000023585228100 imag=0.0000000000000000 real=0.0000037925383400 imag=0.0000000000000000 real=0.0000060377747100 imag=0.0000000000000000 real=0.0000095165796700 imag=0.0000000000000000 real=0.0000148505296000 imag=0.0000000000000000 real=0.0000229435209000 imag=0.0000000000000000 real=0.0000350941882000 imag=0.0000000000000000 real=0.0000531456135000 imag=0.0000000000000000 real=0.0000796813547000 imag=0.0000000000000000 Release the DFTI descriptor Free data array TEST PASSED ... For example, for the 1st element it looks like: Input .: real=0.0000014521359600 imag=0.0000000000000000 Output: real=0.0000014521359600 imag=0.0000000000000000

SergeyKostrov · ‎07-26-2013

The test case I used for all verifications is attached.

mecej4 · ‎07-29-2013

FORUM MODERATOR:

Suggestion: The title of this thread has a mis-spelled word (weired should be weird). There is no edit button for the first post of a thread. Since other users may search for threads using the correct spelling, it would be entirely within your editorial prerogatives to correct the spelling and make such searches work properly.

Gennady_F_Intel · ‎08-28-2013

thanks, the title was corrected.