如何使用英特尔 MKL FFT？有没有一个简单的 C++ 示例？

Question

我需要执行FFT和反FFT转换。输入将是双精度向量和矩阵。理想情况下，输出应该是std :: complex数组，但我可以使用double _Complex。

我没有找到任何简单的示例，所有英特尔示例都在进行很多操作而没有足够的注释。

我只是想要一个简单的C++示例，它以双精度向量（或矩阵）作为输入，并输出经过FFT转换的结果（最好是带有std :: complex）。

- Baptiste Wicht

---

我猜你指的是这些例子，可能更具体地说是"C接口"的那些。如果是这样，请确保阅读顶部附带的傅里叶变换函数链接。 - SleuthEye

@SleuthEye 是的，我指的就是它们。个人而言，我并不觉得这些臃肿的例子有用，但对于其他人来说可能已经足够了。我想找一些更简单的例子，但看起来我可能只能将就着用它们了。 - Baptiste Wicht

我继续检查官方示例，但其中一些甚至无法编译... 这不是一个真正好的开始方式... - Baptiste Wicht

Answer 1

我最终测试了几个东西，最终得到了这三个函数，它们能够完成我的需求，并且我认为它们是简单的示例。

我对一些输入进行了测试，结果良好。尽管如此，我没有进行广泛的测试。

//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;
}

- Baptiste Wicht

Answer 2

虽然Baptiste的答案可行，但通常在将傅里叶变换应用于实数时，人们希望使用更高效的版本。

对于实数的傅里叶变换F，有以下规律:

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

因此，只需要计算大约一半的值。使用mkl的实数傅里叶变换将导致以下代码：

//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;
}

对于逆变换版本，我们需要知道原始向量的大小-这些信息无法从傅里叶变换中恢复。虽然我们知道，如果原始实向量具有偶数个元素，则傅里叶变换中的最后一个元素是实数，但我们不能从傅里叶变换的最后一个元素是实数推断出原始实向量具有偶数个元素！这就是逆函数略微奇怪的签名的原因：

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;
}

注意：相同的描述符用于正向和反向转换。

- ead