本发明提供一种可分离二维FIR滤波器的系数量化方法。本发明不同于普通的量化方法，而是以最大纹波误差最小化(l∞范数)为目标函数，采用分步式整数规划算法求解。首先利用一部分已知的连续函数系数，采用整数规划算法求解其余部分量化整数系数，然后利用求解出的这一部分量化系数再次采用整数规划算法求解原来部分的量化系数，最终获得全部量化系数。本发明的优点是有效减少了可分离二维FIR滤波器的系数量化误差，降低了有限字长效应，提高滤波器性能指标，获得具有最优有限字长属性的可分离二维FIR滤波器。
1.一种可分离二维FIR滤波器的系数量化方法，其特征在于：该方法具体包括以下步骤：步骤一、根据设计要求，确定滤波器的基本参数，包括二维FIR滤波器的类型、阶数N2、频率点数Γ，通带截止频率ωp，阻带截止频率ωs，设计出可分离的二维FIR滤波器原型，称之为原型滤波器；
设计的可分离二维FIR滤波器的脉冲响应系数矩阵描述如式(1)；
其中列向量ri和si是一维FIR子滤波器的脉冲响应，ri尺寸大小为2·N1+1，si尺寸大小为2·N2+1，其中N1＝N2＝(N-1)/2，K是并行结构的数量；为保证滤波器是线性相位，这些列向量满足中心对称或者反中心对称，所以只需要考虑ri和si列向量中的一半元素，分别记为ri_half＝[ri(N1+1) ri(N1+2)…ri(2·N1+1)]T (2)T
si_half＝[si(N2+1) si(N2+1)…si(2·N2+1)]，i＝1,2,...,K (3)所设计具有线性相位的可分离二维FIR滤波器的每个频率点的响应见式(4)；
式(4)中ωl是感兴趣的频率点，l∈{1,2,...,Γ}， 所要设计的可分Τ Τ Τ Τ
离二维FIR滤波器系数向量记为x，定义如下：x＝[r1_half r2_half ...rK_half s1_half s2_halfΤ...sK_halfΤ]Τ；
步骤二、使用信赖域迭代梯度搜索技术优化原型滤波器中的系数向量x；
根据式(4)可分离二维FIR滤波器的整体频率响应为H(ω|x)＝[H(ω1|x) H(ω2|x)...H(ωΓ|x)]Τ，记理想滤波器的频率响应为Hd(ω)＝[Hd(ω1) Hd(ω2)...Hd(ωΓ)]Τ；记w＝[w(ω1) w(ω2)...w(ωΓ)]T为权重向量；采用最大误差最小化设计准则，使用TR-IGS优化系数向量x见式(5)；
Subject to:||w·(H(ω|x)-Hd(ω))||∞≤δ (5.b)式(5)中，||·||∞无穷范数运算，xv指的是x中的非0和非1的元素，“·”指两向量中逐元素相乘，δ指通带、阻带中最大峰值纹波误差；
式(6)中，M为x的长度，若x(e)＝0，e＝1,2,...,M，则令
2-3.求解H(ω|x)的泰勒一阶近似，见式(7)；
式(7)中，xI是x的初始值，▽x是x＝xI时沿G(xI)下降的步长，若x(e)＝0，则令▽x＝0，ε是一个正实数，可根据实际系数情况进行取值；将式(7)代入式(5)中，将问题(5)转化成凸规划问题(8)；
Subject to:||w·(H(ω|xI)+G(xI)·▽x-Hd(ω))||∞≤δ (8.b)||▽x||∞≤ε (8.c)
记j为迭代次数，初始值j＝0，xI(j)＝x，根据式(6)求出此时的G(xI(j))，代入式(8)并进行求解，求出的解记为▽x(j)；然后求出新的x(j+1)＝xI(j)+▽x(j)，并将此时的x(j+1)作为第j+1迭代的初始值xI(j+1)重复上述求解过程，直到x不能进一步优化或不再满足式(8)的约束条件或达到最大迭代次数；
设步骤二得出的x经过量化后的整数系数向量为 量化后的实际频率响应记为设要求的滤波器系数的最大有限字长
为Q比特，其中字长包含符号位，量化目标准则采用最大误差最小化设计准则，则系数量化问题描述为其中δQ表示最大有限字长为Q比特时滤波器的最大峰值纹波误差，也称为Q比特下的量化误差；记 和 分别是ri和si量化后的系数，则记
式(10)中ρ＝1/||x||∞为归一化因子，Ω＝2(Q-1)-1；记则量化后的实际频率响应 又可以表述为式(13)；
将式(13)代入式(9)，可将系数量化问题具体化简为式(14)；
式(14)是一个整数规划问题，由于变量 包含着级联的一维FIR子滤波器的系数变量和 在计算频率响应时，两者是一种相乘的关系，是非线性的，所以无法采用相关的整数规划算法直接求解，提出分步式整数规划算法量化方案；
4.1利用步骤二得出的连续系数ri_half，求解 中的量化系数，目标函数和相应的约束条件为式(15)中，δ′Q是当前量化状态下所获得的最小峰值纹波误差，式(15.b)中除以(Ω·ρ)表示将量化的整数值转化为小数值用于计算实际的频率响应；
式(15.c)表示 的动态取值范围，式(15.d)中，Ψs表示向量xs中0元素的集合，量化过程中，0元素不参与量化；采用Gurobi求解器与yalmip工具箱嵌套使用的方式求解式(15)中的
4.2利用4.1求解出的 再求解 中的量化系数，目标函数和相应的约束条件为式(18)中，
式(18.c)表示 的动态取值范围，式(18.d)中Ψr表示ri_half系数中的0元素的集合；采用Gurobi求解器与yalmip工具箱嵌套使用的方式求解式(15)中的最终获得量化后的整数系数向量
一种可分离二维FIR滤波器的系数量化方法
技术领域
[0001] 本发明属于数字信号处理技术领域，具体涉及一种具有最优有限字长属性的可分离二维FIR滤波器的量化方法。
背景技术
[0002] 二维FIR数字滤波器作为二维数字信号处理的重要组成部分，已经被广泛应用于医学图像处理、卫星图像处理、雷达、声纳以及地震信号处理等很多方面，而拥有较低复杂度的可分离二维FIR数字滤波器也有着广泛的应用。
[0003] 在滤波器设计与硬件实现过程中，滤波器存在系数量化误差和运算误差。实际的数字滤波器是由有限精度的数字和算法实现，因此对于滤波器的系数以及信号的输入输出都需要通过有限数量的比特来量化，并且随着比特数的增加，硬件系统的存储和计算负担会急剧增加，所以需要尽可能减少比特位数，这就造成了滤波器有限字长效应问题。有限位的滤波器系数会影响滤波器零极点的位置，造成滤波的不稳定性，导致实际滤波器的频率响应与所设计的理想频率响应有着细微的差别，这就是所谓的系数量化误差，而这种误差在某种程度上会影响滤波器的滤波性能。
[0004] 在可分离二维FIR滤波器结构中，由多个并联结构的级联一维FIR滤波器构成，这种非线性的滤波器结构更增加了系数量化误差的不确定性。一些适用于一维FIR滤波器的系数量化方法，如果直接使用在这种滤波器结构上，很可能造成不理想的有限字长效应。因此，需要对设计出的可分离二维FIR稀疏滤波器进行系数量化研究，寻找出适用于此滤波器结构的量化方法，得到具有最优有限字长效应的可分离二维FIR滤波器。
发明内容
[0005] 本发明针对现有技术的不足。提供一种可分离二维FIR滤波器的系数量化方法。
[0006] 本发明一种可分离二维FIR滤波器的系数量化方法，该方法具体如下：
[0007] 步骤一、根据设计要求，确定滤波器的基本参数，包括二维FIR滤波器的类型、阶数N2、频率点数Γ，通带截止频率ωp，阻带截止频率ωs，设计出可分离的二维FIR滤波器原型，称之为原型滤波器。
[0008] 设计的可分离二维FIR滤波器的脉冲响应系数矩阵描述如式(1)。
[0009]
[0010] 其中列向量ri和si是一维FIR子滤波器的脉冲响应，ri尺寸大小为2·N1+1，si尺寸大小为2·N2+1，其中N1＝N2＝(N-1)/2，K是并行结构的数量。为保证滤波器是线性相位，这些列向量满足中心对称或者反中心对称，所以只需要考虑ri和si列向量中的一半元素，分别记为
[0011] ri_half＝[ri(N1+1) ri(N1+2) … ri(2·N1+1)]T (2)
[0012] si_half＝[si(N2+1) si(N2+1) … si(2·N2+1)]T，i＝1,2,...,K (3)[0013] 所设计具有线性相位的可分离二维FIR滤波器的每个频率点的响应见式(4)。
[0014]
[0015] 式(4)中ωl是感兴趣的频率点，l∈{1,2,...,Γ}， 所要设计的可分离二维FIR滤波器系数向量记为x，定义如下：x＝[r1_halfT r2_halfT ... rK_halfT s1_halfT T T T
s2_half ... sK_half]。
[0016] 步骤二、使用信赖域迭代梯度搜索(TR-IGS)技术优化原型滤波器中的系数向量x。
[0017] 2-1.确定求解目标函数及约束条件。
[0018] 根据式(4)可分离二维FIR滤波器的整体频率响应为H(ω|x)＝[H(ω1|x) H(ω2|Tx) ... H(ωΓ|x)] ，记理想滤波器的频率响应为Hd(ω)＝[Hd(ω1) Hd(ω2) ... Hd(ωΓ)]T。记w＝[w(ω1) w(ω2) ... w(ωΓ)]T为权重向量。采用最大误差最小化(也称 )设计准则，使用TR-IGS优化系数向量x见式(5)。
[0019]
[0020] Subject to:||w·(H(ω|x)-Hd(ω))||∞≤δ (5.b)
[0021] 式(5)中，||·||∞无穷范数运算，xv指的是x中的非0和非1的元素，“·”指两向量中逐元素相乘，δ指通带、阻带中最大峰值纹波误差。
[0022] 2-2.求解H(ω|x)的雅克比矩阵(梯度矩阵)G(x)。
[0023]
[0024] 式(6)中，M为x的长度，若x(e)＝0，e＝1,2,...,M，则令
[0025] 2-3.求解H(ω|x)的泰勒一阶近似，见式(7)。
[0026]
[0027] 式(7)中，xI是x的初始值， 是x＝xI时沿G(xI)下降的步长，若x(e)＝0，则令ε是一个正实数，可根据实际系数情况进行取值。将式(7)代入式(5)中，将问题(5)转化成凸规划问题(8)。
[0028]
[0029]
[0030]
[0031] 2-4.求解凸规划问题(8)。
[0032] 记j为迭代次数，初始值j＝0，xI(j)＝x，根据式(6)求出此时的G(xI(j))，代入式(8)并进行求解，求出的解记为 然后求出新的 并将此时的x(j+1)作为第j+1迭代的初始值xI(j+1)重复上述求解过程，直到x不能进一步优化或不再满足式(8)的约束条件或达到最大迭代次数。
[0033] 步骤三、确定系数量化的表达式。
[0034] 设步骤二得出的x经过量化后的整数系数向量为 量化后的实际频率响应记为设要求的滤波器系数的最大有限字长为Q比特，其中字长包含符号位，量化目标准则采用最大误差最小化设计准则，则系数量化问题描述为
[0035]
[0036]
[0037]
[0038] 其中δQ表示最大有限字长为Q比特时滤波器的最大峰值纹波误差，也称为Q比特下的量化误差。记 和 分别是ri和si量化后的系数，则记
则
量化后际每个频率点的实际响应见式(10)。
[0039]
[0040] 式(10)中ρ＝1/||x||∞为归一化因子，Ω＝2(Q-1)-1。记
[0041]
[0042]
[0043] 则量化后的实际频率响应 又可以表述为式(13)。
[0044]
[0045] 将式(13)代入式(9)，可将系数量化问题具体化简为式(14)。
[0046]
[0047]
[0048]
[0049] 步骤四、采用分步式整数规划算法求解式(14)。
[0050] 式(14)是一个整数规划问题，由于变量 包含着级联的一维FIR子滤波器的系数变量 和 在计算频率响应时，两者是一种相乘的关系，是非线性的，所以无法采用相关的整数规划算法直接求解，提出分步式整数规划算法量化方案。
[0051] 4.1利用步骤二得出的连续系数ri_half，求解 中的量化系数，目标函数和相应的约束条件为
[0052]
[0053]
[0054]
[0055]
[0056] 式(15)中，δ′Q是当前量化状态下所获得的最小峰值纹波误差，式(15.b)中除以(Ω·ρ)表示将量化的整数值转化为小数值用于计算实际的频率响应。
[0057]
[0058]
[0059] 式(15.c)表示 的动态取值范围，式(15.d)中，Ψs表示向量xs中0元素的集合，量化过程中，0元素不参与量化。采用Gurobi求解器与yalmip工具箱嵌套使用的方式求解式(15)中的
[0060] 4.2利用4.1求解出的 再求解 中的量化系数，目标函数和相应的约束条件为
[0061]
[0062]
[0063]
[0064]
[0065] 式(18)中，
[0066]
[0067]
[0068] 式(18.c)表示 的动态取值范围，式(18.d)中Ψr表示ri_half系数中的0元素的集合。采用Gurobi求解器与yalmip工具箱嵌套使用的方式求解式(15)中的
[0069] 最终获得量化后的整数系数向量
[0070] 本发明具有的有益效果是：
[0071] 本发明可以在硬件实现过程中，不改变可分离二维FIR滤波器结构复杂度，利用分步式整数规划算法降低滤波器系数的量化误差，减少滤波器有限字长效应，提高滤波器性能指标，得到具有最优有限字长属性的可分离二维FIR滤波器。
附图说明
[0072] 图1为本发明的总体设计流程框图。
[0073] 图2为可分离二维FIR滤波器的结构图。
[0074] 图3为设计出的可分离二维FIR圆型滤波器的频率响应图。
[0075] 图4为设计出的可分离二维FIR钻石型滤波器的频率响应图。
[0076] 图5为可分离二维FIR圆型滤波器Q＝5～20时两种量化方法的量化误差曲线图。
[0077] 图6为可分离二维FIR钻石型滤波器Q＝5～20时两种量化方法的量化误差曲线图。
具体实施方式
[0078] 以下结合附图对本发明作进一步说明。
[0079] 如图1所示，具有最优有限字长属性的可分离二维FIR滤波器的设计方法的具体步骤如下：
[0080] 步骤一、根据设计要求，确定滤波器的基本参数，包括二维FIR滤波器的类型、阶数2
N、频率点数Γ，通带截止频率ωp，阻带截止频率ωs，设计出可分离的二维FIR滤波器原型，称之为原型滤波器，可分离二维FIR滤波器的结构如图2所示；。
[0081] 设计的可分离二维FIR滤波器的脉冲响应系数矩阵描述如式(1)。
[0082]
[0083] 其中列向量ri和si是一维FIR子滤波器的脉冲响应，ri尺寸大小为2·N1+1，si尺寸大小为2·N2+1，其中N1＝N2＝N-1/2，K是并行结构的数量。为保证滤波器是线性相位，这些列向量满足中心对称或者反中心对称，所以只需要考虑ri和si列向量中的一半元素，分别记为
[0084] ri_half＝[ri(N1+1) ri(N1+2) … ri(2·N1+1)]T (2)
[0085] si_half＝[si(N2+1) si(N2+1) … si(2·N2+1)]T，i＝1,2,...,K (3)[0086] 所设计具有线性相位的可分离二维FIR滤波器的每个频率点的响应见式(4)。
[0087]
[0088] 式(4)中ωl是感兴趣的频率点，l∈{1,2,...,Γ}， 所要设计的可分离二维FIR滤波器系数向量记为x，定义如下：x＝[r1_halfT r2_halfT ... rK_halfT s1_halfT s2_halfT ... sK_halfT]T。
[0089] 步骤二、使用信赖域迭代梯度搜索(TR-IGS)技术优化原型滤波器中的系数向量x。
[0090] 2-1.确定求解目标函数及约束条件。
[0091] 根据式(4)可分离二维FIR滤波器的整体频率响应为H(ω|x)＝[H(ω1|x) H(ω2|x) ... H(ωΓ|x)]T，记理想滤波器的频率响应为Hd(ω)＝[Hd(ω1) Hd(ω2) ... Hd(ωΓ)]T。记w＝[w(ω1) w(ω2) ... w(ωΓ)]T为权重向量。采用最大误差最小化(也称 )设计准则，使用TR-IGS优化系数向量x见式(5)。
[0092]
[0093] Subject to:||w·(H(ω|x)-Hd(ω))||∞≤δ (5.b)
[0094] 式(5)中，||·||∞无穷范数运算，xv指的是x中的非0和非1的元素，“·”指两向量中逐元素相乘，δ指通带、阻带中最大峰值纹波误差。
[0095] 2-2.求解H(ω|x)的雅克比矩阵(梯度矩阵)G(x)。
[0096]
[0097] 式(6)中，M为x的长度，若x(e)＝0，e＝1,2,...,M，则令
[0098] 2-3.求解H(ω|x)的泰勒一阶近似，见式(7)。
[0099]
[0100] 式(7)中，xI是x的初始值， 是x＝xI时沿G(xI)下降的步长，若x(e)＝0，则令ε是一个正实数，可根据实际系数情况进行取值。将式(7)代入式(5)中，将问题(5)转化成凸规划问题(8)。
[0101]
[0102]
[0103]
[0104] 2-4.求解凸规划问题(8)。
[0105] 记j为迭代次数，初始值j＝0，xI(j)＝x，根据式(6)求出此时的G(xI(j))，代入式(8)并进行求解，求出的解记为 然后求出新的 并将此时的x(j+1)作为第j+1迭代的初始值xI(j+1)重复上述求解过程，直到x不能进一步优化或不再满足式(8)的约束条件或达到最大迭代次数。
[0106] 步骤三、确定系数量化的表达式。
[0107] 设步骤二得出的x经过量化后的整数系数向量为 量化后的实际频率响应记为设要求的滤波器系数的最大有限字长(包含符号位)为Q比特，量化目标准则采用最大误差最小化设计准则，则系数量化问题可以描述为
[0108]
[0109]
[0110]
[0111] 其中δQ表示最大有限字长为Q比特时滤波器的最大峰值纹波误差，也称为Q比特下的量化误差。记 和 分别是ri和si量化后的系数，则记
则 量化后际每个频率点的实际响应见式
(10)。
[0112]
[0113] 式(10)中ρ＝1/||x||∞为归一化因子，Ω＝2(Q-1)-1。记
[0114]
[0115]
[0116] 则量化后的实际频率响应 又可以表述为式(13)。
[0117]
[0118] 将式(13)代入式(9)，可将系数量化问题具体化简为式(14)。
[0119]
[0120]
[0121]
[0122] 步骤四、采用分步式整数规划算法求解式(14)。
[0123] 式(14)是一个整数规划问题，由于变量 包含着级联的一维FIR子滤波器的系数变量 和 在计算频率响应时，两者是一种相乘的关系，是非线性的，所以无法采用相关的整数规划算法直接求解，提出分步式整数规划算法量化方案。
[0124] 4.1利用步骤二得出的连续系数ri_half，求解 中的量化系数，目标函数和相应的约束条件为
[0125]
[0126]
[0127]
[0128]
[0129] 式(15)中，d′Q是当前量化状态下所获得的最小峰值纹波误差，式(15.b)中除以(Ω·ρ)表示将量化的整数值转化为小数值用于计算实际的频率响应。
[0130]
[0131]
[0132] 式(15.c)表示 的动态取值范围，式(15.d)中，Ψs表示向量xs中0元素的集合，量化过程中，0元素不参与量化。采用Gurobi求解器与yalmip工具箱嵌套使用的方式求解式(15)中的
[0133] 4.2利用4.1求解出的 再求解 中的量化系数，目标函数和相应的约束条件为
[0134]
[0135]
[0136]
[0137]
[0138] 式(18)中，
[0139]
[0140]
[0141] 式(18.c)表示 的动态取值范围，式(18.d)中Ψr表示ri_half系数中的0元素的集合。采用Gurobi求解器与yalmip工具箱嵌套使用的方式求解式(15)中的
[0142] 最终获得量化后的整数系数向量
[0143] 为了本发明的有效性，对本发明进行了计算机模拟仿真。
[0144] 实例1.设计一个四分之一对称的圆型二维FIR滤波器，理想的频率响应如下：
[0145]
[0146] 其中ωp＝0.5·π，ωs＝0.7·π，滤波器的阶数为N2＝17×17。频率采样点数Γ取值为1521；K的取值为4。步骤二的参数ε＝0.15，迭代次数为15。系数量化的有限字长Q的取值范围为5～20比特，ρ＝1.59689438141。连续系数状态下的频率响应见图3；表1为可分离二维FIR圆型滤波器舍入法和分步式整数规划算法下的量化误差结果，量化误差曲线图见图5；表2和表3是可分离二维FIR圆型滤波器的连续系数；表4和表5分别是使用舍入法和分步式整数规划法得出的量化系数。
[0147] 表1圆型滤波器舍入法和分步式整数规划算法下的量化误差结果
[0148]
[0149]
[0150] 表2可分离二维FIR圆型滤波器的归一化的连续系数(ρ·ri_half)
[0151]
[0152] 表3可分离二维FIR圆型滤波器的归一化的连续系数(ρ·si_half)
[0153]
[0154] 表4可分离二维FIR圆型滤波器系数在两种量化法下的量化系数( Q＝9)
[0155]
[0156] 表5可分离二维FIR圆型滤波器系数在两种量化法下的量化系数( Q＝9)
[0157]
[0158] 实例2.设计一个四分之一对称的钻石型二维FIR滤波器，理想的频率响应如下：
[0159]
[0160] 其中ωp＝0.6·π，ωs＝π，滤波器的阶数为N2＝17×17阶。对于不同的阶数，频率采样点数Γ取值1521；K的取值依次为5。步骤二的参数ε＝0.15，迭代次数为15。系数量化的有限字长Q的取值范围为5～20比特，ρ＝1.706190735。连续系数状态下的频率响应见图4；表6为可分离二维FIR钻石型滤波器舍入法和分步式整数规划算法下的量化误差结果，量化误差曲线图见图6；表7和表8是可分离二维FIR钻石型滤波器的连续系数；表9和表10分别是使用舍入法和分步式整数规划法得出的量化系数。
[0161] 表6钻石型滤波器舍入法和分步式整数规划算法下的量化误差结果
[0162]
[0163] 表7可分离二维FIR钻石型滤波器的归一化的连续系数(ρ·ri_half)
[0164]
[0165] 表8可分离二维FIR滤波器的归一化的连续系数(ρ·si_half)
[0166]
[0167] 表9可分离二维FIR圆型滤波器系数在两种量化法下的量化系数( Q＝11)
[0168]
[0169] 表10可分离二维FIR圆型滤波器系数在两种量化法下的量化系数( Q＝11)
[0170]
[0171]
[0172] 从表1和表6中可以看出，本发明所设计的可分离二维FIR滤波器的分步式整数规划算法量化方法明显优于舍入量化法，特别是在量化比特Q较小的情况下，差距更为明显。可见，本发明设计的方法能有效降低可分离二维FIR滤波器的系数量化误差，获得最优的有限字长属性。
