为解决高阶数字FIR滤波器并行处理效率的问题，提出了一种高效的、适合GPU体系结构的高阶数字FIR滤波器并行处理算法，该方法采用重叠保留方法结合GPU自身结构特点优化实现高阶数字FIR的频域并行化处理。通过计算FIR频率响应系数，将待处理的输入数据传送给GPU；数据重叠搬移；滤波计算处理；数据合并搬移；将合并搬移结果Y＝{Y0，Y1，....，Yk-1}传送到主机内存等步骤完成高阶数字FIR滤波器频域并行处理。对比在CPU上单线程所实现的FIR频域重叠保留方法，其吞吐率，即每秒处理样点的数量有着极大地提高，典型的加速比在100倍以上。
1.一种基于GPU的高阶数字FIR滤波器频域并行处理实现方法，其特征在于：根据通用图形处理器GPU众核体系结构的特点，将重叠保留方法的处理过程划分为六个核函数：频响计算核函数、重叠搬移核函数、傅里叶变换核函数、乘法核函数、傅里叶逆变换核函数、合并搬移核函数；
当给定的响应系数为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内存中，保存时，第一个字节的地址为存储器位宽的整数倍；
将一块长度为Nblk的待滤波样点数据X＝｛B0，B1，...，Bk-1}从主机内存中传入GPU的内存中，其中Nblk＝k*L，k为整数，L为每个数据块的长度；其中Bi＝{Ci，Di}，0≤i≤k-1，Ci表示在Bi中起点为0，而长度为L-M+1的连续样点数据块，Di表示在Bi中起点为L-M+1，而长度为M-1的连续样点数据块；
启动重叠搬移核函数完成数据重叠搬移操作，即将待滤波样点数据X＝{B0，B1，...，Bk-1}以及输入的长度为M-1的初始状态数据S0，重叠搬移为E＝{E0，E1，...，Ek-1}，其中E0＝{S0，B0}，其中，S0为上一次处理待滤波样点数据X时其中的Dk-1，对于i≠0的数据块Ei＝{Di-1，Bi｝，同时将Dk-1搬移到S0中，作为下一次处理过程的初始状态数据使用；
重叠数据的搬移分为三个步骤，第一步，启动M-1个线程同时工作，线程Ti，
0≤i≤M-2，完成将S0中第i个数据搬移到重叠搬移结果数据E0的S0中的第i的位置，然后再完成将Dk-1中第i个数据搬移到S0中的第i的位置；第二步，启动L个线程同时工作，线程Ti，0≤i≤L-1，完成将待滤波样点数据X的Bj，0≤j≤k-1，块中第i个数据搬移到重叠搬移结果数据E的Bj，0≤j≤k-1，块中的第i的位置；第三步，启动M-1个线程同时工作，线程Ti，0≤i≤M-2，完成将待滤波样点数据X的Dj，0≤j≤k-1，块中第i个数据搬移到重叠搬移结果数据E的Dj，0≤j≤k-1，块中的第i的位置；
启动傅里叶变换核函数，完成对每个Ei做N点的傅里叶变换运算，得到运算结果F＝{F0，F1，...，Fk-1}；然后启动乘法核函数完成Fi与系数H相乘操作，即Ii＝Fi*H；然后启动傅里叶逆变换核函数，完成对每个Ii做N点的傅里叶逆变换运算，并得结果R＝{R0，R1，...，Rk-1｝；
在启动乘法核函数时，启动N个线程同时工作，线程Ti，0≤i≤N-1，完成将Fj，
0≤j≤k-1，中第i个数据与频响系数H中的第i个系数相乘，并将结果保存在Rj，
启动合并搬移核函数完成数据的合并搬移操作，即将结果数据R＝{R0，R1，...，Rk-1}合并搬移为Y＝{Y0，Y1，...，Yk-1｝，其中Ri＝{Zi，Yi}，其中Zi表示在Ri中起点为0，长度为M-1的连续数据点所组成的数据块，Yi表示在Ri中起点为M-1，长度为L的连续数据点所组成的数据块；
在启动合并搬移核函数时，启动L个线程同时工作，线程Ti，0≤i≤L-1，完成将乘积Rj，0≤j≤k-1，中第M-1+i个数据搬移到结果数据Yj，0≤j≤k-1，块中的第i的位置，其中，启动的L个线程在同一时刻搬运长度为L的连续数据块；
步骤六、将合并搬移结果Y＝{Y0，Y1，...，Yk-1｝传送到主机内存，若还有剩余数据需要处理重复步骤二，否则结束处理过程，完成高阶数字FIR滤波器频域并行处理。
一种基于GPU的高阶数字FIR滤波器频域并行处理实现方
法
技术领域
[0001] 本发明涉及一种高阶数字FIR滤波器频域并行处理实现方法，特别是一种基于GPU的高阶数字FIR滤波器频域并行处理实现方法，属于数字信号处理领域。技术背景
[0002] 在数字信号处理系统中，数字有限冲击响应(FIR)滤波器是最为核心和基础的数字信号处理算法之一。由于数字FIR滤波器，特别是高阶的数字FIR滤波器的计算复杂度是相当高的，因而高效的、并行化的FIR实现方法对于加速FIR的处理是极其至关重要的。
[0003] 目前，随着通用图形处理器(GPU)技术的迅猛发展，GPU正在被广泛地应用于众多的应用领域之中，由于GPU所具有的众核体系结构，使其能够提供十分强大的计算能力。
[0004] 因此，对于FIR滤波器这种计算复杂度较高的算法，通过设计合理高效且适合GPU结构的并行算法，就可以使得FIR处理在GPU上得到很高的加速比，从而极大地缩短其处理时间。因此，研发适合于GPU体系结构的数字FIR滤波器并行处理算法，具有很高的实际价值和现实意义。
发明内容
[0005] 本发明所要解决的技术问题是高阶数字FIR滤波器并行处理效率的问题，提出了一种高效的、适合GPU体系结构的高阶数字FIR滤波器并行处理算法，该方法采用重叠保留方法结合GPU自身结构特点优化实现高阶数字FIR的频域并行化处理。
[0006] 为了发挥GPU众核的体系结构，需要将重叠保留方法的计算过程加以分解，并合理地将计算负载分配到在GPU上执行的每个线程I，同时还要优化各个线程对于GPU内存的访问，以最大限度地利用GPU所提供的存储带宽。在GPU的编程模型中，核函数(kernel)是由编程人员定义，并可以在GPU上由众多线程加以并行执行的功能单元。作为CPU加速器的GPU，是典型的fork-jion并行模式。在主机端启动核函数时，通过指定调用配置参数，编程人员可以控制启动执行核函数的线程数量以及线程的组织结构。如何将重叠保留方法中的计算过程分解成为不同的核函数，并确定各个核函数所要完成的处理，对于整个算法的并发执行效率起着至关重要的影响。如果核函数内部控制流比较复杂，则会大大地降低GPU的并行效率；然而如果核函数功能过于简单，将导致核函数的数量增加，从而增大启动核函数的总体时间，因此同样会降低GPU并行效率。
[0007] 因此，根据GPU的自身结构特点，以及重叠保留方法所需的数据处理过程特性，提出了优化的核函数划分方法，即将重叠保留方法的处理划分为六个核函数：频响计算核函数、重叠搬移核函数、傅里叶变换核函数、乘法核函数、傅里叶逆变换核函数、合并搬移核函数，算法的总体处理过程如下：
[0008] 根据通用图形处理器GPU众核体系结构的特点，将重叠保留方法的处理过程划分为六个核函数：频响计算核函数、重叠搬移核函数、傅里叶变换核函数、乘法核函数、傅里叶逆变换核函数、合并搬移核函数；
[0009] 实现该方法的具体步骤如下：
[0010] 步骤一、确定有限冲击响应FIR频率响应系数；
[0011] 当给定的响应系数为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内存中；
[0012] 当给定的响应系数为FIR频率响应系数H＝{H(0)，H(1)，......，H(N-1)}时，将频率响应系数H传给GPU，并保存在GPU内存中，保存时，第一个字节的地址为存储器位宽的整数倍；
[0013] 步骤二、将待处理的输入数据传送给GPU；
[0014] 将一块长度为Nblk的待滤波样点数据X＝{B0，B1，...，Bk-1}从主机内存中传入GPU的内存中，其中Nblk＝k*L，k为整数，L为每个数据块的长度；其中Bi＝{Ci，Di}，0≤i≤k-1，Ci表示在Bi中起点为0，而长度为L-M+1的连续样点数据块，Di表示在Bi中起点为L-M+1，而长度为M-1的连续样点数据块；
[0015] 步骤三、数据的重叠搬移；
[0016] 启动重叠搬移核函数完成数据重叠搬移操作，即将待滤波样点数据X＝{B0，B1，...，Bk-1｝以及输入的长度为M-1的初始状态数据S0，重叠搬移为E＝{E0，B1，...，Ek-1}，其中E0＝{S0，B0}，其中，S0为上一次处理待滤波样点数据X时其中的Dk-1，对于i≠0的数据块Ei＝{Di-1，Bi}，同时将Dk-1搬移到S0中，作为下一次处理过程的初始状态数据使用；
[0017] 重叠数据的搬移分为三个步骤，第一步，启动M-1个线程同时工作，线程Ti(0≤i≤M-2)完成将S0中第i个数据搬移到重叠搬移结果数据E0的S0中的第i的位置，然后再完成将Dk-1中第i个数据搬移到S0中的第i的位置；第二步，启动L个线程同时工作，线程Ti(0≤i≤L-1)完成将待滤波样点数据X的Bj(0≤j≤K-1)块中第i个数据搬移到重叠搬移结果数据E的Bj(0≤j≤k-1)块中的第i的位置；第三步，启动M-1个线程同时工作，线程Ti(0≤i≤M-2)完成将待滤波样点数据X的Dj(0≤j≤k-1)块中第i个数据搬移到重叠搬移结果数据E的Dj(0≤j≤k-1)块中的第i的位置；
[0018] 步骤四、滤波计算处理；
[0019] 启动傅里叶变换核函数，完成对每个Ei做N点的傅里叶变换运算，得到运算结果F＝{F0，F1，...，Fk-1}；然后启动乘法核函数完成Fi与系数H相乘操作，即Ii＝Fi*H；然后启动傅里叶逆变换核函数，完成对每个Ii做N点的傅里叶逆变换运算，并得结果R＝{R0，R1，...，Rk-1}；
[0020] 在启动乘法核函数时，启动N个线程同时工作，线程Ti(0≤i≤N-1)完成将Fj(0≤j≤k-1)中第i个数据与频响系数H中的第i个系数相乘，并将结果保存在Rj(0≤j≤k-1)块中的第i的位置；
[0021] 步骤五、数据的合并搬移；
[0022] 启动合并搬移核函数完成数据的合并搬移操作，即将结果数据R＝{R0，R1，...，Rk-1｝合并搬移为Y＝{Y0，Y1，...，Yk-1}，其中Ri＝{Zi，Yi}，其中Zi表示在Ri中起点为0，长度为M-1的连续数据点所组成的数据块，Yi表示在Ri中起点为M-1，长度为L的连续数据点所组成的数据块；
[0023] 在启动合并搬移核函数时，启动L个线程同时工作，线程Ti(0≤i≤L-1)完成将乘积Rj(0≤j≤k-1)中第M-1+i个数据搬移到结果数据Yj(0≤j≤k-1)块中的第i的位置，其中，启动的L个线程在同一时刻搬运长度为L的连续数据块；
[0024] 步骤六、将合并搬移结果Y＝{Y0，Y1，...，Yk-1}传送到主机内存，若还有剩余数据需要处理重复步骤二，否则结束处理过程，完成高阶数字FIR滤波器频域并行处理。
[0025] 有益效果：
[0026] 通过对于GPU内存访问的优化设计，即所设计的重叠搬移以及合并搬移操作方法，较其他非连续的存储访问方法，能够最大限度地避免内存访问竞争的发生频率，从而可以有效地提高内存访问效率，并提高整体算法的执行效率。另外，通过优化线程组织结构以及启动线程的数量，能够达到最大限度地均衡各个流多处理器上的工作负载，从而提高GPU整体使用率。通过这些优化设计手段，较一般非优化设计，可平均提高执行效率10％左右。同时保存数据时，第一个字节的地址为存储器位宽的整数倍，以确保每次访问内存时能读取到最多的数据。
[0027] 本发明方法对比在CPU上单线程所实现的FIR频域重叠保留方法，其吞吐率，即每秒处理样点的数量有着极大地提高，典型的加速比在100倍以上。如在GeForce580对于复数数据的处理可达到400M样点/秒。
附图说明
[0028] 图1是通用数据重叠搬移操作图；
[0029] 图2是实例数据重叠搬移操作图；
[0030] 图3是通用数据合并搬移操作图；
[0031] 图4是实例数据合并搬移操作图；
[0032] 图5是并行FIR滤波处理流程框图；
具体实施方式
[0033] 本实施例以1025点FIR为例，说明一个具体的实施方法。
[0034] 我们选择傅里叶变换的长度为：211＝2048＝2K，L＝1024，每次处理的数据长度：20 10
Nblk＝2 ＝1M，故k＝Nblk/L＝2 ＝1024。
[0035] 根据前面所述方法，我们创建下列在GPU上执行的核函数：
[0036] 1)频响计算核函数
[0037] 其输入为1025点的FIR时域冲击响应，所进行的处理是将1024个点的时域冲击响应在其尾部填充1023个0点，并进行2048点的傅里叶变换，从而得到相应的频率响应系数，并将其存储到GPU的共享内存中。
[0038] 2)重叠搬移核函数
[0039] 其处理是将输入的220点样本数据，进行重叠搬移。图1为通用的重叠搬移示意图，X＝{B0，B1，...，Bk-1}是输入的待滤波的数据块，其长度为Nblk＝k*L；X是由k个子块Bi组成，0≤i≤k-1，每个子块的长度为L。每个子块Bi则是由长度为L-M+1的数据块Ci和长度为M-1的数据块Di组成，即Bi＝{Ci，Di｝。S0是重叠搬移状态数据块，其长度为M-1，1
初始值为M-1个0。重叠搬移的结果是E＝{E0，E1，...，Ek-1}和S0，其中E0＝{S0，B0｝，对
1
于i≠0的数据块Ei＝{Di-1，Bi}，S0＝Dk-1。而在本具体实例中由于L＝M-1＝1024，因此Ci数据块长度为0，即Di＝Bi，具体搬移操作参见图2。重叠搬移过程分为三个步骤：
[0040] 步骤201：启动1024个线程完成将状态数据块S0搬移到E0中，即将S0(m)搬移到1
e(m)，其中0≤m≤1023，以及将D1023搬移到S0中。在此步骤中，线程根据下列搬移公式进行数据搬移：
[0041] 线程t(i)需要搬移e(i)＝S0(i)，然后需要搬移S0(i)＝x(220-210+i)，其中20 10
0≤i≤1023。例如：线程t(0)将S0(0)搬移到x(0)，将x(2 -2 )搬移到S0(0)；
[0042] 步骤202：启动1024个线程完成将待滤波数据X搬移到E的奇数块中。在此步骤中，线程根据下列搬移公式进行数据搬移：
[0043] 线程t(i)需要搬移e((2*m+1)*210+i)＝x(m*210+i)，其中0≤i≤1023，m＝0，10
1，...，2 -1。例如：线程t(0)在m＝0时，将数据x(0)搬移到e(1024)，在m＝1时，数据x(1024)搬移到e(3072)；
[0044] 步骤203：启动1024个线程完成待滤波数据X搬移到E的偶数块中。在此步骤中，线程根据下列搬移公式进行数据搬移：
[0045] 线程t(i)需要搬移e((2*(m+1)*210+i)＝x(m*210+i)，其中0≤i≤1023，m＝0，10
1，...，2 -1。例如：线程t(0)在m＝0时，将数据x(0)搬移到e(2048)，在m＝1时，数据x(1024)搬移到e(4096)；
[0046] 3)傅里叶变换核函数
[0047] 其处理是将经过重叠搬移后的数据，进行傅里叶变换计算，其中傅里叶变换的长11
度为2 ＝2048。
[0048] 4)乘法核函数
[0049] 其处理是将傅里叶变换的结果数据，与FIR频响系数相乘。在此过程中，启动204811 11 11
个线程，并根据下列公式进行计算：线程t(i)计算p(m＊2 +i)＝e(m＊2 +i)*H(i)/2 ，
10
其中，0≤i≤2047，m＝0，1，...，2 -1，e为重叠搬移操作的结果数据。例如：线程t(0)
11 11 11 11
在m＝0时，计算p(0)＝e(0)*H(0)/2 ，在m＝1时，计算p(2 )＝e(2 )*H(0)/2 ；
[0050] 5)傅里叶逆变换核函数
[0051] 其处理是将与频响系数相乘所得的结果数据，进行傅里叶逆变换计算，其中傅里11
叶逆变换的长度为2 ＝2048。
[0052] 6)合并搬移核函数
[0053] 其处理是将傅里叶逆变换的结果数据，进行合并搬移。图3为通用的合并搬移示意图，将傅里叶逆变换的结果数据R＝{R0，R1，...，Rk-1}合并搬移为Y＝{Y0，Y1，...，Yk-1}，其中Ri＝{Zi，Yi}，其中Zi长度为M-1的连续数据点，Yi则是长度为L的连续数据点。而在10
本具体实例中，L＝M-1＝2 ＝1024，具体搬移操作参见图4。在此过程中，启动1024个线程，并根据下列公式进行数据搬移操作：
[0054] 线程t(i)需要搬移y(m＊210+i)＝r((2*m+1)*210+i)，其中，0≤i≤1023，m＝10 10
0，1，...，2 -1，y为搬移后的结果数据。例如：线程t(0)在m＝0时，将数据r(2 )搬移到y(0)，在m＝1时，将数据r(3072)搬移到y(1024)；
[0055] 上面六个核函数是条件，有了上述的六个可在GPU上并行执行的核函数之后，我们就可以按照下述步骤来完成FIR的并行处理，其中具体的滤波计算处理操作参见图5：
[0056] 步骤501：将1025点的FIR冲击响应数据从主机内存拷贝到GPU内存中，启动频响计算核函数，完成FIR频响系数的计算。
[0057] 步骤502：将一块220＝1M点的数据从主机内存拷贝到GPU内存中，在GPU上启动10
2 个线程，用以执行重叠搬移核函数。
[0058] 步骤503：执行傅里叶变换核函数。
[0059] 步骤504：在GPU上启动211个线程，用以执行乘法核函数。
[0060] 步骤505：执行傅里叶逆变换核函数。
[0061] 步骤506：在GPU上启动210个线程，用以执行合并搬移核函数。
[0062] 步骤507：将合并搬移的结果从GPU内存中拷贝到主机内存中，如果有待处理的数据，则转到步骤二执行。
