英特尔 oneAPI 黑客松 —— 使用oneMKL加速FFT

目录

问题背景

题目描述

运行环境

代码实现

心得与感悟

---

问题背景

傅里叶变换
傅里叶变换（Fourier Transform，简称FT）是⼀种线性积分变换，⽤于函数（应用上称作“信号”）在时域和频域之间的变换。因其基本思想首先由法国学者约瑟夫·傅里叶系统地提出，所以以其名字来命名以示纪念。（连续）傅里叶变换的定义如下：

其中,  为原函数,  是原函数的傅里叶变换, 表示时间，其定义域为时域， 表示频率，其定义域为频域,  是虚数单位。通常情况下,  是一个实函数, 而 则是一个复数值函数。

离散傅里叶变换
离散傅里叶变换（Discrete Fourier Transform，简称DFT），是傅里叶变换在时域和频域上都离散的形式，将信号的时域采样变换为频域采样。对于 点序列 , 它的傅里叶变换（DFT）是一个同样长度为 的序列:

其中 为自然底数, 是虚数单位。通常以符号 表示这一变换, 即:

快速傅里叶变换
快速傅里叶变换（Fast Fourier Transform, FFT），是快速计算序列的离散傅里叶变换（DFT）或其逆变换的方法。FFT能够将计算DFT的复杂度从降低到，其中为数据大小。

题目描述

使用 oneMKL 工具，对 FFT 算法进行加速与优化。

1. 下载和并安装配置 oneMKL。
2. 调用 oneMKL 相应 API 函数， 产生 2048 * 2048 个 随机单精度实数。
3. 根据 2 产生的随机数据作为输入，实现两维 Real to complex FFT 参考代码。
4. 根据 2 产生的随机数据作为输入， 调用 oneMKL API 计算两维 Real to complex FFT。
5. 结果正确性验证，对 3 和 4 计算的两维 FFT 输出数据进行全数据比对（允许适当精度误 差）， 输出 “结果正确”或“结果不正确”信息。
6. 平均性能数据比对（比如运行 1000 次），输出 FFT 参考代码平均运行时间和 oneMKL FFT 平均运行时间。

运行环境

从Intel官网下载OneAPI base以安装OneMKL，使用visual studio2022进行项目配置，在16GB RAM、i5-8300H的笔记本上运行。

代码实现

按照题目表述中的要求，首先生成2048*2048随机数：

VSLStreamStatePtr stream;
vslNewStream(&stream, VSL_BRNG_MT19937, 1);
vsRngUniform(VSL_RNG_METHOD_UNIFORM_STD, stream, rows * columns, inputData, 0.0f, 1.0f);

之后基于fftwf_plan类，实现fftw3的fft计算，使用chrono计算运行时间：

fftwf_complex* fftwResult = (fftwf_complex*)fftwf_malloc((rows / 2 + 1) * columns * sizeof(fftwf_complex));
fftwf_plan fftwPlan = fftwf_plan_dft_r2c_2d(rows, columns, inputData, fftwResult, FFTW_ESTIMATE);
auto fftwStartTime = std::chrono::high_resolution_clock::now();
fftwf_execute(fftwPlan);
auto fftwEndTime = std::chrono::high_resolution_clock::now();
totalFFTWTime += std::chrono::duration<double, std::milli>(fftwEndTime - fftwStartTime).count();

创建OneMKL handle，完成OneMKL计算：

DFTI_DESCRIPTOR_HANDLE handle = NULL;
DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 2, strides0);
DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
DftiSetValue(handle, DFTI_INPUT_STRIDES, strides1);
DftiSetValue(handle, DFTI_OUTPUT_STRIDES, strides2);
DftiCommitDescriptor(handle);
auto mklStartTime = std::chrono::high_resolution_clock::now();
DftiComputeForward(handle, inputData, mklResult);
auto mklEndTime = std::chrono::high_resolution_clock::now();
totalMKLTime += std::chrono::duration<double, std::milli>(mklEndTime - mklStartTime).count();

最后检查fftw3和One MKL的运行结果是否一致，以验证准确性：

for (int j = 0; j < size; ++j) {
     if (fabs(mklResult[j].real - fftwResult[j][0]) > 1e-6 || fabs(mklResult[j].imag - fftwResult[j][1]) > 1e-6) {
          flag = true;
          break;
     }
}

最终测试结果见下图：

心得与感悟

OneMKL提供了简洁的API接口和详细的文档说明，使我能够快速上手并充分利用其功能。OneMKL提供的优化的FFT实现，能够快速高效地处理大规模数据集。通过利用多核处理器和向量化指令集，OneMKL能够充分发挥硬件的潜力，实现快速的傅里叶变换。这对于需要频繁进行信号处理、图像处理、频谱分析等应用而言至关重要。