# [FORTRAN] FFT of delta function, issue w/ MKL & Intel compiler

1. Aug 14, 2015

### avikarto

I am trying to program something using a backwards FFT, and am attempting to feed it a delta function as a test condition since this result is known. However, my results are nonsense compared to what is expected.

It should be the case that if we have $$F(\omega)=\delta(\omega-1.5)=\frac{1}{\sqrt{2 \pi}}\int f(t)\:e^{-i \omega t}$$ the result would look like $$f(t)=\frac{Cos(1.5 t)-i Sin(1.5)t}{\sqrt{2 \pi}}$$ However, my method returns the same exact thing for the real and imag parts of f(t); these parts oscillate with frequency much greater than the expected 1.5 and are immediately on top of each other when plotted.

Numerically, i tried inputing two types of "delta function" arrays: one that was zero everywhere on the grid except at w=1.5, where it was equal to 1/dw (normalization), and another where each grid point was calculated as in a Gaussian, $$F(\omega)=\frac{1}{\sigma \sqrt{2\pi}} exp\left(-\frac{(\omega-1.5)^2}{2 \sigma^2}\right),$$ where I let sigma be small for a narrow spike. Both of these input types give the same result. I even tried a much finer grid, but the result was again the same. I am left to believe that the method is somehow flawed, but it took it almost exactly off of Intel's MKL help page with example FFTs (link), using the 1D out-of-place complex-to-complex transform.

My code follows. Any assistance in locating the issue would be greatly appreciated. Thanks.

Code (Fortran):
use MKL_DFTI
function FFTB(nListIn,nLength)
integer status,nLength
type(DFTI_DESCRIPTOR), POINTER :: My_Desc_Handle
complex*16 nListIn(-nLength:nLength),FFTB(-nLength:nLength)

My_Desc_Handle => null()

status=DftiCreateDescriptor(My_Desc_Handle,DFTI_DOUBLE,DFTI_COMPLEX,1,2*nLength)
status=DftiSetValue(My_Desc_Handle,DFTI_PLACEMENT,DFTI_NOT_INPLACE)
status=DftiCommitDescriptor(My_Desc_Handle)
if (status .ne. 0) then
if (.not. DftiErrorClass(status,DFTI_NO_ERROR)) then
print *, 'Error: ', DftiErrorMessage(status)
endif
endif
status=DftiComputeBackward(My_Desc_Handle,nListIn,FFTB)
status=DftiFreeDescriptor(My_Desc_Handle)
return
end function FFTB

2. Aug 15, 2015

### Staff: Mentor

Why not test it with variations of sin cos functionw? This will test if it can pick out the frequencies of the function.

3. Aug 16, 2015

### avikarto

Because the function I am using is for the backwards transform - starting with F(w) and ending with f(t).

4. Aug 16, 2015

### Staff: Mentor

In your code, the identifier FFTB is the name of a function and also the name of a complex*16 array. It has been many years since I've written any Fortran code, so I'm not sure whether the compiler maintains separate namespaces for functions and arrays.

I would change two lines of code, and see if that makes a difference. I would change this:
Code (Fortran):
complex*16 nListIn(-nLength:nLength),FFTB(-nLength:nLength)
to this:
Code (Fortran):
complex*16 nListIn(-nLength:nLength), nListOut(-nLength:nLength)

And I would also change this line:
Code (Fortran):
status=DftiComputeBackward(My_Desc_Handle,nListIn,FFTB)
to this:
Code (Fortran):
status=DftiComputeBackward(My_Desc_Handle,nListIn, nListOut)
Also, you need to declare your function properly to inform the compiler that it will return an array. This web page (http://www.phy.ornl.gov/csep/pl/node18.html) shows an example of a function definition in Fortran-90 code that returns an array of reals.

5. Aug 16, 2015

### avikarto

Thanks for the idea, but the way that I have defined the function is valid in FORTRAN. Naming the function the same thing as the output variable allows the final line, "return" to return said variable. Then, when calling the function, one can say nListOut=FFTB(input_stuff). This is not the issue.

6. Aug 16, 2015

### Staff: Mentor

You don't show how nListIn is being initialized. Possibly there is a problem with that.

Also, you have what I believe is a typo in post 1 for f(t)
I'm guessing that the argument for the sine function should be 1.5 t, not 1.5, with t multiplying Sin(1.5).

7. Aug 16, 2015

### avikarto

Yep that's a mistype in the OP, it should obviously be Sin(1.5 t).

nListIn is generated as follows:

Code (Fortran):
! for the gaussian
real*8, parameter :: pi=4.d0*atan(1.d0)
real*8, parameter :: sig=1.d-5
real*8, parameter :: const=1.d0/(sig*sqrt(2*pi))

! for the FFT
real*8, parameter :: nmax=3.14*30, dn=1.0d-3
integer, parameter :: nLength=nmax/dn
complex*16 nListIn(-nLength:nLength), nListOut(-nLength:nLength)
integer j

do j=-nLength,nLength
! true delta function
!       if (dn*j.eq.1.5) then
!          nListIn(j)=(1000.d0,0.0d0)
!       else
!          nListIn(j)=(0.0d0,0.0d0)
!       end if

! narrow gaussian
nListIn(j)=(const*exp(-(j*dn-1.5d0)**2/(2.d0*sig**2)),0.d0)

end do

nListOut=FFTB(nListIn,nLength)

8. Aug 16, 2015

### Staff: Mentor

There's nothing that pops out at me. The only thing I can think of is to consider using a smaller value for sigma. 10-5 isn't really all that small, so you're getting a considerably wider spike than would be the case for the actual delta function. You might try a value of 10-10 or possibly even 10-15. With 64-bit reals, the granularity is down around 2.2 X 10-16 or so.

9. Aug 17, 2015

### avikarto

I tried a smaller sigma, with no change in the results. This is expected, because the actual delta input (effectively a super small sigma gaussian) showed these results as well.

I wonder if there is just something about a discrete delta-type array that makes it transform differently than in the known continuous case. I tried extending the length of the input array to many times larger than shown above to simulate more continuous data, but that didn't help either. I don't know what else it could be at this point.