Is there any simple C++ example on how to use Intel MKL FFT?

Question

I need to perform FFT and Inverse-FFT transformations. The input would be vector and matrices of double. Ideally, the output should be an array of std::complex but I can live with double _Complex.

I haven't found any simple example, all Intel examples are doing a lot of things at once without enough comments.

I just want a simple example in C++ taking a vector (or a matrix) of double as input and outputting the FFT-transformed result (ideally with std::complex).

Answer 1

I ended up testing several things and I finally ended up with this three functions that do what I want and that I considered simple examples.

I tested it against some inputs and I had the good results. I haven't done extensive testing though.

//Note after each operation status should be 0 on success 

std::vector<std::complex<float>> fft_complex(std::vector<std::complex<float>>& in){
    std::vector<std::complex<float>> out(in.size());

    DFTI_DESCRIPTOR_HANDLE descriptor;
    MKL_LONG status;

    status = DftiCreateDescriptor(&descriptor, DFTI_SINGLE, DFTI_COMPLEX, 1, in.size()); //Specify size and precision
    status = DftiSetValue(descriptor, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //Out of place FFT
    status = DftiCommitDescriptor(descriptor); //Finalize the descriptor
    status = DftiComputeForward(descriptor, in.data(), out.data()); //Compute the Forward FFT
    status = DftiFreeDescriptor(&descriptor); //Free the descriptor

    return out;
}

std::vector<std::complex<float>> fft_real(std::vector<float>& in_real){
    std::vector<std::complex<float>> in(in_real.size());

    std::copy(in_real.begin(), in_real.end(), in.begin());

    return fft_complex(in);
}

std::vector<float> ifft(std::vector<std::complex<float>>& in){
    std::vector<std::complex<float>> out(in.size());

    DFTI_DESCRIPTOR_HANDLE descriptor;
    MKL_LONG status;

    status = DftiCreateDescriptor(&descriptor, DFTI_SINGLE, DFTI_COMPLEX, 1, in.size()); //Specify size and precision
    status = DftiSetValue(descriptor, DFTI_PLACEMENT, DFTI_NOT_INPLACE); //Out of place FFT
    status = DftiSetValue(descriptor, DFTI_BACKWARD_SCALE, 1.0f / in.size()); //Scale down the output
    status = DftiCommitDescriptor(descriptor); //Finalize the descriptor
    status = DftiComputeBackward(descriptor, in.data(), out.data()); //Compute the Forward FFT
    status = DftiFreeDescriptor(&descriptor); //Free the descriptor

    std::vector<float> output(out.size());

    for(std::size_t i = 0; i < out.size(); ++i){
        output[i] = out[i].real();
    }

    return output;
}

Answer 2

While Baptiste's answer works, usually one would like to use a more efficient version when applying fourier transformation to real values.

For the fourier transformation F of real values, the following holds:

F(k) = conj(F(-k))

and thus only about half of the values must be calculated. Using mkl's real fourier transformation leads to the following code:

//helper function for fft and ifft:
DFTI_DESCRIPTOR* create_descriptor(MKL_LONG length) {
    DFTI_DESCRIPTOR* handle = nullptr;
    // using DFTI_DOUBLE for double precision
    // using DFTI_REAL for using the real version
    bool valid = (DFTI_NO_ERROR == DftiCreateDescriptor(&handle, DFTI_DOUBLE, DFTI_REAL, 1, length)) &&
        // the result should not be inplace:
        (DFTI_NO_ERROR == DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)) &&
        // make clear that the result should be a vector of complex:
        (DFTI_NO_ERROR == DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX));
        // chosen normalization is fft(constant)[0] = constant:
        (DFTI_NO_ERROR == DftiSetValue(handle, DFTI_FORWARD_SCALE, 1. / length)) &&
        (DFTI_NO_ERROR == DftiCommitDescriptor(handle));
    if (!valid) {
        DftiFreeDescriptor(&handle);
        return nullptr; //nullptr means error
    }
    return handle;
}

std::vector<std::complex<double>> real_fft(std::vector<double>& in) {
    size_t out_size = in.size() / 2 + 1; //so many complex numbers needed
    std::vector<std::complex<double>> result(out_size);
    DFTI_DESCRIPTOR* handle = create_descriptor(static_cast<MKL_LONG>(in.size()));
    bool valid = handle &&
        (DFTI_NO_ERROR == DftiComputeForward(handle, in.data(), result.data()));
    if (handle) {
        valid &= (DFTI_NO_ERROR == DftiFreeDescriptor(&handle));
    }
    if (!valid) {
        result.clear(); //empty vector -> error
    }
    return result;
}

For the inverse version we need to know the size of the original vector - this information cannot be restored from the fourier-transform. While we know, that if the original real vector had even number of elements, last element in the fourier transform is real, we cannot follow from last element of the fourier transform being real, that the original real vector had an even number of elements! That is the reason for a slightly strange signature of the inverse function:

std::vector<double> real_fft(std::vector<std::complex<double>> & in, size_t original_size) {
    size_t expected_size = original_size / 2 + 1;
    if (expected_size != in.size()) {
        return {};// empty vector -> error
    }
    std::vector<double> result(original_size);
    DFTI_DESCRIPTOR* handle = create_descriptor(static_cast<MKL_LONG>(original_size));
    bool valid = handle &&
        (DFTI_NO_ERROR == DftiComputeBackward(handle, in.data(), result.data()));
    if (handle) {
        valid &= (DFTI_NO_ERROR == DftiFreeDescriptor(&handle));
    }
    if (!valid) {
        result.clear(); //empty vector -> error
    }
    return result;
}

Note: the same descriptor is used for forward and backward transformations.