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

## FFT Example Program New Contributor I
1,365 Views

This morning I needed an FFT routine, normally I use the Numerical Recipes FFT, but I thought, I would give the MKL routine a whirl.

It was easy to turn the sample program into a module, then I started to look at the code, then I started to smile. I could not decide if I had fallen down some Fortran rabbit hole and someone with a sense of humour was playing mind games.

It is an interesting cross between my Fortran class notes from 1978 and modern Fortran.

It provides a well documented example, thank you, but a student of Fortran would have the book out to look for explanation.  I actually had to get my Fortran book out to look at some of the commands, and a sheet of paper to work out some of the equations.

`````` MODULE FFT

integer, parameter :: WP = selected_real_kind(15,307)
contains

!===============================================================================
!
! This software and the related documents are Intel copyrighted  materials,  and
! your use of  them is  governed by the  express license  under which  they were
! provided to you (License).  Unless the License provides otherwise, you may not
! use, modify, copy, publish, distribute,  disclose or transmit this software or
! the related documents without Intel's prior written permission.
!
! This software and the related documents  are provided as  is,  with no express
! or implied  warranties,  other  than those  that are  expressly stated  in the
!===============================================================================

! Content:
! A simple example of double-precision real-to-complex out-of-place 1D
! FFT using Intel(R) MKL DFTI
!
!*****************************************************************************
subroutine basic_dp_real_dft_1d()

use MKL_DFTI, forget => DFTI_DOUBLE, DFTI_DOUBLE => DFTI_DOUBLE_R

! Size of 1D transform
integer, parameter :: N = 7

! Arbitrary harmonic used to test FFT
integer, parameter :: H = 1

! Need double precision

! Execution status
integer :: status = 0, ignored_status
integer i

! Data arrays
real(WP), allocatable :: x_real (:)
complex(WP), allocatable :: x_cmplx (:)

type(DFTI_DESCRIPTOR), POINTER :: hand

hand => null()

print *,"Example basic_dp_real_dft_1d"
print *,"Forward-Backward double-precision real-to-complex",               &
&      " out-of-place 1D transform"
print *,"Configuration parameters:"
print *,"DFTI_PRECISION              = DFTI_DOUBLE"
print *,"DFTI_FORWARD_DOMAIN         = DFTI_REAL"
print *,"DFTI_DIMENSION              = 1"
print '(" DFTI_LENGTHS                = /"I0"/" )', N
print *,"DFTI_PLACEMENT              =  DFTI_NOT_INPLACE"
print *,"DFTI_CONJUGATE_EVEN_STORAGE = DFTI_COMPLEX_COMPLEX"

print *,"Create DFTI descriptor"
status = DftiCreateDescriptor(hand, DFTI_DOUBLE, DFTI_REAL, 1, N)
if (0 /= status) goto 999

print *,"Set out-of-place"
status = DftiSetValue(hand, DFTI_PLACEMENT, DFTI_NOT_INPLACE)
if (0 /= status) goto 999

print *,"Set CCE storage"
status = DftiSetValue(hand, DFTI_CONJUGATE_EVEN_STORAGE,                   &
&                    DFTI_COMPLEX_COMPLEX)
if (0 /= status) goto 999

print *,"Commit DFTI descriptor"
status = DftiCommitDescriptor(hand)
if (0 /= status) goto 999

print *,"Allocate data arrays"
allocate ( x_real(N), STAT = status)
if (0 /= status) goto 999
allocate ( x_cmplx(INT(N/2.0) + 1), STAT = status)
if (0 /= status) goto 999

print *,"Initialize data for real-to-complex FFT"
call init_r(x_real, N, H)

print *,"Compute forward transform"
status = DftiComputeForward(hand, x_real, x_cmplx)
if (0 /= status) goto 999

print *,"Verify the result"
status = verify_c(x_cmplx, N, H)
if (0 /= status) goto 999

print *,"Initialize data for complex-to-real FFT"
call init_c(x_cmplx, N, H)

do 101 i = 1, 4
write(*,*)i,x_cmplx(i)
101 end do
print *,"Compute backward transform"
status = DftiComputeBackward(hand, x_cmplx, x_real)
if (0 /= status) goto 999

print *,"Verify the result"
status = verify_r(x_real, N, H)
if (0 /= status) goto 999

100 continue

print *,"Release the DFTI descriptor"
ignored_status = DftiFreeDescriptor(hand)

if (allocated(x_real) .or. allocated(x_cmplx)) then
print *,"Deallocate data arrays"
endif
if (allocated(x_real))     deallocate(x_real)
if (allocated(x_cmplx))    deallocate(x_cmplx)

if (status == 0) then
print *, "TEST PASSED"
return
else
print *, "TEST FAILED"
call exit(1)
end if

999 print '("  Error, status = ",I0)', status
return
end subroutine basic_dp_real_dft_1d

! Compute mod(K*L,M) accurately
pure real(WP) function moda(k,l,m)
integer, intent(in) :: k,l,m
integer*8 :: k8
k8 = k
moda = real(mod(k8*l,m),WP)
end function moda

! Initialize x to be inverse of unit peak at H
subroutine init_r(x, N, H)
integer N, H
real(WP) :: x(:)

real(WP), parameter:: TWOPI = 6.2831853071795864769_WP
integer k

forall (k=1:N)
x(k) = 2 * cos( TWOPI * moda(H,k-1,N) / N) / N
end forall
if (mod(2*(N-H),N)==0) x(1:N) =  x(1:N) / 2
end subroutine init_r

! Initialize x to be inverse of unit peak at H
subroutine init_c(x, N, H)
integer N, H
complex(WP) :: x(:)

complex(WP), parameter :: ITWOPI = (0.0_WP,6.2831853071795864769_WP)
integer k

forall (k=1:N/2+1)
x(k) = exp( -ITWOPI * moda(H,k-1,N) / N) / N
end forall
end subroutine init_c

! Verify that x(N) is unit peak at H and N-H
integer function verify_c(x, N, H)
integer N, H
complex(WP) :: x(:)

integer k
real(WP) err, errthr, maxerr
complex(WP) res_exp, res_got

! Note, this simple error bound doesn't take into account error of
! input data
errthr = 2.5 * log(real(N,WP)) / log(2.0_WP) * EPSILON(1.0_WP)
print '("  Check if err is below errthr " G10.3)', errthr

maxerr = 0.0_WP
do k = 1, N/2+1
if (mod(k-1-H,N)==0 .or. mod(-k+1-H,N)==0) then
res_exp = 1.0_WP
else
res_exp = 0.0_WP
end if
res_got = x(k)
err = abs(res_got - res_exp)
maxerr = max(err,maxerr)
if (.not.(err < errthr)) then
print '("  x("I0"):"\$)', k
print '(" expected ("G24.17","G24.17"),"\$)', res_exp
print '(" got ("G24.17","G24.17"),"\$)', res_got
print '(" err "G10.3)', err
print *," Verification FAILED"
verify_c = 1
return
end if
end do
print '("  Verified,  maximum error was " G10.3)', maxerr
verify_c = 0
end function verify_c

! Verify that x is unit peak at H
integer function verify_r(x, N, H)
integer N, H
real(WP) :: x(:)

integer k
real(WP) err, errthr, maxerr, res_exp

! Note, this simple error bound doesn't take into account error of
! input data
errthr = 2.5 * log(real(N,WP)) / log(2.0_WP) * EPSILON(1.0_WP)
print '("  Check if err is below errthr " G10.3)', errthr

maxerr = 0.0_WP
do k = 1, N
if (mod(k-1-H, N)==0) then
res_exp = 1.0_WP
else
res_exp = 0.0_WP
end if
err = abs(x(k) - res_exp)
maxerr = max(err,maxerr)
if (.not.(err < errthr)) then
print '("  x("I0"):"\$)', k
print '(" expected "G24.17","\$)', res_exp
print '(" got "G24.17","\$)', x(k)
print '(" err "G10.3)', err
print *," Verification FAILED"
verify_r = 100
return
end if
end do
print '("  Verified,  maximum error was " G10.3)', maxerr
verify_r = 0
end function verify_r

end module
``````

10 Replies New Contributor I
1,357 Views
``do 101 i = 1, (N+1)/2``

At line 101 in the code, the author hard codes in the size of the complex, which was originally set as a variable.   It is a simple fix. New Contributor I
1,342 Views
``````! Note, this simple error bound doesn't take into account error of
! input data
errthr = 2.5 * log(real(N,WP)) / log(2.0_WP) * EPSILON(1.0_WP)
print '("  Check if err is below errthr " G10.3)', errthr``````

The note about input data is interesting as this error means if I have a 15 unit long number, an error of 1 in the last place will generate an error.  I do engineering stuff, we have COV of 20% on the input data. New Contributor I
1,319 Views
``````print *,"Create DFTI descriptor"
status = DftiCreateDescriptor(hand, DFTI_DOUBLE, DFTI_REAL, 1, NB)
if (0 /= status) goto 999``````

I was looking for the status code meanings in the Intel Fortran MKL manual. I assume because this says that if status == 0 continue,  that any other number is error, but what other numbers are used, there is not a list I can find in the manual.

All the status errors push straight to 999, it is then a interesting exercise to use DEBUG to step through the code, an incorrect operation on a complex array, just makes for some interesting program behaviour. New Contributor I
1,297 Views

I found the status checking program.  It took a while to work out this function existed, using it in the example would be nice, it was not easy to find. Moderator
1,261 Views

Hi,

Thanks for reaching out to us.

>>I started to look at the code

Could you please let us know where did you get the sample code.

In addition, For your reference, latest sample for double-precision real-to-complex out-of-place 1D FFT is available under below path.

\oneAPI\mkl\2021.3.0\examples\examples_core_f.zip\f\dft\source\basic_dp_real_dft_1d.f90

>>I found the status checking program. It took a while to work out this function existed, using it in the example would be nice, it was not easy to find.

Thanks for the suggestion. We'll take it into consideration.

Best Regards,

Shanmukh.SS New Contributor I
1,246 Views

I had a copy of the basic_dp_real_dft_1d.f90 from an older version of Parallel Studio. After I started using it, I stumbled across the zip files.

I used the sample to create a working program for damping calculations.  The problem is with damping you need to estimate the lowest effective frequency and with a graph like this as set out in the standard textbooks it is easy: Standard problem

But the real world of real data is not that nice and one has something like this: 500 milliseconds of data

A car passing a sensor on a bridge.  Normally I use the Numerical Recipes FFT, I have since the early 90's, but I decided what the heck, use the MKL.

I just had to take your sample apart, work out how to feed it the correct data and then make sense of the output.  The graph is from DISLIN.

As I worked through the problem, I commented, mainly as I was bored mostly, but it is a way to communicate with the likes of mecej4 and SL, if I make a mistake, they will tell me. usually kindly, but bluntly.

The data length is 1024 and the rate is 2000 Hz, so one has large FFT steps.

hang on this is getting long. New Contributor I
1,239 Views

So the direct output from the FFT has to be turned into something that visually makes sense, so you can determine the FORTRAN to get the real single number you want to know.  So your output turned into a graph FFT

The red is the amplitude, normalized to match the theta range of 0 to 2pi.  For plotting purposes the ATAN2 function needs to be amended to turn the negative range into pi to 2pi.

A really good Russian Statistician developed the method to determine the required numbers from this mess.  After nearly 25 years of using her method, I can say the answer is obvious, and it is trivial to prove if you understand her method.  The results around 400 Hz are ground noise, this is why when people hear an earthquake wave passing they think it is a train.  The ground can make music up to about the Stuttgart Pitch, A above Middle C, but then it breaks down.

So the MKL  FFT worked a treat, the example is ok, it could provide a bit more output, the status codes are a bit skinny, but they are there. So it was fun.

In the words of the great British Statistician FN David, (a lady), when asked how to bomb effectively in WW2. Her mathematical response, which I have read can be summarized as - try not to stick more than one bomb at one location.  I paid the UK government to turn her report into a pdf, got a copy, lost it, and then when I went back for another copy, as I was told if I payed to digitize it, it would be free after that, found out it was reclassified.  Oxbridge did not have many female professors before the war.

So the real question is:  Are these forums for people who just want to know how to install the methods or are they also for people who use the methods. New Contributor I
1,237 Views

The story of FN David is irrelevant, but interesting, and the last question is rhetorical. Moderator
1,106 Views

Hi,

If you want to talk to any specific person/users, you could add "@" in front of their name so that they would be notified. Do you need any help from us on this?

>>Are these forums for people who just want to know how to install the methods or are they also for people who use the methods.

In this Intel forum, any body is free to talk and hear about their problems which are related to Intel products.

Best Regards,

Shanmukh.SS Moderator
1,063 Views

Hi,

We assume that your issue is resolved. If you need any additional information, please post a new question as this thread will no longer be monitored by Intel.

Best Regards,

Shanmukh.SS 