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

Fast Fourier Transform

JNichols
New Contributor I
1,764 Views

I am using the 

call basic_dp_real_dft_1d(XF, num, numR, RC)

routine.  Your sample code in Fortran has the input matrix with a size 7, ie 1 to 7 in the array, but I thought I picked up that the C sample is 0 to 8 as expected.  

We normally stick to arrays with a size of 16384 -  2000 records per second for about 8 seconds.  This gives us 1000 Hz at 1000/8192 steps.  I normally use the Numerical Recipes FFT Solver, but I thought I would give this a go.  

I am looking at a data length of 960001 records, so it is not even close to a power or 2.  The FFT routine returns me a complex array of half the length as expected.   I just threw in the whole array to see what happened, I expected it to throw an error. But it did not. 

I spent an hour and read the MKL FFT manual chapter, I am none the wiser as to how you get away from the power of 2 requirement.  I checked the output  - the first column is the unit step, so the 9000 is about 10 Hz, for a concrete floor - reasonable, the second of 19000 is about 20 Hz, not unexpected and the 28000 is the electrical cables shaking the building.  We see this a lot. I have not looked for the higher order frequencies 120, 180 etc. yet. 

Am I correct, you do not need a power of 2 input array, in this case I can supply 524288 or 2 to the 19th?

Can I assume that the step is 1000/480000?

 TEST PASSED
     1    ********   3.142  ******
     2      77.753   1.777  ******
     3      43.402   1.754  ******
     4      37.835   1.710  ******
     5      19.610   1.479  ******
     6      18.179   1.167  ******
     7      23.374   1.554  ******
     8      12.224   2.323  ******
     9      17.760   0.986  ******
    18      11.204   1.820  ******
  7399      11.094   4.902  ******
  9070      13.475   3.112  ******
  9248      11.965   3.648  ******
 19147      10.214   2.266  ******
 28196      13.586   1.374  ******
148194      10.736   3.721  ******

 

0 Kudos
1 Solution
Ruqiu_C_Intel
Moderator
1,322 Views

This issue is closing and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only. 


View solution in original post

0 Kudos
10 Replies
JNichols
New Contributor I
1,754 Views

If I graph the frequency output, I get the following for the 960000 elements.  The 15 to 20 Hz is about what one could expect, the power at about 58 Hz means the Houston Supply is less than perfect, which does not surprise me as my old electrical clock changes time at a fairly annoying constant rate.  The 120 is double the 60 and quite good, the 180 is high, which is unusual.  Usually they drop below the 60 multiples. 

image_2021-09-13_153558.png

0 Kudos
ShanmukhS_Intel
Moderator
1,717 Views

Hi,


Thanks for reaching out to us.


We are working on your issue. Could you please share us your environment details(OS details, MKL version etc..)


Best Regards,

Shanmukh.SS


0 Kudos
JNichols
New Contributor I
1,712 Views

The dimension of the transform is a positive integer value represented in an integer scalar of Integer data
type. For a one-dimensional transform, the transform length is specified by a positive integer value
represented in an integer scalar of Integer data type.

I found this in your latest stuff on the FFT, when you are discussing lengths.  

Microsoft Visual Studio Community 2019
Version 16.11.2
VisualStudio.16.Release/16.11.2+31624.102
Microsoft .NET Framework
Version 4.8.04161

Installed Version: Community

Visual C++ 2019   00435-60000-00000-AA912
Microsoft Visual C++ 2019

.NET Core Debugging with WSL   1.0
.NET Core Debugging with WSL

ADL Tools Service Provider   1.0
This package contains services used by Data Lake tools

ASA Service Provider   1.0

ASP.NET and Web Tools 2019   16.11.75.64347
ASP.NET and Web Tools 2019

ASP.NET Web Frameworks and Tools 2019   16.11.75.64347
For additional information, visit https://www.asp.net/

Azure App Service Tools v3.0.0   16.11.75.64347
Azure App Service Tools v3.0.0

Azure Data Lake Node   1.0
This package contains the Data Lake integration nodes for Server Explorer.

Azure Data Lake Tools for Visual Studio   2.6.3000.0
Microsoft Azure Data Lake Tools for Visual Studio

Azure Functions and Web Jobs Tools   16.11.75.64347
Azure Functions and Web Jobs Tools

Azure Stream Analytics Tools for Visual Studio   2.6.3000.0
Microsoft Azure Stream Analytics Tools for Visual Studio

C# Tools   3.11.0-4.21403.6+ae1fff344d46976624e68ae17164e0607ab68b10
C# components used in the IDE. Depending on your project type and settings, a different version of the compiler may be used.

Color Theme Designer   1.0
Designer for creating new color themes

Common Azure Tools   1.10
Provides common services for use by Azure Mobile Services and Microsoft Azure Tools.

Cookiecutter   16.11.21196.2
Provides tools for finding, instantiating and customizing templates in cookiecutter format.

Intel Libraries for oneAPI   Package ID: w_oneAPI_2021.3.0.306
Intel Libraries for oneAPI – toolkit version: 2021.3.0, extension version 21.1.0.16, Package ID: w_oneAPI_2021.3.0.306, Copyright © 2019-2021 Intel Corporation. All rights reserved.
* Other names and brands may be claimed as the property of others.

Intel® Advisor 2021.3   609871
Intel® Advisor 2021.3, (build 609871), Copyright © 2009-2021 Intel Corporation. All rights reserved.

Intel® C++ Compiler   Package ID: w_oneAPI_2021.3.0.306
Intel® C++ Compiler – toolkit version: 2021.3.0, extension version 21.0.16, Package ID: w_oneAPI_2021.3.0.306, Copyright © 2002-2021 Intel Corporation. All rights reserved.
* Other names and brands may be claimed as the property of others.

Intel® C++ Compiler Classic   Package ID: w_oneAPI_2021.3.0.306
Intel® C++ Compiler Classic – toolkit version: 2021.3.0, extension version 19.2.7.16, Package ID: w_oneAPI_2021.3.0.306, Copyright © 2002-2021 Intel Corporation. All rights reserved.
* Other names and brands may be claimed as the property of others.

Intel® Distribution for GDB*   1.0
Provides debugging capabilities for compute shaders running on the Intel GPU.

Intel® DPC++ Compatibility Tool   Package ID: w_oneAPI_2021.3.0.315
Intel® DPC++ Compatibility Tool – toolkit version: 2021.3.0, extension version 2021.1.0.16, Package ID: w_oneAPI_2021.3.0.315, Copyright © 2019-2021 Intel Corporation. All rights reserved.
* Other names and brands may be claimed as the property of others.

Intel® Fortran Compiler   Package ID: w_oneAPI_2021.3.0.306
Intel® Fortran Compiler – toolkit version: 2021.3.0, extension version 19.2.0062.16, Package ID: w_oneAPI_2021.3.0.306, Copyright © 2002-2021 Intel Corporation. All rights reserved.
* Other names and brands may be claimed as the property of others.

Intel® Inspector 2021.3   605451
Intel® Inspector 2021.3, (build 605451), Copyright © 2009-2021 Intel Corporation. All rights reserved.

Intel® oneAPI DPC++ Compiler   Package ID: w_oneAPI_2021.3.0.306
Intel® oneAPI DPC++ Compiler – toolkit version: 2021.3.0, extension version 21.1.0.16, Package ID: w_oneAPI_2021.3.0.306, Copyright © 2019-2021 Intel Corporation. All rights reserved.
* Other names and brands may be claimed as the property of others.

Intel® oneAPI Menu & Samples   1.1.410
Intel® oneAPI Visual Studio Menu & Samples Extension

Intel® VTune™ Profiler 2021   618064
Intel® VTune™ Profiler 2021, (build 618064),  Copyright © 2009-2021 Intel Corporation. All rights reserved.

IntelliCode Extension   1.0
IntelliCode Visual Studio Extension Detailed Info

Microsoft Azure HDInsight Azure Node   2.6.3000.0
HDInsight Node under Azure Node

Microsoft Azure Hive Query Language Service   2.6.3000.0
Language service for Hive query

Microsoft Azure Stream Analytics Language Service   2.6.3000.0
Language service for Azure Stream Analytics

Microsoft Azure Stream Analytics Node   1.0
Azure Stream Analytics Node under Azure Node

Microsoft Azure Tools for Visual Studio   2.9
Support for Azure Cloud Services projects

Microsoft Continuous Delivery Tools for Visual Studio   0.4
Simplifying the configuration of Azure DevOps pipelines from within the Visual Studio IDE.

Microsoft JVM Debugger   1.0
Provides support for connecting the Visual Studio debugger to JDWP compatible Java Virtual Machines

Microsoft Library Manager   2.1.113+g422d40002e.RR
Install client-side libraries easily to any web project

Microsoft MI-Based Debugger   1.0
Provides support for connecting Visual Studio to MI compatible debuggers

Microsoft Visual C++ Wizards   1.0
Microsoft Visual C++ Wizards

Microsoft Visual Studio Tools for Containers   1.2
Develop, run, validate your ASP.NET Core applications in the target environment. F5 your application directly into a container with debugging, or CTRL + F5 to edit & refresh your app without having to rebuild the container.

Microsoft Visual Studio VC Package   1.0
Microsoft Visual Studio VC Package

MySQL for Visual Studio   1.2.9
Data design and management tools for MySQL. Copyright (c) 2007, 2019, Oracle and/or its affiliates. All rights reserved.

NuGet Package Manager   5.11.0
NuGet Package Manager in Visual Studio. For more information about NuGet, visit https://docs.nuget.org/

ProjectServicesPackage Extension   1.0
ProjectServicesPackage Visual Studio Extension Detailed Info

Python   16.11.21196.2
Provides IntelliSense, projects, templates, debugging, interactive windows, and other support for Python developers.

Python - Conda support   16.11.21196.2
Conda support for Python projects.

Python - Django support   16.11.21196.2
Provides templates and integration for the Django web framework.

Python - Profiling support   16.11.21196.2
Profiling support for Python projects.

Razor (ASP.NET Core)   16.1.0.2122504+13c05c96ea6bdbe550bd88b0bf6cdddf8cde1725
Provides languages services for ASP.NET Core Razor.

SQL Server Data Tools   16.0.62107.28140
Microsoft SQL Server Data Tools

Test Adapter for Boost.Test   1.0
Enables Visual Studio's testing tools with unit tests written for Boost.Test.  The use terms and Third Party Notices are available in the extension installation directory.

Test Adapter for Google Test   1.0
Enables Visual Studio's testing tools with unit tests written for Google Test.  The use terms and Third Party Notices are available in the extension installation directory.

ToolWindowHostedEditor   1.0
Hosting json editor into a tool window

TypeScript Tools   16.0.30526.2002
TypeScript Tools for Microsoft Visual Studio

Visual Basic Tools   3.11.0-4.21403.6+ae1fff344d46976624e68ae17164e0607ab68b10
Visual Basic components used in the IDE. Depending on your project type and settings, a different version of the compiler may be used.

Visual F# Tools   16.11.0-beta.21322.6+488cc578cafcd261d90d748d8aaa7b8b091232dc
Microsoft Visual F# Tools

Visual Studio Code Debug Adapter Host Package   1.0
Interop layer for hosting Visual Studio Code debug adapters in Visual Studio

Visual Studio Container Tools Extensions   1.0
View, manage, and diagnose containers within Visual Studio.

Visual Studio Tools for CMake   1.0
Visual Studio Tools for CMake

Visual Studio Tools for Containers   1.0
Visual Studio Tools for Containers

 

Screenshot 2021-09-14 104915.png

I think this is all you want.  

The status error messages from the FFT that are predefined do not include an error message for incorrect length. 

0 Kudos
JNichols
New Contributor I
1,709 Views

This is the real sample from your files.  It has been around for a long time as I have older versions in my backup files.  The length of the real array is 7, I thought at one stage I went and had a look at the C version and there was a 0 to 7 array as I expected, I have not rechecked any my memory may be flawed.  

The first time I used this method, I stuck to 1048 as the input array, but I was lazy and just threw in the whole array, as your 7 appears to make it ok?

Your references include the standard such as Brigham, which is definitely 2 to the nth power.  

    !===============================================================================
    ! 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
    !
    !*****************************************************************************
    program 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
    integer, parameter :: WP = selected_real_kind(15,307)

    ! 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"
        call exit(0)
    else
        print *, "TEST FAILED"
        call exit(1)
    end if

999 print '("  Error, status = ",I0)', status
    goto 100


    contains

    ! 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 program basic_dp_real_dft_1d

 

0 Kudos
JNichols
New Contributor I
1,714 Views

The other thing that is an interesting annoyance is the KIND=WP, on the Fortran Forum, the usual is KIND=DP.  when you mix modules you need to fix the problem. 

I always wondered what WP stood for, I know what it means.  

0 Kudos
JohnNichols
Valued Contributor III
1,682 Views

I reran the problem with a 2**19 data set, about 70% of the original and aside from the scaling factor it gives a similar answer. 

0 Kudos
Ruqiu_C_Intel
Moderator
1,545 Views

Hi,


MKL DFT supports both the input length is a power of 2, for example 2^19 = 524288, or not a power of 2, for example 1048.


I tested the sample code basic_dp_real_dft_1d for both input length 1048 and 524288. Both of them are passed as below:


 Example basic_dp_real_dft_1d

 Forward-Backward double-precision real-to-complex out-of-place 1D transform

 Configuration parameters:

 DFTI_PRECISION       = DFTI_DOUBLE

 DFTI_FORWARD_DOMAIN     = DFTI_REAL

 DFTI_DIMENSION       = 1

 DFTI_LENGTHS        = /1048/

 DFTI_PLACEMENT       = DFTI_NOT_INPLACE

 DFTI_CONJUGATE_EVEN_STORAGE = DFTI_COMPLEX_COMPLEX

 Create DFTI descriptor

 Set out-of-place

 Set CCE storage

 Commit DFTI descriptor

 Allocate data arrays

 Initialize data for real-to-complex FFT

 Compute forward transform

 Verify the result

 Check if err is below errthr 0.557E-14

 Verified, maximum error was 0.722E-15

 Initialize data for complex-to-real FFT

 Compute backward transform

 Verify the result

 Check if err is below errthr 0.557E-14

 Verified, maximum error was 0.111E-14


and


 Example basic_dp_real_dft_1d

 Forward-Backward double-precision real-to-complex out-of-place 1D transform

 Configuration parameters:

 DFTI_PRECISION       = DFTI_DOUBLE

 DFTI_FORWARD_DOMAIN     = DFTI_REAL

 DFTI_DIMENSION       = 1

 DFTI_LENGTHS        = /524288/

 DFTI_PLACEMENT       = DFTI_NOT_INPLACE

 DFTI_CONJUGATE_EVEN_STORAGE = DFTI_COMPLEX_COMPLEX

 Create DFTI descriptor

 Set out-of-place

 Set CCE storage

 Commit DFTI descriptor

 Allocate data arrays

 Initialize data for real-to-complex FFT

 Compute forward transform

 Verify the result

 Check if err is below errthr 0.105E-13

 Verified, maximum error was 0.884E-16

 Initialize data for complex-to-real FFT

 Compute backward transform

 Verify the result

 Check if err is below errthr 0.105E-13

 Verified, maximum error was 0.390E-16


​Not sure whether you have read this article below about the CCS output for MKL sample code before

https://software.intel.com/content/www/us/en/develop/articles/unpack-result-of-intel-mkl-fft-to-align-with-matlab.html



0 Kudos
JohnNichols
Valued Contributor III
1,519 Views

By default, the FFT size is the first equal or superior power of 2 of the window size.

So I am assuming and this is a question, that the MKL FFT packs the 1048 array to 2048 with 1000 extra zeros. 

This should be true, is it?

I would have thought this should trigger a warning to say you are packing the data with zeros?

0 Kudos
Ruqiu_C_Intel
Moderator
1,503 Views

No, for 1048 size length which is not power of 2, MKL doesn't packet it to 2048 with extra 1000 zeros. Actually MKL will perform DFT 1048 API.


0 Kudos
Ruqiu_C_Intel
Moderator
1,323 Views

This issue is closing and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only. 


0 Kudos
Reply