高阶有限冲击响应数字滤波器设计方法

摘要

本发明公开了一种高阶有限冲击响应数字滤波器设计方法，包括以下步骤：利用神经网络方法，通过建立待设计的滤波器与理想滤波器的幅频响应误差平方和作为神经网络的计算能量函数，采用梯度下降学习算法训练神经网络的权值使待设计滤波器与理想滤波器的幅频响应误差平方和最小化，当神经网络稳定时即可获得有限冲击响应数字滤波器各项参数，即完成有限冲击响应数字滤波器的设计。本发明方法所设计的高阶有限冲击响应数字滤波器具有通带波动小、阻带衰减大、边界频率可控、精度高等特点，特别是设计中不需要进行矩阵求逆运算、计算速度快，在数据传输、高精度电视、雷达和声纳系统、语音和图像处理等工程领域具有广泛的应用前景。

说明书

技术领域

本发明涉及一种滤波器的设计方法，特别涉及一种高阶有限冲击响应数字滤波器设计方法。

背景技术

有限冲击响应数字滤波器具有无限冲击响应数字滤波器难以实现的线性相位特性，在数据传输、高精度电视、雷达和声纳系统、语音和图像处理等工程领域具有广泛的应用前景。

采用窗口设计法与频率采样设计法来实现的有限冲击响应数字滤波器，在实际应用中不易精确控制通带与阻带边界频率。一些设计方法，例如利用最大误差最小化准则的Remez交换算法与线性规划算法，需要计算一个矩阵的逆，当滤波器的阶数很高时，矩阵的求逆会很困难。

发明内容

为了解决现有有限冲击响应数字滤波器设计存在的上述技术问题，本发明提供一种高阶有限冲击响应数字滤波器设计方法。采用本发明方法设计的有限冲击响应数字滤波器具有高精度和可控制性。

本发明解决上述的技术问题的技术方案包括以下步骤：

利用神经网络方法，通过建立待设计的滤波器与理想滤波器的幅频响应误差平方和作为神经网络的计算能量函数，采用梯度下降学习算法训练神经网络的权值使待设计滤波器与理想滤波器的幅频响应误差平方和最小化，当神经网络稳定时即可获得有限冲击响应数字滤波器各项参数，即完成有限冲击响应数字滤波器的设计。

本发明的技术效果在于：本发明方法设计的高阶有限冲击响应数字滤波器具有通带波动小、阻带衰减大、边界频率可控、精度高等特点，特别是设计中不需要进行矩阵求逆运算、计算速度快，在数据传输、高精度电视、雷达和声纳系统、语音和图像处理等工程领域具有广泛的应用前景。

下面结合附图及实施例对本发明作详细的说明。

附图说明

图1为本发明设计有限冲击响应数字滤波器所应用的神经网络模型。

图2为本发明应用示例1的有限冲击响应数字陷波器幅频响应。

图3为本发明应用示例2的有限冲击响应数字低通滤波器幅频响应及幅频响应误差。

图4为本发明应用示例3的有限冲击响应数字低通微分器幅频响应及幅频响应误差。

图5为本发明应用示例4的有限冲击响应数字窄带滤波器幅频响应。

具体实施方式

本发明中首先考虑所设计的N-1阶有限冲击响应数字滤波器的系统函数为$H\left(z\right)=\underset{n=0}{\overset{N-1}{\Sigma }}h\left(n\right)z^{-n},$其中h(n)为滤波器的冲激响应，N为奇数，且h(n)＝h(N-1-n)。

此时$H\left(e^{\mathrm{j\Omega }}\right)=e^{\mathrm{j\theta }\left(\Omega \right)}H\left(\Omega \right)=e^{-j\frac{N-1}{2}\Omega }\underset{n=0}{\overset{\frac{N-1}{2}}{\Sigma }}a\left(n\right)\cos \left(\mathrm{n\Omega }\right),$线性相位为$\theta \left(\Omega \right)=-\frac{N-1}{2}\Omega ,$幅频响应为$H\left(\Omega \right)=\underset{n=0}{\overset{\frac{N-1}{2}}{\Sigma }}a\left(n\right)\cos \left(\mathrm{n\Omega }\right),$其中$a\left(0\right)=h\left(\frac{N-1}{2}\right),$$a\left(n\right)=2h(\frac{N-1}{2}-n),n=\mathrm{1,2},...,\frac{N-1}{2}.$

该有限冲击响应数字滤波器系数a(n)和h(n)采用下述神经网络算法获取：

神经网络结构参见图1，该神经网络为三层：输入层、隐层和输出层，其中每个输入层和第i个隐层神经元的权值分别为1和w_i，隐层第i个神经元的激励函数为f_i＝cos[iΩ(m)]，i＝0，1，...，n。

神经网络能量函数为$V\left(e\right)=\frac{1}{2}\underset{m=1}{\overset{M}{\Sigma }}{e}^2\left(m\right),$式中e＝M_d-H_d，$H_d={\mathrm{ct}}_1^{T}W,$${\mathrm{ct}}_1={[c_1^T\left(\Omega \left(1\right)\right);c_1^T\left(\Omega \left(2\right)\right);...;c_1^T\left(\Omega \left(M\right)\right)]}^T,$M_d、H_d分别为待求的与理想的滤波器幅频响应，$H_d\left(m\right)=\underset{i=0}{\overset{n}{\Sigma }}W\left(i\right)\cos \left[\mathrm{i\Omega }\left(m\right)\right]=W^T{c}_1\left[\Omega \left(m\right)\right].$

神经网络的权值调整方式W(k+1)＝W(k)+ηct₁e，其中η为学习率，$0<\eta <\frac{2}{n+1},$n+1为隐层神经元个数。

当神经网络稳定时即V(e)最小，此时由神经网络权值即可获得滤波器系数：W＝[a(0)，a(1)，...，a(n)]^T，c₁(Ω)＝[1，cos(Ω)，...，cos(nΩ)]^T，$n=\frac{N-1}{2}.$

应用示例1：

为便于比较和验证本发明方法的有效性，本发明以文献[Pavel Zahradnikand Miroslav Fast analytical design algorithms for FIR notch filters]，IEEETrans.Circuits and Syst.I，vol.51(3)，pp.608-623，2004]中的例1为例，设计一个80阶的陷波频率Ω_m＝0.35π的有限冲击响应陷波器。

首先设N＝81，ε＝+10^-10，η＝0.0244，并在Ω∈[0，π]内均匀采样以获得41个训练样本，例如${\Omega }_l=\frac{\pi }{40}l,$l＝0，1，2，...，40，并设初始权值为随机的。为使通带和阻带无过冲和波动现象，在过渡带中分别取两个样本，其幅值分别为0.78和0.25。按上步骤，经21次训练，神经网络收敛，所设计的有限冲击响应陷波器幅频响应见图2。实际陷波器的参数为Ω_m＝0.3499π和ΔΩ＝0.0924π(相对于a＝-3.0103dB)。而文献[IEEE Trans.Circuits and Syst.I，vol.51(3)，pp.608-623，2004]中的例1用解析算法设计指标相同的89阶滤波器，其实际滤波器参数为Ω_m＝0.3498π和ΔΩ＝0.1496π(相对于a＝-3.0103dB)。显见，利用本发明方法设计的有限冲击响应陷波器更逼近理想情况。

其次，本发明考虑所设计的N-1阶有限冲击响应数字滤波器的系统函数为$H\left(z\right)=\underset{n=0}{\overset{N-1}{\Sigma }}h\left(n\right)z^{-n},$其中h(n)为滤波器的冲激响应，N为偶数，且h(n)＝h(N-1-n)。

此时$H\left(e^{\mathrm{j\Omega }}\right)=e^{\mathrm{j\theta }\left(\Omega \right)}H\left(\Omega \right)=e^{-j\frac{N-1}{2}\Omega }\underset{n=1}{\overset{\frac{N}{2}}{\Sigma }}b\left(n\right)\cos \left[(n-\frac{1}{2})\Omega \right],$线性相位为$\theta \left(\Omega \right)=-\frac{N-1}{2}\Omega ,$幅频响应为$H\left(\Omega \right)=\underset{n=1}{\overset{\frac{N}{2}}{\Sigma }}b\left(n\right)\cos \left[(n-\frac{1}{2})\Omega \right],$其中$b\left(n\right)=2h(\frac{N}{2}-n),n=\mathrm{1,2},...,\frac{N}{2}..$

该有限冲击响应数字滤波器系数b(n)和h(n)采用下述神经网络(结构参见图1)算法获取：

神经网络能量函数为$V\left(e\right)=\frac{1}{2}\underset{m=1}{\overset{M}{\Sigma }}{e}^2\left(m\right),$式中e＝M_d-H_d，$H_d={\mathrm{ct}}_2^{T}W,$${\mathrm{ct}}_2={[c_2^T\left(\Omega \left(1\right)\right);c_2^T\left(\Omega \left(2\right)\right);...;c_2^T\left(\Omega \left(M\right)\right)]}^T,$M_d、H_d分别为待求的与理想的滤波器幅频响应。

神经网络的权值调整方式W(k+1)＝W(k)+ηct₂e，其中η为学习率，$0<\eta <\frac{2}{n},$n为神经网络隐层神经元的个数。

当神经网络稳定时即V(e)最小，此时由神经网络权值即可获得滤波器系数：W＝[b(1)，a(2)，...，b(n)]^T，c₂(Ω)＝[cos(0.5ω)，...，cos((n-0.5)Ω)]^T，$n=\frac{N}{2};$

应用示例2：

为便于比较和验证本发明方法的有效性，本发明以文献[Soo-Chang Peiand Peng-Hua Wang.Design of equiripple FIR filters with constraint using amultiple exchange algorithm，IEEE Trans.Circuits and Syst.I，vol.49(1)，pp.113-116，2002]中的例1为例，设计一个77阶、通带为[0，0.4π]、阻带为[0.45π，π]的有限冲击响应低通滤波器。

设N＝78，ε＝+10^-5，η＝0.02564，并在Ω∈[0，π]内均匀采样以获得40个训练样本，例如${\Omega }_l=\frac{\pi }{39}l,$l＝0，1，2，...39，设神经网络初始权值为随机的，并在过渡带中分别取两个样本，其幅值分别为0.78和0.25。经12次训练，神经网络收敛，所设计的有限冲击响应滤波器幅频响应及幅频响应误差见图3。从图中可看出，采用本发明方法设计的滤波器在低通带和高阻带的波动与文献[IEEE Trans.Circuits and Syst.I，vol.49(1)，pp.113-116，2002]用多变换算法设计的78阶滤波器类似。

第三，本发明考虑所设计的N-1阶有限冲击响应数字滤波器的系统函数为$H\left(z\right)=\underset{n=0}{\overset{N-1}{\Sigma }}h\left(n\right)z^{-n},$其中h(n)为滤波器的冲激响应，N为奇数，且h(n)＝-h(N-1-n)。

此时$H\left(e^{\mathrm{j\Omega }}\right)=e^{\mathrm{j\theta }\left(\Omega \right)}H\left(\Omega \right)=e^{j(\frac{\pi }{2}-\frac{N-1}{2}\Omega )}\underset{n=1}{\overset{\frac{N-1}{2}}{\Sigma }}c\left(n\right)\sin \left(\mathrm{n\Omega }\right),$线性相位为$\theta \left(\Omega \right)=\frac{\pi }{2}-\frac{N-1}{2}\Omega ,$幅频响应为$H\left(\Omega \right)=\underset{n=1}{\overset{\frac{N-1}{2}}{\Sigma }}c\left(n\right)\sin \left(\mathrm{n\Omega }\right),$其中$c\left(n\right)=2h(\frac{N-1}{2}-n),n=\mathrm{1,2},...,\frac{N-1}{2}.$

该有限冲击响应数字滤波器系数c(n)和h(n)采用下述神经网络(结构参见图1)算法获取：

神经网络能量函数为$V\left(e\right)=\frac{1}{2}\underset{m=1}{\overset{M}{\Sigma }}{e}^2\left(m\right),$式中e＝M_d-H_d，$H_d=c{t}_3^{T}W,$${\mathrm{ct}}_3={[c_3^T\left(\Omega \left(1\right)\right);c_3^T\left(\Omega \left(2\right)\right);...;c_3^T\left(\Omega \left(M\right)\right)]}^T,$M_d、H_d分别为待求的与理想的滤波器幅频响应。

神经网络的权值调整方式W(k+1)＝W(k)+ηct₃e，其中η为学习率，$0<\eta <\frac{4}{n+1},$n为神经网络隐层神经元的个数。

当神经网络稳定时即V(e)最小，此时由神经网络权值即可获得滤波器系数：W＝[c(1)，c(2)，...，c(n)]^T，c₃(Ω)＝[sin(Ω)，...，sin(nΩ)]^T，$n=\frac{N-1}{2}.$

应用示例3：

为便于比较和验证本发明方法的有效性，本发明以文献[IEEE Trans.Circuits and Syst.I，vol.49(1)，pp.113-116，2002]中的例2为例，设计一个78阶低通频率响应的线性相位数字微分器。该低通微分器在低频带具有理想的频率响应jΩ，在高频带具有零频率响应。且低频带边界频率为0.8π，高频带边界频率为0.85π。

设N＝79，ε＝+10^-10，η＝0.05，并在Ω∈[0，π]内均匀采样以获得40个训练样本，例如${\Omega }_l=\frac{\pi }{39}l,$l＝0，1，2，...39。设神经网络初始权值为随机的，并在过渡带中分别取两个样本，其幅值分别为2.0和0.6。经5次训练，神经网络收敛，所设计的FIR滤波器幅频响应及幅频响应误差见图4。从图中可看出，所设计的滤波器在通带和阻带的波动小于0.015，而文献[IEEE Trans.Circuits and Syst.I，vol.49(1)，pp.113-116，2002]用多变换算法设计的78阶滤波器在通带和阻带的波动大于0.02。表明本发明方法所设计的滤波器精度高，更有效。

最后，本发明考虑所设计的N-1阶有限冲击响应数字滤波器的系统函数为$H\left(z\right)=\underset{n=0}{\overset{N-1}{\Sigma }}h\left(n\right)z^{-n},$其中h(n)为滤波器的冲激响应，N为偶数，且h(n)＝-h(N-1-n)。

此时，$H\left(e^{\mathrm{j\Omega }}\right)=e^{\mathrm{j\theta }\left(\Omega \right)}H\left(\Omega \right)=e^{j(\frac{\pi }{2}-\frac{N-1}{2}\Omega )}\underset{n=1}{\overset{\frac{N}{2}}{\Sigma }}d\left(n\right)\sin \left[(n-\frac{1}{2})\Omega \right],$线性相位为$\theta \left(\Omega \right)=\frac{\pi }{2}-\frac{N-1}{2}\Omega ,$幅频响应为$H\left(\Omega \right)=\underset{n=1}{\overset{\frac{N}{2}}{\Sigma }}d\left(n\right)\sin \left[(n-\frac{1}{2})\Omega \right],$其中$d\left(n\right)=2h(\frac{N}{2}-n),n=\mathrm{1,2},...,\frac{N}{2}..$

该有限冲击响应数字滤波器系数d(n)和h(n)采用下述神经网络(结构参见图1)算法获取：

经网络能量函数为$V\left(e\right)=\frac{1}{2}\underset{m=1}{\overset{M}{\Sigma }}{e}^2\left(m\right),$式中e＝M_d-H_d，$H_d={\mathrm{ct}}_4^{T}W,$${\mathrm{ct}}_4=[c_4^T\left(\Omega \left(1\right)\right);c_4^T\left(\Omega \left(2\right)\right);...;c_4^T\left(\Omega \left(M\right)\right){]}^T,$M_d、H_d分别为待求的与理想的滤波器幅频响应。

神经网络的权值调整方式W(k+1)＝W(k)+ηct₄e，其中η为学习率，$0<\eta \le \frac{2}{n},$n为神经网络隐层神经元的个数。

当神经网络稳定时即V(e)最小，此时由神经网络权值即可获得滤波器系数：W＝[d(1)，d(2)，...，d(n)]^T，c₄(Ω)＝[sin(0.5Ω)，...，sin((n-0.5)Ω)]^T，

应用示例4：

为便于比较和验证本发明方法的有效性，本发明设计一个79阶窄带通线性相位有限冲击响应数字滤波器，其中心频率为Ω_m＝0.35π。

设N＝80，ε＝+10^-5，η＝0.025，并在Ω∈[0，π]内均匀采样以获得41个训练样本，例如${\Omega }_l=\frac{\pi }{40}l,$l＝0，1，2...，40。设神经网络初始权值为随机的，并在过渡带中分别取两个样本，其幅值分别为0.5999和0.015。经13次训练，神经网络收敛，所设计的有限冲击响应数字滤波器幅频响应见图5。实际滤波器参数为Ω_m＝0.3500π和ΔΩ＝0.0580π(a＝-6.0206dB)。可见本发明方法设计的精度高，几乎与理想情况一致。