FFT Example Program

JNichols · ‎09-05-2021

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
    
    !===============================================================================
    ! Copyright 2011-2019 Intel Corporation.
    !
    ! 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
    ! License.
    !===============================================================================

    ! 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

JNichols · ‎09-05-2021

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.

JNichols · ‎09-05-2021

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

JNichols · ‎09-05-2021

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.

JNichols · ‎09-05-2021

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.

ShanmukhS_Intel · ‎09-06-2021

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

JNichols · ‎09-06-2021

Thank you for your response.

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.

JNichols · ‎09-06-2021

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.

JNichols · ‎09-06-2021

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

ShanmukhS_Intel · ‎09-14-2021

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

ShanmukhS_Intel · ‎09-21-2021

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