Intel® oneAPI Math Kernel Library Cookbook

ID 758503
Date 3/31/2023
Public

A newer version of this document is available. Customers should click here to go to the newest version.

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.