CN105429610A - 一种基于spt系数的fir滤波器优化方法 - Google Patents

Abstract

本发明公开了一种基于SPT系数的FIR滤波器优化方法，(1)首先得到具有无限字长的FIR滤波器系数向量，对无限字长滤波器系数向量中的元素按照绝对值从大到小进行排序；(2)设定滤波器系数增益取值空间，对无限字长滤波器系数向量进行截断处理，截断步骤在每次循环中同时分配多个SPT项，从而大幅度提高了分配效率；(3)选出使得所需SPT项最少的情况下所对应的滤波器系数增益；(4)利用贪婪算法进行优化，在截断步骤基础上进一步改进滤波器的性能，降低Minimax误差，从而充分利用SPT项提升滤波器设计性能。本发明能够取得良好的滤波器性能及较高的计算效率。

Description

一种基于SPT系数的FIR滤波器优化方法

技术领域

本发明涉及一种FIR滤波器优化方法，具体涉及一种基于SPT系数的FIR滤波器优化方法。

背景技术

有限冲激响应(FiniteImpulseResponse，FIR)滤波器是广泛运用于通信与信号处理的滤波器之一，它具有很多优点。

首先FIR滤波器具有固有的稳定性，与无限冲激响应(InfiniteImpulseResponse，IIR)滤波器相比，其极点均位于Z平面原点；其次可方便的利用对称结构，设计者容易获得具有精确线性相位特性的滤波器，避免信号发生相位失真；以及其结构相对简单，便于实现。

传统的FIR滤波器设计基于无限字长这一基本假设。此时，许多的FIR滤波器设计可以表述为凸优化问题，利用较为成熟的凸优化算法，从而获得设计问题的全局最优解。基于集成电路技术，数字FIR滤波器的系数可采用SPT形式构成，但是设计问题也发生了本质的变化，变成了离散空间中的组合优化问题。然而，在离散空间中寻找最优解是极度耗时的，因此必须要在优化质量与优化效率上达到折中。通常使用增益q来优化滤波器性能，这样处理既不影响SPT分配，又能保证滤波器系数的动态范围在设定的动态范围内，尽可能的减少所需要的SPT项的数量以达到减小滤波器实现成本的目的。然而传统方法由于每次只寻找一个分配位置进行SPT项的分配，因此传统算法SPT项计算效率不高。

发明内容

针对现有技术存在的不足，本发明目的是提供一种基于SPT系数的FIR滤波器优化方法，在给定系数表示中SPT项总数的情况下，在保证精确的前提下，将每个SPT项快速且合理地分配给滤波器系数，能够取得良好的滤波器性能及较高的计算效率。

为了实现上述目的，本发明是通过如下的技术方案来实现：

本发明的一种基于SPT系数的FIR滤波器优化方法，其包括以下步骤：

(1)首先得到具有无限字长的FIR滤波器系数向量g，对g中的元素按照绝对值从大到小的顺序进行排序，得到排好序的向量h；

(2)对向量h进行截断处理得到截断后的滤波器系数向量h'，设定滤波器系数增益q取值空间□＝{q₁,q₂,q₃...q_P}，对于计算增益后的滤波器系数向量e＝qh；截断步骤在每次循环中同时分配多个SPT项；

截断后的滤波器系数向量h'可表示为：

h'＝Cs,C∈C(M)(1)

其中，C为N×2K矩阵，其元素为0或1，N为滤波器阶数，K为除去符号位的字长，s＝[-2^-1-2^-2...-2^-K2^-12^-2...2^-K]^T，M为给定SPT总数；

(3)选出使得所需SPT项最少的情况下所对应的滤波器系数增益q₀以及

选择标准为：

ρ(C_qs)表示成本函数，即所需要的SPT项总数；

(4)利用贪婪算法对进行优化，通过试探性的对中其中的一项进行取反，判断此时最大误差是否得到改进，从而决定是否对这一项进行调整，逐项遍历中的所有可调项；

优化模型为：

$max\left|\right|s^T{C_h}^{T}v\left(\mathrm{\ω}\right)|-q_0{H}_d\left(\mathrm{\ω}\right)|=min\[\underset{\mathrm{\ω}\∈{\mathrm{\Ω}}_I}{max}\left|\right|s^T{C_{q_0}}^{T}v\left(\mathrm{\ω}\right)|-q_0{H}_d\left(\mathrm{\ω}\right)|\]---\left(3\right)$

其中，v(ω)＝[1e^-jωe^-j2ω...e^-jNω]，N为滤波器阶数，H_d(ω)代表滤波器的理想幅度响应，C_h由经过贪婪算法优化最终得到，Ω_I表示所研究的频率段。

上述步骤(1)通过Remezexchange算法得到具有无限字长的FIR滤波器系数向量g。

步骤(1)中，对g中的元素按照绝对值从大到小的顺序进行排序时，保存排序前与排序后滤波器系数的对应性，方便后续步骤中进行相应的SPT分配。

步骤(2)中，对向量h进行截断处理的方法如下,p用来指定滤波器系数增益空间□＝{q₁,q₂,q₃...q_P}中的q值，m用来表示所用的SPT项数量，k用来指定2的幂，集合□₁与□₂用来存储待分配SPT项的系数的序号，向量e中任一元素可用e(i)表示，i∈{1,2,...}，集合□₁与□₂中的元素可分别用j₁、j₂来指定；

步骤2-1：初始化p＝0；

步骤2-2：p＝p+1，计算e＝q_ph，初始化m＝0，k＝0以及矩阵C＝[0]_N×2K；

步骤2-3：k＝k+1，计算n只是一个中间变量，用来表示从k+1到K的取值，并初始化和

步骤2-4：对于

(4)如果e(i)＞2^-k或|2^-k-e(i)|≤|sum-e(i)|，□₁＝□₁∪{i}；

(5)如果e(i)＜-2^-k或|-2^-k-e(i)|≤|-sum-e(i)|，□₂＝□₂∪{i}；

步骤2-5：

(6)如果M-m≥|□₁|+|□₂|，对于C(j₁,k+K)＝1，e(j₁)＝e(j₁)-2^-k；

对于C(j₂,k)＝1，e(j₂)＝e(j₂)+2^-k，m＝m+|□₁|+□₂|；

(7)如果M-m＜|□₁|+|□₂|，对e中元素按照绝对值大小排序得e'，在e'中从前到后选出M-m个系数，分配SPT项，若此时p＜P，转到步骤2-2，否则终止算法；

步骤2-6：如果k＜K，转至步骤2-3；

步骤2-7：如果p＜P，转至步骤2-2；否则，终止算法。

步骤(4)中，基于C的对称性，只考虑中可调部分，记为A，表示前A行，式(3)可转换为以下求解步骤：

步骤4-1：计算幅度响应以及最大幅度响应误差 ${\mathrm{\Δ}}_M=\underset{\mathrm{\ω}\∈{\mathrm{\Ω}}_I}{max}|H\left(\mathrm{\ω}\right)-q_0{H}_d\left(\mathrm{\ω}\right)|$ 及其对应的频率点ω_M；

步骤4-2：初始化k＝0，r＝1，k表示取矩阵C的特定列，r表示取矩阵C的特定行；

步骤4-3：k＝k+1，若k＝2K+1，则令k＝1，r＝r+1；首先对取反，并根据对称性修改m；

步骤4-4：计算此时的幅度响应误差和对应的频率点如果

且m≤M，则 ${\mathrm{\Δ}}_M={\widetilde{\mathrm{\Δ}}}_M,{\mathrm{\ω}}_M={\widetilde{\mathrm{\ω}}}_M;$ 否则，复原

步骤4-5：如果k＝2K且r＝A，终止算法；否则，转到步骤4-3。

本发明的有益效果如下：

本发明将具有SPT系数的滤波器SPT项计算效率与设计误差进行联合考虑，滤波器系数截断步骤大幅度提高了SPT项分配的效率，在每次循环中并不只分配一个SPT项，而是寻找多个位置，同时分配多个SPT项；而贪婪算法步骤降低设计误差，在截断步骤基础上进一步改进滤波器的性能，降低Minimax误差，从而充分利用SPT项提升滤波器设计性能，因此本发明能达到较高的计算效率与良好的FIR滤波器设计效果。

附图说明

图1为本发明的工作流程图；

图2为本发明应用示例1中的线性相位低通滤波器幅度响应；

图3为本发明应用示例2中的线性相位高通滤波器幅度响应。

具体实施方式

为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。

一、模型建立

N阶FIR系统的传递函数由下式给出：

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

其中，为FIR滤波器的冲激响应，即FIR滤波器系数。如果待设计FIR滤波器具有线性相位，那么其幅度响应可由下式计算所得

其中{ν_n(ω)}_n表示三角函数，例如设计Ⅰ型线性相位FIR滤波器时，{ν_n(ω)}_n为：

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

待设计滤波器需要满足以下约束条件：

其中δ为给定的误差门限，H_d(ω)代表滤波器的理想幅度响应。假设C(M)表示N×2K矩阵集合，其元素为0或1，K为不包括符号位的字长，且设为1的比特位总数应不超过M，若s＝[-2^-1-2^-2...-2^-K2^-12^-2...2^-K]^T，则FIR滤波器系数可表示为：h'＝Cs，C∈C(M)。可以看到，FIR滤波器系数可由2的整数次幂之和计算所得，非常便于硬件实现。在给定字长情况下，SPT项分配的问题就成为了关于矩阵C的优化问题。设g为利用传统设计算法所获得的具有无限字长的最优滤波器系数向量，q为滤波器增益系数，q的取值空间为□。对于□中每个q，截断过程得到的C_q应与g具有最短的切比雪夫距离，即满足以下条件：

||C_qs-qg||_∞＝min||Cs-qg||_∞,C∈C(M)(5)

其中||.||_∞表示向量的无穷范数。然后，在预先设定的q的取值空间内，寻找最优的q₀∈□，使其满足：

ρ(C_qS)表示成本函数，即所需要的SPT项总数。截断过程得到的是我们应该运用贪婪算法来进一步处理的目标，通过调节来降低Minimax误差，因此，经最小化最大误差方法优化后得到的C_h应满足以下条件：

$max\left|\right|s^T{C_h}^{T}v\left(\mathrm{\ω}\right)|-q_0{H}_d\left(\mathrm{\ω}\right)|=min\[\underset{\mathrm{\ω}\∈{\mathrm{\Ω}}_I}{max}\left|\right|s^T{C_{q_0}}^{H}v\left(\mathrm{\ω}\right)|-q_0{H}_d\left(\mathrm{\ω}\right)|\]---\left(7\right)$

其中v(ω)＝[1e^-jωe^-j2ω...e^-jNω]，N为滤波器阶数，H_d(ω)代表滤波器的理想幅度响应，C_h由C_q0经过优化最终得到，Ω_I表示我们感兴趣的频率段。

二、问题求解

首先得到具有无限字长的FIR滤波器系数向量g，对g中元素按照绝对值从大到小排序，得到排好序的向量h。

对无限字长滤波器系数向量g按照绝对值从大到小排序得到h时，保存排序前与排序后滤波器系数的对应性，方便后续步骤中进行相应的SPT分配。

为了提高截断步骤的效率，截断步骤在每次循环中并不只分配一个SPT项，而是寻找多个位置，同时分配多个SPT项，从而大幅度提高了分配效率。在对h进行截断处理时，分为以下步骤。

步骤2-1：初始化p＝0。

步骤2-2：p＝p+1，计算e＝q_ph，初始化m＝0，k＝0以及矩阵C＝[0]_N×2K。

步骤2-3：k＝k+1，计算并初始化和

步骤2-4：对于

(8)如果e(i)＞2^-k或|2^-k-e(i)|≤|sum-e(i)|，□₁＝□₁∪{i}；

(9)如果e(i)＜-2^-k或|-2^-k-e(i)|≤|-sum-e(i)|，□₂＝□₂∪{i}。

步骤2-5：

(10)如果M-m≥|□₁|+|□₂|，对于C(j₁,k+K)＝1，e(j₁)＝e(j₁)-2^-k；

对于C(j₂,k)＝1，e(j₂)＝e(j₂)+2^-k，m＝m+|□₁|+|□₂|。

(11)如果M-m＜|□₁|+|□₂|，对e中元素按照绝对值大小排序得e'，在e'中从前到后选出M-m个系数，分配SPT项，若此时p＜P，转到步骤2-2，否则终止算法。

步骤2-6：如果k＜K，转至步骤2-3。

步骤2-7：如果p＜P，转至步骤2-2；否则，终止算法。

依据式(12)选出使得所需SPT项最少的情况下所对应的增益q₀以及

为了对设计模型式(13)进行方便求解，将其转换为以下求解方式：

利用贪婪算法对进行优化，在截断步骤基础上进一步改进滤波器的性能，降低Minimax误差，从而充分利用SPT项提升滤波器设计性能。可通过试探性的对中某一项进行取反(从0变为1或从1变为0)判断此时最大误差是否得到改进来决定是否对这一项进行调整，逐项遍历中的所有可调项。其可按以下步骤进行：

步骤4-1：计算幅度响应以及最大幅度响应误差

${\mathrm{\Δ}}_M=\underset{\mathrm{\ω}\∈{\mathrm{\Ω}}_I}{max}|H\left(\mathrm{\ω}\right)-q_0{H}_d\left(\mathrm{\ω}\right)|$ 及其对应的频率点ω_M。

步骤4-2：初始化k＝0，r＝1。

步骤4-3：k＝k+1。若k＝2K+1，则令k＝1，r＝r+1。首先对取反，并根据对称性修改m。

步骤4-4：计算此时的幅度响应误差和对应的频率点如果

且m≤M，则 ${\mathrm{\Δ}}_M={\widetilde{\mathrm{\Δ}}}_M,{\mathrm{\ω}}_M={\widetilde{\mathrm{\ω}}}_M;$ 否则，复原

步骤4-5：如果k＝2K且r＝A，终止算法；否则，转到步骤4-3。

三、设计示例

设计方法的可行性可用以下两个示例进行验证。所有设计均采用归一化频率，即二分之一采样频率对应于ω＝π。

示例1：示例1为一线性相位低通滤波器设计，表1为示例1所用的各项参数，图2所示为线性相位低通滤波器的幅度响应，可以看到，在滤波器阶数较高情况下，该方法仍能取得良好的滤波器性能及较高的计算效率。

表1

表2

示例2：示例2为一线性相位高通滤波器设计，表2为示例2所用的各项参数，图3所示为线性相位高通滤波器的幅度响应，可以看到，在滤波器阶数较高情况下，该方法仍能取得良好的滤波器性能及较高的计算效率。

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (5)

1.一种基于SPT系数的FIR滤波器优化方法，其特征在于，其包括以下步骤：

(1)首先得到具有无限字长的FIR滤波器系数向量g，对g中的元素按照绝对值从大到小的顺序进行排序，得到排好序的向量h；

(2)对向量h进行截断处理得到截断后的滤波器系数向量h'，设定滤波器系数增益q取值空间对于计算增益后的滤波器系数向量e＝qh；截断步骤在每次循环中同时分配多个SPT项；

截断后的滤波器系数向量h'可表示为：

h'＝Cs,C∈C(M)(1)

其中，C为N×2K矩阵，其元素为0或1，N为滤波器阶数，K为除去符号位的字长，s＝[-2^-1-2^-2...-2^-K2^-12^-2...2^-K]^T，M为给定SPT总数；

(3)选出使得所需SPT项最少的情况下所对应的滤波器系数增益q₀以及选择标准为：

ρ(C_qs)表示成本函数，即所需要的SPT项总数；

(4)利用贪婪算法对进行优化，通过试探性的对中其中的一项进行取反，判断此时最大误差是否得到改进，从而决定是否对这一项进行调整，逐项遍历中的所有可调项；

优化模型为：

$max\left|\right|s^T{C_h}^{T}v\left(\mathrm{\ω}\right)|-q_0{H}_d\left(\mathrm{\ω}\right)|=min\[\underset{\mathrm{\ω}\∈{\mathrm{\Ω}}_I}{max}\left|\right|s^T{C_{q_0}}^{T}v\left(\mathrm{\ω}\right)|-q_0{H}_d\left(\mathrm{\ω}\right)|\]---\left(3\right)$

其中，v(ω)＝[1e^-jωe^-j2ω...e^-jNω]，N为滤波器阶数，H_d(ω)代表滤波器的理想幅度响应，C_h由经过贪婪算法优化最终得到，Ω_I表示所研究的频率段。

2.根据权利要求1所述的基于SPT系数的FIR滤波器优化方法，其特征在于，步骤(1)中，通过Remezexchange算法得到具有无限字长的FIR滤波器系数向量g。

3.根据权利要求1所述的基于SPT系数的FIR滤波器优化方法，其特征在于，

步骤(1)中，对g中的元素按照绝对值从大到小的顺序进行排序时，保存排序前与排序后滤波器系数的对应性，方便后续步骤中进行相应的SPT分配。

4.根据权利要求1所述的基于SPT系数的FIR滤波器优化方法，其特征在于，

步骤(2)中，对向量h进行截断处理的方法如下,p用来指定滤波器系数增益空间中的q值，m用来表示所用的SPT项数量，k用来指定2的幂，集合与用来存储待分配SPT项的系数的序号，向量e中任一元素可用e(i)表示，i∈{1,2,...}，集合与中的元素可分别用j₁、j₂来指定；

步骤2-1：初始化p＝0；

步骤2-2：p＝p+1，计算e＝q_ph，初始化m＝0，k＝0以及矩阵C＝[0]_N×2K；

步骤2-3：k＝k+1，计算n只是一个中间变量，用来表示从k+1到K的取值，并初始化和

步骤2-4：对于

(4)如果e(i)＞2^-k或|2^-k-e(i)|≤|sum-e(i)|，

(5)如果e(i)＜-2^-k或|-2^-k-e(i)|≤|-sum-e(i)|，

步骤2-5：

(6)如果对于C(j₁,k+K)＝1，e(j₁)＝e(j₁)-2^-k；

对于C(j₂,k)＝1，e(j₂)＝e(j₂)+2^-k，

(7)如果对e中元素按照绝对值大小排序得e'，在e'中从前到后选出M-m个系数，分配SPT项，若此时p＜P，转到步骤2-2，否则终止算法；

步骤2-6：如果k＜K，转至步骤2-3；

步骤2-7：如果p＜P，转至步骤2-2；否则，终止算法。

5.根据权利要求3所述的基于SPT系数的FIR滤波器优化方法，其特征在于，

步骤(4)中，基于C的对称性，只考虑中可调部分，记为A，表示前A行，式(3)可转换为以下求解步骤：

步骤4-1：计算幅度响应以及最大幅度响应误差 ${\mathrm{\Δ}}_M=\underset{\mathrm{\ω}\∈{\mathrm{\Ω}}_I}{max}|H\left(\mathrm{\ω}\right)-q_0{H}_d\left(\mathrm{\ω}\right)|$ 及其对应的频率点ω_M；

步骤4-2：初始化k＝0，r＝1，k表示取矩阵C的特定列，r表示取矩阵C的特定行；

步骤4-3：k＝k+1，若k＝2K+1，则令k＝1，r＝r+1；首先对取反，并根据对称性修改m；

步骤4-4：计算此时的幅度响应误差和对应的频率点如果且m≤M，则 ${\mathrm{\Δ}}_M={\widetilde{\mathrm{\Δ}}}_M,{\mathrm{\ω}}_M={\widetilde{\mathrm{\ω}}}_M;$ 否则，复原

步骤4-5：如果k＝2K且r＝A，终止算法；否则，转到步骤4-3。