CN106372342A - 一种基于遗传算法的高阶数字微分器设计方法 - Google Patents

Abstract

本发明公开了一种基于遗传算法的高阶数字微分器设计方法，该方法通过不断进行遗传算法的选择、交叉和变异等操作，在有效频段内最小化误差的平方值，从而求解出对应的线性相位FIR微分器系数。该方法简单易行，鲁棒性强，结果令人满意。

Description

一种基于遗传算法的高阶数字微分器设计方法

技术领域

本发明涉及数字滤波器设计领域，更具体地，涉及一种基于遗传算法的高阶数字微分器设计方法。

背景技术

信号导数的求取是一个普遍存在的问题，并且很多工程应用当中都会遇到该问题，如防御系统中的导弹拦截系统以及时下热点的高速时间交织模数转换器(TI-ADC)校正系统等。用解析数学表达式表示的信号，如基础函数等，能够用数学方法求导。但在大多数情况下，信号是没有数学表达式的，所以不能用数学方法直接求导，通常采用差分方法来近似估计信号的导数。通常情况下，由于系统中存在失配或信号中存在噪声，所以通过差分方法不能正确的估计信号的导数。因此，构造微分器是不可避免的，而高阶微分器的设计更一直是一个比较困难的问题。

McClellan提出过一个微分器设计的计算机程序，但只能适用于一阶微分器设计。虽然之后研究人员对此进行多次改进，但计算过程往往非常复杂，有时甚至无法收敛导致错误。

本发明采用遗传算法对数字微分器的系数进行寻优。遗传算法是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。该算法是一类借鉴生物界的进化规律(适者生存，优胜劣汰遗传机制)演化而来的随机化搜索方法。它是由美国的J.Holland教授1975年首先提出，其主要特点是直接对结构对象进行操作，不存在求导和函数连续性的限定；具有内在的隐并行性和更好的全局寻优能力；采用概率化的寻优方法，能自动获取和指导优化的搜索空间，自适应地调整搜索方向，不需要确定的规则。遗传算法的这些性质，已被人们广泛地应用于组合优化、机器学习、信号处理、自适应控制和人工生命等领域。它是现代有关智能计算中的关键技术。

发明内容

本发明为克服上述现有技术所述的至少一种缺陷，提供一种基于遗传算法的高阶数字微分器设计方法。该方法简单易行，鲁棒性强，结果令人满意。

为解决上述技术问题，本发明的技术方案如下：

一种基于遗传算法的高阶数字微分器设计方法，包括以下步骤：

S1：确定高阶数字微分器的阶数k、长度N以及截止频率ω_p，高阶数字微分器的期望传递函数为：

$D\left(\mathrm{\ω}\right)=\left\{\begin{array}{cc}{\left(\frac{j\mathrm{\ω}}{2\mathrm{\π}}\right)}^k,& 0\≤\mathrm{\ω}\≤{\mathrm{\ω}}_p\\ -{\[\frac{j(2\mathrm{\π}-\mathrm{\ω})}{2\mathrm{\π}}\]}^k,& 2\mathrm{\π}-{\mathrm{\ω}}_p\≤\mathrm{\ω}\≤2\mathrm{\π}\end{array}\right.$

其中，k代表高阶数字微分器的阶数，当k为偶数时，为偶数阶数字微分器，当k为奇数时为奇数阶数字微分器；ω_p为高阶数字微分器需要作用的截至频率，而在ω_p之外的频带，高阶数字微分器的特性不需要关心；

高阶数字微分器的系统函数表示为：

$H\left(z\right)=\underset{n=0}{\overset{N-1}{\Σ}}h\left(n\right)z^{-n}$

其中h(n)为对应的冲击响应；

S2：根据k、N的奇偶性选择合适类型的FIR滤波器模型：

当k为偶数时、N为奇数时，设计为第一类FIR滤波器；

当k为偶数时、N为偶数时，设计为第二类FIR滤波器；

当k为奇数时、N为奇数时，设计为第三类FIR滤波器；

当k为奇数时、N为偶数时，设计为第四类FIR滤波器；

S3：建立适应度函数表达式；

S4：遗传算法初始化，建立初始种群；

S5：采用遗传算法获得寻优结果向量B，寻优结果向量B为使适应度函数取最小的B向量；

S6：根据向量B计算为微分器系数。

在一种优选的方案中，步骤S2中，当k为偶数时、N为奇数时，设计为第一类FIR滤波器， 表示向下取整；此时：

H(z)的频率响应函数表示为：

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{H}_e\left(\mathrm{\ω}\right);$

$H_e\left(\mathrm{\ω}\right)=\underset{n=0}{\overset{(N-1)/2}{\Σ}}b\left(n\right)cos\left(n\mathrm{\ω}\right)$

其中：

$b\left(n\right)=\left\{\begin{array}{cc}h\left(\frac{N-1}{2}\right)& n=0\\ 2h(\frac{N-1}{2}-n)& n=1,2,...,\frac{N-1}{2}\end{array}\right.$

定义：

$B={\[b\left(1\right),b\left(2\right),...,b\left(\frac{N}{2}\right)\]}^T$

$C\left(\mathrm{\ω}\right)={\[1,cos\left(\mathrm{\ω}\right),...,cos\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_e(ω)可表示为H_e(ω)＝B^TC(ω)。

在一种优选的方案中，步骤S2中，当k为偶数时、N为偶数时，设计为第二类FIR滤波器， 表示向下取整；此时：

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{H}_e\left(\mathrm{\ω}\right);$

$H_e\left(\mathrm{\ω}\right)=\underset{n=1}{\overset{N/2}{\Σ}}b\left(n\right)cos\left(\frac{2n-1}{2}\mathrm{\ω}\right)$

其中：

$b\left(n\right)=2h(\frac{N}{2}-n),n=1,2,...,\frac{N}{2}$

定义：

$B={\[b\left(1\right),b\left(2\right),...,b\left(\frac{N}{2}\right)\]}^T$

$C\left(\mathrm{\ω}\right)={\[cos\left(\frac{\mathrm{\ω}}{2}\right),cos\left(\frac{3}{2}\mathrm{\ω}\right),...,cos\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_e(ω)可表示为H_e(ω)＝B^TC(ω)。

在一种优选的方案中，步骤S2中，当k为奇数时、N为奇数时，设计为第三类FIR滤波器，此时：

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{e}^{\frac{j\mathrm{\π}}{2}}{H}_o\left(\mathrm{\ω}\right);$

$H_o\left(\mathrm{\ω}\right)=\underset{n=1}{\overset{(N-1)/2}{\Σ}}b\left(n\right)sin\left(n\mathrm{\ω}\right)$

其中：

$b\left(n\right)=2h\left(\frac{N-1}{2}-n\right),n=1,2,\mathrm{...},\frac{N-1}{2}$

定义：

$B={\[b\left(1\right),b\left(2\right),...,b\left(\frac{N-1}{2}\right)\]}^T$

$S\left(\mathrm{\ω}\right)={\[sin\left(\mathrm{\ω}\right),sin\left(2\mathrm{\ω}\right),...,sin\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_O(ω)可表示为H_O(ω)＝B^TS(ω)。

在一种优选的方案中，步骤S2中，当k为奇数时、N为偶数时，设计为第四类FIR滤波器，此时：

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{e}^{\frac{j\mathrm{\π}}{2}}{H}_o\left(\mathrm{\ω}\right);$

$H_O\left(\mathrm{\ω}\right)=\underset{n=1}{\overset{N/2}{\Σ}}b\left(n\right)\sin \left(\frac{2n-1}{2}\mathrm{\ω}\right)$

其中：

$b\left(n\right)=2h(\frac{N}{2}-n),n=1,2,...,\frac{N}{2}$

定义：

$B={\[b\left(1\right),b\left(2\right),...,b\left(\frac{N}{2}\right)\]}^T$

$S\left(\mathrm{\ω}\right)={\[sin\left(\frac{\mathrm{\ω}}{2}\right),sin\left(\frac{3}{2}\mathrm{\ω}\right),...,sin\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_o(ω)可表示为H_o(ω)＝B^TS(ω)。

在一种优选的方案中，步骤S3中：

k为偶数时，适应函数为：

$e^2={\∫}_0^{{\mathrm{\ω}}_p}|D_e\left(\mathrm{\ω}\right)-H_e\left(\mathrm{\ω}\right){|}^{2}d\mathrm{\ω};$

k为奇数时，适应函数为：

$e^2={\∫}_0^{{\mathrm{\ω}}_p}|D_o\left(\mathrm{\ω}\right)-H_o\left(\mathrm{\ω}\right){|}^{2}d\mathrm{\ω}.$

在一种优选的方案中，步骤S5中，采用遗传算法获得寻优结果向量B包括以下步骤：

S5.1：执行交叉操作，产生后代种群；

S5.2：执行变异操作，随机产生变异个体；

S5.3：计算个体适应度函数；

S5.4：执行选择操作，对优秀个体选择出来产生下一代；

S5.5：重复步骤S5.1到S5.4，直到不能搜索到更好的后代，此时即可获得寻优结果向量B。

与现有技术相比，本发明技术方案的有益效果是：本发明公开了一种基于遗传算法的高阶数字微分器设计方法，该方法通过不断进行遗传算法的选择、交叉和变异等操作，在有效频段内最小化误差的平方值，从而求解出对应的线性相位FIR微分器系数。该方法简单易行，鲁棒性强，结果令人满意。

附图说明

图1为本发明基于遗传算法的高阶数字微分器设计方法的流程图。

图2为二阶微分器示意图(N＝32，偶阶，偶长度)，其中：(a)冲击响应；(b)幅度响应；(c)最优与平均适应度函数。

图3为三阶微分器示意图(N＝32，奇阶，偶长度)，其中：(a)冲击响应；(b)幅度响应(c)最优与平均适应度函数。

图4为三阶微分器示意图(N＝27，奇阶，奇长度)，其中：(a)冲击响应；(b)幅度响应(c)最优与平均适应度函数

图5为四阶微分器示意图(N＝34，偶阶，奇长度)，其中：(a)冲击响应；(b)幅度响应；(c)最优与平均适应度函数。

具体实施方式

下面结合附图和实施例对本发明的技术方案做进一步的说明。

实施例1

如图1所示，本实施例提供一种基于遗传算法的高阶数字微分器设计方法，包括以下步骤：

1、确定高阶微分器的期望传递函数：

$D\left(\mathrm{\ω}\right)=\left\{\begin{array}{cc}{\left(\frac{j\mathrm{\ω}}{2\mathrm{\π}}\right)}^k,& 0\≤\mathrm{\ω}\≤{\mathrm{\ω}}_p\\ -{\[\frac{j(2\mathrm{\π}-\mathrm{\ω})}{2\mathrm{\π}}\]}^k,& 2\mathrm{\π}-{\mathrm{\ω}}_p\≤\mathrm{\ω}\≤2\mathrm{\π}\end{array}\right.$

其中k代表微分器的阶数。当k为偶数时，为偶数阶微分器，当k为奇数时为奇数阶微分器。ω_p为微分器需要作用的截至频率，而在ω_p之外的频带，滤波器的特性我们并不需要关心。

2、当k为偶数时，即待设计微分器为偶数阶微分器时，可设计为第一类或第二类FIR滤波器。

假设待设计的FIR滤波器H(z)长度为N，则有：

$H\left(z\right)=\underset{n=0}{\overset{N-1}{\Σ}}h\left(n\right)z^{-n}$

其中h(n)＝h(N-1-n)为对应的冲击响应， 表示向下取整

显然，这满足对称性质，当N为奇数时，为第一类FIR滤波器；当N为偶数时，为第二类FIR滤波器。

3、H(z)的频率响应函数就可以表示为

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{H}_e\left(\mathrm{\ω}\right)$

4、当N为奇数时(第一类FIR)

$H_e\left(\mathrm{\ω}\right)=\underset{n=0}{\overset{(N-1)/2}{\Σ}}b\left(n\right)cos\left(n\mathrm{\ω}\right)$

其中

$b\left(n\right)=\left\{\begin{array}{cc}h\left(\frac{N-1}{2}\right)& n=0\\ 2h(\frac{N-1}{2}-n)& n=1,2,...,\frac{N-1}{2}\end{array}\right.$

定义

$B={\[b\left(0\right),b\left(1\right),...,b\left(\frac{N-1}{2}\right)\]}^T$

$C\left(\mathrm{\ω}\right)={\[1,cos\left(\mathrm{\ω}\right),...,cos\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_e(ω)可表示为H_e(ω)＝B^TC(ω)。

5、当N为偶数时(第二类FIR)

$H_e\left(\mathrm{\ω}\right)=\underset{n=1}{\overset{N/2}{\Σ}}b\left(n\right)cos\left(\frac{2n-1}{2}\mathrm{\ω}\right)$

其中

$b\left(n\right)=2h(\frac{N}{2}-n),n=1,2,...,\frac{N}{2}$

定义

$B={\[b\left(1\right),b\left(2\right),...,b\left(\frac{N}{2}\right)\]}^T$

$C\left(\mathrm{\ω}\right)={\[cos\left(\frac{\mathrm{\ω}}{2}\right),cos\left(\frac{3}{2}\mathrm{\ω}\right),...,cos\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_e(ω)可表示为H_e(ω)＝B^TC(ω)

6、所谓数字微分器设计，即在0～ω_p频段通过确定b(n)来逼近D(ω)，也就是最小化

其中D_e(ω)＝D(ω)。

7、为了应用遗传算法寻找符合要求的最优的B向量。以H_e(ω)的值为染色体，以为适应函数，对B进行多变量寻优。

遗传算法作如下设定：

种群数量定义为100；

在交叉操作中，应用Heuristic准则，即child＝parent2+1.2*(parent1-parent2)，其中parent1和parent2为两父代，且适应度parent1优于parent2，child为子代；

在变异操作中，采用随机产生新染色体的方法，变异率设为0.01

在选择操作中，从父代、子代以及变异的染色体当中，选择适应度最优的个体进行下一次遗传的父代产生新的子代，以逐步逼近理想微分器。

8、当适应度函数e²取得最小值时，得出对应的B向量，即可计算出偶次阶微分器的系数h(n)，n＝0,…,N-1，完成偶次阶微分器的设计。

9、当k为奇数时，即待设计微分器为奇数阶微分器时，可设计为第三类或第四类FIR滤波器。

假设待设计的FIR滤波器H(z)长度为N，则有

$H\left(z\right)=\underset{n=0}{\overset{N-1}{\Σ}}h\left(n\right)z^{-n}$

当N为奇数时，为第三类FIR滤波器；当N为偶数时，为第四类FIR滤波器，

对于第三类FIR滤波器，

对于第四类FIR滤波器，

显然，这满足反对称性质，对应的频率响应函数就可以表示为

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{e}^{\frac{j\mathrm{\π}}{2}}{H}_o\left(\mathrm{\ω}\right)$

10、当N为奇数时(第三类FIR)

$H_o\left(\mathrm{\ω}\right)=\underset{n=1}{\overset{(N-1)/2}{\Σ}}b\left(n\right)sin\left(n\mathrm{\ω}\right)$

其中

$b\left(n\right)=2h(\frac{N-1}{2}-n),n=1,2,...,\frac{N-1}{2}$

定义

$B={\[b\left(1\right),b\left(2\right),...,b\left(\frac{N-1}{2}\right)\]}^T$

$S\left(\mathrm{\ω}\right)={\[sin\left(\mathrm{\ω}\right),sin\left(2\mathrm{\ω}\right),...,sin\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_O(ω)可表示为H_O(ω)＝B^TS(ω)。

11、当N为偶数时(第四类FIR)

$H_O\left(\mathrm{\ω}\right)=\underset{n=1}{\overset{N/2}{\Σ}}b\left(n\right)sin\left(\frac{2n-1}{2}\mathrm{\ω}\right)$

其中

$b\left(n\right)=2h(\frac{N}{2}-n),n=1,2,...,\frac{N}{2}$

定义

$B={\[b\left(1\right),b\left(2\right),...,b\left(\frac{N}{2}\right)\]}^T$

$S\left(\mathrm{\ω}\right)={\[sin\left(\frac{\mathrm{\ω}}{2}\right),sin\left(\frac{3}{2}\mathrm{\ω}\right),...,sin\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_o(ω)可表示为H_o(ω)＝B^TS(ω)。

仿照步骤7，应用遗传算法，寻找使最小的B向量，其中D_o(ω)取D(ω)的虚部。当适应度函数e²取得最小值时，得出对应的B向量，即可计算出奇次阶微分器的系数h(n)，n＝0,…,N-1，完成奇次阶微分器的设计。

实施例2

对于图2(偶数阶，偶长度)所实现的二阶微分器。

先确定要实现的微分器阶数k＝2，要设计的FIR滤波器长度N＝32，以及滤波器截至频率ω_p＝0.92π。

待寻优向量B为1×16行向量，即B＝[B(1)，B(2)，…，B(16)]

染色体的值被积函数为f(ω)＝Deω-He(ω)2,对上述被积函数在0～ωp频段内进行Gauss-Kronrod数值积分。所得积分的值作为遗传算法的适应度函数fitnessfcn(B)。

对上述条件代入遗传算法中，并设定估计值的上下界，本设计的下界和上界分别设置为-2和2，应用Heuristic交叉操作准则，对B向量进行寻优。寻优结果如图2(c)所示。得到B向量的值为

B＝[-0.0600380314773367 0.0949230407895234 -0.0585146692182696

0.0409091719995725 -0.0304305776845123 0.0230406656188536

-0.0186201074932732 0.0143041725561841 -0.0108589682873630

0.00882499852915601 -0.00692292034036513 0.00524604238067849

-0.00397217080385697 0.00293827042906896 -0.00170079087374337

0.00110692084969835]

最终该B向量对应的适应度函数的值(误差的平方值)为:3.0382647694387307×10^-6，利用关系h(16-n)＝0.5×B(n),n＝1,…,16，可确定微分器系数h(0)～h(15)，利用对称性质h(32-1-n)＝h(n)，n＝0,…,15，可确定微分器系数h(16)～h(31)。由此求得滤波器系数向量

h＝[0.000553460424849174 -0.000850395436871684 0.00146913521453448

-0.00198608540192849 0.00262302119033925 -0.00346146017018256

0.00441249926457801 -0.00542948414368151 0.00715208627809206

-0.00931005374663661 0.0115203328094268 -0.0152152888422562

0.0204545859997863 -0.0292573346091348 0.0474615203947617

-0.0300190157386683 -0.0300190157386683 0.0474615203947617

-0.0292573346091348 0.0204545859997863 -0.0152152888422562

0.0115203328094268 -0.00931005374663661 0.00715208627809206

-0.00542948414368151 0.00441249926457801 -0.00346146017018256

0.00262302119033925 -0.00198608540192849 0.00146913521453448

-0.000850395436871684 0.000553460424849174]

该微分器的冲击响应和频率响应如图2(a)、(b)所示。

对于图3(奇数阶，偶长度)所实现的微分器阶数k＝3，要设计的FIR滤波器长度N＝32，以及滤波器截至频率ω_p＝π。

待寻优向量B为1×16行向量，即B＝[B(1)，B(2)，…，B(16)]

染色体的值被积函数为f(ω)＝Doω-Ho(ω)2，对上述被积函数在0～ωp频段内进行Gauss-Kronrod数值积分。所得积分的值作为遗传算法的适应度函数，寻优得到的向量B，寻优过程如图3(c)所示。得到B向量的值为

B＝[-0.05723829329724478 0.030906005212783753 -0.012267970604854156

0.005992112473784427 -0.0037494751967100995 0.0026949416412901366

-0.0016117315686561785 0.001156580412998436 -6.018754603643743×10^-4

9.86342488373556×10^-4 -6.870418563593486×10^-4 3.290447088439463×10^-4

-5.39421333833814×10^-4 1.6218959029910924×10^-4 -4.43980181941184×10^-4

5.846863404313215×10^-4]

最终该B向量对应的适应度函数的值(误差的平方值)为:2.2264872295677593×10^-6，利用关系h(16-n)＝0.5×B(n),n＝1,…,16，可确定微分器系数h(0)～h(15)，利用对称性质h(32-1-n)＝-h(n)，n＝0,…,15，可确定微分器系数h(16)～h(31)。由此求得滤波器系数向量

h＝[0 -0.000221990090970592 8.10947951495546×10^-5

-0.000269710666916907 0.000164522354421973 -0.000343520928179674

0.000493171244186778 -0.000300937730182187 0.000578290206499218

-0.000805865784328089 0.00134747082064507 -0.00187473759835505

0.00299605623689221 -0.00613398530242708 0.0154530026063919

-0.0286191466486224 0.0286191466486224 -0.0154530026063919

0.00613398530242708 -0.00299605623689221 0.00187473759835505

-0.00134747082064507 0.000805865784328089 -0.000578290206499218

0.000300937730182187 -0.000493171244186778 0.000343520928179674

-0.000164522354421973 0.000269710666916907 -8.10947951495546×10^-5

0.000221990090970592 0]

该微分器的冲击响应和频率响应如图3(a)、(b)所示。

对于图4(奇数阶，奇长度)所实现的微分器阶数k＝3，要设计的FIR滤波器长度N＝27，以及滤波器截至频率ω_p＝0.88π。

待寻优向量B为1×13行向量，即B＝[B(1)，B(2)，…，B(13)]。染色体的值被积函数为f(ω)＝|D_o(ω)-H_o(ω)|²,对上述被积函数在0～ω_p频段内进行Gauss-Kronrod数值积分。所得积分的值作为遗传算法的适应度函数，寻优得到的向量B，寻优过程如图4(c)所示。得到B向量的值为

B＝[-0.024083536674745287 0.018910141128054402

-0.004742738495578499 -0.00545961436928196 0.012812098527400258

-0.01788529220610157 0.020185936513756147 -0.021450090753923512

0.01966465417017256 -0.01747981740243351 0.01528034502076689

-0.01218138939772307 0.008602450370204556]

最终该B向量对应的适应度函数的值(误差的平方值)为:1.3228584269256425×10^-4，利用关系h(13)＝0，h(13-n)＝0.5×B(n),n＝1,…,13，可确定微分器系数h(0)～h(12)，利用对称性质h(27-1-n)＝-h(n)，n＝0,…,13，可确定微分器系数h(14)～h(26)。由此求得滤波器系数向量

h＝[0.00430122518510228 -0.00609069469886154 0.00764017251038345

-0.00873990870121675 0.00983232708508628 -0.0107250453769618

0.0100929682568781 -0.00894264610305079 0.00640604926370013

-0.00272980718464098 -0.00237136924778925 0.00945507056402720

-0.0120417683373726 0 0.0120417683373726

-0.00945507056402720 0.00237136924778925 0.00272980718464098

-0.00640604926370013 0.00894264610305079-0.0100929682568781

0.0107250453769618 -0.00983232708508628 0.00873990870121675

-0.00764017251038345 0.00609069469886154 -0.00430122518510228]

该微分器的冲击响应和频率响应如图4(a)、(b)所示。

对于图5(偶数阶，奇长度)所实现的微分器阶数k＝4，要设计的FIR滤波器长度N＝25，以及滤波器截至频率ω_p＝π。

待寻优向量B为1×13行向量，即B＝[B(1)，B(2)，…，B(13)]。染色体的值被积函数为f(ω)＝|D_e(ω)-H_e(ω)|²,对上述被积函数在0～ω_p频段内进行Gauss-Kronrod数值积分。所得积分的值作为遗传算法的适应度函数，寻优得到的向量B，寻优过程如图5(c)所示。得到B向量的值为

B＝[0.012488175693597886 -0.019796484314248068

0.010498435892450246 -0.005565438371338432 0.00326215711081659

-0.0020947541809127617 0.0013384845435775621 -0.00107433724452206

7.047398283569596×10^-4 -6.701172202185925×10^-4 5.473029250253703×10^-4

-5.094238501774658×10^-4 1.035418511600765×10^-4]

最终该B向量对应的适应度函数的值(误差的平方值)为:1.1680956317584083×10^-4，利用关系h(12)＝B(1)，h(12-n)＝0.5×B(n+1),n＝1,…,12，可确定微分器系数h(0)～h(11)，利用对称性质h(25-1-n)＝-h(n)，n＝0,…,11，可确定微分器系数h(13)～h(24)。由此求得滤波器系数向量

h＝[5.17709255800383×10^-5 -0.000254711925088733 0.000273651462512685

-0.000335058610109296 0.000352369914178480 -0.000537168622261030

0.000669242271788781 -0.00104737709045638 0.00163107855540830

-0.00278271918566922 0.00524921794622512 -0.00989824215712403

0.0124881756935979 -0.00989824215712403 0.00524921794622512

-0.00278271918566922 0.00163107855540830 -0.00104737709045638

0.000669242271788781 -0.000537168622261030 0.000352369914178480

-0.000335058610109296 0.000273651462512685 -0.000254711925088733

5.17709255800383×10^-5]

该微分器的冲击响应和频率响应如图5(a)、(b)所示。

总的来说，本发明为一种基于遗传算法的高阶数字微分器设计方法。从以上的实验结果可以看出本发明对高阶数字微分器的实现，具有简单有效、鲁棒性强的优点。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (7)

1.一种基于遗传算法的高阶数字微分器设计方法，其特征在于，包括以下步骤：

S1：确定高阶数字微分器的阶数k、长度N以及截止频率ω_p，高阶数字微分器的期望传递函数为：

$D\left(\mathrm{\ω}\right)=\left\{\begin{array}{cc}{\left(\frac{j\mathrm{\ω}}{2\mathrm{\π}}\right)}^k,& 0\≤\mathrm{\ω}\≤{\mathrm{\ω}}_p\\ -{\[\frac{j(2\mathrm{\π}-\mathrm{\ω})}{2\mathrm{\π}}\]}^k,& 2\mathrm{\π}-{\mathrm{\ω}}_p\≤\mathrm{\ω}\≤2\mathrm{\π}\end{array}\right.$

其中，k代表高阶数字微分器的阶数，当k为偶数时，为偶数阶数字微分器，当k为奇数时为奇数阶数字微分器；ω_p为高阶数字微分器需要作用的截至频率；

高阶数字微分器的系统函数表示为：

$H\left(z\right)=\underset{n=0}{\overset{N-1}{\Σ}}h\left(n\right)z^{-n}$

其中h(n)为对应的冲击响应；

S2：根据k、N的奇偶性选择合适类型的FIR滤波器模型：

S3：建立适应度函数表达式；

S4：遗传算法初始化，建立初始种群；

S5：采用遗传算法获得寻优结果向量B，寻优结果向量B为使适应度函数取最小的B向量；

S6：根据向量B计算为微分器系数。

2.根据权利要求1所述的基于遗传算法的高阶数字微分器设计方法，其特征在于，步骤S2中，当k为偶数时、N为奇数时，设计为第一类FIR滤波器，h(n)＝h(N-1-n)， 表示向下取整；此时：

H(z)的频率响应函数表示为：

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{H}_e\left(\mathrm{\ω}\right);$

$H_e\left(\mathrm{\ω}\right)=\underset{n=0}{\overset{(N-1)/2}{\Σ}}b\left(n\right)cos\left(n\mathrm{\ω}\right)$

其中：

$b\left(n\right)=\left\{\begin{array}{cc}h\left(\frac{N-1}{2}\right)& n=0\\ 2h(\frac{N-1}{2}-n)& n=1,2,\mathrm{...},\frac{N-1}{2}\end{array}\right.$

定义：

$B={\[b\left(0\right),b\left(1\right),\mathrm{...},b\left(\frac{N-1}{2}\right)\]}^T$

$C\left(\mathrm{\ω}\right)={\[1,cos\left(\mathrm{\ω}\right),\mathrm{...},cos\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_e(ω)可表示为H_e(ω)＝B^TC(ω)。

3.根据权利要求1所述的基于遗传算法的高阶数字微分器设计方法，其特征在于，步骤S2中，当k为偶数时、N为偶数时，设计为第二类FIR滤波器，h(n)＝h(N-1-n)， 表示向下取整；此时：

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{H}_e\left(\mathrm{\ω}\right);$

$H_e\left(\mathrm{\ω}\right)=\underset{n=1}{\overset{N/2}{\Σ}}b\left(n\right)cos\left(\frac{2n-1}{2}\mathrm{\ω}\right)$

其中：

$b\left(n\right)=2h(\frac{N}{2}-n),n=1,2,\mathrm{...},\frac{N}{2}$

定义：

$B={\[b\left(1\right),b\left(2\right),\mathrm{...},b\left(\frac{N}{2}\right)\]}^T$

$C\left(\mathrm{\ω}\right)={\[cos\left(\frac{\mathrm{\ω}}{2}\right),cos\left(\frac{3}{2}\mathrm{\ω}\right),\mathrm{...},cos\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_e(ω)可表示为H_e(ω)＝B^TC(ω)。

4.根据权利要求1所述的基于遗传算法的高阶数字微分器设计方法，其特征在于，步骤S2中，当k为奇数时、N为奇数时，设计为第三类FIR滤波器，h(n)＝-h(N-1-n)，此时：

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{e}^{\frac{j\mathrm{\π}}{2}}{H}_o\left(\mathrm{\ω}\right);$

$H_o\left(\mathrm{\ω}\right)=\underset{n=1}{\overset{(N-1)/2}{\Σ}}b\left(n\right)sin\left(n\mathrm{\ω}\right)$

其中：

$b\left(n\right)=2h(\frac{N-1}{2}-n),n=1,2,\mathrm{...},\frac{N-1}{2}$

定义：

$B={\[b\left(1\right),b\left(2\right),\mathrm{...},b\left(\frac{N-1}{2}\right)\]}^T$

$S\left(\mathrm{\ω}\right)={\[sin\left(\mathrm{\ω}\right),sin\left(2\mathrm{\ω}\right),\mathrm{...},sin\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_O(ω)可表示为H_O(ω)＝B^TS(ω)。

5.根据权利要求1所述的基于遗传算法的高阶数字微分器设计方法，其特征在于，步骤S2中，当k为奇数时、N为偶数时，设计为第四类FIR滤波器，h(n)＝-h(N-1-n)，此时：

$H\left(e^{-j\mathrm{\ω}}\right)=e^{-\frac{j\mathrm{\ω}(N-1)}{2}}{e}^{\frac{j\mathrm{\π}}{2}}{H}_o\left(\mathrm{\ω}\right);$

$H_O\left(\mathrm{\ω}\right)=\underset{n=1}{\overset{N/2}{\Σ}}b\left(n\right)sin\left(\frac{2n-1}{2}\mathrm{\ω}\right)$

其中：

$b\left(n\right)=2h(\frac{N}{2}-n),n=1,2,\mathrm{...},\frac{N}{2}$

定义：

$B={\[b\left(1\right),b\left(2\right),\mathrm{...},b\left(\frac{N}{2}\right)\]}^T$

$S\left(\mathrm{\ω}\right)={\[sin\left(\frac{\mathrm{\ω}}{2}\right),sin\left(\frac{3}{2}\mathrm{\ω}\right),\mathrm{...},sin\left(\frac{N-1}{2}\mathrm{\ω}\right)\]}^T$

则H_o(ω)可表示为H_o(ω)＝B^TS(ω)。

6.根据权利要求1所述的基于遗传算法的高阶数字微分器设计方法，其特征在于，步骤S3中：

k为偶数时，适应函数为：

$e^2={\∫}_0^{{\mathrm{\ω}}_p}|D_e\left(\mathrm{\ω}\right)-H_e\left(\mathrm{\ω}\right){|}^{2}d\mathrm{\ω};$

k为奇数时，适应函数为：

$e^2={\∫}_0^{{\mathrm{\ω}}_p}|D_o\left(\mathrm{\ω}\right)-H_o\left(\mathrm{\ω}\right){|}^{2}d\mathrm{\ω}.$

7.根据权利要求1所述的基于遗传算法的高阶数字微分器设计方法，其特征在于，步骤S5中，采用遗传算法获得寻优结果向量B包括以下步骤：

S5.1：执行交叉操作，产生后代种群；

S5.2：执行变异操作，随机产生变异个体；

S5.3：计算个体适应度函数；

S5.4：执行选择操作，对优秀个体选择出来产生下一代；

S5.5：重复步骤S5.1到S5.4，直到不能搜索到更好的后代，此时即可获得寻优结果向量B。