专利摘要

本发明公开了一种基于粒子群算法的FIR滤波器设计方法，以FIR滤波器的单位冲激响应作为粒子群算法的粒子位置向量，基于FIR滤波器的滤波后信号的信噪比和FIR滤波器的均方误差来设计粒子群算法的适应度函数，迭代得到FIR滤波器。本发明基于信噪比和均方误差最小化准则设计粒子群算法的适应度函数，采用粒子群算法对FIR滤波器进行设计，从而得到兼顾信噪比性能的FIR滤波器。

权利要求

1.一种基于粒子群算法的FIR滤波器设计方法，其特征在于，包括以下步骤：

S1：设置粒子群算法中的参数，包括粒子群大小M、最大迭代次数Gmax和适应度阈值Fmin，其中粒子位置向量为FIR滤波器的单位冲激响应序列；

S2：初始化粒子群中M个粒子的位置和速度，以及FIR滤波器的均方误差参考值

S3：令迭代次数k＝1；

S4：采用每个粒子所对应的FIR滤波器对输入信号进行滤波，得到其滤波后信号，根据以下公式计算粒子群中每个粒子的适应度值F：

F＝SNR+E

其中，SNR表示滤波后信号的信噪比，表示粒子位置对应的单位冲激响应序列与理想滤波器的单位冲击响应序列的均方误差；

S5：更新每个粒子的局部最优位置和全局最优位置，以全局最优位置的均方误差对均方差参考值进行更新；

S6：如果k＝Gmax或Fbest＞Fmin，其中Fbest表示当前粒子群中全局最优位置的适应度值，则根据全局最优位置得到FIR滤波器的单位冲激响应序列，FIR滤波器设计完毕；否则更新粒子群中每个粒子的位置和速度，令k＝k+1，返回步骤S4。

2.根据权利要求1所述的FIR滤波器设计方法，其特征在于，所述步骤S2中粒子位置的初始化方法为：采用常规FIR滤波器设计方法对FIR滤波器进行设计，得到初始FIR滤波器，记其单位冲激响应序列为[h0(0),h0(1),h0(2)…h0(N-1)]，N表示FIR滤波器的阶数，每个粒子从[h0(0),h0(1),h0(2)…h0(N-1)]中随机选择Ki个单位冲激响应作为自身位置向量的对应元素，其余元素随机设置，Ki∈(0,N-1]，i＝1,2,…,M。

3.根据权利要求2所述的FIR滤波器设计方法，其特征在于，所述步骤S2中均方误差参考值的初始值为初始FIR滤波器对应的均方误差。

4.根据权利要求2所述的FIR滤波器设计方法，其特征在于，所述步骤S6中，粒子群算法迭代过程中各个粒子的速度上限vmax＜max(h0(n))。

说明书

技术领域

本发明属于滤波器设计技术领域，更为具体地讲，涉及一种基于粒子群算法的FIR滤波器设计方法。

背景技术

FIR(Finite Impulse Response，有限长单位冲激响应)滤波器又称为非递归型滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而该滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。

目前常用的FIR滤波器的设计方法有窗函数法、频率取样法和切比雪夫等波纹逼近的最优化设计方法等。在实际的工程应用中，FIR滤波器的设计实现，通常需要同时满足多个技术指标或达到较高的精度，设计工作比较复杂，而且只能是逼近工程的指标要求。那么此时FIR滤波器的设计其实就是寻找滤波器系数使频率响应不断的逼近理想的滤波器，可以看作一个最优化问题。求解最优化问题常用的算法有遗传算法、神经网络算法、粒子群算法等等。

在现有的基于最优化问题的FIR滤波器设计方法中，常用的最优化准则包括均方误差最小准则和最大误差最小化准则。其中均方误差最小准则是使FIR滤波器的频率响应误差的平方和最小，最大误差最小化准则是使设计的滤波器幅频响应最大逼近误差达到最小。这两种准则都是考虑的仅关注了要设计滤波器与理想低通滤波器的逼近程度，所考虑的因素较为单一，并未考虑FIR滤波器的其他性能，例如信噪比等，当输入信号所受干扰变化时，难以保证滤波效果。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于粒子群算法的FIR滤波器设计方法，基于信噪比和均方误差最小化准则设计粒子群算法的适应度函数，采用粒子群算法对FIR滤波器进行设计，从而得到兼顾信噪比性能的FIR 滤波器。

为实现上述发明目的，本发明基于粒子群算法的FIR滤波器设计方法包括以下步骤：

S1：设置粒子群算法中的参数，包括粒子群大小M、最大迭代次数Gmax和适应度阈值Fmin，其中粒子位置向量为FIR滤波器的单位冲激响应序列；

S2：初始化粒子群中的M个粒子的位置和速度，以及FIR滤波器的均方误差参考值

S3：令迭代次数k＝1；

S4：采用每个粒子所对应的FIR滤波器对输入信号进行滤波，得到其滤波后信号，根据以下公式计算粒子群中每个粒子的适应度值F：

F＝SNR+E

其中，SNR表示滤波后信号的信噪比， 表示粒子位置对应的单位冲激响应序列与理想滤波器的单位冲激响应序列的均方误差；

S5：更新每个粒子的局部最优位置和全局最优位置，以全局最优位置的均方误差对均方差参考值 进行更新；

S6：如果k＝Gmax或Fbest＞Fmin，其中Fbest表示当前粒子群中全局最优位置的适应度值，则根据全局最优位置得到FIR滤波器的单位冲激响应序列，FIR滤波器设计完毕；否则更新粒子群中每个粒子的位置和速度，令k＝k+1，返回步骤S4。

本发明基于粒子群算法的FIR滤波器设计方法，以FIR滤波器的单位冲激响应作为粒子群算法的粒子位置向量，基于FIR滤波器的滤波后信号的信噪比和FIR滤波器的均方误差来设计粒子群算法的适应度函数，迭代得到FIR滤波器。

本发明具有以下有益效果：

1)本发明基于FIR滤波器的滤波后信号的信噪比和FIR滤波器的均方误差来设计粒子群算法的适应度函数，可以很好地适应外界干扰，保证滤波后信号的信噪比；

2)在对粒子群算法中的粒子进行初始化时，可以采用常规设计方法设计得到初始FIR滤波器，根据初始FIR滤波器对粒子进行初始化，在已存在较好位置基础上寻找更佳位置，这样不仅可以加快搜索速度，还可以大大减少陷入局 部最优位置的概率。

附图说明

图1是本发明基于粒子群算法的FIR滤波器设计方法的流程图；

图2是本实施例中初始FIR低通滤波器的幅频特征曲线；

图3是本实施例中初始FIR低通滤波器与理想低通滤波器的冲击响应对比图；

图4是本实施例中输入信号波形图；

图5是图4所示输入信号经初始FIR低通滤波器滤波后的信号；

图6是本实施例中设计得到的FIR低通滤波器与初始FIR低通滤波器和理想低通滤波器的冲击响应对比图；

图7是图4所示输入信号经设计得到的FIR低通滤波器滤波后的信号。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

由于FIR滤波器设计可以看作一个最优化问题，因此本发明采用粒子群算法来实现FIR滤波器的设计，并基于对滤波后信号信噪比的考虑，对适应度函数进行了改进。图1是本发明基于粒子群算法的FIR滤波器设计方法的流程图。如图1所示，本发明基于粒子群算法的FIR滤波器设计方法包括以下步骤：

S101：设置粒子群算法参数：

首先需要对粒子群算法中的参数进行设置，包括粒子群大小M、最大迭代次数Gmax和适应度阈值Fmin。FIR滤波器的系统函数可以记为：

其中h(n)表示FIR滤波器的单位冲激响应，N表示FIR滤波器的阶数。因此粒子群中的每个粒子的位置向量为FIR滤波器的单位冲激响应序列。

本实施例中，采用N阶线性相位FIR低通滤波器，其单位冲激响应h(n)为实数，且以对称中心n＝(N-1)/2对称，因此有h(n)＝h(N-1-n)。因此本实施例中粒子群中的每个粒子，其位置向量的格式为[h(0),h(1),h(2)…h((N-1)/2)]，即只取单位冲激响应序列中的前一半。

根据粒子群算法的流程可知，在算法执行过程中每个粒子的位置和速度都会进行更新，自然还需要对粒子群算法自身的一些参数进行设置，包括加速度常数c1和c2，惯性权重ω。本实施例中粒子群大小M＝100，最大迭代次数Gmax＝1000，适应度阀值Fmin＝110dB，加速度常数c1＝c2＝0.2，惯性权重ω＝0.4。

S102：初始化粒子：

初始化粒子群中的M个粒子的位置和速度，以及FIR滤波器的均方误差参考值

在粒子群算法中，粒子的位置和速度的初始值通常是随机获取的。为了提高粒子群算法的效率，需要尽可能设置一个合适的初始值。本实施例中，采用常规FIR滤波器设计方法对FIR滤波器进行设计，得到初始FIR滤波器，记其单位冲激响应序列为[h0(0),h0(1),h0(2)…h0(N-1)]，每个粒子从[h0(0),h0(1),h0(2)…h0(N-1)]中随机选择Ki个单位冲激响应作为自身位置向量的对应元素，其余元素随机设置，Ki∈(0,N-1]，i＝1,2,…,M。每个粒子对应的Ki可以相同，也可以不同，显然当Ki＝N-1时，所有初始粒子的位置都是相同的。基于初始FIR滤波器，均方误差参考值 的初始值为初始FIR滤波器对应的均方误差，粒子群算法迭代过程中各个粒子的速度上限vmax＜max(h0(n))。

FIR滤波器的常规FIR滤波器设计方法通常有三种主要的方法：窗函数法、频率取样法和切比雪夫等波纹逼近的最优化设计方法。本实施例中FIR滤波器以FIR低通滤波器为例，由于窗函数中的凯泽窗可以调节主半宽度与旁瓣衰减，适应性较强，而且其阻带衰减大，通带衰减很小，因此本实施例采用窗函数中的凯泽窗来进行初始FIR滤波器的设计。

理想线性低通滤波器的传输函数为：

其中，e表示自然常数，ω表示频率，ωc表示截止频率，τ为群延。

其单位冲击响应为：

其时域的单位hd(n)为无限长。窗函数设计法的基本思路是用一定宽度窗函数截取理想情况下的脉冲响应，忽略窗函数宽度以外的序列，这样得到的脉冲响便是有限长的。设窗函数序列为ω(n)，将hd(n)加以截断后的序列为h(n)：

h(n)＝hd(n)×ω(n) 0≤n≤N-1 (4)

加窗后的实际频率响应H(ejω)为：

窗函数法设计FIR滤波器就是寻找N和窗的形状ω(n)使加窗后的实际频率响应H(ejω)满足设计要求。

凯泽窗函数形式为：

其中，I0(·)是第一类零阶变型贝塞尔函数，RN(n)表示Lagrange型余项。β是一个可调的参数，用来调节主瓣宽度和旁瓣衰减，一般说来，β可以调节窗的形状，β越大，窗函数主瓣宽度增加，旁瓣衰减也就越大。β可采用下面公式来计算：

其中As为阻带衰减。

窗长度N的计算公式为：

其中过渡带角频率Δω＝2πΔf/fs，Δf表示过滤带频率，fs表示输入信号的采样频率。

本实施例中输入信号的采样频率fs＝1220Hz，FIR滤波器的截止频率fc＝1Hz，此处选择Δf＝5Hz，设置阻带衰减As＝80dB。则凯泽窗参数计算结果如下：

β＝0.1102×(As-8.7)＝0.1102×(80-8.7)＝7.857 (9)

为了使粒子群优化后的滤波器满足线性相位条件取N＝1227，于是凯泽窗函数为：

图2是本实施例中初始FIR低通滤波器的幅频特征曲线。如图2所示，可以看出初始FIR滤波器的阻带衰减比较大。图3是本实施例中初始FIR低通滤波器与理想低通滤波器的冲击响应对比图。如图3所示，本实施例中采用凯泽窗设计得到的初始FIR低通滤波器，与理想的低通滤波器仍然存在一定误差，其均方误差为：

因此本实施例中

图4是本实施例中输入信号波形图。如图4所示，该输入信号的信噪比为66.76dB。图5是图4所示输入信号经初始FIR低通滤波器滤波后的信号。如图5所示，经初始FIR低通滤波器滤波后，信号的信噪比提高到106.27dB。

S103：令迭代次数k＝1。

S104：计算适应度值：

滤波器设计有两种最优化准则，一种是均方误差最小准则，另一种是最大误差最小化准则。由于最大误差最小化准则比较复杂，本发明基于滤波后信号的信噪比，参考均方误差最小准则，修改适应度函数。下面先分别对于信噪比 和均方误差最小准则进行说明。

●信噪比

信噪比用来衡量信号中有用信号功率与噪声功率的比例，其计算公式为：

其中Vs为信号有效值，Vn为噪声有效值。Vs、Vn难以计算，对于直流信号，信噪比一般用下面公式计算：

其中，x(n)表示FIR滤波器滤波后信号，滤波后的信号信噪比越高，表明滤波方法越有效。

●均方误差最小准则

均方误差最小准则原则上是使FIR滤波器与理想滤波器频率响应的误差的平方和最小。以Hd(ejω)表示所要设计的滤波器理想的频率响应，H(ejω)表示设计得到的滤波器频率响应，Ed(ejω)表示所要设计的滤波器理想的频率响应与实际响应的误差，则：

Ed(ejω)＝Hd(ejω)-H(ejω) (15)

均方误差(误差的平方和的平均值)为：

均方误差最小准则就是找到一组h(n)使e2最小。将 代入式(16)得：

其中

根据帕赛瓦 其中X(eiω)为x(n)的傅里叶变换，均方误差可以用下式表示：

只与理想滤波器特性有关，与设计的实际滤波器特性无关，因此在实际计算均方误差时，其计算公式采用：

因此，本发明基于信噪比和均方误差最小准则，设计得到粒子群算法的适应度函数表示如下：

F＝SNR+E (20)

其中，F表示适应度值，SNR表示滤波后信号的信噪比， 表示粒子位置对应的单位冲激响应序列与理想滤波器的单位冲击响应序列的均方误差， 表示FIR滤波器的均方误差参考值。可见，在本发明中，适应度值越大越好。

从适应度函数F的表达式可以看出，如果粒子新位置的均方差满足 则是以新位置作为滤波器系数滤波后的信号的信噪比作为评价粒子位置的优劣，当粒子新位置的均方差不满足 时，则该粒子的适应度为无穷小，即忽略该位置，不作为个体最优位置、全局最优位置的考虑。在实际应用中，-∞一般以一个绝对小值来代替即可。

S105：更新最优位置和均方误差参考值：

更新每个粒子的局部最优位置和全局最优位置，以全局最优位置的均方误差对均方差参考值 进行更新，即以当前粒子群中全局最优位置的均方误差作 为新的均方差参考值

局部最优位置是指粒子在历史迭代过程中适应度值最大的粒子位置，全局最优位置是指历史迭代过程中所有粒子适应度值最大的粒子位置。

S106：判断是否k＝Gmax或Fbest＞Fmin，其中Fbest表示当前粒子群中全局最优位置的适应度值，如果是，说明已经达到迭代结束条件，进入步骤S108，否则进入步骤S107。

S107：更新粒子：

更新粒子群中每个粒子的位置和速度，令k＝k+1，返回步骤S104。

在粒子群算法中，记第i个粒子在第k次迭代的位置 其中 表示第i个粒子第k次迭代位置的第q个元素，q＝1,2,…,Q。记第i个粒子的局部最优位置pbesti＝[pbesti1,pbesti2,pbesti3…,pbestiQ]，其中pbestiq表示第i个粒子局部最优位置的第q个元素。记全局最优位置gbest＝[gbest1,gbest2,gbest3…,gbestQ]，其中gbestq表示粒子群最优位置的第q个元素。记第i个粒子在第k次迭代的速度 其中 表示第i个粒子第k次迭代速度的第q个元素。那么粒子的位置和速度的更新公式如下：

在粒子群算法中，粒子的位置和速度一般会存在一个限定范围。

S108：得到FIR滤波器：

根据全局最优位置得到FIR滤波器的单位冲激响应序列，FIR滤波器设计完毕。

图6是本实施例中设计得到的FIR低通滤波器与初始FIR低通滤波器和理想低通滤波器的冲击响应对比图。如图6所示，采用本发明优化后的FIR低通滤波器，其冲击响应更接近于理想低通滤波器的冲击响应，明显优于采用凯泽窗所设计得到的初始低通滤波器。

图7是图4所示输入信号经设计得到的FIR低通滤波器滤波后的信号。对比图5和图7，可知采用本发明设计得到的FIR低通滤波器，较采用凯泽窗设计得到的FIR低通滤波器得到的滤波后信号波动范围更小，其信噪比提高到 114.72dB，明显优于初始FIR低通滤波器。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

基于粒子群算法的FIR滤波器设计方法专利购买费用说明

Q：办理专利转让的流程及所需资料

A：专利权人变更需要办理著录项目变更手续，有代理机构的，变更手续应当由代理机构办理。

1：专利变更应当使用专利局统一制作的“著录项目变更申报书”提出。

2：按规定缴纳著录项目变更手续费。

3：同时提交相关证明文件原件。

4：专利权转移的，变更后的专利权人委托新专利代理机构的，应当提交变更后的全体专利申请人签字或者盖章的委托书。

Q:专利著录项目变更费用如何缴交

A：（1）直接到国家知识产权局受理大厅收费窗口缴纳,（2）通过代办处缴纳，（3）通过邮局或者银行汇款,更多缴纳方式

Q：专利转让变更，多久能出结果

A：著录项目变更请求书递交后，一般1-2个月左右就会收到通知，国家知识产权局会下达《转让手续合格通知书》。