r如何gpu加速 gpu加速fft

软硬件FFT性能测试

  FFT在很多算法中都有应用，M6678 DSP支持软件FFT（调用DSP库），和硬件FFT（有一个独立的FFT硬件加速模块）。

测试条件

• 操作系统 Win11
• CCS 6.2.0
• CGT-Tools 7.4.4
• XDCTools 3.25.5.94
• SYS/BIOS 6.33.6.50
• DSPLIB C66x 3.4.0.4
• MATHLIB C66x 3.1.2.4
• 256kB L2 Cache, 32kB L1P/D Cache
• 堆区放在DDR，待处理数据与结果数据都在堆区开辟；
• 待处理的数据均为单精度浮点数，实部虚部交叉存储；
• 处理器主时钟 1GHz

测试结果

1、 单次运算
  软件FFT以调用FFT函数开始计时，Cache写回之后结束计时。硬件FFT以启动FFT计算开始，收到FFT完成中断后结束计时。

FFT点数	软件FFT	硬件FFT
32	1327 ns	5292 ns
64	1840 ns	5828 ns
128	3078 ns	6372 ns
256	3956 ns	6918 ns

  可以看到单次FFT运算，软件FFT更快，我猜测可能是中断开销比较大，或者是硬件FFT在数据传输的时候有一个固有的开销。

2、连续多次运算
  我目前有一个二维FFT计算的需求，对MxN的二维数据，需要连续对M行数据分别进行N点FFT计算再转置，再进行M点N行FFT计算和转置。以同样的起始和结束条件进行测试，得到以下平均每次FFT计算的时间开销：

测试条件	软件FFT平均每次耗时	硬件FFT平均每次耗时
32次32点	696 ns	180 ns
64次64点	806 ns	163 ns
128次128点	3219 ns	271 ns
256次256点	5113 ns	563 ns

相关测试工程

二维FFT计算分解

  对二维FFT运算的使用，最普适的就是从复数到复数的变换。即使我想对一个二维的实数矩阵进行FFT，也可以将实数转换成复数表示然后做FFT。

  但这样做的效率不太高，结合实际的应用场景，我最近想做二维FFT是为了加速卷积运算，“时域卷积等于频域相乘”。因此FFT的源数据是实数，而且在这种应用场景下IFFT的结果数据也是实数。

  综合了一些考虑之后，我设计了一个用于调用硬件FFT的类CFftManager。它主要支持下面这些功能：

FFT2d_R2R

  有些实数数据我们先验地知道它是中心对称的，那么它做FFT的结果也会是实数，这个在我之前KCF算法推导的第三篇文章中有提到，因为旋转因子是中心对称的，只要实数也满足同样的对称性，那么计算中的虚部就会都抵消只剩下实数。

  “FFT_R2R_T”表示从实数到实数的变换，并且会进行转置。只要连续调用该函数两次就能实现二维的实数到实数的FFT变换。虽然这里命名是“R2R”，但实际上中间的计算结果是有可能有虚数存在的。只有输入数据不光在二维上中心对称，而且在一维上也中心对称结果才只有实数。若输入数据仅仅满足二维的中心对称，第一次FFT之后得到的结果中是有虚数的，然后在完成第二次FFT之后才只剩实数。

  用户给定需要进行FFT的实数数据SRC，结果可以输出到SRC中，也可以存放到另外的DST空间。

  由CFftManager类管理两块临时空间TempR和TempI用于存放FFT的结果，再由DMA将FFT结果搬运到目的空间，并进行转置。TempR空间不仅仅是用来存放FFT结果中的实数，也起到了存放输入的虚数的作用。因此一个完整的二维FFT2d_R2R的计算，第一步需要先将TempR清零。然后连续调用两次FFT_R2R_T即可得到正确结果。

FFT2d_R2C

  由CFftManager类管理的TempR和TempI是连续的两块存储空间。一般的实数到复数的二维FFT的计算过程如下图所示：

  首先计算由一维的实数到复数的FFT，此时输入的数据的实部和虚部是分开存储的，可以在第一次转置的时候顺便将存储方式改为交叉存储。而后需要进行第二次FFT，然后再进行第二次转置。
  为什么要交叉存储？交叉存储能够方便处理器从存储器中取操作数。比如取一个复数，这个复数如果是连续存放的，很可能就只需要一条指令就能够完成取操作数。这个在我的前面的文章里由提到关于复数乘法运算时，如果复数存放方式是[imag, real, imag, …]，那么每次就只需要一条LDDW指令就可以完成取操作数，并且接着就可以进行复数乘法计算。
  但是硬件FFT所支持的复数交叉存放方式是[real, imag, real, …]的形式。为了应对不同的需求我写了两个FFT_C2C_T的函数，它们都是完成FFT复数到复数的变换，只是最后一步的转置有一点点区别。“FFT_C2C_T”输出结果是按照[imag, real, imag, …]的方式存放的，更适合后续需要进行复数乘法的应用场景；“FFT_C2C_TN”结尾的“N”取自Normal，结果是按照[real, imag, real, …]的方式存放的，更加符合一般的用法。

IFFT2d_C2R

  IFFT没有什么特殊的地方，在确定输出一定为实数的情况下，我只是将最后计算的结果中的实数部分进行了提取并转置后存放到指定的DST空间中。

CFftManager类定义

  每次计算二维的FFT都需要从两个维度分别计算一维的FFT，每次都需要等待一维的FFT计算完成，即接收到完成中断才能进行下一步。这部分前后衔接的工作涉及中断，所以不适合放到类里去实现。我把前面提到的计算步骤分别写成几个独立的成员函数，在实际应用过程中进行组合即可。

/**
 * @brief Hardware FFT caculation manager
 * 
 */
class CFftManager{
private:
    CSL_Edma3ChannelHandle m_hCha;
    CSL_Edma3ChannelObj m_ChaObj;
    float *m_pTempR;
    float *m_pTempI;
    float *m_pMid1;
    float *m_pMid2;
    Bool m_bInitDone;
    int m_nCicId;

    uint8_t getPowOfLen(uint16_t nLength);

public:
    CFftManager();
    ~CFftManager();

    /**
     * @brief allocate temp space， TempR and TempI.
     *
     * @param maxWidth max FFT length
     * @param maxHeight max FFT lines
     * @return int
     */
    int init(int maxWidth, int maxHeight);

    /**
     * @brief attach FFT finish interrupt to EDMA channel and enable EDMA interrupt
     * user should setup a Hwi for the event
     * 
     * @param coreId core ID
     * @return int
     */
    int HwInit(int coreId);

    /**
     * @brief do FFT and transpose. 
     * source data is central symmetry real data
     * pSrc and pDst can point to the same space
     * using m_pTempR as input imag data and output imag data
     * 
     * @param pSrc input only real part data
     * @param pDst output only real part data
     * @param nLen FFT points
     * @param nRows number of rows
     * @return int 
     */
    int FFT_R2R_T(float *pSrc, float *pDst, int nLen, int nRows);

    /**
     * @brief do FFT and transpose
     * source data is half size of dst data(dst data is complex)
     * when pSrc and pDst is same, carefully check if source data space is large enough
     * 
     * @param pSrc input only real part data
     * @param pDst output real part and imaginary part data, cross storeage stratagy
     * @param nLen FFT points
     * @param nRows number of rows
     * @return int 
     */
    int FFT_R2C_T(float *pSrc, float *pDst, int nLen, int nRows);

    /**
     * @brief do FFT and transpose
     *  real is in high address, imag is in low address
     * [imag, real, imag, real , ...]
     * 
     * @param pSrc 
     * @param pDst 
     * @param nLen FFT points
     * @param nRows number of rows
     * @return int 
     */
    int FFT_C2C_T(float *pSrc, float *pDst, int nLen, int nRows);

    /**
     * @brief do FFT and transpose normal mode
     * real is in low address, imag is in high address
     * [real, imag, real, imag , ...]
     * 
     * @param pSrc 
     * @param pDst 
     * @param nLen FFT points
     * @param nRows number of rows
     * @return int 
     */
    int FFT_C2C_TN(float *pSrc, float *pDst, int nLen, int nRows);

    int IFFT_C2C_T(float *pSrc, float *pDst, int nLen, int nRows);
    int IFFT_C2R_T(float *pSrc, float *pDst, int nLen, int nRows);

    void ClearInterruptFlag();
    void ClearTempR(size_t nSize);
};

具体应用方式

  创建一个专门用来执行FFT的Task，其他Task可以通过命令队列将要执行的FFT计算发送到队列里，同时需要发送一个启动信号。FFTTask收到启动信号后读取命令队列，然后执行，完成之后发送FFT完成的信号量。

  在调用FFT硬件进行计算的过程中，MainTask仍然可以进行一些其它的计算。MainTask也可以一股脑儿地将要进行计算地任务都发送到命令队列里，处理完其它事情后再回过头来检查FFT完成的信号量的个数是否达到预期，然后进行后续处理，这样可以最大限度地利用硬件FFT的加速功能。

typedef enum FFT2dType
{
    FFT2d_R2R = 0,
    FFT2d_R2C,
    FFT2d_R2C_N,
    IFFT2d_C2R
} FFT2dType;

typedef struct FFTCmd{
    Queue_Elem _elem;
    FFT2dType eType;
    float *pSrc;
    float *pDst;
    int nWidth;
    int nHeight;
} FFTCmd;

Void FFTTskFxn(UArg a0, UArg a1)
{
    FFTCmd *pCmd = NULL;
    int nWidth, nHeight;

    // init in main
    // FftMgr.init(FFT_POINT, FFT_POINT);
    // FftMgr.HwInit(0);

    for (;;){
        Semaphore_pend(hFFTStartSem, BIOS_WAIT_FOREVER);

        pCmd = (FFTCmd *)Queue_dequeue(hFFTQue);
        nWidth = pCmd->nWidth;
        nHeight = pCmd->nHeight;

        switch (pCmd->eType)
        {
        case FFT2d_R2R:
            FftMgr.ClearTempR(sizeof(float) * pCmd->nWidth * pCmd->nHeight);
            FftMgr.FFT_R2R_T(pCmd->pSrc, pCmd->pDst, nWidth, nHeight);
            Semaphore_pend(hFFTHwiSem, BIOS_WAIT_FOREVER);
            FftMgr.FFT_R2R_T(pCmd->pDst, pCmd->pDst, nHeight, nWidth);
            Semaphore_pend(hFFTHwiSem, BIOS_WAIT_FOREVER);
            // TODO check result
            Cache_inv(pCmd->pDst, sizeof(float) * nWidth * nHeight, Cache_Type_ALLD, 1);
            break;
        case FFT2d_R2C:
            FftMgr.FFT_R2C_T(pCmd->pSrc, pCmd->pDst, nWidth, nHeight);
            Semaphore_pend(hFFTHwiSem, BIOS_WAIT_FOREVER);
            FftMgr.FFT_C2C_T(pCmd->pDst, pCmd->pDst, nHeight, nWidth);
            Semaphore_pend(hFFTHwiSem, BIOS_WAIT_FOREVER);

            Cache_inv(pCmd->pDst, sizeof(float) * nWidth * nHeight * 2, Cache_Type_ALLD, 1);
            break;
        case FFT2d_R2C_N:
            FftMgr.FFT_R2C_T(pCmd->pSrc, pCmd->pDst, nWidth, nHeight);
            Semaphore_pend(hFFTHwiSem, BIOS_WAIT_FOREVER);
            FftMgr.FFT_C2C_TN(pCmd->pDst, pCmd->pDst, nHeight, nWidth);
            Semaphore_pend(hFFTHwiSem, BIOS_WAIT_FOREVER);

            Cache_inv(pCmd->pDst, sizeof(float) * nWidth * nHeight * 2, Cache_Type_ALLD, 1);
            break;
        case IFFT2d_C2R:
            FftMgr.IFFT_C2C_T(pCmd->pSrc, pCmd->pSrc, nWidth, nHeight);
            Semaphore_pend(hFFTHwiSem, BIOS_WAIT_FOREVER);
            FftMgr.IFFT_C2R_T(pCmd->pSrc, pCmd->pDst, nHeight, nWidth);
            Semaphore_pend(hFFTHwiSem, BIOS_WAIT_FOREVER);

            Cache_inv(pCmd->pDst, sizeof(float) * nWidth * nHeight, Cache_Type_ALLD, 1);
            break;
        default:
            assert(0); // FFT type error
            break;
        }
        Semaphore_post(hFFTDoneSem);
    }
}