• 中文
  • 3.9.1

    • Fast Fourier Transform

    快速傅里叶变换（Fast Fourier Transform，以下简称FFT）是离散傅里叶变换（Discrete Fourier Transform，以下简称DFT）的快速算法，是时域一频域变换分析中最基本的方法之一。该算法可以实现数据从时域到频域之间的转换，为数据提供了另一维度的处理方法支持。在工程实践中，DFT算法由于计算量过大，无法在实际中广泛实现，FFT算法通过改进DFT的计算方法，将计算量优化到了可以实践的量级，离散傅里叶变换才真正在工程领域广泛应用了起来。

    算子效果

    时域输入数据参数频域输出数据
    imagepSize = HB_FFT_POINT_SIZE_16image

    原理

    FFT是基于DFT算法的改进方法，在计算结果上与DFT等效，优化了计算过程。若F(n)为f(n)的离散傅里叶变换，DFT计算公式表达如下：

    F(n)=∑k=0N−1f(k)WNkn,n=0,1,⋯ ,N−1F(n)=\sum_{k=0}^{N-1}f(k)W_N^{kn},\qquad n=0,1,\cdots,N-1

    其中，旋转因子为 WNkn=e−j2πNknW_N^{kn}=e^{-j\frac{2\pi}{N}kn} ，N为输入序列长度。

    IDFT为DFT的逆运算，其计算公式表达如下：

    f(n)=1N∑k=0N−1F(k)WN−kn,n=0,1,⋯ ,N−1f(n)=\frac{1}{N}\sum_{k=0}^{N-1}F(k)W_N^{-kn},\qquad n=0,1,\cdots,N-1

    另一个基础公式为欧拉公式：

    eix=cos(x)+i∗sin(x)e^{ix}=cos(x)+i*sin(x)

    由欧拉公式可以推导出旋转因子具有以下特性：

    WN0=e−j2πN0=e−j0=cos(0)+jsin(0)=1W_N^{0}=e^{-j\frac{2\pi}{N}0}=e^{-j0}=cos(0)+jsin(0)=1 WNkN2=(cos(−π)+jsin(−π))k=(−1)kW_N^{k\frac{N}{2}}=(cos(-\pi)+jsin(-\pi))^k=(-1)^k WN2kn=ej(2πN2)kn=WN2knW_N^{2kn}=e^{j(\frac{2\pi}{\frac{N}{2}})kn}=W_{\frac{N}{2}}^{kn}

    对于N点DFT公式，可以按照奇偶顺序将公式分解如下：

    F(n)=∑k=0N−1f(k)WNkn=E+OF(n)=\sum_{k=0}^{N-1}f(k)W_N^{kn}=E+O E=∑k=0N2−1f(2k)e−j2πknN2E=\sum_{k=0}^{\frac{N}{2}-1}f(2k)e^{\frac{-j2\pi kn}{\frac{N}{2}}} O=WNn∑k=0N2−1f(2k+1)e−j2πknN2O=W_N^n\sum_{k=0}^{\frac{N}{2}-1}f(2k+1)e^{\frac{-j2\pi kn}{\frac{N}{2}}}

    分解效果如图，将整体的计算任务逐层二分。

    HPL5FFT_split

    由上可知，N为2的幂次时，FFT任务都可以被分解为单点的计算任务。

    由欧拉公式， WNnW_N^{n} 可以写成如下形式，可以看出 WNnW_N^{n} 具有随N相关的周期性和对称性。

    WNn=e−j2πNn=cos(−2πNn)+i∗sin(−2πNn)=cos(2πNn)−i∗sin(2πNn)W_N^{n}=e^{-j\frac{2\pi}{N}n} =cos(-\frac{2\pi}{N}n)+i*sin(-\frac{2\pi}{N}n) =cos(\frac{2\pi}{N}n)-i*sin(\frac{2\pi}{N}n)

    根据此性质就可以得到 WNn=WNN2+nW_N^{n}=W_N^{\frac{N}{2}+n} ，其中 nn 取值范围是 [0,N/2)[0, N/2) 。

    以N为8的FFT运算为例，具体分解过程如下：

    1. 将8点拆解为4点
    HPL5FFT_bufferfly1
    1. 4点拆解为2点
    HPL5FFT_bufferfly2
    1. 计算单点FFT
    HPL5FFT_bufferfly3

    其中，节点两两交错的计算过程被称为蝶形运算，形式如下：

    HPL5FFT_bufferfly

    具体计算公式为：

    FN/4′(0)=fN(0)+WN/40∗fN(4)F'_{N/4}(0) = f_N(0)+W^0_{N/4}*f_N(4) FN/4′(1)=fN(0)−WN/40∗fN(4)F'_{N/4}(1) = f_N(0)-W^0_{N/4}*f_N(4)

    上图中可以看出，F(n)和f(n)之间的n存在比特反的关系。即FFT运算中，输入输出参数之间的索引值是互相的比特相反数，为优化计算效率提供了更好的支持。

    IFFT的计算过程和FFT类似，根据IDFT公式推导即可。

    API接口

    int32_t hbFFT1D(hbUCPTaskHandle_t *taskHandle, hbHPLImaginaryData *dst,
                hbHPLImaginaryData const *src, hbFFTParam const *param);

    详细接口信息请查看 hbFFT1D。

    使用方法

    // Include the header
#include "hobot/hb_ucp.h"
#include "hobot/hpl/hb_hpl.h"
#include "hobot/hpl/hb_fft.h"

// init, allocate memory for data
src_length = 1024 * 1024 * 5;
hbUCPMallocCached(&src_re_mem, src_length, 0);
hbUCPMallocCached(&src_im_mem, src_length, 0);
hbHPLImaginaryData src;
src.realDataVirAddr = src_re_mem.virAddr;
src.realDataPhyAddr = src_re_mem.phyAddr;
src.imDataVirAddr = src_im_mem.virAddr;
src.imDataPhyAddr = src_im_mem.phyAddr;
src.numDimensionSize = 1;
src.dataType = HB_HPL_DATA_TYPE_I16;
src.imFormat = HB_IM_FORMAT_SEPARATE;
src.dimensionSize[0] = 16 * 11;

hbUCPMallocCached(&dst_re_mem, src_length, 0);
hbUCPMallocCached(&dst_im_mem, src_length, 0);
hbHPLImaginaryData dst;
dst.realDataVirAddr = dst_re_mem.virAddr;
dst.realDataPhyAddr = dst_re_mem.phyAddr;
dst.imDataVirAddr = dst_im_mem.virAddr;
dst.imDataPhyAddr = dst_im_mem.phyAddr;
dst.numDimensionSize = 1;
dst.dataType = HB_HPL_DATA_TYPE_I16;
dst.imFormat = HB_IM_FORMAT_SEPARATE;
dst.dimensionSize[0] = 16 * 11;

// init param
hbFFTParam param;
param.pSize = HB_FFT_POINT_SIZE_16;

// init task handle and schedule param
hbUCPTaskHandle_t task_handle{nullptr};
hbUCPSchedParam sched_param;
HB_UCP_INITIALIZE_SCHED_PARAM(&sched_param);
sched_param.backend = HB_UCP_DSP_CORE_0;

// create task
hbFFT1D(&task_handle, &dst, &src, &param);

// submit task
hbUCPSubmitTask(task_handle, &sched_param);

// wait for task done 
hbUCPWaitTaskDone(task_handle, 0);

// release task handle
hbUCPReleaseTask(task_handle);

// release memory
hbUCPFree(&src_re_mem);
hbUCPFree(&src_im_mem);
hbUCPFree(&dst_re_mem);
hbUCPFree(&dst_im_mem);