weird results from MKL FFT

weird results from MKL FFT

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? 

12 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

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?

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?

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.

>>...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?

Quote:

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

Quote:

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

>>>>...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.

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

The test case I used for all verifications is attached.

Attachments: 

AttachmentSize
Download test29.cpp5.59 KB

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.

thanks, the title was corrected.

Leave a Comment

Please sign in to add a comment. Not a member? Join today