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

## Fast Fourier Transform

New Contributor I
2,218 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  ******``````

1 Solution
Moderator
1,776 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.

10 Replies
New Contributor I
2,208 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.

Moderator
2,171 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

New Contributor I
2,166 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

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

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.

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

Intel Libraries for oneAPI   Package ID: w_oneAPI_2021.3.0.306
* Other names and brands may be claimed as the property of others.

Intel® C++ Compiler   Package ID: w_oneAPI_2021.3.0.306
* Other names and brands may be claimed as the property of others.

Intel® C++ Compiler Classic   Package ID: w_oneAPI_2021.3.0.306
* 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
* Other names and brands may be claimed as the property of others.

Intel® Fortran Compiler   Package ID: w_oneAPI_2021.3.0.306
* Other names and brands may be claimed as the property of others.

Intel® Inspector 2021.3   605451

Intel® oneAPI DPC++ Compiler   Package ID: w_oneAPI_2021.3.0.306
* 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

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

NuGet Package Manager   5.11.0

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

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.

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

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.

New Contributor I
2,163 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.

``````    !===============================================================================
!
! 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
!
!*****************************************************************************
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
``````

New Contributor I
2,168 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.

Valued Contributor III
2,136 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.

Moderator
1,999 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

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

Valued Contributor III
1,973 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?

Moderator
1,957 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.

Moderator
1,777 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.