Intel® oneAPI Math Kernel Library Cookbook

Download PDF
ID 758503
Date 11/07/2023
Public
Document Table of Contents
Document Table of Contents x
Intel® oneAPI Math Kernel Library Cookbook
Intel® oneAPI Math Kernel Library Cookbook x
Notices and Disclaimers Getting Technical Support What's New Notational Conventions Related Information Finding an approximate solution to a stationary nonlinear heat equation Factoring general block tridiagonal matrices Solving a system of linear equations with an LU-factored block tridiagonal coefficient matrix Factoring block tridiagonal symmetric positive definite matrices Solving a system of linear equations with a block tridiagonal symmetric positive definite coefficient matrix Computing principal angles between two subspaces Computing principal angles between invariant subspaces of block triangular matrices Evaluating a Fourier integral Using Fast Fourier Transforms for computer tomography image reconstruction Noise filtering in financial market data streams Using the Monte Carlo method for simulating European options pricing Using the Black-Scholes formula for European options pricing Multiple simple random sampling without replacement Using a histospline technique to scale images Speeding up Python* scientific computations Bibliography
Intel® oneAPI Math Kernel Library Cookbook

Evaluating a Fourier integral

Goal

Use a fast Fourier transform (FFT) to numerically evaluate the continuous Fourier transform integral

Solution

Let’s assume that the real-valued function f(x) is zero outside the interval [a, b] and is sampled at N equidistant points xn = a + nT/N, where T = |b - a| and n = 0, 1, ... , N-1. An FFT will be used to evaluate the integral at points ξk = k2π/T, where k = 0, 1, ... , N/2.

Using Intel® oneAPI Math Kernel Library FFT Interface in C/C++

float *f;   // input: f[n] = f(a + n*T/N), n=0...N-1
complex *F; // output: F[k] = F(2*k*PI/T), k=0...N/2
DFTI_DESCRIPTOR_HANDLE h = NULL;
DftiCreateDescriptor(&h,DFTI_SINGLE,DFTI_REAL,1,(MKL_LONG)N);
DftiSetValue(h,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX);
DftiSetValue(h,DFTI_PLACEMENT,DFTI_NOT_INPLACE);
DftiCommitDescriptor(h);
DftiComputeForward(h,f,F);
for (int k = 0; k <= N/2; ++k)
{
  F[k] *= (T/N)*complex( cos(2*PI*a*k/T), -sin(2*PI*a*k/T) );
}

Discussion

The evaluation follows this derivation, based on step-function approximation of the integral:



The sum in the last line is an FFT by definition. When the support of the function f extends symmetrically around zero, that is, [a, b] = [-T/2, T/2], the factor before the sum turns into (T/N)(-1)k.

When the function f is real-valued, F(ξk) = conj(F(ξN-k)). The first N/2 + 1 complex values of the real-to-complex FFT occupy approximately the same memory as the real input, and they suffice to compute the whole result by conjugation. If the FFT computation is configured to perform a real-to-complex transform, it also takes approximately half as much time as a complex-to-complex FFT.

Parent topic: Intel® oneAPI Math Kernel Library Cookbook
Level Two Title