Problem calling the MKL FFT functions DftiComputeForward and DftiComputeBackward

[CSA_O]

Hello,

My ultimate goal is to convert a Matlab code that has calls to the ifft function into C.

Firstly, I'm just trying to see if I can call the MKL correctly and even with a simple test I can't seem to be able to make it work.

My idea is to create an array [1, 2, 3, 4, 5, 6, 7, 8] pass it through the FFT, then its backward function and I should get back my simple array.

Here is the code I wrote:

#define _USE_MATH_DEFINES

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

#include "mkl_dfti.h"

#define SIZEFFT 8

void verySimpleTest()
{
    MKL_LONG status;

    //Create input array
    float inputRe[SIZEFFT];
    for (int i = 0; i < SIZEFFT; i++) 
    {
        inputRe[i] = i+1;
    }

    //Create output array
    float _Complex outputComp[SIZEFFT];

    // Create a descriptor handle
    DFTI_DESCRIPTOR_HANDLE desc;

    // Create the descriptor
    status = DftiCreateDescriptor(&desc, DFTI_DOUBLE, DFTI_COMPLEX, 1, SIZEFFT);

    // Initialize the placement to not in place
    status = DftiSetValue(desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE);

    // Commit the descriptor
    status = DftiCommitDescriptor(desc);

    // Perform the forward Fourier transform
    status = DftiComputeForward(desc, inputRe, outputComp);

    // Erase data from the input array to make sure the backward function writes in it
    for (int i = 0; i < SIZEFFT; i++) 
    {
        inputRe[i] = 0;
    }
    
    // Perform the backward Fourier transform
    status = DftiComputeBackward(desc, outputComp, inputRe);

    // Print the data
    for (int i = 0; i < SIZEFFT; i++) 
    {
        printf("inputRe[%d] = %f\n", i, inputRe[i]);
    }

    // Free the descriptor
    DftiFreeDescriptor(&desc);
}

int main(void)
{
    verySimpleTest();

    return 0;
}

The result I get is:

inputRe[0] = 1.000000
inputRe[1] = 2.750000
inputRe[2] = 3.000000
inputRe[3] = 5.500000
inputRe[4] = 5.000000
inputRe[5] = 7.500000
inputRe[6] = 7.000000
inputRe[7] = 11.000000

What am I missing to make it work?

Thank you.

[CSA_O · Accepted solution]

I found the problem myself. If anyone ever find that question and has the same issue, the problem was that the FFT is symetrical and one needs to complete the output array as the DftiComputeForward function only fills half the result .

void verySimpleTest()
{
    MKL_LONG status;

    //Create input array
    double inputRe[SIZEFFT];
    for (int i = 0; i < SIZEFFT; i++) 
    {
        inputRe[i] = i+1;
    }

    //Create output array
    double _Complex outputComp[SIZEFFT];

    // Create a descriptor handle
    DFTI_DESCRIPTOR_HANDLE desc;

    // Create the descriptor
    status = DftiCreateDescriptor(&desc, DFTI_DOUBLE, DFTI_REAL, 1, SIZEFFT);

    // Initialize the placement to not in place
    status = DftiSetValue(desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE);

    status = DftiSetValue(desc, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);

    // Commit the descriptor
    status = DftiCommitDescriptor(desc);

    // Perform the forward Fourier transform
    status = DftiComputeForward(desc, inputRe, outputComp);

    // Erase data from the input array to make sure the backward function writes in it
    for (int i = 0; i < SIZEFFT; i++) 
    {
        inputRe[i] = 0;
    }
    
    // Complete the output with its symetry
    for (int i = SIZEFFT/2+1; i < SIZEFFT; i++)
    {
        __real__ outputComp[i] = __real__ outputComp[SIZEFFT-i];
        __imag__ outputComp[i] = (-1) * __imag__ outputComp[SIZEFFT-i];
    }

    // Perform the backward Fourier transform
    status = DftiComputeBackward(desc, outputComp, inputRe);

    // Print the data
    for (int i = 0; i < SIZEFFT; i++) 
    {
        printf("inputRe[%d] = %f\n", i, inputRe[i]);
    }

    for (int i = 0; i < SIZEFFT; i++) 
    {
        printf("outputComp[%d] = %f%+fi\n", i, __real__ outputComp[i], __imag__ outputComp[i]);
    }

    // Free the descriptor
    DftiFreeDescriptor(&desc);
}

[IntelSupport]

Hi,

Thank you for posting in Intel Communities.

We're glad to hear that the issue was resolved. If you ever have any further queries or concerns, please raise a new thread. Thank you.

Regards,

Jilani