[发明专利]一种基于混合粒子群算法的数字全通滤波器设计方法

权利要求书

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

(1)、构造基于二阶节级联的数字全通滤波器的相频响应；

(1.1)、根据实际需求，设置数字全通滤波器的阶数N，N为偶数，对应的数字全通滤波器的二阶节个数为N/2，那么级联的数字全通滤波器的Z域表达式H(z)为：

其中，H_n(z)表示第n阶数字全通滤波器的Z域表达式，a_n1、a_n2是第n阶数字全通滤波器的系数；

(1.2)、设第n个二阶节的极点为M_n、θ_n为第n个二阶节极点的模值与相角，则N阶数字全通滤波器的极点表示为

定义变量U＝[M₁,θ₁,M₂,θ₂,…,M_n,θ_n,…,M_N/2,θ_N/2]，则数字全通滤波器的相频响应表示为：

其中，ω＝[ω₁,ω₂,…,ω_l,…,ω_L]，L为采样的频点个数；

(2)、构造数字全通滤波器的代价函数；

设数字全通滤波器的目标相频响应为φ_goal(ω)，则由数字全通滤波器的相频响应与目标相频响应为φ_goal(ω)的误差值构造而成的代价函数φ_error(U)表示为：

(3)、利用粒子群算法寻找代价函数φ_error(U)的全局最优解；

(3.1)、设定粒子群算法的种群数S以及粒子群算法的最大迭代次数K；

(3.2)、初始化粒子群为了保证数字全通滤波器每个二阶节的稳定性，则种群中每个粒子初始化的过程需要满足ρ为略小于1的常数；

(3.3)、粒子群更新及迭代；

(3.3.1)、计算每个粒子的更新速度；

其中，为第s个粒子第k次更新后的速度，C₁和C₂分别对应加速因子，r₁和r₂为[0,1]内的随机数，zb_s和gb分别为第s个粒子的历史最优解以及S个粒子组成的粒子群的全局最优解；

(3.3.2)、更新粒子位置；

其中，χ(·)为限定函数；

(3.3.3)、计算迭代后的粒子的代价函数并根据下式更新第s个粒子的历史最优解zb_s以及S个粒子组成的粒子群的全局最优解gb；

(3.4)、重复步骤(3.3)，共计迭代K次，得到最终的全局最优解gb；

(4)、利用列文伯格马夸尔特算法进行局部最优求解；

(4.1)、设置列文伯格马夸尔特算法的总迭代次数K_max，k＝1,2,…,K_max，初始化k＝1；

(4.2)、计算列文伯格马夸尔特算法的迭代起始点X₀＝[x₁,θ₁,x₂,θ₂,…,x_n,θ_n,…,x_N/2,θ_N/2]，其中，

(4.3)、计算第k次迭代后X_k对应的数字全通滤波器的相频响应φ_design(ω,X)

(4.4)、计算第k次迭代后X_k对应的残差向量R_k；

R_k＝φ_design(ω,X_k)-φ_goal(ω) (9)

(4.5)、计算第k次迭代后X_k对应的代价函数值；

φ_error(X_k)＝R_k·R_k^T (10)

(4.6)、计算第k次迭代后X_k对应的雅阁比矩阵J_k，黑塞矩阵H_k＝J_k^TJ_k以及对角阵D_k＝diag{H_k}；

(4.7)、计算第k次迭代后的更新向量；

Δ_k＝(H_k+λ·D_k)^-1·J_k·R_k (11)

(4.8)、计算第k次迭代后的新向量X_new＝X_k+Δ_k，以及对应的代价函数值φ_error(X_new)，然后比较φ_error(X_new)与φ_error(X_k)的大小，如果φ_error(X_new)＜φ_error(X_k)，则令X_k+1＝X_new，且λ＝λ/v；如果φ_error(X_new)≥φ_error(X_k)，则令X_k+1＝X_k，且λ＝λ×v；其中，λ、v均为整数，λ＞＞v；

(4.9)、循环执行步骤(4.3)～(4.8)共计K_max次，最终得到迭代终点X_final；

(5)、逆映射；

将X_final带入到下式中，得到最终的极点向量U_final；

(6)、根据极点向量U_final计算数字全通滤波器系数；

完成全通滤波器的设计。

2.根据权利要求1所述的一种基于混合粒子群算法的数字全通滤波器设计方法，其特征在于，所述限定函数χ(·)的表达式为：