CN105680825A - 数字滤波器解析设计法及其滤波器 - Google Patents

数字滤波器解析设计法及其滤波器 Download PDF

Info

Publication number
CN105680825A
CN105680825A CN201610086020.4A CN201610086020A CN105680825A CN 105680825 A CN105680825 A CN 105680825A CN 201610086020 A CN201610086020 A CN 201610086020A CN 105680825 A CN105680825 A CN 105680825A
Authority
CN
China
Prior art keywords
filter
formula
sigma
rsqb
lsqb
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201610086020.4A
Other languages
English (en)
Inventor
黄翔东
张博
马欣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tianjin University
Original Assignee
Tianjin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tianjin University filed Critical Tianjin University
Priority to CN201610086020.4A priority Critical patent/CN105680825A/zh
Publication of CN105680825A publication Critical patent/CN105680825A/zh
Pending legal-status Critical Current

Links

Classifications

    • HELECTRICITY
    • H03ELECTRONIC CIRCUITRY
    • H03HIMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00Networks using digital techniques
    • H03H17/0063R, L, C, simulating networks
    • HELECTRICITY
    • H03ELECTRONIC CIRCUITRY
    • H03HIMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00Networks using digital techniques
    • H03H2017/0072Theoretical filter design
    • H03H2017/0081Theoretical filter design of FIR filters

Landscapes

  • Complex Calculations (AREA)

Abstract

本发明涉及数字信号处理技术领域,为实现不引入迭代优化的措施情况下,自动消除全相位滤波器的过冲,最终生成全过程完全实现解析设计的高效率、高性能的FIR滤波器设计法,并给予数字信号处理器实现。本发明采用的技术方案是,数字滤波器解析设计法及其滤波器,步骤如下:Step?1.确定边界频率整数参数值;Step?2.得到一长度为2N-1的卷积窗wc(n);Step?3.将参数M、N和卷积窗wc(n)直接代入如下解析式,获取滤波器系数g(n);Step?4.计算Lichtenberg比率频点处的幅频响应;Step?5.将T=1-(Q-1)=2-Q代入如下改进的解析式中,算出最终的可消除过冲的滤波器。本发明主要应用于数字信号处理场合。

Description

数字滤波器解析设计法及其滤波器
技术领域
本发明涉及数字信号处理技术领域,尤其涉及一种基于Lichtenberg比率的FIR滤波器解析设计法。
背景技术
在有限冲击响应滤波器(FiniteImpulseResponseFilter,FIRfilter)的设计中,兼顾滤波器良好传输性能(即保证通带波纹足够小和阻带衰减足够大)和滤波器的设计效率一直是个技术难题。无论是经典滤波器设计法、经典优化设计法,还是现代滤波器进化设计法,这个问题都很突出。
经典滤波器设计法,例如窗函数法,可以将边界频带参数ωc直接代入理想滤波器公式得到滤波器系数,但是由于理想滤波器系数是无限长的,因而只能对理想滤波器进行截断,在截断过程中会引入吉布斯(Gibbs)效应[1]而导致滤波器在边界频带附近处的通带、阻带传输曲线出现很大的振荡。加窗虽然可以减轻传输曲线的振荡,但是会导致滤波器过渡带的加宽和边界频带的模糊。再如频率采样法也存在同样的问题,该方法是通过对频率响应向量H直接作傅里叶反变换而得到滤波器系数,虽然可以通过在H的不同位置处设置相应的0、1值来控制边界频带,但是这同样会导致滤波器传输曲线的通带和阻带出现很大的振荡。加过渡点可以减轻这些振荡,但是这是以加宽过渡带、模糊边界频带位置作为代价的。
经典优化设计法通常是在某个数学准则(如最小均方误差准则、切比雪夫等波纹逼近准则)下,通过迭代优化滤波器幅频曲线,使之逼近理想传输特性来实现的。常见的有Parks-McClellan方法、WLS方法[2,3]、神经网络法[4,5]等,这些方法在设计滤波器的优秀传输性能方面具有优势,但是,由于通过进行多次参数迭代直至收敛的全局优化过程,他们往往很难达到较高的效率。例如,Parks-McClellan方法需要对多个频点进行大量迭代才能获得一个等波纹的逼近。往往这些优化过程是针对滤波器的所有系数进行,其阶数越高,需要的计算复杂度越高,这样无疑会耗费大量的资源,在一些需要快速设计的场合如软件无线电、多速信号处理等,经典优化算法是不适用的。
对于近些年出现的现代滤波器进化优化算法(如GA[6,7]、PSO、DE[8,9]、CSO[10,11]等),其核心思想是模拟自然界的生物选择和进化过程,按照“优胜劣汰、适者生存”的法则开发出启发式的搜索算法。为了在进化过程中寻找到全局最优路径,进化算法需要建立大量的粒子种群(代表滤波器系数),从而能够使进化过程快速的跳出局部最优并且获得全局最优。由于类似于缓慢的自然生物进化,这些进化算法耗费大量的迭代致使其同样计算缓慢,并且对资源的耗费量也较大。因此,在需要高阶滤波器以及快速响应的场合,现代的进化算法也是不适用的。
全相位滤波器设计法[12],在优化滤波器的传输性能和设计效率这两方面均具有较突出的优势[13]。全相位滤波器内含了N个子滤波器的叠加过程,这些子滤波器的频率响应在叠加中正负相互抵消,实现了幅度互补,从而保证了最终设计的滤波器传输曲线的通带波纹足够小和阻带衰减足够大[14],因而从优化角度看,全相位滤波器设计过程其实等效于滤波器全局幅频响应的优化过程,它是内在的自然优化过程,不需要像优化算法的设计那样循环迭代,因此拥有很高的设计效率。
但是现有的全相位滤波器设计法有一个缺陷,就是在传输曲线的通带边界和阻带边界附近各存在一个过冲(即存在一定程度的Gibbs效应),这两个过冲若不去除,会影响音视频信号处理、软件无线电、多速信号处理等的性能(即因通带的边界频率成分的幅值过大和阻带的边界频率成分的幅值过小,会引起波形的失真)。
参考文献
[1]高西全,丁玉美.数字信号处理(第三版)[M].西安:西安电子科技大学出版社,2008.
[2]AlgaziVR,SukM,RimC-S.DesignofalmostminimaxFIRfiltersinoneandtwodimensionsbyWLStechniques[J].CircuitsandSystems,IEEETransactionson,1986,33(6):590-596.
[3]NguyenTQ.ThedesignofarbitraryFIRdigitalfiltersusingtheeigenfiltermethod[J].SignalProcessing,IEEETransactionson,1993,41(3):1128-1139.
[4]BhattacharyaD,AntoniouA.Real-timedesignofFIRfiltersbyfeedbackneuralnetworks[J].SignalProcessingLetters,IEEE,1996,3(5):158-161.
[5]WangX,HeY,LiuM.Fouroptimaldesignapproachesofhigh-orderfinite-impulseresponsefiltersbasedonneuralnetwork[J].JournalofCentralSouthUniversityofTechnology,2007,14:94-99.
[6]AbabnehJI,BatainehMH.LinearphaseFIRfilterdesignusingparticleswarmoptimizationandgeneticalgorithms[J].DigitalSignalProcessing,2008,18(4):657-668.
[7]YeWB,YuYJ.Single-stageandcascadedesignofhighordermultiplierlesslinearphaseFIRfiltersusinggeneticalgorithm[J].CircuitsandSystemsI:RegularPapers,IEEETransactionson,2014,60(11):2987-2997.
[8]AbbassHA,SarkerR,NewtonC.PDE:aPareto-frontierdifferentialevolutionapproachformulti-objectiveoptimizationproblems[C].EvolutionaryComputation,2001.Proceedingsofthe2001Congresson.IEEE,2001,2:971-978.
[9]KarabogaN,CetinkayaB.DesignofdigitalFIRfiltersusingdifferentialevolutionalgorithm[J].Circuits,SystemsandSignalProcessing,2006,25(5):649-660.
[10]TsaiP-W,PanJ-S,ChenS-M,etal.Parallelcatswarmoptimization[C].MachineLearningandCybernetics,2008InternationalConferenceon.IEEE,2008,6:3328-3333.
[11]SahaSK,GhoshalSP,KarR,etal.CatSwarmOptimizationalgorithmforoptimallinearphaseFIRfilterdesign[J].ISAtransactions,2013,52(6):781-794.
[12]王兆华,侯正信,苏飞.全相位数字滤波[J].信号处理,2003,19(增刊):1-4.
[13]苏飞,王兆华.DFT域全相位数字滤波器的设计与实现[J].信号处理,2004,20(3):231-235.
[14]王兆华,黄翔东.数字信号全相位谱分析与滤波技术[M].北京:电子工业出版社,2009.
[15]ReddyKS,SahooSK.AnapproachforFIRfiltercoefficientoptimizationusingdifferentialevolutionalgorithm[J].AEU-InternationalJournalofElectronicsandCommunications,2015,69(1):101-108.
发明内容
为克服现有技术的不足,本发明旨在实现不引入迭代优化的措施情况下,自动消除全相位滤波器的过冲,最终生成全过程完全实现解析设计的高效率、高性能的FIR滤波器设计法,并给予数字信号处理器实现。本发明采用的技术方案是,数字滤波器解析设计法及其滤波器,指定一个满足传统奇对称H(k)=H(N-k),k=0,...,N-1的频率采样向量H=[H(0),H(1),...,H(N-1)],设置为如下形式
直接对H进行IDFT得到
采用如下定义域延拓的方法来构造子滤波器:
把式(29)中n的IDFT定义域进行延拓,使得n∈[-m,-m+N-1],即
这样获得N个子滤波器hm=[h(-m),h(-m+1),...,h(-m+N-1)],m=0,...,N-1,其傅立叶变换为
如果将这N个子滤波器做算数平均,那么得到合成后的频响G(jω)为
使用一个长度为N的常见的归一化对称窗{f(m),m=0,...,N-1},将式(5)中对Hm(jω)的简单算术平均替换为加权平均,表示为
将式(31)代入式(33)中,并且交换m和n的求和次序,得到
为了简化上式,定义一个长度为2N-1的卷积窗{wc(n)},由长度为N的对称窗{f(n)}和长度为N的反转矩形窗{RN(-n)}构成如下
wc(n)=f(n)*RN(-n),n=-N+1,...,0,...,N-1(8)
上式进一步表示为
因为{f(n)}和{RN(-n)}的非零元素都定义在区间[0,N-1]中,所以m满足
因此对式(36)分为两种情况,进一步推导可得
相应的,wc(-n)表示为
将式(39)代入式(34),得G(jω)为
由于f(n)是对称的,将f(m)=f(N-1-m)代入式(39)得
对比式(41)和式(39)可看出,wc(n)也是对称的,因此式(40)进一步推导为
与式(33)相比,全相位滤波器传输特性G(jω)得到进一步精简表示,将式(29)代入式(42)得
其中,Wc(jω)定义为
式(42)可以看出,全相位FIR滤波器的系数长度为2N-1,其滤波器系数g(n)为
g(n)=h(n)wc(n),-N+1≤n≤N-1(18)
又因为exp[j2π/N·(n+N)]=exp[j2πn/N],所以根据式(30)得出
h(n)=h(n+N)(19)
因此,定义域延拓的IDFT得到的2N-1长的滤波器系数向量h′表示为
h′=[h(-N+1),...,h(-1),h(0),h(1),...,h(N-1)]
(20)
=[h(1),...,h(N-1),h(0),h(1),...h(N-1)]
基于式(35)、式(43)和式(45),全相位FIR滤波器的设计可以总结为以下三步:
①对指定的频率采样向量H进行IDFT变换得到h=[h(0),h(1),...,h(N-1)],然后插入复制h(1)~h(N-1)为h左半部分得h′;
②选择一长度为N的常用窗f和长度为N的矩形窗b进行卷积并归一化,得到长度为2N-1的卷积窗wc
③将h′和wc点乘即得最终的全相位FIR滤波器系数g(n)。
对H(k)进行定义域延拓的IDFT得到h(n)为
将式(28)代入式(48)得到
使用欧拉公式,上式进一步推导为
因此,结合式(50)和式(43)可得最终滤波器系数g(n)为
上式中包含了n=0的情况,是通过将n=0直接代入式(48)和式(43)得到g(0)。
如果将k=M-1处的采样值1替换为一适当值T(T<1),则k=M处的采样值为1-T,更改后的频率采样向量HT具有以下的形式
明显上式只引入了一个参数T需要被决定;
类似于式(51)中滤波器系数g(n)的推导过程,推导出由改进的频率采样向量HT决定的全相位低通FIR滤波器系数gT(n)的解析式如下
由于G(jω)曲线上的凸起部分非常小,汉明窗情况下只有5.5%的幅度,因此最优T值的取值区间为T∈(0.9,1);通过上述加入过渡点优化的解析式设计,只要使用高效的方法来确定T值,同样可以快速的获得滤波器系数。
频点ωL处的凸起波纹G(ωL)-1近似认为是过渡带两侧凸起的最大值,将其用来补偿ω=(M-1)Δω频点处的采样设置值,并通过以下的三步来完成基于Lichtenberg比率的设计法,获得较优的参数T以及改进的滤波器系数gT(n):
①使用系数解析式(51)设计初始滤波器g(n),-N+1≤n≤N-1;
②计算g(n)在频点处的幅频响应Q,即
③将T=1-(Q-1)=2-Q代入改进的系数解析式(70)中,便得到最终改进的FIR滤波器gT(n),-N+1≤n≤N-1。
前述步骤进一步优化为:
Step1.给定期望的滤波器的6dB边界频率ωc和奇数长度L,算出对应的全相位滤波器阶数N=(L+1)/2和频率单位Δω=2π/N,进而确定边界频率整数参数值M=[ωc/Δω+0.5]
Step2.选择一长度为N的常用窗{f(n)},并与长度为N的反转矩形窗{RN(-n)}做卷积而得到一长度为2N-1的卷积窗wc(n),即
wc(n)=f(n)*RN(-n),n=-N+1,...,0,...,N-1
Step3.将参数M、N和卷积窗wc(n)直接代入如下解析式,获取滤波器系数g(n);
Step4.计算Lichtenberg比率频点处的幅频响应,即
Step5.将T=1-(Q-1)=2-Q代入如下改进的解析式中,算出最终的可消除过冲的滤波器
系数gT(n),-N+1≤n≤N-1,即
本发明的特点及有益效果是:
1、自动消除全相位滤波器的通带边界和阻带边界附近过冲,形成全频带无过冲的传输曲线;
2、只需通过设置相关参数就可以实现滤波器的快速设计,无需迭代优化过程,具有高效性。
3、滤波器全频带的传输曲线通带波纹小,阻带衰减大。
附图说明:
图1子滤波器及其算术平均后的传输曲线。
图2加权平均后的传输曲线|G(jω)|。
图3频率采样法传输曲线|H(jω)|。
图4RN(ω)的频响曲线。
图5Wc(jω)的幅频曲线。
图6传输曲线G(jω)和Gc(jω)。
图7不同窗情况下的传输曲线(N=12,M=4)。
图8不同窗设计下的优化前后滤波器频响曲线和衰减曲线(N=9,M=3)。
图9不同方法设计的滤波器频响曲线和衰减曲线。
图10大阶数下不同方法设计的滤波器频响曲线和衰减曲线。
图11不同阶数下三种方法设计滤波器耗费时间对比。
图12本发明的硬件实现图。
图13DSP内部程序流图。
具体实施方式
本发明所要解决的,就是在不须引入迭代优化的措施情况下,自动消除全相位滤波器的过冲,最终生成全过程完全实现解析设计的高效率、高性能的FIR滤波器设计法,并给予数字信号处理器实现。
本发明旨在实现:
1、自动消除全相位滤波器的通带边界和阻带边界附近过冲,形成全频带无过冲传输曲线;
2、只需通过设置相关参数就可以实现滤波器的快速设计,无需迭代优化过程,具高效性。
3、滤波器全频带的传输曲线通带波纹小,阻带衰减大。
1.本发明总体实现步骤
本发明提出的基于Lichtenberg比率的数字滤波器解析设计法步骤如下:
Step1.给定期望的滤波器的6dB边界频率ωc和奇数长度L,算出对应的全相位滤波器阶数N=(L+1)/2和频率单位Δω=2π/N,进而确定边界频率整数参数值M=[ωc/Δω+0.5]
Step2.选择一长度为N的常用窗{f(n)},并与长度为N的反转矩形窗{RN(-n)}做卷积而得到一长度为2N-1的卷积窗wc(n),即
wc(n)=f(n)*RN(-n),n=-N+1,...,0,...,N-1
Step3.将参数M、N和卷积窗wc(n)直接代入如下解析式,获取滤波器系数g(n);
Step4.计算Lichtenberg比率频点处的幅频响应,即
Step5.将T=1-(Q-1)=2-Q代入如下改进的解析式中,算出最终的可消除过冲的滤波器系数gT(n),-N+1≤n≤N-1,即
2.本发明技术实现原理
2.1.全相位FIR滤波器的基本原理
2.1.1.全相位FIR滤波器的的衍生
全相位FIR滤波器设计法的核心思想来源于频率采样法,通过引入性能更好的卷积窗谱内插函数,设计出的FIR滤波器的通带波纹更小、阻带衰减更大。经典的频率采样法指定一个满足传统奇对称H(k)=H(N-k),k=0,...,N-1的频率采样向量H=[H(0),H(1),...,H(N-1)],可以设置为如下形式
直接对H进行IDFT得到
经典全相位滤波理论中的子滤波器是通过对h(n)进行循环移位得到的,基于此构造循环移位图,导出全相位滤波器系数。本发明中采用新的定义域延拓的方法来构造子滤波器。
我们把式(29)中n的IDFT定义域进行延拓,使得n∈[-m,-m+N-1],即
这样获得N个子滤波器hm=[h(-m),h(-m+1),...,h(-m+N-1)],m=0,...,N-1,其傅立叶变换为
如果将这N个子滤波器做算数平均,那么可以得到合成后的频响G(jω)为
以N=9,M=3为例,各子滤波器及其算数平均后的滤波器传输曲线如图1(a)~(j)所示。
从图1(a)~(i)可以看出,在9个子滤波器的幅频曲线上,分布着较大的正负波纹,但是,对他们做算数平均后,每个子滤波器上的波纹通过正负抵消,合成的幅频曲线具有较小的波纹,正如图1(j)所示。并且,从图中可以看出,每一个子滤波器的幅频曲线|H0(jω)|~|H8(jω)|在ω=k2π/N,k=0,...,N-1频点处精确通过频率采样向量H=[111000011]的设置点,因此其算数平均后的合成频响|G(jω)|也通过这些设置点。基于以上两方面的原因,合成后的滤波器频响|G(jω)|具有更好的传输性能。
为了进一步提高合成后的滤波器G(jω)的传输性能,可以使用一个长度为N的常见的归一化对称窗{f(m),m=0,...,N-1},将式(5)中对Hm(jω)的简单算术平均替换为加权平均,表示为
使用与图1相同的参数,指定式(33)中的加权窗{f(m)}为归一化的汉明窗,绘制传输曲线|G(jω)|如下图2所示,可与图3所示的经典频率采样法获得的传输曲线|H(jω)|对比:
对比图2和图3可看出,相比于传统的频率采样法传输曲线|H(jω)|,对N个子滤波器加权平均后的传输曲线|G(jω)|有更小的通带和阻带波纹,同时,|G(jω)|仍然准确通过设置的每一个频率采样点,其过渡带也没有被加宽。这种N个子滤波器加权平均的方法就是全相位FIR滤波器设计的核心思想,这种加权平均其实质就是对N个子滤波器在叠加过程中正负相互抵消,实现了幅度互补,从而保证了最终设计的滤波器传输曲线的通带波纹足够小和阻带衰减足够大,而这一过程的实现是内在的自然优化过程,不需要像优化算法的设计那样循环迭代,因此拥有很高的设计效率。从图2中可看出G(jω)在通带和阻带处的两个大约幅值为5.5%的凸起,本发明提出的Lichtenberg比率方法可消除该凸起。
2.1.2.全相位FIR滤波器的传输特性与卷积窗的关系
将式(31)代入式(33)中,并且交换m和n的求和次序,可得
为了简化上式,定义一个长度为2N-1的卷积窗{wc(n)},由长度为N的对称窗{f(n)}和长度为N的反转矩形窗{RN(-n)}构成如下
wc(n)=f(n)*RN(-n),n=-N+1,...,0,...,N-1(35)
上式可进一步表示为
因为{f(n)}和{RN(-n)}的非零元素都定义在区间[0,N-1]中,所以m满足
因此对式(36)分为两种情况,进一步推导可得
相应的,wc(-n)可表示为
将式(39)代入式(34),可得G(jω)为
由于f(n)是对称的,将f(m)=f(N-1-m)代入式(39)可得
对比式(41)和式(39)可看出,wc(n)也是对称的,因此式(40)可进一步推导为
与式(33)相比,全相位滤波器传输特性G(jω)得到进一步精简表示。将式(29)代入式(42)可得
其中,Wc(jω)定义为
对比式(43)和式(30)可以看出,全相位FIR滤波器的传输特性G(jω)同样可以表示为类似于频率采样法的内插形式,其实质上是通过新的卷积窗谱对频率采样向量H进行内插,最终表示为H(k)和归一化的长度为2N-1的卷积窗{wc(n)}的傅立叶谱Wc(jω)的离散卷积。
2.2.传输特性存在凸起的全相位滤波器系数解析表达式
2.2.1.等效三步设计法
式(42)可以看出,全相位FIR滤波器的系数长度为2N-1,其滤波器系数g(n)为
g(n)=h(n)wc(n),-N+1≤n≤N-1(45)
又因为exp[j2π/N·(n+N)]=exp[j2πn/N],所以根据式(30)可得出
h(n)=h(n+N)(46)
因此,定义域延拓的IDFT得到的2N-1长的滤波器系数向量h′可表示为
h′=[h(-N+1),...,h(-1),h(0),h(1),...,h(N-1)]
(47)
=[h(1),...,h(N-1),h(0),h(1),...,h(N-1)]
基于式(35)、式(43)和式(45),全相位FIR滤波器的设计可以总结为以下三步:
①对指定的频率采样向量H进行IDFT变换得到h=[h(0),h(1),...,h(N-1)],然后插入复制h(1)~h(N-1)为h左半部分得h′。
②选择一长度为N的常用窗f和长度为N的矩形窗b进行卷积并归一化,得到长度为2N-1的卷积窗wc
③将h′和wc点乘即得最终的全相位FIR滤波器系数g(n)。
显然以上三步设计法绕过了前面复杂的N个子滤波器加权平均过程,对FIR滤波器的设计来说,十分精简。
2.2.2.解析形式滤波器系数表达式的推导
为了进一步简化滤波器的设计步骤,本发明进一步推导全相位FIR滤波器设计的解析式,以便在给出设计参数时便能快速得到滤波器系数。
对H(k)进行定义域延拓的IDFT得到h(n)为
将式(28)代入式(48)得到
使用欧拉公式,上式可进一步推导为
因此,结合式(50)和式(43)可得最终滤波器系数g(n)为
上式中包含了n=0的情况,是通过将n=0直接代入式(48)和式(43)得到g(0)。
总之,通过将参数M(决定滤波器通带宽度)、N(决定滤波器长度)和卷积窗wc(n)(决定滤波器性能)直接代入系数解析式(51),便可以直接得到滤波器系数g(n)。通过上面的解析式直接进行滤波器设计,这一过程是十分高效的。
然而,对于单窗情况,解析式(51)仍有不足,就是没有考虑到边界频带的特性,由于通带和阻带附近的频率设置点为1,因而会导致过渡带两侧频带区间内传输曲线出现凸起,所以要对这一问题进行优化设计。
2.3.基于Lichtenberg比率的优化设计
2.3.1.全相位FIR滤波器的新性质
性质一:近似内插特性
对于任意的由整数部分和小数部分组成的频点ω,即ω=(p+β)Δω,p∈Z+,0≤β<1,其传输特性G(jω)可被近似看成由相邻的4个傅立叶传输函数Wc[j(ω-kΔω)],k=p-1,p,p+1,p+2组成,即
证明:对比图4和图5可发现,相比于矩形窗谱RN(jω),全相位FIR滤波器设计法引入的卷积窗谱Wc(jω)具有更窄、更小的旁瓣,更具体的说,Wc(jω)的旁瓣窄于Δω而RN(jω)的旁瓣大于3Δω。显然,频点ω=(p+β)Δω,p∈Z+,0≤β<1位于函数Wc[j(ω-pΔω)]的主瓣和函数Wc[j(ω-(p+1)Δω)]的主瓣上,同时位于函数Wc[j(ω-(p-1)Δω)]的第一右旁瓣和函数Wc[j(ω-(p+2)Δω)]的第一左旁瓣上,而对于其他的傅立叶传输函数Wc[j(ω-kΔω)],均远离于频点ω=(p+β)Δω,因此对传输特性G(jω)的影响可以忽略不计。因此,结合式(43),G(jω)可以近似表示为4个相邻内插函数的和。
近似内插特性对于进一步分析全相位FIR滤波器的频带组成有重要意义,尤其是在优化边界频带方面,提供了理论依据。
性质二:全通幅度互补特性
该性质在原有文献[14]中已经证明过,不属于新的性质,但为了下一个新性质证明的方便,本发明将该性质进行说明如下
对于两个长度为N的互补的频率采样向量H和Hc
他们对应的传输特性G(jω)和Gc(jω)满足以下的幅度互补关系:
G(jω)+Gc(jω)=1(54)
证明:由式(53)可知
H(k)+Hc(k)=1,k=0,...,N-1(55)
由于全相位FIR滤波器系数g(n)等于卷积窗wc(n)和h(n)定义域延拓的IDFT变换的乘积,如式(43)所示,因此这两个滤波器系数g(n)和gc(n)可表示为
对g(n)和gc(n)求和,得
代入式(55)得
显然,式(58)中的求和部分是对N个均匀分布在单位圆上的旋转因子求和,其和为一个Dirac函数,即
因此,式(58)可以进一步推导为
g(n)+gc(n)=wc(n)δ(n)=wc(0)δ(n)=δ(n)(60)
对上式两边进行傅立叶变换即可得G(jω)+Gc(jω)=1,证明了全相位FIR滤波器的全通幅度互补特性。
性质三:关于过渡带近似对称特性
对任意的频率对ω1=(M-1-λ)Δω,ω2=(M+λ)Δω,0≤λ<1,二者关于过渡带中心ωt=(M-0.5)Δω中心对称,如图6所示,全相位FIR滤波器传输曲线G(jω)位于频点ω1处的通带波纹与频点ω2处的阻带波纹近似相等。即实数的传输函数G(jω)满足
G(jω1)-1≈-G(jω2)(61)
证明:根据全通幅度互补特性,对于频点ω2=(M+λ)Δω,有
再根据近似内插特性,对于频点ω1=(M-1-λ)Δω,将p=M-2,β=1-λ,及H(M-3)=H(M-2)=H(M-1)=1,H(M)=0代入式(52)得
类似的,对于频点ω2=(M+λ)Δω,将p=M,β=λ,Hc(M-1)=0,Hc(M)=Hc(M+1)=Hc(M+2)=1代入式(52)得
因为Wc(jω)为偶函数,如图5所示,对比式(63)和式(64)可得
G(jω1)≈Gc(jω2)(65)
进一步将式(65)代入式(62)得
G(jω1)-1≈-G(jω2),(66)
对于位于频点(M-1)Δω和频点MΔω之间的过渡带中心频点ωt=(M-0.5)Δω,代入式(66)得
G[j(ωt-(0.5+λ)Δω)]-0.5≈-{G[j(ωt+(0.5+λ)Δω)]-0.5}(67)
上式表明,在ω∈[(M-2),(M+1)]Δω区域内,全相位低通FIR滤波器的传输曲线的通带部分是和阻带部分关于过渡带中心频点G(jωt)=0.5处近似奇对称的,其也准确通过6-dB截止频点,从如图6也可看出,式(67)也解释了G(jω)在通带和阻带处的两个大约幅值为5.5%的凸起是近似相等的。
以上三个特性是在低通滤波器的情况下证明的,显然也适用于高通、带通等滤波器,全相位FIR滤波器的这些特性,保证了其高性能和广泛应用的潜力,也为滤波器性能的进一步优化、特殊情况下的设计提供了理论依据。
2.3.2.加入过渡点优化的全相位滤波器系数解析式推导
为了抑制全相位FIR滤波器传输曲线过渡带两侧出现的小的凸起,有效的方法便是修改频率采样向量的过渡点设置值,在将通带和阻带凸起抑制的同时,保证过渡带最小的加宽量和引入最少的参数。
对比式(63)、式(64)和式(67),可以发现G(jω)关于过渡带中心的6dB频点ωt=(M-0.5)Δω对称的凸起部分,主要由过渡带左右相邻的两个边界频带采样点值H(M-1)=1和H(M)=0决定,并且其值满足
H(M-1)+H(M)=1(68)
因此过渡带左右的两个凸起可以通过修改这两个频点的值来有效抑制。利用奇对称属性,如果将k=M-1处的采样值1替换为一适当值T(T<1),则k=M处的采样值为1-T,更改后的频率采样向量HT具有以下的形式
明显上式只引入了一个参数T需要被决定,基于全相位滤波器的三个新性质,若使用适当的T值,便可有效抑制两个凸起,并且获得一个均衡的改进。
类似于式(51)中滤波器系数g(n)的推导过程,可以推导出由改进的频率采样向量HT决定的全相位低通FIR滤波器系数gT(n)的解析式如下
由于G(jω)曲线上的凸起部分非常小,汉明窗情况下只有5.5%的幅度,因此最优T值的选择在一个小的区间里,实际上,经过大量实验,其取值区间一般为T∈(0.9,1)。
通过上述加入过渡点优化的解析式设计,只要使用高效的方法来确定T值,同样可以快速的获得滤波器系数。对比于经典的优化方法和现代的进化优化方法,具有非常低的运算复杂度。因为最终的长度为2N-1的滤波器系数可以直接通过将参数N、M、T、wc(n)代入解析式(41)中而获得,其主要的运算复杂度在wc(n)的卷积操作上,并且可以通过使用FFT和IFFT的结合或者快速卷积(例如使用重叠相加或重叠保存法)来代替卷积的运算,因而使运算量进一步降低。因此,加入过渡点优化的全相位解析设计仍具有较高的效率,非常适用于高阶滤波器的设计,并且进一步提升了全相位FIR滤波器的性能。
2.3.3.最终的基于Lichtenberg比率的全相位滤波器解析式推导
一个适合的在(0.9,1)区间内的T值需要被确定,以便有效抑制过渡带左右的凸起。事实上,这两个凸起是由卷积窗wc(n)=f(n)*RN(-n)的第一旁瓣振荡引起的,选择不同的常用窗f(n)会引起不同程度的凸起,也就是需要确定不同的T值,来适用于不同的应用场景。若要确定全局最优的T值,那么使用优化算法是普遍的选择,但优化算法的引入一定会大大降低设计的效率,同时需要引入冗余参数,所以在一些需要高效滤波器的设计场合,引入优化算法不是一个好的解决方法。本发明通过大量实验,提出一种基于Lichtenberg比率的解析设计法,可以避免降低设计效率,同时选择较优的T值。
显然最优的T值跟过渡带两侧的凸起的幅度密切相关,通过在不同窗情况下构成全相位FIR滤波器的大量实验,如图7举例来说明(N=12,M=4时,选择常用的三角窗、汉宁窗、汉明窗和布莱克曼窗),发现滤波器的过渡带左右的凸起的最大值集中出现在一个特殊位置处,也就是说,Δω与|ωL-(M-2)Δω|(ωL与其相邻的左频率采样点的距离)的比例等于一个常用的Lichtenberg比率因此,频点ωL处的凸起波纹G(ωL)-1可以近似认为是过渡带两侧凸起的最大值,可以将其用来补偿ω=(M-1)Δω频点处的采样设置值。作为总结,通过以下的三步来完成基于Lichtenberg比率的设计法,获得较优的参数T以及改进的滤波器系数gT(n)。
①使用系数解析式(51)设计初始滤波器g(n),-N+1≤n≤N-1;
②计算g(n)在频点处的幅频响应Q,即
③将T=1-(Q-1)=2-Q代入改进的系数解析式(70)中,便得到最终改进的FIR滤波器gT(n),-N+1≤n≤N-1。
3.实验
实验一、加入过渡点优化设计后的滤波器性能改进效果
使用N=9,M=3的参数设计,使用长度为N的4种常用窗:三角窗、汉宁窗、汉明窗、布莱克曼窗来生成4个卷积窗wc(n),然后将这些参数代入式(51)得到四个原始的滤波器g(n),再使用上节中的基于Lichtenberg比率的解析设计法三步骤生成四个过渡点优化的滤波器gT(n),每一种窗设计下的原始滤波器和优化了的滤波器的幅频曲线和衰减曲线对比如图8所示,图中也列举了基于Lichtenberg比率的设计法中自动配置的参数T的值依次为0.9819、0.9554、0.9456、0.9091,表1进一步列举了所有这8个滤波器的性能指标(包括通带最大波纹、边界频带凸起、第1旁瓣衰减、0~1间的过渡带宽度)。
从图8(a)~(d)可看出,原始的滤波器g(n)的传输曲线具有在大部分通带良好的传输性能,除过在通带边沿有小的凸起(四种窗下分别达到2.6%、4.63%、5.47%、9.22%),而基于Lichtenberg比率的设计法优化后的滤波器gT(n)的传输曲线有效的抑制了这些波纹,其过渡带的加宽量特别小。从表1可看出,这4种窗下过渡点优化的设计出的滤波器频响曲线有不同的性能:三角窗设计的滤波器拥有最窄的过渡带带宽(1.16Δω),而其他窗的过渡带带宽分别为汉宁窗1.35Δω、汉明窗1.42Δω,布莱克曼窗1.73Δω,相应的三角窗下滤波器具有通带最大波纹1.55%;汉宁窗和汉明窗设计的滤波器在压缩通带波纹(小于1%)和控制过渡带宽度上有最好的均衡;布莱克曼窗设计的滤波器具有最小的通带最大波纹(0.13%),对比于其它窗的三角窗1.55%、汉宁窗0.76%、汉明窗0.44%,同时其具有最宽的过渡带(1.73Δω)。
表1不同窗设计下的滤波器性能对比
图8(a)~(d)也进一步证明了不管选择哪个窗、不管使用哪个过渡点的值T,全相位FIR滤波器的传输曲线的6-dB截止频点基本位于ωt=(M-0.5)2π/N=2.5Δω不变,这也证明了全相位FIR滤波器的控制6-dB截止频点位置上具有高的鲁棒性。
实验二、与经典设计法的滤波器性能对比
使用本节提出基于Lichtenberg比率的设计法与经典的eigenfilter方法(在WLS算法准则下的设计)[3]、Remez算法和传统差分演进算法[15]来设计阶数为30(滤波器长度L为31)、通带截止频率为ωp=0.47π、过渡带宽度为0.17π的FIR滤波器。为满足这些要求,全相位FIR滤波器的参数设置为N=(L+1)/2=16,wc(n)由长度N的汉明窗和长度N的矩形窗卷积而成;Eigenfilter的的加权因子α设置为0.5;Remez算法通过调用Matlab函数‘firpm.m’生成;DE算法的参数设置为:种群大小200,变异因子0.4,交叉概率0.8。
图9(a)、(b)分别展示了这四个滤波器的频响曲线和衰减曲线,包括局部放大曲线。表2列举了各个滤波器的性能指标,包括边界频带凸起、第1旁瓣衰减、耗费时间(使用AMD奔腾II2.0GHzCPU,4GBRAM)、0~1间的过渡带带宽和RMSE,其中RMSE表示为滤波器通带波纹的最小均方误差值,定义为
RMSE={E[|G(jω)|-1]2}1/2,0<ω≤ωp(72)
表2不同方法设计的滤波器性能对比
从表2可看出,在相同的过渡带0.1712π设计下,对比于其它三种(eigenfilter方法、Remez算法、DE算法)经典设计法,本节提出的基于Lichtenberg比率的设计法设计出的FIR滤波器具有与其相当的性能。其中,eigenfilter具有最小的RMSE值(0.0020),其次是Remez滤波器(0.0027),而本节方法设计的滤波器RMSE值为0.0028。由于等波纹的设计,Remez算法的优秀性能主要表现在对滤波器边界波纹的控制上,如图9所示,其设计的滤波器边界凸起最小,为0.38%,而本文设计法设计的滤波器边界波纹为0.69%。但是在不考虑边界凸起的通带范围内,本节方法设计的滤波器具有等同于eigenfilter方法设计的非常小的通带波纹,小于Remez滤波器和DE算法设计的滤波器。在耗费时间上面,本节基于Lichtenberg比率的设计法最小只有0.004s,而eigenfilter为0.011s、Remez为0.242s、DE算法为98.792s,这是基于本节方法的滤波器系数是通过两个解析式直接求取,不同与Remez的大量迭代、eigenfilter的矩阵分解以及DE的种群进化,解析的设计拥有更高的设计效率。
实验三、滤波器设计效率分析
增加FIR滤波器的阶数至198(滤波器长度为L=199,N=100),指定通带截止频率ωp=0.5π,过渡带宽度为0.02π(M=26),wc(n)由长度N的汉明窗和长度N的矩形窗卷积而成。在这种情况下,DE算法因为所需资源(存储器和处理单元)太大而难以进行设计,因此,使用eigenfilter,Remez算法和基于Lichtenberg比率的设计法三种方法来进行对比分析,如图10(a)、(b)分别展示三种方法设计出的滤波器的频响曲线和衰减曲线(包括局部放大后的曲线),相应地表3列举了这三个滤波器的性能参数值。
表3大阶数下不同方法设计的滤波器性能对比
从图10和表3可以看出,当滤波器阶数较大时,在通带整体波纹上,eigenfilter的最小,RMSE值为0.0009,本节提出的基于Lichtenberg比率的设计法次之,RMSE值为0.0013,Remez算法的最高为0.0026。相反的是,Remez滤波器具有最大的第1旁瓣衰减(-48.66dB),本节提出的基于Lichtenberg比率的设计法次之(-44.12dB),eigenfilter最小(-43.19dB)。此外,从放大的曲线可以看出,三种方法设计的滤波器在不同的通带区域具有不同的波纹:在通带边沿,本节设计法的凸起较大为0.62%,仍然不如Remez滤波器的0.37%,但要优于eigenfilter的0.69%;在通带内部,eigenfilter具有最小的波纹,本节设计法次之,而Remez滤波器最大。
再从表3可以得出,在198阶的设计下,本节提出的基于Lichtenberg比率的设计法具有最小的执行时间,为0.008s,仅为eigenfilter(0.025s)的1/3,Remez滤波器(0.482s)的1/60,换句话说,本节提出的方法在高阶数FIR滤波器的设计上比这两种方法具有更高的效率,为了证明这个优势,我们在滤波器阶数从8到400时将三种方法耗费的时间进行分析对比(固定滤波器通带截止频率ωp=0.5π),其结果如图11所示。
图11表示随着滤波器阶数的增大,eigenfilter方法和Remez算法的滤波器设计时间也大幅度增加,而本节提出的基于Lichtenberg比率的方法设计滤波器所耗费的时间基本不变。具体说来,eigenfilter方法的增长速度(0.01s~0.08s)小于Remez算法(0.01s~0.8s),而本节方法的执行时间基本维持在10-2s不变,图11中也可看到其去曲线基本为直线。并且,随着滤波器阶数的增大,频率分辨率Δω=2π/N也相应的变小,这会增加本节方法对滤波器边界频带的控制。因此,本节基于Lichtenberg比率的解析设计法是非常适合设计高阶FIR滤波器的。
通过以上仿真实验证明,可以得出,使用基于Lichtenberg比率的设计法,T值的确定是一个内部的确定过程,并没有额外增加需要的设计参数,只是在初始的解析设计基础上增加了内部计算过程和基于改进解析式的系数计算,其设计复杂度的增加是很少的。因此相比于经典的FIR滤波器设计法,在滤波器性能有所提高的情况下,基于Lichtenberg比率的全相位设计法还具有更高的设计效率,尤其体现在高阶滤波器的设计上,非常适用于更高性能、高效率的设计场合上。
下面对实施本发明的硬件予以简单说明。
在图12中,首先将所需的滤波器6dB边界频率ωc、奇数长度L及卷积窗wc(n)存入外部RAM中,再将它们实时输入到DSP中,经过DSP内部核心算法,得到全相位滤波器的系数及其传输曲线;根据全相位滤波器传输曲线计算出Lichtenberg比率频点ωL处的幅频响应Q,进而得到优化参数T,返还并存储在外部RAM,由外部RAM将滤波器设计所需的全部参数再次输入DSP,得到所要求的滤波器系数及其传输曲线,由输出驱动及显示电路将其实时显示出来。
其中,图12的DSP(DigitalSignalProcessor,数字信号处理器)为核心器件,在计算滤波器系数的过程中,完成如下主要功能:
(1)调用内部核心算法,对本发明提出的解析式进行构建,计算出所需的滤波器系数,对滤波器系数进行傅里叶变换,得到滤波器传输曲线;
(2)控制滤波器参数输入时间,并根据需要实时调整所需要的参数值;
(3)将滤波器设计结果实时输出至驱动和显示模块。
DSP器件的内部程序流程如图13所示。
本发明将所提出的“基于Lichtenberg比率的数字滤波器解析设计法”的核心算法植入DSP器件内,基于此完成高精度、低复杂度、高效的滤波器设计。
本发明流程分为如下几个步骤:
(1)首先根据具体需要的滤波器截止频率和过渡带带宽计算滤波器设计所需的参数N、m、Δω,并根据Lichtenberg比率计算出幅频响应Q及优化参数T;
(2)然后,CPU主控器从I/O端口读取滤波器参数,进入内部RAM;
(3)根据推导出的优化解析式(70)进行滤波器设计是DSP算法最核心的部分,运行该算法后,即可得到目标滤波器系数及其滤波器传输曲线;
(4)判断本发明方法是否满足实际需求,若不满足,程序返回,重新根据要求设定滤波器参数;
(5)直至设计结果符合实际要求,然后通过DSP的输出总线输出至外部显示驱动设备,将滤波器设计结果进行数码显示。
需指出,由于采用了DSP实现,使得整个滤波器设计变得更为灵活快捷,可根据滤波器设计过程中的实际需要,灵活变换滤波器参数,使之最终符合工程需要。

Claims (6)

1.一种数字滤波器解析设计法及其滤波器,其特征是,指定一个满足传统奇对称H(k)=H(N-k),k=0,...,N-1的频率采样向量H=[H(0),H(1),...,H(N-1)],设置为如下形式:
直接对H进行IDFT得到
h ( n ) = 1 N &Sigma; k = 0 N - 1 H ( k ) e j 2 &pi; n k / N , n = 0 , ... , N - 1 - - - ( 2 )
采用如下定义域延拓的方法来构造子滤波器:
把式(2)中n的IDFT定义域进行延拓,使得n∈[-m,-m+N-1],即
h ( n ) = 1 N &Sigma; k = 0 N - 1 H ( k ) e j 2 &pi; n k / N , n = - m , ... , - m + N - 1 - - - ( 3 )
这样获得N个子滤波器hm=[h(-m),h(-m+1),...,h(-m+N-1)],m=0,...,N-1,其傅立叶变换为
H m ( j &omega; ) = &Sigma; n = - m - m + N - 1 h ( n ) e - j n &omega; , m = 0 , ... , N - 1 - - - ( 4 )
如果将这N个子滤波器做算数平均,那么得到合成后的频响G(jω)为
G ( j &omega; ) = 1 N &Sigma; m = 0 N - 1 H m ( j &omega; ) , m = 0 , ... , N - 1 - - - ( 5 )
使用一个长度为N的常见的归一化对称窗{f(m),m=0,...,N-1},将式(5)中对Hm(jω)的简单算术平均替换为加权平均,表示为
G ( j &omega; ) = &Sigma; m = 0 N - 1 f ( m ) H m ( j &omega; ) , m = 0 , ... , N - 1 - - - ( 6 )
将式(4)代入式(6)中,并且交换m和n的求和次序,得到
G ( j &omega; ) = &Sigma; m = 0 N - 1 f ( m ) &Sigma; n = - m - m + N - 1 h ( n ) e - j n &omega; = &Sigma; n = - N + 1 - 1 h ( n ) e - j n &omega; &Sigma; m = - n N - 1 f ( m ) + &Sigma; n = 0 N - 1 h ( n ) e - j n &omega; &Sigma; m = 0 - n + N - 1 f ( m ) - - - ( 7 )
为了简化上式,定义一个长度为2N-1的卷积窗{wc(n)},由长度为N的对称窗{f(n)}和长度为N的反转矩形窗{RN(-n)}构成如下
wc(n)=f(n)*RN(-n),n=-N+1,...,0,...,N-1(8)
上式进一步表示为
w c ( n ) = &Sigma; m f ( m ) &CenterDot; R N ( m - n ) - - - ( 9 )
因为{f(n)}和{RN(-n)}的非零元素都定义在区间[0,N-1]中,所以m满足
0 &le; m &le; N - 1 0 &le; m - n &le; N - 1 &DoubleRightArrow; 0 &le; m &le; N - 1 n &le; m &le; n + N - 1 - - - ( 10 )
因此对式(9)分为两种情况,进一步推导可得
w c ( n ) = &Sigma; m = 0 n + N - 1 f ( m ) - N + 1 &le; n < 0 &Sigma; m = n N - 1 f ( m ) 0 &le; n &le; N - 1 - - - ( 11 )
相应的,wc(-n)表示为
w c ( - n ) = &Sigma; m = - n N - 1 f ( m ) , - N + 1 &le; n &le; 0 &Sigma; m = 0 - n + N - 1 f ( m ) , 0 < n &le; N - 1 - - - ( 12 )
将式(12)代入式(7),得G(jω)为
G ( j &omega; ) = &Sigma; n = - N + 1 N - 1 h ( n ) w c ( - n ) e - j n &omega; - - - ( 13 )
由于f(n)是对称的,将f(m)=f(N-1-m)代入式(12)得
对比式(14)和式(12)可看出,wc(n)也是对称的,因此式(13)进一步推导为
G ( j &omega; ) = &Sigma; n = - N + 1 N - 1 h ( n ) w c ( n ) e - j n &omega; - - - ( 15 )
与式(6)相比,全相位滤波器传输特性G(jω)得到进一步精简表示,将式(2)代入式(15)得
G ( j &omega; ) = 1 N &Sigma; n = - N + 1 N - 1 &Sigma; k = 0 N - 1 H ( k ) e j 2 &pi; n k / N w c ( n ) e - j n &omega; = &Sigma; k = 0 N - 1 H ( k ) 1 N &Sigma; n = - N + 1 N - 1 w c ( n ) e - j n ( &omega; - 2 &pi; k / N ) = &Sigma; k = 0 N - 1 H ( k ) W c &lsqb; j ( &omega; - 2 &pi; k / N ) &rsqb; - - - ( 16 )
其中,Wc(jω)定义为
W c ( j &omega; ) = &Sigma; n = - N + 1 N - 1 w c ( n ) e - j n &omega; - - - ( 17 )
2.如权利要求1所述的数字滤波器解析设计法及其滤波器,其特征是,由式(15),全相位FIR滤波器的系数长度为2N-1,其滤波器系数g(n)为:
g(n)=h(n)wc(n),-N+1≤n≤N-1(18)
又因为exp[j2π/N·(n+N)]=exp[j2πn/N],所以根据式(3)得出
h(n)=h(n+N)(19)
因此,定义域延拓的IDFT得到的2N-1长的滤波器系数向量h′表示为
h &prime; = &lsqb; h ( - N + 1 ) , ... , h ( - 1 ) , h ( 0 ) , h ( 1 ) , ... , h ( N - 1 ) &rsqb; = &lsqb; h ( 1 ) , ... , h ( N - 1 ) , h ( 0 ) , h ( 1 ) , ... , h ( N - 1 ) &rsqb; - - - ( 20 )
基于式(8)、式(16)和式(18),全相位FIR滤波器的设计可以总结为以下三步:
①对指定的频率采样向量H进行离散傅里叶逆变换IDFT变换得到h=[h(0),h(1),...,h(N-1)],然后插入复制h(1)~h(N-1)为h左半部分得h′;
②选择一长度为N的常用窗f和长度为N的矩形窗b进行卷积并归一化,得到长度为2N-1的卷积窗wc
③将h′和wc点乘即得最终的全相位FIR滤波器系数g(n)。
3.如权利要求2所述的数字滤波器解析设计法及其滤波器,其特征是,对H(k)进行定义域延拓的IDFT得到h(n)为:
h ( n ) = 1 N &Sigma; k = 0 N - 1 H ( k ) e j 2 &pi; n k / N , - N + 1 &le; n &le; N - 1 - - - ( 21 )
将式(1)代入式(21)得到
h ( n ) = 1 N &lsqb; &Sigma; k = 0 M - 1 e j 2 &pi; n k / N + &Sigma; k = N - M + 1 N - 1 e j 2 &pi; n k / N &rsqb; = 1 N &lsqb; 1 - e j 2 &pi; n M / N 1 - e j 2 &pi; n / N + e j 2 &pi; n ( N - M + 1 ) / N ( 1 - e j 2 &pi; n ( M - 1 ) / N ) 1 - e j 2 &pi; n / N &rsqb; = 1 N &lsqb; e j &pi; n ( 2 M - 1 ) / N - e - j &pi; n ( 2 M - 1 ) / N e j &pi; n / N - e - j &pi; n / N &rsqb; , n &Element; &lsqb; - N + 1 , - 1 &rsqb; &cup; &lsqb; 1 , N - 1 &rsqb; - - - ( 22 )
使用欧拉公式,上式进一步推导为
h ( n ) = s i n &lsqb; n ( 2 M - 1 ) &pi; / N &rsqb; N sin ( &pi; n / N ) , n &Element; &lsqb; - N + 1 , - 1 &rsqb; &cup; &lsqb; 1 , N - 1 &rsqb; - - - ( 23 )
因此,结合式(23)和式(16)可得最终滤波器系数g(n)为
g ( n ) = w c ( n ) N s i n &lsqb; n ( 2 M - 1 ) &pi; / N &rsqb; s i n ( &pi; n / N ) , n &Element; &lsqb; - N + 1 , - 1 &rsqb; &cup; &lsqb; 1 , N - 1 &rsqb; 2 M - 1 N , n = 0 - - - ( 24 )
上式中包含了n=0的情况,是通过将n=0直接代入式(21)和式(16)得到g(0)。
4.如权利要求3所述的数字滤波器解析设计法及其滤波器,其特征是,如果将k=M-1处的采样值1替换为一适当值T(T<1),则k=M处的采样值为1-T,更改后的频率采样向量HT具有以下的形式:
上式只引入了一个参数T需要被决定;
类似于式(24)中滤波器系数g(n)的推导过程,推导出由改进的频率采样向量HT决定的全相位低通FIR滤波器系数gT(n)的解析式如下
g T ( n ) = w c ( n ) N sin &lsqb; &pi; n ( 2 M - 3 ) / N &rsqb; sin ( &pi; n / N ) + 2 T cos &lsqb; 2 &pi; n ( M - 1 ) / N &rsqb; + 2 ( 1 - T ) cos ( 2 &pi; n M / N ) , n &Element; &lsqb; - N + 1 , - 1 &rsqb; &cup; &lsqb; 1 , N - 1 &rsqb; ( 2 M - 1 ) / N , n = 0 - - - ( 43 )
由于G(jω)曲线上的凸起部分非常小,汉明窗情况下只有5.5%的幅度,因此最优T值的取值区间为T∈(0.9,1);通过上述加入过渡点优化的解析式设计,只要使用高效的方法来确定T值,同样可以快速的获得滤波器系数。
5.如权利要求3所述的数字滤波器解析设计法及其滤波器,其特征是,频点ωL处的凸起波纹G(ωL)-1近似认为是过渡带两侧凸起的最大值,将其用来补偿ω=(M-1)△ω频点处的采样设置值,并通过以下的三步来完成基于Lichtenberg比率的设计法,获得较优的参数T以及改进的滤波器系数gT(n):
①使用系数解析式(24)设计初始滤波器g(n),-N+1≤n≤N-1;
②计算g(n)在频点处的幅频响应Q,即
Q = &Sigma; n = - N + 1 N - 1 g ( n ) e - jn&omega; L - - - ( 44 )
③将T=1-(Q-1)=2-Q代入改进的系数解析式(43)中,便得到最终改进的FIR滤波器gT(n),-N+1≤n≤N-1。
6.如权利要求所述的数字滤波器解析设计法及其滤波器,其特征是,设计法及其滤波器进一步优化为:
Step1.给定期望的滤波器的6dB边界频率ωc和奇数长度L,算出对应的全相位滤波器阶数N=(L+1)/2和频率单位△ω=2π/N,进而确定边界频率整数参数值M=[ωc/△ω+0.5]
Step2.选择一长度为N的常用窗{f(n)},并与长度为N的反转矩形窗{RN(-n)}做卷积而得到一长度为2N-1的卷积窗wc(n),即
wc(n)=f(n)*RN(-n),n=-N+1,...,0,...,N-1
Step3.将参数M、N和卷积窗wc(n)直接代入如下解析式,获取滤波器系数g(n);
g ( n ) = w c ( n ) N s i n &lsqb; n ( 2 M - 1 ) &pi; / N &rsqb; sin ( &pi; n / N ) , n &Element; &lsqb; - N + 1 , - 1 &rsqb; &cup; &lsqb; 1 , N - 1 &rsqb; 2 M - 1 N , n = 0
Step4.计算Lichtenberg比率频点处的幅频响应,即
Q = &Sigma; n = - N + 1 N - 1 g ( n ) e - jn&omega; L
Step5.将T=1-(Q-1)=2-Q代入如下改进的解析式中,算出最终的可消除过冲的滤波器系数gT(n),-N+1≤n≤N-1,即
g T ( n ) = w c ( n ) N sin &lsqb; &pi; n ( 2 M - 3 ) / N &rsqb; sin ( &pi; n / N ) + 2 T cos &lsqb; 2 &pi; n ( M - 1 ) / N &rsqb; + 2 ( 1 - T ) cos ( 2 &pi; n M / N ) , n &Element; &lsqb; - N + 1 , - 1 &rsqb; &cup; &lsqb; 1 , N - 1 &rsqb; ( 2 M - 1 ) / N , n = 0 .
CN201610086020.4A 2016-02-16 2016-02-16 数字滤波器解析设计法及其滤波器 Pending CN105680825A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610086020.4A CN105680825A (zh) 2016-02-16 2016-02-16 数字滤波器解析设计法及其滤波器

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610086020.4A CN105680825A (zh) 2016-02-16 2016-02-16 数字滤波器解析设计法及其滤波器

Publications (1)

Publication Number Publication Date
CN105680825A true CN105680825A (zh) 2016-06-15

Family

ID=56304441

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610086020.4A Pending CN105680825A (zh) 2016-02-16 2016-02-16 数字滤波器解析设计法及其滤波器

Country Status (1)

Country Link
CN (1) CN105680825A (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106452389A (zh) * 2016-09-30 2017-02-22 中国人民解放军63908部队 基于指数再生窗时域调制滤波器设计方法及滤波器
CN106849909A (zh) * 2017-01-04 2017-06-13 天津大学 一种基于两层次优化的fir滤波器设计方法及其装置
CN108205701A (zh) * 2016-12-20 2018-06-26 联发科技股份有限公司 一种执行卷积计算的系统及方法
CN108390698A (zh) * 2018-03-16 2018-08-10 贵州电网有限责任公司 一种基于插值fft算法的电力线载波参数测量方法
CN109921763A (zh) * 2019-02-26 2019-06-21 华南理工大学 一种用于减少乘法器的fir滤波器及其输出计算方法
CN113676156A (zh) * 2021-08-09 2021-11-19 成都玖锦科技有限公司 一种基于lms的任意幅频响应fir滤波器设计方法
CN113766239A (zh) * 2020-06-05 2021-12-07 于江鸿 数据处理的方法和系统

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106452389A (zh) * 2016-09-30 2017-02-22 中国人民解放军63908部队 基于指数再生窗时域调制滤波器设计方法及滤波器
CN106452389B (zh) * 2016-09-30 2018-10-26 中国人民解放军63908部队 基于指数再生窗时域调制滤波器设计方法及滤波器
CN108205701A (zh) * 2016-12-20 2018-06-26 联发科技股份有限公司 一种执行卷积计算的系统及方法
CN108205701B (zh) * 2016-12-20 2021-12-28 联发科技股份有限公司 一种执行卷积计算的系统及方法
CN106849909A (zh) * 2017-01-04 2017-06-13 天津大学 一种基于两层次优化的fir滤波器设计方法及其装置
CN108390698A (zh) * 2018-03-16 2018-08-10 贵州电网有限责任公司 一种基于插值fft算法的电力线载波参数测量方法
CN108390698B (zh) * 2018-03-16 2021-08-10 贵州电网有限责任公司 一种基于插值fft算法的电力线载波参数测量方法
CN109921763A (zh) * 2019-02-26 2019-06-21 华南理工大学 一种用于减少乘法器的fir滤波器及其输出计算方法
CN113766239A (zh) * 2020-06-05 2021-12-07 于江鸿 数据处理的方法和系统
CN113676156A (zh) * 2021-08-09 2021-11-19 成都玖锦科技有限公司 一种基于lms的任意幅频响应fir滤波器设计方法
CN113676156B (zh) * 2021-08-09 2024-01-26 成都玖锦科技有限公司 一种基于lms的任意幅频响应fir滤波器设计方法

Similar Documents

Publication Publication Date Title
CN105680825A (zh) 数字滤波器解析设计法及其滤波器
CN107294511B (zh) 一种低复杂度的可变分数时延滤波方法及滤波器
CN104283527B (zh) 一种边界频带可快速配置的高效滤波器方法及其装置
CN103888104B (zh) Fir数字滤波器设计方法和系统
CN103199822B (zh) 一种带宽可变低通数字滤波器的设计方法
CN111222088B (zh) 一种改进的平顶自卷积窗加权电力谐波幅值估计方法
CN103197300A (zh) 一种基于gpu的外辐射源雷达直达波杂波对消实时处理方法
Goel et al. Design of FIR filter using FCSD representation
Kisler et al. Plate silencers for broadband low frequency sound attenuation
Saini et al. Modification of chebyshev pseudospectral method to minimize the gibbs oscillatory behaviour in resynthesizing process
Budunova et al. On a new method for approximation of squares of atomic functions ha (x) by nonnegative rational fractions
Ferdi Impulse invariance-based method for the computation of fractional integral of order 0< α< 1
CN104240230A (zh) 一种提高相位相关算法匹配精度的方法
CN106849909A (zh) 一种基于两层次优化的fir滤波器设计方法及其装置
Chen et al. Pyramid-shaped grid for elastic wave propagation
CN112332809B (zh) 一种幅度非衰减均衡相位的分数阶椭圆滤波器设计方法
CN105677957B (zh) 一种近似精确重构余弦调制滤波器组的设计方法与装置
Khoo et al. General formulation of shift and delta operator based 2-D VLSI filter structures without global broadcast and incorporation of the symmetry
Gronczyński Recursive Fourier transform algorithms with integrated windowing
CN106338743B (zh) 基于ccd解算的空频联合自适应调零算法
CN107947760B (zh) 一种稀疏fir陷波器的设计方法
Babic et al. Optimum low-order windows for discrete Fourier transform systems
CN104300938B (zh) 均衡由耦合损耗引起的失真的方法及利用其生产的滤波器
Meyer-Baese et al. Infinite impulse response (IIR) digital filters
Gorinevsky et al. Optimization-based design and implementation of multi-dimensional zero-phase IIR filters

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20160615