胡仕兵, 聶 喜, 陳子為

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

0 引言

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

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

在確定FIR數(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)解)和程序運行時間較長等缺陷。差分進(jìn)化(differential evolution,DE)算法是由Rainer Storn等[19-21]于20世紀(jì)末提出的一種基于模擬自然界生物的進(jìn)化過程、求解連續(xù)空間上優(yōu)化問題的高效、強大的隨機種群搜索算法,由于其控制參數(shù)少、計算簡單、易于實現(xià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è)計FIR數(shù)字濾波器時過渡帶樣點值的優(yōu)化中,取得了良好的實驗效果。

1 MDE算法

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

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

1.1 種群初始化

設(shè)進(jìn)化代數(shù)為 G∈{0,1,…,Gmax}時的種群為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個個體(染色體),xi,j(G)表示第 G代中第 i個個體的第 j個分量(基因)。每個個體代表解空間內(nèi)的某一個解,而基因代表解的各個分量。種群初始化就是在解空間范圍內(nèi)隨機均勻地產(chǎn)生M個個體,每個個體由D維向量組成,取值方式如下:

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

1.2 變異運算

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

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

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

1.3 交叉運算

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

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

常數(shù)CR值較大時,試驗個體中的信息大多來自于變異個體,使種群的多樣性越來越好,但容易忽略掉全局最優(yōu)解,算法容易出現(xiàn)早熟收斂現(xiàn)象;CR值較小時,試驗個體中的信息大多來自目標(biāo)個體,會使種群性降低,導(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 選擇運算

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

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

1.5 邊界檢查和處理

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

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

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

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

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

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

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

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

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

4 實驗與結(jié)果分析

4.1 FIR高通數(shù)字濾波器設(shè)計實驗

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

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

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

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

在安裝有Windows10 64位操作系統(tǒng)的Intel(R)Celeron(R)N3350@1.10 GHz CPU、4.00GB RAM 的PC機上用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。實驗結(jié)果如圖3和圖4所示,可見采用LUT法和MDE算法設(shè)計的FIR高通濾波器都達(dá)到了設(shè)計指標(biāo),但LUT法設(shè)計的濾波器通帶最大波動值(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;程序運行的平均時間(取10次運行的平均值)為53.3141秒。

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

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

4.2 FIR低通數(shù)字濾波器設(shè)計實驗

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

解:因為阻帶最小衰減As=60 dB,過渡帶樣點數(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)的抽樣點值為

在用MATLAB7.0編程實現(xiàn)該算法時,計算個體適應(yīng)度值函數(shù)與式(23)相同,只是min函數(shù)操作的對象是設(shè)計出濾波器幅度函數(shù)Hg(ω)的前512個樣點值(因為第一類線性相位FIR數(shù)字濾波器在單位脈沖響應(yīng)長度N為偶數(shù)時,Hg(ω)關(guān)于ω=π奇對稱);算法參數(shù)與4.1節(jié)完全相同,不同的是解空間維數(shù)D=2,尋優(yōu)范圍上界值U1,2=1,下界值L1,2=0。實驗結(jié)果如圖5和圖6所示,可見MDE算法設(shè)計出的FIR低通濾波器通帶最大波動值(0.3006 dB)小于LUT法(0.3603 dB),阻帶最小衰減值(73.7692 dB)比LUT法(59.8452 dB)大得多,并且LUT法設(shè)計的濾波器阻帶最小衰減值沒有達(dá)到設(shè)計指標(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;程序的平均運行時間為52.8375秒。

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

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

4.3 FIR帶通數(shù)字濾波器設(shè)計實驗

試?yán)芒蛐皖l率采樣法設(shè)計一個理想頻率響應(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。

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

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

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

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

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

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

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

5 結(jié)束語

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

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