专利名称：一种基于gpu的高阶数字fir滤波器频域并行处理实现方法
技术领域：
本发明涉及一种高阶数字FIR滤波器频域并行处理实现方法，特别是一种基于 GPU的高阶数字HR滤波器频域并行处理实现方法，属于数字信号处理领域。
背景技术：
在数字信号处理系统中，数字有限冲击响应(FIR)滤波器是最为核心和基础的数字信号处理算法之一。由于数字FIR滤波器，特别是高阶的数字HR滤波器的计算复杂度是相当高的，因而高效的、并行化的FIR实现方法对于加速FIR的处理是极其至关重要的。目前，随着通用图形处理器(GPU)技术的迅猛发展，GPU正在被广泛地应用于众多的应用领域之中，由于GPU所具有的众核体系结构，使其能够提供十分强大的计算能力。因此，对于FIR滤波器这种计算复杂度较高的算法，通过设计合理高效且适合GPU 结构的并行算法，就可以使得HR处理在GPU上得到很高的加速比，从而极大地缩短其处理时间。因此，研发适合于GPU体系结构的数字HR滤波器并行处理算法，具有很高的实际价值和现实意义。
发明内容
本发明所要解决的技术问题是高阶数字FIR滤波器并行处理效率的问题，提出了一种高效的、适合GPU体系结构的高阶数字HR滤波器并行处理算法，该方法采用重叠保留方法结合GPU自身结构特点优化实现高阶数字FIR的频域并行化处理。为了发挥GPU众核的体系结构，需要将重叠保留方法的计算过程加以分解，并合理地将计算负载分配到在GPU上执行的每个线程中，同时还要优化各个线程对于GPU内存的访问，以最大限度地利用GPU所提供的存储带宽。在GPU的编程模型中，核函数(kernel) 是由编程人员定义，并可以在GPU上由众多线程加以并行执行的功能单元。作为CPU加速器的GPU，是典型的fork-jion并行模式。在主机端启动核函数时，通过指定调用配置参数， 编程人员可以控制启动执行核函数的线程数量以及线程的组织结构。如何将重叠保留方法中的计算过程分解成为不同的核函数，并确定各个核函数所要完成的处理，对于整个算法的并发执行效率起着至关重要的影响。如果核函数内部控制流比较复杂，则会大大地降低 GPU的并行效率；然而如果核函数功能过于简单，将导致核函数的数量增加，从而增大启动核函数的总体时间，因此同样会降低GPU并行效率。因此，根据GPU的自身结构特点，以及重叠保留方法所需的数据处理过程特性，提出了优化的核函数划分方法，即将重叠保留方法的处理划分为六个核函数频响计算核函数、重叠搬移核函数、傅里叶变换核函数、乘法核函数、傅里叶逆变换核函数、合并搬移核函数，算法的总体处理过程如下根据通用图形处理器GPU众核体系结构的特点，将重叠保留方法的处理过程划分为六个核函数频响计算核函数、重叠搬移核函数、傅里叶变换核函数、乘法核函数、傅里叶逆变换核函数、合并搬移核函数；实现该方法的具体步骤如下步骤一、确定有限冲击响应FIR频率响应系数；当给定的响应系数为FIR的冲击响应系数h = {h(0)，h(l)，......，h (M-I)}时，
将FIR的冲击响应系数h = {h (0)，h (1)，......，h (M-I)}，经过尾部填0扩展至长度为N，
将扩展后的系数传送给GPU，启动频响计算核函数对FIR的冲击响应系数h进行N点的傅里
叶变换，结果为有限冲击响应FIR频率响应系数H = {H (0)，H (1)，......，H (N-I)}，保存在
GPU内存中；当给定的响应系数为HR频率响应系数H= IH(O)，H(I)，......，H(N-I)}时，将
频率响应系数H传给GPU，并保存在GPU内存中，保存时，第一个字节的地址为存储器位宽的整数倍；步骤二、将待处理的输入数据传送给GPU ；将一块长度为Nblk的待滤波样点数据X = {B0, B1, ...，Bk_J从主机内存中传入GPU的内存中，其中Nblk = k*L，k为整数，L为每个数据块的长度；其中Bi = ICi, Dj, O^ I^k-LCi表示在Bi中起点为0，而长度为L-M+1的连续样点数据块，Di表示在Bi中起点为L-M+1，而长度为M-I的连续样点数据块；步骤三、数据的重叠搬移；启动重叠搬移核函数完成数据重叠搬移操作，即将待滤波样点数据X = {B0, B^.^Bk—J以及输入的长度为M-I的初始状态数据&，重叠搬移为E = (EojE1,... ,Ek_J, 其中= {&，BJ，其中，&为上一次处理待滤波样点数据X时其中的Dlri，对于i兴0的数据块Ei = {Dh，BJ，同时将01;_1搬移到中&，作为下一次处理过程的初始状态数据使用；重叠数据的搬移分为三个步骤，第一步，启动M-I个线程同时工作，线程 Ti(0彡i彡M-2)完成将&中第i个数据搬移到重叠搬移结果数据的&中的第i的位置，然后再完成将Dlri中第i个数据搬移到&中的第i的位置；第二步，启动L个线程同时工作，线程Ti (0彡i彡L-1)完成将待滤波样点数据X的Bj (0彡j彡k-Ι)块中第i个数据搬移到重叠搬移结果数据E的Bj (0彡j彡k-Ι)块中的第i的位置；第三步，启动M-I个线程同时工作，线程Ti (0彡i彡M-2)完成将待滤波样点数据X的Dj (0彡j彡k-Ι)块中第i 个数据搬移到重叠搬移结果数据E的DdO ^ j ^ k-1)块中的第i的位置；步骤四、滤波计算处理；启动傅里叶变换核函数，完成对每个EJiA N点的傅里叶变换运算，得到运算结果 F= {F0, F1, ... , Fk_J ；然后启动乘法核函数完成Fi与系数H相乘操作，即Ii = FjH;然后启动傅里叶逆变换核函数，完成对每个Ii做N点的傅里叶逆变换运算，并得结果R = {R0,
Ri，· · ·，Rk—工}；在启动乘法核函数核函数时，启动N个线程同时工作，线程Ti (0彡i ^ N-1)完成将h (0 ^ j ^ k-1)中第i个数据与频响系数H中的第i个系数相乘，并将结果保存在 Rj (0 ( j ( k-1)块中的第i的位置；步骤五、数据的合并搬移；启动合并搬移核函数完成数据的合并搬移操作，即将结果数据R = {R0, R1,...， Rk-J合并搬移为Y = {Y0，Y1, ...，Yk-J，其中Ri = IZi, YJ，其中Zi表示在Ri中起点为0，长度为M-I的连续数据点所组成的数据块，Yi表示在氏中起点为M-I，长度为L的连续数据点所组成的数据块；在启动合并搬移核函数时，启动L个线程同时工作，线程Ti (0彡i彡L-1)完成将乘积民(0彡j彡k-Ι)中第M-1+i个数据搬移到结果数据YjO彡j彡k-Ι)块中的第i的位置，其中，启动的L个线程在同一时刻搬运长度为L的连续数据块；步骤六、将合并搬移结果Y = {Y0, Y1, ... , Yk_J传送到主机内存，若还有剩余数据需要处理重复步骤二，否则结束处理过程，完成高阶数字HR滤波器频域并行处理。有益效果通过对于GPU内存访问的优化设计，即所设计的重叠搬移以及合并搬移操作方法，较其他非连续的存储访问方法，能够最大限度地避免内存访问竞争的发生频率，从而可以有效地提高内存访问效率，并提高整体算法的执行效率。另外，通过优化线程组织结构以及启动线程的数量，能够达到最大限度地均衡各个流多处理器上的工作负载，从而提高GPU 整体使用率。通过这些优化设计手段，较一般非优化设计，可平均提高执行效率10%左右。 同时保存数据时，第一个字节的地址为存储器位宽的整数倍，以确保每次访问内存时能读取到最多的数据。本发明方法对比在CPU上单线程所实现的HR频域重叠保留方法，其吞吐率，即每秒处理样点的数量有着极大地提高，典型的加速比在100倍以上。如在GeR)rce580对于复数数据的处理可达到400M样点/秒。
图1是通用 ft据重叠搬移操作图2是实例ift据重叠搬移操作图3是通用ift据合并搬移操作图4是实例ift据合并搬移操作图5是并行FIR滤波处理流程框图
具体实施例方式本实施例以1025点FIR为例，说明一个具体的实施方法。我们选择傅里叶变换的长度为2^ = 2048 = 14= IOM，每次处理的数据长度 Nblk = 220 = 1M，故 k = Nblk/L = 210 = 1024。根据前面所述方法，我们创建下列在GPU上执行的核函数1)频响计算核函数其输入为1025点的HR时域冲击响应，所进行的处理是将IOM个点的时域冲击响应在其尾部填充1023个0点，并进行2048点的傅里叶变换，从而得到相应的频率响应系数，并将其存储到GPU的共享内存中。2)重叠搬移核函数其处理是将输入的22°点样本数据，进行重叠搬移。图1为通用的重叠搬移示意图，X = {B0, B1, ...，Bk_J是输入的待滤波的数据块，其长度为Nblk = k*L ；X是由k个子块 Bi组成，0 ( i ( k-Ι，每个子块的长度为L。每个子块Bi则是由长度为L-M+1的数据块Ci
6和长度为M-I的数据块DJi成，即Bi= (Ci5DJ0 &是重叠搬移状态数据块，其长度为M-I， 初始值为M-I个0。重叠搬移的结果是E = {E0, E1, ... , Ek_J和S1tl，其中= {S0, Bj，对于i乒0的数据块Ei = {D^，BJ，S10 = Dh。而在本具体实例中由于L = M-I = 1024，因此Ci数据块长度为0，即Di =Bi,具体搬移操作参见图2。重叠搬移过程分为三个步骤步骤201 启动IOM个线程完成将状态数据块&搬移到中，即将搬移到 e(m)，其中0彡m彡1023，以及将Dlte3搬移到S1tl中。在此步骤中，线程根据下列搬移公式进行数据搬移线程t(i)需要搬移e(i) = S^i)，然后需要搬移= H22°-21Q+i)，其中 0彡i彡1023。例如线程t(0)将S0(O)搬移到X(O)，将X(22°-21Q)搬移到S0(O)；步骤202 启动IOM个线程完成将待滤波数据X搬移到E的奇数块中。在此步骤中，线程根据下列搬移公式进行数据搬移线程t ⑴需要搬移 e ((2*m+l) *21Q+i) = χ (m*210+i)，其中 0 彡 i 彡 1023，m = 0， 1，. . .，210-1。例如线程t (0)在m = 0时，将数据χ (0)搬移到e (1024)，在m = 1时，数据 χ (1024)搬移到 e (3072)；步骤203 启动IOM个线程完成待滤波数据X搬移到E的偶数块中。在此步骤中， 线程根据下列搬移公式进行数据搬移线程t ⑴需要搬移 e((2*(m+l)*21Q+i) = χ (m*2lcl+i)，其中 0 彡 i 彡 1023，m = 0， 1，. . .，210-1。例如线程t (0)在m = 0时，将数据χ (0)搬移到e (2048)，在m = 1时，数据 χ (1024)搬移到 e ^)96)；3)傅里叶变换核函数其处理是将经过重叠搬移后的数据，进行傅里叶变换计算，其中傅里叶变换的长度为 211 = 2048。4)乘法核函数其处理是将傅里叶变换的结果数据，与HR频响系数相乘。在此过程中，启动2048 个线程，并根据下列公式进行计算线程t(i)计算 p(m * 2n+i) = e(m * 2^+1)^(1)/211,其中，0 彡 i 彡 2047,m = 0，l，...，21(l-l，e为重叠搬移操作的结果数据。例如线程t(0)在m = 0时，计算p(0)= θ(0)*Η(0)/2"，在 m = 1 时，计算 W211) = e (211) *H(0)/211 ；5)傅里叶逆变换核函数其处理是将与频响系数相乘所得的结果数据，进行傅里叶逆变换计算，其中傅里叶逆变换的长度为211 = 2048。6)合并搬移核函数其处理是将傅里叶逆变换的结果数据，进行合并搬移。图3为通用的合并搬移示意图，将傅里叶逆变换的结果数据R= IR0, R1,..., Rk-J合并搬移为Y= {、，Y1,... ,YkJ， 其中氏={ZpYj，其中长度为M-I的连续数据点,Yi则是长度为L的连续数据点。而在本具体实例中，L = M-I = 210 = 1024，具体搬移操作参见图4。在此过程中，启动IOM个线程，并根据下列公式进行数据搬移操作线程t(i)需要搬移 y(m * 2lcl+i) = r ((2*m+l) *2lcl+i)，其中，0 彡 i 彡 1023，m = 0,1, ...,210-1, y为搬移后的结果数据。例如线程t(0)在m = 0时，将数据H21(l)搬移到y (0)，在m = 1时，将数据r (3072)搬移到y (1024)；上面六个核函数是条件，有了上述的六个可在GPU上并行执行的核函数之后，我们就可以按照下述步骤来完成FIR的并行处理，其中具体的滤波计算处理操作参见图5 步骤501 将1025点的HR冲击响应数据从主机内存拷贝到GPU内存中，启动频响计算核函数，完成FIR频响系数的计算。步骤502 将一块22° = IM点的数据从主机内存拷贝到GPU内存中，在GPU上启动 210个线程，用以执行重叠搬移核函数。步骤503 执行傅里叶变换核函数。步骤504 在GPU上启动211个线程，用以执行乘法核函数。步骤505 执行傅里叶逆变换核函数。步骤506 在GPU上启动21°个线程，用以执行合并搬移核函数。步骤507 将合并搬移的结果从GPU内存中拷贝到主机内存中，如果有待处理的数据，则转到步骤二执行。
权利要求
1. 一种基于GPU的高阶数字FIR滤波器频域并行处理实现方法，其特征在于根据通用图形处理器GPU众核体系结构的特点，将重叠保留方法的处理过程划分为六个核函数 频响计算核函数、重叠搬移核函数、傅里叶变换核函数、乘法核函数、傅里叶逆变换核函数、 合并搬移核函数；实现该方法的具体步骤如下步骤一、确定有限冲击响应FIR频率响应系数；当给定的响应系数为HR的冲击响应系数h= {h(0)，h(l)，......，h (M-I)}时，将FIR的冲击响应系数h = {h (0)，h (1)，......，h (M-I)}，经过尾部填0扩展至长度为N,将扩展后的系数传送给GPU，启动频响计算核函数对FIR的冲击响应系数h进行N点的傅里叶变换，结果为有限冲击响应FIR频率响应系数H = IH(O)，H(I)，......，H(N-I)}，保存在GPU内存中；当给定的响应系数为HR频率响应系数H= IH(O)，H(I)，......，H(N-I)}时，将频率响应系数H传给GPU，并保存在GPU内存中，保存时，第一个字节的地址为存储器位宽的整数倍；步骤二、将待处理的输入数据传送给GPU ；将一块长度为Nblk的待滤波样点数据X= {Β。』”...，Bk_J从主机内存中传入GPU的内存中，其中Nblk = k*L，k为整数，L为每个数据块的长度；其中Bi = ICi, DJ，0彡i彡k-1， Ci表示在Bi中起点为0，而长度为L-M+1的连续样点数据块，Di表示在Bi中起点为L-M+1， 而长度为M-I的连续样点数据块；步骤三、数据的重叠搬移；启动重叠搬移核函数完成数据重叠搬移操作，即将待滤波样点数据X = {B0, B1,...， Bk_J以及输入的长度为M-I的初始状态数据Stl，重叠搬移为E = {E0, E1, ...，Ek_J，其中 ={S0, BJ，其中，&为上一次处理待滤波样点数据X时其中的Dlri，对于i Φ 0的数据块Ei ={Dg，BJ，同时将01;_1搬移到中Stl，作为下一次处理过程的初始状态数据使用；重叠数据的搬移分为三个步骤，第一步，启动M-I个线程同时工作，线程 \(0彡i彡M-2)完成将&中第i个数据搬移到重叠搬移结果数据的&中的第i的位置，然后再完成将Dlri中第i个数据搬移到&中的第i的位置；第二步，启动L个线程同时工作，线程Ti (0彡i彡L-1)完成将待滤波样点数据X的Bj (0彡j彡k-Ι)块中第i个数据搬移到重叠搬移结果数据E的Bj (0彡j彡k-Ι)块中的第i的位置；第三步，启动M-I个线程同时工作，线程Ti (0彡i彡M-2)完成将待滤波样点数据X的Dj (0彡j彡k-Ι)块中第i 个数据搬移到重叠搬移结果数据E的DdO ^ j ^ k-1)块中的第i的位置；步骤四、滤波计算处理；启动傅里叶变换核函数，完成对每个Ei做N点的傅里叶变换运算，得到运算结果F = (FojF1,.. ·，Fk_J ；然后启动乘法核函数完成Fi与系数H相乘操作，即Ii = _ ；然后启动傅里叶逆变换核函数，完成对每个Ii做N点的傅里叶逆变换运算，并得结果R = {R0, R1,...，Rk-J ；在启动乘法核函数核函数时，启动N个线程同时工作，线程1\(0彡i彡N-1)完成将。(0 ^ j ^ k-1)中第i个数据与频响系数H中的第i个系数相乘，并将结果保存在 Rj (0 ( j ( k-1)块中的第i的位置；步骤五、数据的合并搬移；启动合并搬移核函数完成数据的合并搬移操作，即将结果数据R = {礼，R1, ...，Rk-J 合并搬移为Y = H · ·，Yk-J，其中Ri = IZi, YJ，其中Zi表示在Ri中起点为0，长度为 M-I的连续数据点所组成的数据块，Yi表示在氏中起点为M-1，长度为L的连续数据点所组成的数据块；在启动合并搬移核函数时，启动L个线程同时工作，线程Ti (0 ^ i ^ L-1)完成将乘积 Rj (0彡j彡k-Ι)中第M-1+i个数据搬移到结果数据Yj (0彡j彡k-Ι)块中的第i的位置， 其中，启动的L个线程在同一时刻搬运长度为L的连续数据块；步骤六、将合并搬移结果Y = {Yo，Y1,...，Yk-J传送到主机内存，若还有剩余数据需要处理重复步骤二，否则结束处理过程，完成高阶数字HR滤波器频域并行处理。
全文摘要
为解决高阶数字FIR滤波器并行处理效率的问题，提出了一种高效的、适合GPU体系结构的高阶数字FIR滤波器并行处理算法，该方法采用重叠保留方法结合GPU自身结构特点优化实现高阶数字FIR的频域并行化处理。通过计算FIR频率响应系数，将待处理的输入数据传送给GPU；数据重叠搬移；滤波计算处理；数据合并搬移；将合并搬移结果Y＝{Y0，Y1，....，Yk-1}传送到主机内存等步骤完成高阶数字FIR滤波器频域并行处理。对比在CPU上单线程所实现的FIR频域重叠保留方法，其吞吐率，即每秒处理样点的数量有着极大地提高，典型的加速比在100倍以上。
文档编号H03H17/02GK102340296SQ20111020494
公开日2012年2月1日 申请日期2011年7月21日 优先权日2011年7月21日
发明者宋昕, 张春宏, 汪晋宽, 韩英华, 高静 申请人:东北大学秦皇岛分校