基于改進(jìn)差分進(jìn)化算法的Ⅱ型頻率采樣法FIR數(shù)字濾波器設(shè)計(jì)

胡仕兵, 聶 喜, 陳子為

(成都信息工程大學(xué)電子工程學(xué)院,四川 成都 610225)

0 引言

頻率采樣法是設(shè)計(jì)FIR數(shù)字濾波器的常用方法之一[1-2],分為Ⅰ型和Ⅱ型兩種,對應(yīng)的濾波器頻率響應(yīng)Hd(ejω)的 N 點(diǎn)采樣值分別為和。對Ⅰ型頻率采樣法應(yīng)用和研究較廣泛,而對Ⅱ型頻率采樣法很少述及。實(shí)際上在對FIR濾波器頻率響應(yīng)擬合程度要求較高的場合,可以將Ⅰ型頻率采樣法和Ⅱ型頻率采樣法交替使用,以獲得更精準(zhǔn)的邊界頻率,兩種采樣法的邊界頻率擬合誤差曲線如圖1所示。

圖1 Ⅰ型和Ⅱ型頻率采樣法FIR低通濾波器邊界頻率擬合誤差曲線(過渡帶寬0.15 π、過渡帶內(nèi)1個(gè)樣點(diǎn))

在確定FIR數(shù)字濾波器的過渡帶樣點(diǎn)值時(shí),采用了查表(look-up table,LUT)法[3]、遺傳算法(genetic algorithm,GA)[4-7]、模擬退火算法(simulated annealing algorithm,SAA)[8]、蟻群優(yōu)化算法(ant colony optimization algorithm,ACOA)[9]、粒子群優(yōu)化算法(particle swarm optimization algorithm,PSOA)[5,10-12]、免疫算法(immune algorithm,IA)[13]、克隆選擇算法(clonal selection algorithm,CSA)[14-15]、自由搜索算法(free search algorithm,FSA)[16]、水循環(huán)算法(water cycle algorithm,WCA)[17]、人工魚群算法(artificial fish swarm algorithm,AFSA)[18]等,這些方法在一定程度上存在求得最優(yōu)解的偏差過大、算法群體多樣性較小、算法控制變量多、結(jié)構(gòu)較復(fù)雜、尋優(yōu)精度不高、收斂速度較慢、易于出現(xiàn)早熟收斂(易陷入局部最優(yōu)解)和程序運(yùn)行時(shí)間較長等缺陷。差分進(jìn)化(differential evolution,DE)算法是由Rainer Storn等[19-21]于20世紀(jì)末提出的一種基于模擬自然界生物的進(jìn)化過程、求解連續(xù)空間上優(yōu)化問題的高效、強(qiáng)大的隨機(jī)種群搜索算法,由于其控制參數(shù)少、計(jì)算簡單、易于實(shí)現(xiàn)等特點(diǎn),引起廣泛關(guān)注和研究。但是對于特定的優(yōu)化問題,DE算法也表現(xiàn)出對控制參數(shù)的設(shè)置敏感、收斂速度慢、種群搜索停滯和早熟收斂等問題[21]。本文對DE算法進(jìn)行改進(jìn),增加縮放因子自適應(yīng)生成策略、交叉概率拋物線式動態(tài)產(chǎn)生策略和基因邊界檢查和處理方法,提出一種新的改進(jìn)差分進(jìn)化(modified differential evolution,MDE)算法并且應(yīng)用于采用Ⅱ型頻率采樣法設(shè)計(jì)FIR數(shù)字濾波器時(shí)過渡帶樣點(diǎn)值的優(yōu)化中,取得了良好的實(shí)驗(yàn)效果。

1 MDE算法

對于待優(yōu)化的問題函數(shù)f,尋找參數(shù)x1,x2,…,xD使f達(dá)到最小值:

其中:D是解空間的維數(shù),Uj和Lj分別是第j個(gè)分量xj取值范圍的上界值和下界值。MDE算法包括4個(gè)基本操作:初始化種群、變異、交叉和選擇。在初始化種群后,算法循環(huán)執(zhí)行變異、交叉和選擇操作,直到計(jì)算結(jié)果在預(yù)期的誤差精度內(nèi)或達(dá)到最大進(jìn)化代數(shù);最后獲得最優(yōu)解。

1.1 種群初始化

設(shè)進(jìn)化代數(shù)為 G∈{0,1,…,Gmax}時(shí)的種群為P(G)={X1(G),X2(G),…,XM(G)},其中:M為種群規(guī)模大小,Xi(G)={xi,1(G),xi,2(G),…,xi,D(G)}(i=1,2,…,M)表示第G代中的第i個(gè)個(gè)體(染色體),xi,j(G)表示第 G代中第 i個(gè)個(gè)體的第 j個(gè)分量(基因)。每個(gè)個(gè)體代表解空間內(nèi)的某一個(gè)解,而基因代表解的各個(gè)分量。種群初始化就是在解空間范圍內(nèi)隨機(jī)均勻地產(chǎn)生M個(gè)個(gè)體,每個(gè)個(gè)體由D維向量組成,取值方式如下:

其中,rand(0,1)表示在(0,1)區(qū)間均勻分布的隨機(jī)數(shù)。

1.2 變異運(yùn)算

傳統(tǒng)差分進(jìn)化算法通過差分策略實(shí)現(xiàn)個(gè)體變異。對于每個(gè)個(gè)體Xi(G),從種群中隨機(jī)選擇3個(gè)不同的個(gè)體 Xp1(G)、Xp2(G)和 Xp3(G)(i≠p1≠p2≠p3)生成變異向量Hi(G+1):

其中:F∈[0,2]為縮放因子(常數(shù)),用于控制差分向量Xp2(G)-Xp3(G)的影響力。

F值較大時(shí)能使算法在全局范圍內(nèi)搜索到有效解,但會降低收斂速度;F值較小時(shí)會降低種群的多樣性,導(dǎo)致陷入局部最優(yōu)解,收斂性能差。本文提出改進(jìn)方法:將Xp1(G)、Xp2(G)、Xp3(G)按照適應(yīng)度值優(yōu)劣進(jìn)行排序得到個(gè)體向量Xbest(G)、Xmiddle(G)和Xworst(G),且有為計(jì)算個(gè)體適應(yīng)度值函數(shù),采用以下的自適應(yīng)調(diào)整策略生成變異向量Hi(G+1):Fi的取值根據(jù)生成差分向量的兩個(gè)個(gè)體適應(yīng)度值作自適應(yīng)變化。式(5)中,FL和FU分別是縮放因子的下限和上限,取FL=0.1,FU=0.9。自適應(yīng)變異可以加快算法的收斂速度,避免有時(shí)陷入局部最優(yōu)點(diǎn)的缺陷。

1.3 交叉運(yùn)算

引入交叉運(yùn)算可以增強(qiáng)種群的多樣性。對于每個(gè)目標(biāo)個(gè)體Xi(G)及其生成的變異個(gè)體Hi(G+1)進(jìn)行個(gè)體間的交叉操作生成實(shí)驗(yàn)向量Vi(G+1),即按照一定的概率CR選擇Hi(G+1)或Xi(G)的等位基因作為Vi(G+1)的等位基因。為保證變異中間體{Hi(G+1)}的每個(gè)染色體至少有一個(gè)基因遺傳給下一代,需要隨機(jī)取出Hi(G+1)中的第jrand位基因作為Vi(G+1)的第jrand位等位基因。Vi(G+1)的基因計(jì)算方式如下:

其中:CR∈[0,1]為交叉概率(常數(shù)),jrand為{1,2,…,D}間的隨機(jī)整數(shù)。

常數(shù)CR值較大時(shí),試驗(yàn)個(gè)體中的信息大多來自于變異個(gè)體,使種群的多樣性越來越好,但容易忽略掉全局最優(yōu)解,算法容易出現(xiàn)早熟收斂現(xiàn)象;CR值較小時(shí),試驗(yàn)個(gè)體中的信息大多來自目標(biāo)個(gè)體,會使種群性降低,導(dǎo)致種群進(jìn)化出現(xiàn)停滯現(xiàn)象。故提出如下拋物線式動態(tài)CR策略:

其中,CR(G)為第G代的交叉概率,CRU和CRL分別是交叉概率的上限值和下限值,一般取CRU=1.0、CRL=0.9。這樣CR在0.9處緩慢增長,可以減小固定參數(shù)帶來的不足。

1.4 選擇運(yùn)算

采用貪婪算法從實(shí)驗(yàn)向量Vi(G+1)和原向量Xi(G)中選擇出適應(yīng)度值更高的個(gè)體作為下一代種群的個(gè)體,計(jì)算方式為

選擇操作保證了個(gè)體Xi(G+1)一定優(yōu)于或持平于Xi(G),使算法最終一定會收斂到某個(gè)最優(yōu)點(diǎn)(可能是局部最優(yōu))上;而變異、交叉操作有助于解跳出局部最優(yōu)達(dá)到全局最優(yōu)。

1.5 邊界檢查和處理

為確保在進(jìn)化過程中解的有效性,必須判斷染色體中各基因是否超出邊界,如果超出,則基因用隨機(jī)方法重新生成或取邊界值。假設(shè)是經(jīng)過變異或者交叉操作后所得某個(gè)新個(gè)體的第j個(gè)基因,邊界檢查和處理算法為

2 Ⅱ型頻率采樣法設(shè)計(jì)FIR數(shù)字濾波器的原理

設(shè)待設(shè)計(jì)的理想濾波器頻率特性為Hd(ejω),對其在ω的一個(gè)周期[0,2π)等間隔地抽樣N個(gè)點(diǎn)(第一個(gè)抽樣點(diǎn)在ω=π/N處),得

以上便是Ⅱ型頻率采樣設(shè)計(jì)法的原理,H(ejω)即為設(shè)計(jì)出的濾波器頻率響應(yīng)。

3 線性相位FIR數(shù)字濾波器的Ⅱ型頻率采樣設(shè)計(jì)公式

FIR數(shù)字濾波器具有線性相位的約束條件是h(n)是實(shí)序列且滿足h(n)=±h(N-1-n)。將Hd(ejω)表示成幅度函數(shù)Hdg(ω)(可正可負(fù)的實(shí)函數(shù))和相位函數(shù)θ(ω)的形式,即 Hd(ejω)=Hdg(ω)ejθ(ω)。 從而有

將式(18)和式(19)分別代入式(13)中,經(jīng)化簡后得頻率響應(yīng)H(ejω)的表示式為

與理想頻響 Hd(ejω)嚴(yán)格相等;在抽樣點(diǎn)之間,H(ejω)為 Hd(k)的加權(quán)內(nèi)插函數(shù)的延伸疊加;Hd(ejω)變化越平緩,內(nèi)插越接近理想值,逼近誤差越小。圖2的仿真曲線證實(shí)了這一結(jié)論,其中Hg(ω)是實(shí)際濾波器頻響函數(shù)H(ejω)去掉相位特性后的幅度響應(yīng)函數(shù)。通常設(shè)計(jì)出的濾波器的通帶最大波動Rp為1.2 dB左右,阻帶最小衰減As為18 dB左右,不能滿足一般工程上的要求[1-2];可以在濾波器過渡帶內(nèi)增加一些樣點(diǎn),以達(dá)到提高阻帶最小衰減As和減小通帶最大波動Rp的目的。本文就是采用MDE優(yōu)化算法尋找出這些樣點(diǎn)的最佳值。

圖2 Ⅱ型頻率抽樣法設(shè)計(jì)的FIR高通數(shù)字濾波器幅度響應(yīng)(N=25、ωp=0.55π)

4 實(shí)驗(yàn)與結(jié)果分析

4.1 FIR高通數(shù)字濾波器設(shè)計(jì)實(shí)驗(yàn)

用Ⅱ型頻率采樣法設(shè)計(jì)一個(gè)第一類線性相位FIR高通數(shù)字濾波器,要求理想頻率響應(yīng)為陡變矩形,通帶截止頻率ωp=0.43π,阻帶最小衰減As≥40 dB,過渡帶寬度Δω≤0.08π。

解:由于要求As≥40 dB,確定濾波器過渡帶樣點(diǎn)數(shù)為m=1。FIR濾波器長度為N≥2π(m+1)/Δω=50,因?yàn)樵O(shè)計(jì)的是高通濾波器,N必須為奇數(shù),取N=51。通帶范圍 k的取值為「ωpN/(2π)-1/2」≤k≤N-1-「ωpN/(2π)-1/2」(表示對 d朝+∞方向取最小整數(shù)),即11≤k≤39,故過渡帶的樣點(diǎn)值D1在k=10、40處,全部頻率響應(yīng)的樣點(diǎn)值如下:

本文以獲得最大的濾波器阻帶最小衰減值為優(yōu)化目標(biāo),因此計(jì)算個(gè)體適應(yīng)度值函數(shù)為

其中:DFT(·)和IDFT(·)分別為離散傅里葉變換和離散傅里葉逆變換運(yùn)算,為節(jié)省計(jì)算時(shí)間,采用快速的FFT和IFFT算法。為了只保留濾波器頻譜函數(shù)中幅度特性部分,式中乘以以抵消相位特性信息,min(·)為取最小值運(yùn)算。

在安裝有Windows10 64位操作系統(tǒng)的Intel(R)Celeron(R)N3350@1.10 GHz CPU、4.00GB RAM 的PC機(jī)上用MATLAB7.0編寫程序?qū)崿F(xiàn)本例的MDE算法。其中算法參數(shù)設(shè)置為:種群規(guī)模大小M=250,解空間維數(shù)為D=1,尋優(yōu)范圍上界值U1=1和下界值L1=0,最大進(jìn)化代數(shù)Gmax=100。實(shí)驗(yàn)結(jié)果如圖3和圖4所示,可見采用LUT法和MDE算法設(shè)計(jì)的FIR高通濾波器都達(dá)到了設(shè)計(jì)指標(biāo),但LUT法設(shè)計(jì)的濾波器通帶最大波動值(0.5202 dB)大于MDE算法(0.4721 dB),且阻帶最小衰減值(40.1331 dB)比MDE算法(43.0613 dB)小,所以MDE算法效果優(yōu)越。該例MDE算法在經(jīng)過22代的進(jìn)化搜索后獲得最佳解D1=0.38324,對應(yīng)的適應(yīng)度值為1.4226×102;程序運(yùn)行的平均時(shí)間(取10次運(yùn)行的平均值)為53.3141秒。

圖3 Ⅱ型頻率抽樣法設(shè)計(jì)的FIR高通數(shù)字濾波器幅頻響應(yīng)曲線(N=51、ωp=0.43π、D=1)

圖4 設(shè)計(jì)FIR高通數(shù)字濾波器時(shí)MDE算法適應(yīng)度進(jìn)化曲線

4.2 FIR低通數(shù)字濾波器設(shè)計(jì)實(shí)驗(yàn)

利用Ⅱ型頻率采樣法設(shè)計(jì)一個(gè)第一類線性相位FIR低通濾波器,理想通帶截止頻率ωp=0.35π,允許過渡帶寬Δω=0.12π,阻帶最小衰減As=60 dB。

解:因?yàn)樽鑾ё钚∷pAs=60 dB,過渡帶樣點(diǎn)數(shù)取為m=2,FIR濾波器長度為N=2π(m+1)/Δω=50。通帶范圍0～ωp相當(dāng)于k的范圍為0≤k≤表示對 d 朝 0 方向取最大整數(shù)),即0≤k≤8,故過渡帶的抽樣值D1、D2應(yīng)分別選在k=9、40和k=10、39處,寫出頻率響應(yīng)的抽樣點(diǎn)值為

在用MATLAB7.0編程實(shí)現(xiàn)該算法時(shí),計(jì)算個(gè)體適應(yīng)度值函數(shù)與式(23)相同,只是min函數(shù)操作的對象是設(shè)計(jì)出濾波器幅度函數(shù)Hg(ω)的前512個(gè)樣點(diǎn)值(因?yàn)榈谝活惥€性相位FIR數(shù)字濾波器在單位脈沖響應(yīng)長度N為偶數(shù)時(shí),Hg(ω)關(guān)于ω=π奇對稱);算法參數(shù)與4.1節(jié)完全相同,不同的是解空間維數(shù)D=2,尋優(yōu)范圍上界值U1,2=1,下界值L1,2=0。實(shí)驗(yàn)結(jié)果如圖5和圖6所示,可見MDE算法設(shè)計(jì)出的FIR低通濾波器通帶最大波動值(0.3006 dB)小于LUT法(0.3603 dB),阻帶最小衰減值(73.7692 dB)比LUT法(59.8452 dB)大得多,并且LUT法設(shè)計(jì)的濾波器阻帶最小衰減值沒有達(dá)到設(shè)計(jì)指標(biāo),因而MDE算法效果明顯優(yōu)于傳統(tǒng)LUT法。該例MDE算法在經(jīng)過77代的進(jìn)化搜索后獲得最佳解 D1=0.54297、D2=0.084864,對應(yīng)的適應(yīng)度值為4.88139×103;程序的平均運(yùn)行時(shí)間為52.8375秒。

圖5 Ⅱ型頻率抽樣法設(shè)計(jì)的FIR低通數(shù)字濾波器幅頻響應(yīng)曲線(N=50、ωp=0.35π、D=2)

圖6 設(shè)計(jì)FIR低通數(shù)字濾波器時(shí)MDE算法適應(yīng)度進(jìn)化曲線

4.3 FIR帶通數(shù)字濾波器設(shè)計(jì)實(shí)驗(yàn)

試?yán)芒蛐皖l率采樣法設(shè)計(jì)一個(gè)理想頻率響應(yīng)為矩形的線性相位FIR帶通濾波器,技術(shù)指標(biāo)為:通帶下邊界頻率和上邊界頻率分別為 ωp1=0.35π和 ωp2=0.7π,阻帶下邊界頻率和上邊界頻率分別為 ωs1=0.25π和 ωs2=0.8π,通帶波紋 Rp=0.5 dB,阻帶最小衰減 As≥60 dB。

解:需要過渡帶抽樣點(diǎn)數(shù)為m=2,過渡帶寬度Δω=ωp1-ωs1=ωs2-ωp2=0.1π,于是濾波器的長度點(diǎn)數(shù)為 N≥2π(m+1)/Δω=60,取 N=60,各邊界頻率的抽樣點(diǎn)值 k分別為:。故過渡帶抽樣值D1、D2應(yīng)分別位于 k=8、22、37、51 和 k=9、21、38、50處,于是有

在用MATLAB7.0編程實(shí)現(xiàn)時(shí),算法參數(shù)以及計(jì)算個(gè)體適應(yīng)度值函數(shù)與4.2節(jié)完全相同。實(shí)驗(yàn)結(jié)果如圖7和圖8所示,由結(jié)果可以看出:采用MDE算法優(yōu)化設(shè)計(jì)的FIR帶通濾波器通帶最大波動Rp和阻帶最小衰減As均滿足設(shè)計(jì)要求,且具有比傳統(tǒng)LUT法更好的幅頻響應(yīng)性能,即濾波器的通帶波動較小、阻帶衰減更大。該例MDE算法在經(jīng)過52代的進(jìn)化搜索后獲得最佳解D1=0.11173、D2=0.59742,對應(yīng)的適應(yīng)度值為1.21487×103;程序的平均運(yùn)行時(shí)間為54.1125秒。

圖7 Ⅱ型頻率抽樣法設(shè)計(jì)的FIR帶通數(shù)字濾波器幅頻響應(yīng)曲線(N=60、ωs1=0.25π、ωp1=0.35π、ωp2=0.7π、ωs2=0.8π、D=2)

圖8 設(shè)計(jì)FIR帶通數(shù)字濾波器時(shí)MDE算法適應(yīng)度進(jìn)化曲線

4.4 實(shí)驗(yàn)結(jié)果數(shù)據(jù)對比分析

表1為實(shí)驗(yàn)結(jié)果數(shù)據(jù)對比,其中GA法數(shù)據(jù)為根據(jù)3個(gè)不同實(shí)驗(yàn)的適應(yīng)度函數(shù)調(diào)用MATLAB7.0優(yōu)化工具箱內(nèi)置函數(shù)[x,fval]=ga(fitnessfun,nvars,options)[22]求得最優(yōu)解的結(jié)果。可以清楚地看出,本文提出的MDE算法設(shè)計(jì)出的FIR高通、低通、帶通數(shù)字濾波器各項(xiàng)指標(biāo)均達(dá)到設(shè)計(jì)要求,具有比傳統(tǒng)LUT法和GA法更為優(yōu)越的性能(通帶最大波動Rp更小、阻帶最小衰減As更大)。由于MDE算法保留最好的個(gè)體到下一代中,所以在進(jìn)化的過程中,適應(yīng)度值是一隨進(jìn)化代數(shù)增加而單調(diào)上升的過程,算法逐漸收斂且收斂精度高。多次運(yùn)算結(jié)果表明:MDE算法能以最短的時(shí)間、最為精確地搜索到全局最優(yōu)解,且算法的健壯性更強(qiáng)。

表1 MDE算法與LUT法、GA法的實(shí)驗(yàn)結(jié)果比較

5 結(jié)束語

本文采用Ⅱ型頻率采樣法設(shè)計(jì)FIR數(shù)字濾波器的基本原理,推導(dǎo)第一類線性相位FIR數(shù)字濾波器的Ⅱ型頻率采樣設(shè)計(jì)公式;針對傳統(tǒng)差分進(jìn)化算法的特點(diǎn),提出一種改進(jìn)的差分進(jìn)化算法并設(shè)計(jì)相應(yīng)的適應(yīng)度函數(shù),應(yīng)用于采用Ⅱ型頻率采樣法設(shè)計(jì)FIR數(shù)字濾波器的過渡帶樣點(diǎn)值優(yōu)化中。實(shí)驗(yàn)驗(yàn)證了改進(jìn)差分進(jìn)化算法控制參數(shù)少、結(jié)構(gòu)簡單、收斂速度更快、尋優(yōu)精度更高、穩(wěn)定性和魯棒性更好,設(shè)計(jì)出的FIR數(shù)字濾波器性能顯著優(yōu)于傳統(tǒng)查表法和遺傳算法。若更改優(yōu)化的目標(biāo)函數(shù),該改進(jìn)差分進(jìn)化算法也可以實(shí)現(xiàn)常規(guī)計(jì)算方法難以求解的高維、多目標(biāo)、多模型、多峰值、非線性、非凸形等復(fù)雜環(huán)境中的函數(shù)優(yōu)化問題和組合優(yōu)化問題等。

致謝:感謝2022年成都信息工程大學(xué)大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練計(jì)劃項(xiàng)目(202210621118);2021年成都信息工程大學(xué)電子工程學(xué)院教學(xué)改革項(xiàng)目對本文的資助