Search code examples
c++fftintel-mkl

Using Intel MKL to convert real time domain data to complex frequency and phase data and the inverse as well


I'm writing an audio plugin and I would like to convert a series of real numbers representing sample values into a complex array representing frequency and phase. Then I want to be able to do the opposite, turning a complex array of frequency and phases to a series of real numbers, reconstructing the original data.

I'm using Intel MKL and I see only the possibility to perform real->real conversions or complex->complex conversions. Here's the reference I'm using: Intel MKL FFT Functions.

In that reference, there are two overloaded functions: DftiComputeForward and DftiComputeBackward. So, I would like to use these to do the conversions. In the reference for DftiCreateDescriptor, the only options available for DFTI_FORWARD_DOMAIN are DFTI_COMPLEX and DFTI_REAL, but no option for mixed conversions.

Edit I found that the phase is actually equal to atan(imaginary/real). I don't want to mislead anyone getting information from questions.

Edit I just learned that it's best to use atan2(imaginary,real). More information is in the comments.


Solution

  • Every real number is a complex number: ℝ ⊂ ℤ. So going forward from float or double in time domain to complex is trivial. The language does that for you!

    #include <complex>
    #include <iostream>
    int main() {
        double d = 42.0;
        std::complex<double> z = d;
        std::cout << d << " = " << z << '\n';
    }
    

    Output: 42 = (42,0)

    And the C++ standard library also does everything else. It's quite simple, in fact. For once, the library does pretty much what it says on the box.

    Even better: std::complex offers array access. You can reinterpret-cast std::complex<T> to T[2], whether through a reference or a pointer. And thus, std::complex can be "stripped" of its identity and passed into any lower-level API that requires pairs of floats or pairs of doubles.

    The complex frequency domain data can be converted to magnitude and phase, and back, as follows:

    #include <complex>
    #include <iostream>
    int main() {
        std::complex<double> z{0.7071, 0.7071};
        double magnitude = abs(z);
        double phase = arg(z); // in radians
    
        std::cout << z << " ≈ (" << magnitude << "∠" << phase*180.0/M_PI << "°)\n";
    
        std::complex<double> z2 = std::polar(magnitude, phase);
        std::cout << " ≈ " << z2 << '\n';
    }
    

    Output:

    (0.7071,0.7071) ≈ (0.99999∠45°)
     ≈ (0.7071,0.7071)
    

    Once you get the "real" data back, it's not likely that the imaginary part of the time domain data will be zero - it depends on what processing you'll do with the frequency domain data. What you want to convert back is the magnitude of each complex time sample, using the abs function.

    There's stuff in the C++ library that's mind-bogglingly overcomplicated, to the point where you have to have the reference open or you won't remember how to use it. See e.g. the mess known as the random number support. Ugh. But complex number support is relatively sane and even follows the notation used in teaching complex number arithmetic :)