瞧瞧Inter的MKL数学库是怎样做的FFT和IFFT

MKL数学库

Intel数学核心函数库（MKL）是一套高度优化、线程安全的数学例程、函数，面向高性能的工程、科学与财务应用。英特尔 MKL 的集群版本包括 ScaLAPACK 与分布式内存快速傅立叶转换，并提供了线性代数 (BLAS、LAPACK 和Sparse Solver)、快速傅立叶转换、矢量数学 (Vector Math) 与随机号码生成器支持。

1、FFT

1.1、创建FFT句柄

DFTI_DESCRIPTOR_HANDLE data_hand_ ;

1.2、创建一个FFT接收标志

MKL_LONG status;

1.3、创建描述符

status = DftiCreateDescriptor( &data_hand_ , precision, forward_domain, dimension, length );

参数1: FFT句柄

参数2: 转换精度，DFTI_SINGLE或DFTI_DOUBLE。

参数3: 转换域，DFTI_COMPLEX或DFTI_REAL。

参数4: 变换的维数

参数5: FFT阶数

1.4、设置参数

status = DftiSetValue(data_hand_ , config_param, config_val )

参数1: FFT句柄

参数2: 配置参数(可设值请参考官网帮助文档)，DFTI_PLACEMENT指输出数据存储位置

参数3: 配置参数的值(可设值请参考官网帮助文档)，DFTI_NOT_INPLACE指输出数据不存储于输入数据内存中。

1.5、确认描述符

status = DftiCommitDescriptor(data_hand_);

函数使描述符为计算做好准备。一旦描述被确认，之前所设置的参数将不可修改！！

1.6、开始FFT计算

status = DftiComputeForward(data_hand_, datain, dataout)

参数1: FFT句柄

参数2: 输入数据

参数3: 输出数据

1.7、释放描述符

status = DftiFreeDescriptor(&data_hand_)

释放内存，以上所有函数操作成功时，返回0。

#include <iostream>
#include "mkl_dfti.h"
#include "mkl.h"

void fft_complex(MKL_Complex8 *datain, MKL_Complex8 *dataout, uint16_t fft_jie)
{
    DFTI_DESCRIPTOR_HANDLE data_hand_;
    MKL_LONG status;
    status = DftiCreateDescriptor(&data_hand_, DFTI_SINGLE, DFTI_COMPLEX, 1, fft_jie); 
    status = DftiSetValue(data_hand_, DFTI_PLACEMENT, DFTI_NOT_INPLACE); 
    status = DftiCommitDescriptor(data_hand_); 
    status = DftiComputeForward(data_hand_, datain, dataout); //计算fft
    status = DftiFreeDescriptor(&data_hand_); 
}
int main()
{
    MKL_Complex8 *datain = (MKL_Complex8 *)mkl_calloc(1024, sizeof(MKL_Complex), 64);
    for (int i=0;i<1024;i++)
    {
        datain[i].real = i+1;
        datain[i].imag = i+1;
    }
    MKL_Complex8 *dataout = (MKL_Complex8 *)mkl_calloc(1024, sizeof( MKL_Complex8), 64);
    uint16_t fft_jie = 1024;
    fft_complex(datain, dataout, fft_jie);
    cout << "FFT输出结果: " << endl;
    for (int i=0;i<1024;i++)
    {
        cout << dataout[i] << "\t";
    }
    return 0;
}

2、IFFT

IFFT与FFT操作流程大致相同,不同之处在于第四步设置完数据存储位置之后还需要设置一下计算比例因子σ；以及二者计算函数不同。

2.1、设置比例因子

status = DftiSetValue(data_hand_, DFTI_BACKWARD_SCALE, 1.0f / fft_jie)

参数1：IFFT句柄

参数2：计算方式，可选DFTI_FORWARD_SCALE和DFTI_BACKWARD_SCALE，前者表示FFT比例因子，后者表示IFFT比例因子。该配置参数值不设置默认为1，若计算IFFT需要设置DFTI_BACKWARD_SCALE为1 / fft_jie。

2.2、IFFT计算

status = DftiComputeBackward(data_hand_, datain, dataout)

参数1：IFFT句柄

参数2：输入数据

参数3：输出数据

void ifft_complex(MKL_Complex8 *datain, MKL_Complex8 *dataout, uint16_t fft_jie)
{
	DFTI_DESCRIPTOR_HANDLE data_hand_;
    MKL_LONG status;
    status = DftiCreateDescriptor(&data_hand_, DFTI_SINGLE, DFTI_COMPLEX, 1, fft_jie); 
    status = DftiSetValue(data_hand_, DFTI_PLACEMENT, DFTI_NOT_INPLACE); 
    status = DftiSetValue(data_hand_, DFTI_BACKWARD_SCALE, 1.0f / fft_jie); //设置ifft比例因子
    status = DftiCommitDescriptor(data_hand_); 
    status = DftiComputeBackward(data_hand_, datain, dataout); //计算ifft
    status = DftiFreeDescriptor(&data_hand_); 
}

int main()
{
    MKL_Complex8 *datain = (MKL_Complex8 *)mkl_calloc(1024, sizeof(MKL_Complex), 64);
    for (int i=0;i<1024;i++)
    {
        datain[i].real = i+1;
        datain[i].imag = i+1;
    }
    MKL_Complex8 *dataout = (MKL_Complex8 *)mkl_calloc(1024, sizeof( MKL_Complex8), 64);
    uint16_t fft_jie = 1024;
    ifft_complex(datain, dataout, fft_jie);
    cout << "IFFT输出结果: " << endl;
    for (int i=0;i<1024;i++)
    {
        cout << dataout[i] << "\t";
    }
    return 0;
}