Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

weird results from MKL FFT

MH_A_
Beginner
688 Views

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? 








0 Kudos
11 Replies
Gennady_F_Intel
Moderator
688 Views

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?

0 Kudos
SergeyKostrov
Valued Contributor II
688 Views
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?
0 Kudos
MH_A_
Beginner
688 Views

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.

0 Kudos
SergeyKostrov
Valued Contributor II
688 Views
>>...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?
0 Kudos
Dmitry_B_Intel
Employee
688 Views

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

0 Kudos
Dmitry_B_Intel
Employee
688 Views

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

0 Kudos
SergeyKostrov
Valued Contributor II
688 Views
>>>>...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.
0 Kudos
SergeyKostrov
Valued Contributor II
688 Views
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
0 Kudos
SergeyKostrov
Valued Contributor II
688 Views
The test case I used for all verifications is attached.
0 Kudos
mecej4
Honored Contributor III
688 Views

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.

0 Kudos
Gennady_F_Intel
Moderator
688 Views

thanks, the title was corrected.

0 Kudos
Reply