CN102340296A - 一种基于gpu的高阶数字fir滤波器频域并行处理实现方法 - Google Patents

Abstract

为解决高阶数字FIR滤波器并行处理效率的问题，提出了一种高效的、适合GPU体系结构的高阶数字FIR滤波器并行处理算法，该方法采用重叠保留方法结合GPU自身结构特点优化实现高阶数字FIR的频域并行化处理。通过计算FIR频率响应系数，将待处理的输入数据传送给GPU；数据重叠搬移；滤波计算处理；数据合并搬移；将合并搬移结果Y＝{Y₀，Y₁，....，Y_k-1}传送到主机内存等步骤完成高阶数字FIR滤波器频域并行处理。对比在CPU上单线程所实现的FIR频域重叠保留方法，其吞吐率，即每秒处理样点的数量有着极大地提高，典型的加速比在100倍以上。

Description

一种基于GPU的高阶数字FIR滤波器频域并行处理实现方法

技术领域

本发明涉及一种高阶数字FIR滤波器频域并行处理实现方法，特别是一种基于GPU的高阶数字FIR滤波器频域并行处理实现方法，属于数字信号处理领域。

技术背景

在数字信号处理系统中，数字有限冲击响应(FIR)滤波器是最为核心和基础的数字信号处理算法之一。由于数字FIR滤波器，特别是高阶的数字FIR滤波器的计算复杂度是相当高的，因而高效的、并行化的FIR实现方法对于加速FIR的处理是极其至关重要的。

目前，随着通用图形处理器(GPU)技术的迅猛发展，GPU正在被广泛地应用于众多的应用领域之中，由于GPU所具有的众核体系结构，使其能够提供十分强大的计算能力。

因此，对于FIR滤波器这种计算复杂度较高的算法，通过设计合理高效且适合GPU结构的并行算法，就可以使得FIR处理在GPU上得到很高的加速比，从而极大地缩短其处理时间。因此，研发适合于GPU体系结构的数字FIR滤波器并行处理算法，具有很高的实际价值和现实意义。

发明内容

本发明所要解决的技术问题是高阶数字FIR滤波器并行处理效率的问题，提出了一种高效的、适合GPU体系结构的高阶数字FIR滤波器并行处理算法，该方法采用重叠保留方法结合GPU自身结构特点优化实现高阶数字FIR的频域并行化处理。

为了发挥GPU众核的体系结构，需要将重叠保留方法的计算过程加以分解，并合理地将计算负载分配到在GPU上执行的每个线程中，同时还要优化各个线程对于GPU内存的访问，以最大限度地利用GPU所提供的存储带宽。在GPU的编程模型中，核函数(kernel)是由编程人员定义，并可以在GPU上由众多线程加以并行执行的功能单元。作为CPU加速器的GPU，是典型的fork-jion并行模式。在主机端启动核函数时，通过指定调用配置参数，编程人员可以控制启动执行核函数的线程数量以及线程的组织结构。如何将重叠保留方法中的计算过程分解成为不同的核函数，并确定各个核函数所要完成的处理，对于整个算法的并发执行效率起着至关重要的影响。如果核函数内部控制流比较复杂，则会大大地降低GPU的并行效率；然而如果核函数功能过于简单，将导致核函数的数量增加，从而增大启动核函数的总体时间，因此同样会降低GPU并行效率。

因此，根据GPU的自身结构特点，以及重叠保留方法所需的数据处理过程特性，提出了优化的核函数划分方法，即将重叠保留方法的处理划分为六个核函数：频响计算核函数、重叠搬移核函数、傅里叶变换核函数、乘法核函数、傅里叶逆变换核函数、合并搬移核函数，算法的总体处理过程如下：

根据通用图形处理器GPU众核体系结构的特点，将重叠保留方法的处理过程划分为六个核函数：频响计算核函数、重叠搬移核函数、傅里叶变换核函数、乘法核函数、傅里叶逆变换核函数、合并搬移核函数；

实现该方法的具体步骤如下：

步骤一、确定有限冲击响应FIR频率响应系数；

当给定的响应系数为FIR的冲击响应系数h＝{h(0)，h(1)，......，h(M-1)}时，将FIR的冲击响应系数h＝{h(0)，h(1)，......，h(M-1)}，经过尾部填0扩展至长度为N，将扩展后的系数传送给GPU，启动频响计算核函数对FIR的冲击响应系数h进行N点的傅里叶变换，结果为有限冲击响应FIR频率响应系数H＝{H(0)，H(1)，......，H(N-1)}，保存在GPU内存中；

当给定的响应系数为FIR频率响应系数H＝{H(0)，H(1)，......，H(N-1)}时，将频率响应系数H传给GPU，并保存在GPU内存中，保存时，第一个字节的地址为存储器位宽的整数倍；

步骤二、将待处理的输入数据传送给GPU；

将一块长度为N_blk的待滤波样点数据X＝{B₀，B₁，...，B_k-1}从主机内存中传入GPU的内存中，其中N_blk＝k*L，k为整数，L为每个数据块的长度；其中B_i＝{C_i，D_i}，0≤i≤k-1，C_i表示在B_i中起点为0，而长度为L-M+1的连续样点数据块，D_i表示在B_i中起点为L-M+1，而长度为M-1的连续样点数据块；

步骤三、数据的重叠搬移；

启动重叠搬移核函数完成数据重叠搬移操作，即将待滤波样点数据X＝{B₀，B₁，...，B_k-1}以及输入的长度为M-1的初始状态数据S₀，重叠搬移为E＝{E₀，E₁，...，E_k-1}，其中E₀＝{S₀，B₀}，其中，S₀为上一次处理待滤波样点数据X时其中的D_k-1，对于i≠0的数据块E_i＝{D_i-1，B_i}，同时将D_k-1搬移到中S₀，作为下一次处理过程的初始状态数据使用；

重叠数据的搬移分为三个步骤，第一步，启动M-1个线程同时工作，线程T_i(0≤i≤M-2)完成将S₀中第i个数据搬移到重叠搬移结果数据E₀的S₀中的第i的位置，然后再完成将D_k-1中第i个数据搬移到S₀中的第i的位置；第二步，启动L个线程同时工作，线程T_i(0≤i≤L-1)完成将待滤波样点数据X的B_j(0≤j≤k-1)块中第i个数据搬移到重叠搬移结果数据E的B_j(0≤j≤k-1)块中的第i的位置；第三步，启动M-1个线程同时工作，线程T_i(0≤i≤M-2)完成将待滤波样点数据X的D_j(0≤j≤k-1)块中第i个数据搬移到重叠搬移结果数据E的D_j(0≤j≤k-1)块中的第i的位置；

步骤四、滤波计算处理；

启动傅里叶变换核函数，完成对每个E_i做N点的傅里叶变换运算，得到运算结果F＝{F₀，F₁，...，F_k-1}；然后启动乘法核函数完成F_i与系数H相乘操作，即I_i＝F_i*H；然后启动傅里叶逆变换核函数，完成对每个I_i做N点的傅里叶逆变换运算，并得结果R＝{R₀，R₁，...，R_k-1}；

在启动乘法核函数核函数时，启动N个线程同时工作，线程T_i(0≤i≤N-1)完成将F_j(0≤j≤k-1)中第i个数据与频响系数H中的第i个系数相乘，并将结果保存在R_j(0≤j≤k-1)块中的第i的位置；

步骤五、数据的合并搬移；

启动合并搬移核函数完成数据的合并搬移操作，即将结果数据R＝{R₀，R₁，...，R_k-1}合并搬移为Y＝{Y₀，Y₁，...，Y_k-1}，其中R_i＝{Z_i，Y_i}，其中Z_i表示在R_i中起点为0，长度为M-1的连续数据点所组成的数据块，Y_i表示在R_i中起点为M-1，长度为L的连续数据点所组成的数据块；

在启动合并搬移核函数时，启动L个线程同时工作，线程T_i(0≤i≤L-1)完成将乘积R_j(0≤j≤k-1)中第M-1+i个数据搬移到结果数据Y_j(0≤j≤k-1)块中的第i的位置，其中，启动的L个线程在同一时刻搬运长度为L的连续数据块；

步骤六、将合并搬移结果Y＝{Y₀，Y₁，...，Y_k-1}传送到主机内存，若还有剩余数据需要处理重复步骤二，否则结束处理过程，完成高阶数字FIR滤波器频域并行处理。

有益效果：

通过对于GPU内存访问的优化设计，即所设计的重叠搬移以及合并搬移操作方法，较其他非连续的存储访问方法，能够最大限度地避免内存访问竞争的发生频率，从而可以有效地提高内存访问效率，并提高整体算法的执行效率。另外，通过优化线程组织结构以及启动线程的数量，能够达到最大限度地均衡各个流多处理器上的工作负载，从而提高GPU整体使用率。通过这些优化设计手段，较一般非优化设计，可平均提高执行效率10％左右。同时保存数据时，第一个字节的地址为存储器位宽的整数倍，以确保每次访问内存时能读取到最多的数据。

本发明方法对比在CPU上单线程所实现的FIR频域重叠保留方法，其吞吐率，即每秒处理样点的数量有着极大地提高，典型的加速比在100倍以上。如在GeForce580对于复数数据的处理可达到400M样点/秒。

附图说明

图1是通用数据重叠搬移操作图；

图2是实例数据重叠搬移操作图；

图3是通用数据合并搬移操作图；

图4是实例数据合并搬移操作图；

图5是并行FIR滤波处理流程框图；

具体实施方式

本实施例以1025点FIR为例，说明一个具体的实施方法。

我们选择傅里叶变换的长度为：2¹¹＝2048＝2K，L＝1024，每次处理的数据长度：N_blk＝2²⁰＝1M，故k＝N_blk/L＝2¹⁰＝1024。

根据前面所述方法，我们创建下列在GPU上执行的核函数：

1)频响计算核函数

其输入为1025点的FIR时域冲击响应，所进行的处理是将1024个点的时域冲击响应在其尾部填充1023个0点，并进行2048点的傅里叶变换，从而得到相应的频率响应系数，并将其存储到GPU的共享内存中。

2)重叠搬移核函数

其处理是将输入的2²⁰点样本数据，进行重叠搬移。图1为通用的重叠搬移示意图，X＝{B₀，B₁，...，B_k-1}是输入的待滤波的数据块，其长度为N_blk＝k*L；X是由k个子块B_i组成，0≤i≤k-1，每个子块的长度为L。每个子块B_i则是由长度为L-M+1的数据块C_i和长度为M-1的数据块D_i组成，即B_i＝{C_i，D_i}。S₀是重叠搬移状态数据块，其长度为M-1，初始值为M-1个0。重叠搬移的结果是E＝{E₀，E₁，...，E_k-1}和S¹ ₀，其中E₀＝{S₀，B₀}，对于i≠0的数据块E_i＝{D_i-1，B_i}，S¹ ₀＝D_k-1。而在本具体实例中由于L＝M-1＝1024，因此C_i数据块长度为0，即D_i＝B_i，具体搬移操作参见图2。重叠搬移过程分为三个步骤：

步骤201：启动1024个线程完成将状态数据块S₀搬移到E₀中，即将S₀(m)搬移到e(m)，其中0≤m≤1023，以及将D₁₀₂₃搬移到S¹ ₀中。在此步骤中，线程根据下列搬移公式进行数据搬移：

线程t(i)需要搬移e(i)＝S₀(i)，然后需要搬移S₀(i)＝x(2²⁰-2¹⁰+i)，其中0≤i≤1023。例如：线程t(0)将S₀(0)搬移到x(0)，将x(2²⁰-2¹⁰)搬移到S₀(0)；

步骤202：启动1024个线程完成将待滤波数据X搬移到E的奇数块中。在此步骤中，线程根据下列搬移公式进行数据搬移：

线程t(i)需要搬移e((2*m+1)*2¹⁰+i)＝x(m*2¹⁰+i)，其中0≤i≤1023，m＝0，1，...，2¹⁰-1。例如：线程t(0)在m＝0时，将数据x(0)搬移到e(1024)，在m＝1时，数据x(1024)搬移到e(3072)；

步骤203：启动1024个线程完成待滤波数据X搬移到E的偶数块中。在此步骤中，线程根据下列搬移公式进行数据搬移：

线程t(i)需要搬移e((2*(m+1)*2¹⁰+i)＝x(m*2¹⁰+i)，其中0≤i≤1023，m＝0，1，...，2¹⁰-1。例如：线程t(0)在m＝0时，将数据x(0)搬移到e(2048)，在m＝1时，数据x(1024)搬移到e(4096)；

3)傅里叶变换核函数

其处理是将经过重叠搬移后的数据，进行傅里叶变换计算，其中傅里叶变换的长度为2¹¹＝2048。

4)乘法核函数

其处理是将傅里叶变换的结果数据，与FIR频响系数相乘。在此过程中，启动2048个线程，并根据下列公式进行计算：

线程t(i)计算p(m＊2¹¹+i)＝e(m＊2¹¹+i)*H(i)/2¹¹，其中，0≤i≤2047，m＝0，1，...，2¹⁰-1，e为重叠搬移操作的结果数据。例如：线程t(0)在m＝0时，计算p(0)＝e(0)*H(0)/2¹¹，在m＝1时，计算p(2¹¹)＝e(2¹¹)*H(0)/2¹¹；

5)傅里叶逆变换核函数

其处理是将与频响系数相乘所得的结果数据，进行傅里叶逆变换计算，其中傅里叶逆变换的长度为2¹¹＝2048。

6)合并搬移核函数

其处理是将傅里叶逆变换的结果数据，进行合并搬移。图3为通用的合并搬移示意图，将傅里叶逆变换的结果数据R＝{R₀，R₁，...，R_k-1}合并搬移为Y＝{Y₀，Y₁，...，Y_k-1}，其中R_i＝{Z_i，Y_i}，其中Z_i长度为M-1的连续数据点，Y_i则是长度为L的连续数据点。而在本具体实例中，L＝M-1＝2¹⁰＝1024，具体搬移操作参见图4。在此过程中，启动1024个线程，并根据下列公式进行数据搬移操作：

线程t(i)需要搬移y(m＊2¹⁰+i)＝r((2*m+1)*2¹⁰+i)，其中，0≤i≤1023，m＝0，1，...，2¹⁰-1，y为搬移后的结果数据。例如：线程t(0)在m＝0时，将数据r(2¹⁰)搬移到y(0)，在m＝1时，将数据r(3072)搬移到y(1024)；

上面六个核函数是条件，有了上述的六个可在GPU上并行执行的核函数之后，我们就可以按照下述步骤来完成FIR的并行处理，其中具体的滤波计算处理操作参见图5：

步骤501：将1025点的FIR冲击响应数据从主机内存拷贝到GPU内存中，启动频响计算核函数，完成FIR频响系数的计算。

步骤502：将一块2²⁰＝1M点的数据从主机内存拷贝到GPU内存中，在GPU上启动2¹⁰个线程，用以执行重叠搬移核函数。

步骤503：执行傅里叶变换核函数。

步骤504：在GPU上启动2¹¹个线程，用以执行乘法核函数。

步骤505：执行傅里叶逆变换核函数。

步骤506：在GPU上启动2¹⁰个线程，用以执行合并搬移核函数。

步骤507：将合并搬移的结果从GPU内存中拷贝到主机内存中，如果有待处理的数据，则转到步骤二执行。

Claims (1)

1.一种基于GPU的高阶数字FIR滤波器频域并行处理实现方法，其特征在于：根据通用图形处理器GPU众核体系结构的特点，将重叠保留方法的处理过程划分为六个核函数：频响计算核函数、重叠搬移核函数、傅里叶变换核函数、乘法核函数、傅里叶逆变换核函数、合并搬移核函数；

实现该方法的具体步骤如下：

步骤一、确定有限冲击响应FIR频率响应系数；

当给定的响应系数为FIR的冲击响应系数h＝{h(0)，h(1)，......，h(M-1)}时，将FIR的冲击响应系数h＝{h(0)，h(1)，......，h(M-1)}，经过尾部填0扩展至长度为N，将扩展后的系数传送给GPU，启动频响计算核函数对FIR的冲击响应系数h进行N点的傅里叶变换，结果为有限冲击响应FIR频率响应系数H＝{H(0)，H(1)，......，H(N-1)}，保存在GPU内存中；

当给定的响应系数为FIR频率响应系数H＝{H(0)，H(1)，......，H(N-1)}时，将频率响应系数H传给GPU，并保存在GPU内存中，保存时，第一个字节的地址为存储器位宽的整数倍；

步骤二、将待处理的输入数据传送给GPU；

将一块长度为N_blk的待滤波样点数据X＝{B₀，B₁，...，B_k-1}从主机内存中传入GPU的内存中，其中N_blk＝k*L，k为整数，L为每个数据块的长度；其中B_i＝{C_i，D_i}，0≤i≤k-1，C_i表示在B_i中起点为0，而长度为L-M+1的连续样点数据块，D_i表示在B_i中起点为L-M+1，而长度为M-1的连续样点数据块；

步骤三、数据的重叠搬移；

启动重叠搬移核函数完成数据重叠搬移操作，即将待滤波样点数据X＝{B₀，B₁，...，B_k-1}以及输入的长度为M-1的初始状态数据S₀，重叠搬移为E＝{E₀，E₁，...，E_k-1}，其中E₀＝{S₀，B₀}，其中，S₀为上一次处理待滤波样点数据X时其中的D_k-1，对于i≠0的数据块E_i＝{D_i-1，B_i}，同时将D_k-1搬移到中S₀，作为下一次处理过程的初始状态数据使用；

重叠数据的搬移分为三个步骤，第一步，启动M-1个线程同时工作，线程T_i(0≤i≤M-2)完成将S₀中第i个数据搬移到重叠搬移结果数据E₀的S₀中的第i的位置，然后再完成将D_k-1中第i个数据搬移到S₀中的第i的位置；第二步，启动L个线程同时工作，线程T_i(0≤i≤L-1)完成将待滤波样点数据X的B_j(0≤j≤k-1)块中第i个数据搬移到重叠搬移结果数据E的B_j(0≤j≤k-1)块中的第i的位置；第三步，启动M-1个线程同时工作，线程T_i(0≤i≤M-2)完成将待滤波样点数据X的D_j(0≤j≤k-1)块中第i个数据搬移到重叠搬移结果数据E的D_j(0≤j≤k-1)块中的第i的位置；

步骤四、滤波计算处理；

启动傅里叶变换核函数，完成对每个E_i做N点的傅里叶变换运算，得到运算结果F＝{F₀，F₁，...，F_k-1}；然后启动乘法核函数完成F_i与系数H相乘操作，即I_i＝F_i*H；然后启动傅里叶逆变换核函数，完成对每个I_i做N点的傅里叶逆变换运算，并得结果R＝{R₀，R₁，...，R_k-1}；

在启动乘法核函数核函数时，启动N个线程同时工作，线程T_i(0≤i≤N-1)完成将F_j(0≤j≤k-1)中第i个数据与频响系数H中的第i个系数相乘，并将结果保存在R_j(0≤j≤k-1)块中的第i的位置；

步骤五、数据的合并搬移；

启动合并搬移核函数完成数据的合并搬移操作，即将结果数据R＝{R₀，R₁，...，R_k-1}合并搬移为Y＝{Y₀，Y₁，...，Y_k-1}，其中R_i＝{Z_i，Y_i}，其中Z_i表示在R_i中起点为0，长度为M-1的连续数据点所组成的数据块，Y_i表示在R_i中起点为M-1，长度为L的连续数据点所组成的数据块；

在启动合并搬移核函数时，启动L个线程同时工作，线程T_i(0≤i≤L-1)完成将乘积R_j(0≤j≤k-1)中第M-1+i个数据搬移到结果数据Y_j(0≤j≤k-1)块中的第i的位置，其中，启动的L个线程在同一时刻搬运长度为L的连续数据块；

步骤六、将合并搬移结果Y＝{Y₀，Y₁，...，Y_k-1}传送到主机内存，若还有剩余数据需要处理重复步骤二，否则结束处理过程，完成高阶数字FIR滤波器频域并行处理。

Images