本发明属于数字信号处理领域，具体涉及一种基于软件查表法的fir滤波器实现方法。
背景技术：
在数字通信理论中，数字基带信号需要通过基带内插成型滤波后才能够进行调制输出。基带成型滤波器一般采用fir滤波器，常用的fir滤波器的实现方法分为硬件和软件两种。其中，硬件实现包括利用dsp或fpga来实现，软件实现是指借助通用cpu来实现。
在硬件实现方面，dsp芯片实现的fir滤波器具有稳定性好，精确度高以及不受环境影响的特点，并且dsp支持流水线操作，能够实现并行处理，进而能够提高fir滤波器的运算速度。fpga能够并行处理数据并具有可编程逻辑性，其运算速度快、实时性高、程序移植性强、可靠性强、实现成本较低，因此目前是fir滤波器设计的主要方法。
在软件实现方面，通用cpu具有远高于fpga的主频，能够进行高速运算，但是因为通用cpu采用串行方式执行指令，fir滤波器的大量乘累加计算会使得整体运算速度下降。目前还没有能够在通用cpu上快速实现fir滤波器的较好算法。
技术实现要素：
本发明的目的在于针对上述现有fir滤波器cpu实现技术中的问题，提供一种基于软件查表法的fir滤波器实现方法，利用计算机的大容量内部存储器来存储fir滤波器的查找表，并通过查表运算代替卷积运算，以计算机自身的ram资源换取整体计算速度的提高。
为了实现上述目的，本发明采用的技术方案包括以下步骤：
步骤一、计算查找表的大小并与设定阈值进行比较，阈值指的是实际应用中所能接受的最大查找表大小，若查找表的大小不超过设定阈值即生成查找表，否则修改滤波器参数；
步骤二、生成查找表，查找表采用静态生成法，特定fir滤波器对应一个固定的查找表；
2.1)根据设计指标得到滤波器的各抽头系数h(n)和抽头数总数nt；
2.2)根据滤波器的阶数n，确定在每次卷积过程中有效的码元位数k；
2.3)码元序列s(n)的长度设为k位，s(n)依次按照集合s＝{0,1,2,…,2^k-1}取值；
2.4)码元序列s(n)根据内插倍数m在两个码元之间插入m-1个0，生成输入序列x(n)；
2.5)根据离散信号的卷积公式：
2.6)将每个s(n)对应的y(n)存入查找表，回到步骤2.3)，直至s(n)取完s中的所有值；
步骤三、计算机根据实时输入码元流计算索引值，读取查找表得到对应的滤波器输出值。
所述查找表的大小根据以下公式进行计算：size＝m*2p*n*l(byte)；
式中，p为输入码元分辨率，l为表格数据精度，n为滤波器阶数，m为内插倍数，size为查找表大小，查找表大小的单位为字节；所述的输入码元分辨率p指的是每个码元的比特数，表格数据精度l指的是查找表中每个数据以l个字节进行存储。
步骤2.1)所述的设计指标包括滤波器的阶数n、内插倍数m、采样频率、通带和阻带截止频率、通带和阻带衰减值等；所述的各抽头系数值h(n)通过滤波器设计软件filterlab、filtercad获得，或者通过matlab的fdatool工具得到。
步骤2.2)确定有效的码元位数k时，分有效码元与滤波器抽头中心点对齐和不对齐两种情况进行计算：
2.2a)当有效码元与滤波器抽头中心点对齐时，输入序列的首尾两点为有效码元，则码元有效位数k'＝n；
2.2b)当有效码元与滤波器抽头中心点不对齐时，滤波器的抽头中心点与内插0值对齐，则码元的有效位数k＝n-1。
步骤2.3)针对有效的码元位数k分别计算两种情况下的s(n)取值：
2.3a)当有效码元与滤波器抽头中心点对齐时，集合s'＝{0,1,2,…,2^k'-1}，k'＝n；s(n)在集合s'中依次取值，计算出的各输出值存入查找表的第一行；
2.3b)当滤波器的抽头中心点与内插0值对齐时，集合s＝{0,1,2,…,2^k-1}，k＝n-1；s(n)在集合s中依次取值，计算出的各输出值存入查找表的剩余m-1行。
步骤2.4)实现内插的方法为：码元序列s(n)的每一比特依次和间隔m位的滤波器抽头系数h(n+m)相乘后累加；若内插倍数为1，则x(n)＝s(n)。
步骤2.5)针对有效的码元位数k分别计算两种情况下的滤波器输出值：
2.5a)当有效码元与滤波器抽头中心点对齐时：
2.5a-1)将滤波器抽头系数h(n)按照内插倍数m间隔取值，存入数组h_center中，即h_center＝{h(0),h(m),…,h(nt-m-1),h(nt-1)}；
2.5a-2)s(n)从零开始取值，按位与h_center数组中的每个元素相乘并累加；
2.5a-3)计算出对应的y(n)后存入查找表的第一行第一列中，然后s(n)加1，并计算对应的y(n)，直到s(n)取完集合s'中的所有值，填满查找表的第一行；
2.5b)当滤波器的抽头中心点与内插0值对齐时：
2.5b-1)由滤波器的对称性得到h(0)＝h(nt-1)，仅取前nt-1个h(n)，按照m间隔取值，分成m个有k个元素的h_sample数组，只取后m-1个数组；
2.5b-2)设h_sample数组集合为{h_sample(m),1≤m≤m-1}，数组内每个元素的取值如下：h_sample(m)＝{h(m),h(m+m),…,h(nt-m-1+m)}；
2.5b-3)m首先取值为1，则h_sample(1)＝{h(1),h(1+m),…,h(nt-m)}；
2.5b-4)s(n)从零开始取值，按位与h_sample(m)数组中的每个元素相乘并累加；
2.5b-5)计算出对应的y(n)存入查找表中第m+1行第1列，然后s(n)加1，并计算对应y(n)，直到s(n)取完集合s中的所有值，填满查找表的第m+1行；
2.5b-6)m加1，返回到步骤2.5b-4)，如此循环，直至m＝m-1，计算完m-1个数组对应的y(n)，填满查找表的剩余m-1行。
步骤2.6)生成的查找表为m行2^k'列，第一行是有效码元与滤波器抽头中心点对齐时，输入码元所有情况的输出值；其余m-1行是滤波器抽头中心点与内插0值对应时，输入码元所有情况的输出值，但此m-1行只有2^k列的数据，根据滤波器系数的对称性，前2^k列数值与后2^k列数值相等，故将前2^k列数值翻转并复制到2^k+1至2^k’列，补充完查找表。
所述的步骤三读取查找表得到对应滤波器输出值的方法为：将输入码元序列按照2^k'取余加1作为索引值index，读取查找表中第index列的所有行的数值作为滤波器的输出值。
与现有技术相比，本发明基于查表法，以查表取代卷积运算，通过计算机生成并读取查找表完成了滤波功能。用内存资源换取滤波速度，克服了通用cpu在实现fir滤波器时，由于串行执行导致运算次数大，计算速率低的缺点。利用了现有计算机大内存的优势，通过计算机生成并存储fir滤波器的大查找表，配合cpu的高速运算，整体提高了运算速度。fir滤波器的抽头数主要由滤波器的阶数和码元序列的内插倍数决定，阶数越大，内插倍数越高，滤波器的抽头数越多。高阶fir滤波器与低阶fir滤波器相比，具有更窄的过渡带，阻带衰减快，滤波效果好。高内插fir滤波器和普通fir滤波器相比，具有更多的波形采样点，有利于接收端的波形重构，降低误码率；同时采样速率的提高，压缩了频谱，能够进行频分复用，提高了信道利用率。高阶和高内插fir滤波器的优点使其应用广泛，但由于运算量大，普通的软件法可实现的滤波器阶数和码元内插倍数有限。一个20阶的2倍内插fir滤波器，若采用卷积计算，每次滤波需要进行20次的乘法运算和19次的加法运算，在本发明中只需要进行一个计算索引值的运算，然后根据索引值查表即可完成滤波功能。本发明基于查表法，以查表运算替代卷积运算，同时利用计算机的大容量内存存储fir滤波器的大查找表，以内存换取速度，能够实现普通fir滤波器的功能。本发明充分考虑了当前计算机内存容量不断增大、cpu速度逐步提高的发展趋势，具有强大的可延展性和广阔的发展前景。
进一步的，在本发明中，实现内插的方法为，码元序列s(n)的每一比特依次和间隔m位的滤波器抽头系数h(n+m)相乘并累加，若内插倍数为1，则x(n)＝s(n)。如此即可不再另外生成输入序列x(n)，达到减少数据存储量以及减少乘法和累加次数的效果，从而提高运算速度。
具体实施方式
下面结合具体实施例对本发明做进一步的详细说明。
本发明是fir滤波器的查表法实现，重点是用查找表取代卷积运算，提高运算速度。查找表的大小与滤波器阶数以及码元内插倍数两个因素都有关，但因为滤波器阶数的增加会使得查找表大小成指数性增大，而内插倍数和查找表大小是线性关系，所以相比之下，本发明更适合高内插fir滤波器的实现。下面对本发明采用的技术方案进行详细描述：
a.计算查找表大小
根据输入码元分辨率p、表格数据精度l、滤波器阶数n和内插倍数m，计算该滤波器的查找表大小size(单位：字节)。输入码元的分辨率p指的是每个码元的比特数，表格数据精度l指的是查找表中每个数据以l个字节存储，阈值thr根据实际应用设定，指的是所能接受的最大查找表大小。查找表大小根据以下公式进行计算：
size＝m*2p*n*l(byte)
若size≤thr，进入下一步生成查找表；否则，修改滤波器参数。
b.生成查找表
2.1)根据设计要求，如滤波器的阶数n、内插倍数m、采样频率、通带和阻带截止频率、通带和阻带衰减值等设计指标，得到指定滤波器的各抽头系数h(n)和抽头总数nt。可利用filterlab、filtercad等滤波器设计软件，或者matlab的fdatool工具得到各抽头系数h(n)。
2.2)根据滤波器阶数n，确定在每次卷积过程中有效的码元位数k。计算方法为：分析计算卷积的乘累加过程，由于输入序列采取m倍内插，考虑有效码元与滤波器抽头中心点对齐和不对齐的两种情况：
2.2a)当有效码元与滤波器抽头中心点对齐时，输入序列的首尾两点正好是有效码元，则码元有效位数k'＝n；
2.2b)当有效码元与滤波器抽头中心点不对齐时，滤波器的抽头中心点与内插0值对齐，则码元的有效位数k＝n-1。
2.3)码元序列s(n)的长度设为k位，s(n)依次按照集合s取值。根据步骤2.2)中两种情况下k的取值，分别计算对应的s(n)取值：
2.3a)当有效码元与滤波器抽头中心点对齐时，集合s'＝{0,1,2,…,2^k'-1}，k'＝n；s(n)在集合s'中依次取值，计算出的各输出值存入查找表的第一行；
2.3b)当滤波器的抽头中心点与内插0值对齐时，集合s＝{0,1,2,…,2^k-1}，k＝n-1；s(n)在集合s中依次取值，计算出的各输出值存入查找表的剩余m-1行。
2.4)码元序列s(n)根据内插倍数m在两个码元之间插入m-1个0，生成输入序列x(n)。实现内插的方法为：s(n)的每一比特依次和间隔m位的滤波器抽头系数h(n+m)相乘后累加；若内插倍数为1，则x(n)＝s(n)。如此可不再另外生成输入序列x(n)，达到减少数据的存储量以及减少乘法和累加次数的效果，从而提高运算速度。
2.5)根据离散信号的卷积公式：
2.5a)当有效码元与滤波器抽头中心点对齐时：
2.5a-1)将滤波器抽头系数h(n)按照内插倍数m间隔取值，存入数组h_center中，即h_center＝{h(0),h(m),…,h(nt-m-1),h(nt-1)}；
2.5a-2)s(n)从零开始取值，按位与h_center数组中的每个元素相乘并累加；
2.5a-3)计算出对应的y(n)后存入查找表的第一行第一列中；然后s(n)加1，并计算对应的y(n)，直到s(n)取完集合s'中的所有值，填满查找表的第一行。
2.5b)当滤波器的抽头中心点与内插0值对齐时：
2.5b-1)由滤波器的对称性知，h(0)＝h(nt-1)。仅取前nt-1个h(n)，按照m间隔取值，分成m个有k个元素的h_sample数组，只取后m-1个数组；
2.5b-2)设h_sample数组集合为{h_sample(m),1≤m≤m-1}，数组内每个元素的取值如下：h_sample(m)＝{h(m),h(m+m),…,h(nt-m-1+m)}；
2.5b-3)m首先取值为1，则h_sample(1)＝{h(1),h(1+m),…,h(nt-m)}；
2.5b-4)s(n)从零开始取值，按位与h_sample(m)数组中的每个元素相乘累加；
2.5b-5)计算出对应的y(n)存入查找表中第m+1行第1列，然后s(n)加1，并计算对应的y(n)，直到s(n)取完集合s中的所有值，填满查找表的第m+1行；
2.5b-6)m加1，回到步骤2.5b-4)如此循环，直至m＝m-1，计算完m-1个数组对应的y(n)，填满查找表的剩余m-1行。
2.6)生成的查找表为m行2^k'列，第一行是有效码元与滤波器抽头中心点对齐时，输入码元所有情况的输出值；其余m-1行是滤波器抽头中心点与内插0值对应时，输入码元所有情况的输出值，但此m-1行只有2^k列的数据，由滤波器系数的对称性知，前2^k列数值与后2^k列数值相等，故将前2^k列数值翻转并复制到2^k+1至2^k'列，补充完查找表。
c.读取查找表
计算机根据实时输入码元流计算索引值，读取查找表得到对应滤波器输出值，完成滤波功能。根据查找表读取滤波器输出值的方法是：将输入码元序列按照2^k'取余加1作为索引值index，读取查找表中第index列的所有行的数值作为滤波器的输出值。
下面以一个实例，对码元序列采取m倍内插情况下的整个发明过程做进一步的详细说明：
设计要求：码元采用qpsk调制，码元速率为19kbps，内插倍数为64倍；采用滚降系数为0.4的根升余弦滤波器，滤波器阶数为17阶，通带衰减在1db以内，阻带衰减为40db，查找表数据采用short类型存储，查找表大小在1gb内。
1.根据查找表大小计算公式，计算出查找表大小size＝16.78mb，满足设计要求。
2.根据设计要求，计算可知滤波器的采样频率为64*19khz。查找表中每个数据以2个字节存储，即数据精度l＝2。虽然码元采用qpsk调制，但由于qpsk的两个码元是正交的，所以可将两个码元单独分析，认为每个码元采用bpsk调制，则码元分辨率p＝1。
3.为加快运算速度，简化程序，选择在matlab平台下通过矩阵运算完成本次滤波器设计；调用matlab的rcosfir函数，得到指定根升余弦滤波器抽头系数矩阵h，抽头总数为1025。
4.由于h是一个1*1025的矩阵，为了便于进行卷积运算，首先按照内插倍数64的间隔取出滤波器抽头系数，组成1*17的矩阵h_center。然后，取矩阵h前1024个数，再依次按照内插倍数64间隔取值，生成64*16的矩阵h_sample_1。接着，根据卷积过程，将整个h_sample_1矩阵进行上下翻转生成矩阵64*16的矩阵h_sample_2，即h_sample_1的第64行是h_sample_2的第1行，h_sample_1的第63行是h_sample_2的第2行，以此类推，最后，仅取h_sample_2的前63行生成63*16的矩阵h_sample；
5.将64*2^17的查找表矩阵置零，按照有效码元与滤波器抽头中心点对齐、内插0值与滤波器抽头中心点对齐两种情况分别计算滤波器输出值。
注：有效码元采用bpsk调制，即-1→0，1→1。
5a)当有效码元与滤波器抽头中心点对齐时：
5a-1)有效码元位数k'＝n＝17，集合s'＝{0,1,2,…,2^17-1}，令s(n)＝0；
5a-2)按位取出s(n)的17位数值，该位若为0，对应的h_center值乘上-1加到查找表矩阵对应位，该位若为1，对应的h_center值乘上1加到查找表矩阵对应位。即，若s(1)为0，减去h_center(1)的值，为1则加上h_center(1)的值；若s(2)为0，减去h_center(2)的值，为1则加上h_center(2)的值，以此类推；
5a-3)计算s(n)每一位和h_center矩阵各数值相乘的结果，累加后作为查找表矩阵第一行第一列的值；
5a-4)s(n)加1，回到步骤5a-2)，按位依次计算和h_center矩阵各数值的乘积并累加，填入查找表第一行的下一列；以此循环，直至s(n)取完集合s'的值，填完查找表矩阵的第一行。
5b)当内插0值与滤波器抽头中心点对齐时：
5b-1)有效码元位数k＝n-1＝16，集合s＝{0,1,2,…,2^16-1}，令s(n)＝0，取出h_sample的第一行，即h_sample(1,x)；
5b-2)按位取出s(n)的16位的数值，该位若为0，对应的h_sample(1,x)值乘上-1加到查找表矩阵对应位，该位若为1，对应的h_sample(1,x)值乘上1加到查找表矩阵对应位置。即，若s(1)为0，减去h_sample(1,1)的值，为1则加上h_sample(1,1)的值；s(2)为0，减去h_sample(1,2)的值，为1则加上h_sample(1,2)的值，以此类推；
5b-3)计算s(n)每一位和h_sample(1,x)各数值相乘的结果，累加后作为查找表矩阵第二行第一列的值；
5b-4)s(n)加1，回到步骤5b-2)，按位依次计算和h_sample(1,x)各数值的乘积并累加，填入查找表矩阵第二行的下一列；以此循环，直至s(n)取完集合s的值，填满查找表矩阵的第二行前2^16列，进入下一步；
5b-5)取出h_sample的下一行，即h_sample(2,x)，回到步骤5b-1)，直至取完h_sample的所有63行，进入下一步；
5b-6)至此，查找表矩阵的第2至64行的前2^16列都已填满。根据滤波器的对称性，将前2^16行数据翻转复制至后2^16行，填满整个查找表矩阵。
6.将查找表矩阵写入指定文件中，计算机根据实时码流计算查表索引值，读取查找表得到滤波器输出值，将输入码元流按照2^17取余加1作为索引值index，读取查找表中第index列的64个值作为滤波器的输出值。
以一个20阶的2倍内插fir滤波器为例，若采用卷积计算，每次滤波时需要进行20次的乘法运算和19次的加法运算。在本发明中，只需要进行一次计算索引值的运算，然后根据索引值查表即可完成滤波功能。本发明充分考虑了当前计算机内存容量不断增大、cpu速度逐步提高的发展趋势，具有强大的可延展性和广阔的发展前景。