第11章線性相位濾波器與自適應(yīng)濾波器11.1線性相位濾波器11.2自適應(yīng)濾波器11.1線性相位濾波器

11.1.1線性相位與線性相位濾波器

1.線性相位圖例

設(shè)濾波前初始信號y(t)有三個頻率分量y1(t)、y2(t)和y3(t)，y(t)=y1(t)+y2(t)+y3(t)，三個分量的頻率依次為ω1=2πf1，ω2=2πf2和ω3=2πf3，且f2=2f1，f3=3f1。于是，令ω=ω1，有

y(t)=sinω1t+sinω2t+sinω3t=sinωt+sin2ωt+sin3ωt

y(t)及其分量的波形如圖11-1(a)所示。圖11-1線性相位圖示

(a)濾波前的三個頻率分量；(b)濾波后的三個頻率分量經(jīng)過濾波后，基波分量y1(t)=sinωt產(chǎn)生滯后時間t0，折合成相位為

j1=ω1t0=ωt0(11-1)

對于二次諧波分量y2(t)=sin2ωt與三次諧波分量y3(t)=sin3ωt而言，它們也必須滯后相同的時間t0，才能保證合成的信號y(t)的波形不會發(fā)生畸變,如圖11-1(b)所示。同樣的時間t0，對于不同的頻率ω2和ω3有不同的相位。二次諧波y2(t)的相位滯后j2為

j2=ω2t0=2ωt0(11-2)三次諧波y3(t)的相位滯后j3為

j3=ω3t0=3ωt0(11-3)

推論：若將一信號以采樣間隔T0進(jìn)行采樣，共采樣N點，經(jīng)FFT(快速傅里葉變換)得到N條譜線，其中有效譜線N/2條，序號k為0,1,…,N/2－1；則基波k=1和k≥2的諧波必須有線性相移

jk=ω3t0=kωt0

諧波次數(shù)(11-4)

于是，濾波前后信號波形保持不變。

2.線性相位濾波器概念

線性濾波器對各次諧波分量產(chǎn)生相同時間t0的滯后，各次諧波k相移有如式(11-4)所示的與頻率ω的線性比例關(guān)系。11.1.2線性相位有限沖擊響應(yīng)濾波器的數(shù)學(xué)模型

FIR脈沖傳遞函數(shù)表達(dá)式及其對應(yīng)差分方程已由4.5節(jié)

式(4-67)和式(4-68)給出，脈沖傳遞函數(shù)及對應(yīng)差分方程重寫如下：H(z)=

=C0+C1z－1+···+CN－1z－(N－1)=(11-5a)y(n)=c0x(n)+c1x(n－1)+c2x(n－2)+…+cN－1x(n－N+1)=(11-5b)根據(jù)4.5節(jié)內(nèi)容可知，由一個等效模擬濾波器H(s)經(jīng)采用不同的變換法，如后向差分法、雙線性變換法等，就得到了Z域的脈沖傳遞函數(shù)H(z)：H(z)=

=

(11-6)但求得的是無限沖擊響應(yīng)(IIR)濾波器的脈沖傳遞函數(shù)，對式(11-6)進(jìn)行多項式除法，保留前N項，即可獲得如式(11-5)所示的N－1階FIR濾波器。

1.一階低通模擬濾波器生成的FIR濾波器

一階低通模擬濾波器的S域傳遞函數(shù)Hd(s)為(11-7)式中a=ωcd=100rad/s，為截止角頻率，相應(yīng)的截止頻率fcd=ωcd/2π=15.915Hz。采用雙線性變換法求與Hd(s)等效的Z域脈沖傳遞函數(shù)

Hd(z)，得到經(jīng)整理可得(11-8)(11-9)對式(11-9)進(jìn)行多項式直接相除，其商的表達(dá)式即為FIR的一階低通濾波器的脈沖傳遞函數(shù)

Hd(z)=c0+c1z－1+c2z－2+…+cN－1z－(N－1)(11-10)

式中c0=aT/(2+aT)，c1=4aT/(2+aT)2，c2=－4aT(aT－2)/

(2+aT)3，…，cN－1=－4aT(aT－2)N－2/(2+aT)N。把ci(i=0,1,···,N－1)的值代入式(11-5)得：

2.一階高通模擬濾波器生成的FIR濾波器

一階高通模擬濾波器的S域傳遞函數(shù)Hg(s)為(11-11)式中a=ωcd=100rad/s，為截止角頻率，相應(yīng)的截止頻率fcd=ωcd/2π=15.915Hz。采用雙線性變換法，也就是把代入到式(11-11)中，得到與Hg(s)等效的Z域脈沖傳遞函數(shù)Hg(z)：(11-12)與一階低通濾波器FIR系數(shù)的求取方法一樣，對式(11-12)進(jìn)行多項式直接相除，其商的表達(dá)式即為FIR一階高通濾波器的脈沖傳遞函數(shù)，其表達(dá)式的形式也與式(11-10)完全一樣，只是此時的ci(i=0,1,…,N－1)分別為c0=2/(2+aT)，c1=－4aT/(2－aT)2，c2=－4aT(2－aT)/(2+aT)3，…，cN－1=－4aT(2－aT)N－2/(2+aT)N。當(dāng)采樣周期為T=1.0ms時，aT=0.1，一階高通FIR濾波器的系數(shù)分別為

3.有限沖激響應(yīng)濾波器(FIR)的頻率特性H(ω)由Z變換知識z=ejωT，將其代入式(11-5a)可得(11-13)而時間序列h(n)的離散傅里葉變換式為(11-14)式(11-14)與式(11-13)形式相同，故式(11-14)中的傅里葉系數(shù)h(n)就是FIR系統(tǒng)的濾波系數(shù)Cn，即

h(0)=C0,h(1)=C1,…,h(N－1)=CN－1;N→∞

4.生成線性相位FIR濾波器需要考慮的問題

1)無限長時間序列的不可實現(xiàn)與有限截取

由式(11-5b)差分方程可見，當(dāng)前時刻(nT)的輸出y(n)是由當(dāng)前時刻的輸入值x(n)乘以系數(shù)C0，以及無限長歷史時刻的輸入值x(n－1)，x(n－2)，…，x(n－N+1)乘以相應(yīng)的系數(shù)組成的無限長序列。實際上，若N趨于無窮大，這是不可能實現(xiàn)的，故必須對輸入x(nT)進(jìn)行有限截取，截取長度N應(yīng)為有限值。

2)具有線性相位FIR濾波器的形式

當(dāng)FIR系統(tǒng)的系數(shù)Cn(h(n))對稱時，濾波器將具有線性相位。對稱有偶對稱和奇對稱兩種情況。對于偶對稱，即FIR濾波器系數(shù)h滿足關(guān)系式h(N－1－n)=h(n)；而對于奇對稱，即FIR濾波器系數(shù)h滿足關(guān)系式h(N－1－n)=－h(huán)(n)。另外，F(xiàn)IR系統(tǒng)系數(shù)的長度有奇數(shù)個和偶數(shù)個之分，因此需要分為四種情況來論證。（1）奇數(shù)長度偶對稱。

奇數(shù)長度偶對稱的濾波器有奇數(shù)個參數(shù)，可以分成三部分，前(N－3)/2個參數(shù)為一部分，后(N－3)/2個參數(shù)為一部分，中間的參數(shù)h((N－1)/2)為一部分。于是，該濾波器的離散傅里葉變換可寫為(11-15)令n=m+(N－1)/2，則h(n)e－jωn=h((N－1)/2+m)e－jω((N－1)/2+m),于是把m改為n，即得到對于一個類似1+2+…+100的求和計算，我們知道，它等同于100+99+…+1的求和計算。因此，有再考慮到濾波器是偶對稱的，h(N－1－n)=h(n)，因此，把上式中的h((N－1)/2－n)替換為

h(N－1－［(N－1)/2－n］)=h((N－1)/2+n)得到：于是，式(11-15)可以重寫為從上式右邊三項中提取e－jω（N－1）/2得到：(11-16)考慮到ejωn+e－jωn=2cos（ωn），式(11-16)可以寫作(11-17)由式(11-17)中可以看出，H(jω)的幅值為對于給定的ω，它是一個確定的數(shù)值。H(jω)的相位θ=－ω(N－1)/2，負(fù)號的存在表明相位是滯后的。當(dāng)y和x滿足形如y=ax+b的關(guān)系式時，我們稱y和x是線性的，其中

x為自變量，y為因變量，a為比例系數(shù)，b為截長。明顯地，H(jω)的相位θ與頻率ω滿足此條件，自變量為ω，因變量為θ，比例系數(shù)為－(N－1)/2，截長為0，因此，相位與頻率滿足線性關(guān)系，也就是說H(jω)具有線性相位。在采用形如式(11-17)的濾波器濾波后的信號中，所有的頻率分量的滯后時間為Δt=θ/ω=(N－1)/2。注意到離散傅里葉變換中，所選時間窗長度為N-1，忽略采樣時間參數(shù)Ts，實為(N－1)Ts，因此，經(jīng)濾波后的信號中，各分量的滯后時間均為Ts(N－1)/2。對于選定長度的濾波器，它是一個常數(shù)，經(jīng)濾波后的信號中，各頻率分量合成的信號等同于原信號的通帶范圍內(nèi)各頻率分量分別滯后Ts(N－1)/2后重新合成的信號，也等同于原信號的通帶范圍內(nèi)各頻率分量的合成信號整體滯后Ts(N－1)/2，因此不會發(fā)生畸變。

（2）偶數(shù)長度偶對稱。

偶數(shù)長度的濾波器系數(shù)可分成前半部分和后半部分兩部分，因此其傅里葉變換可寫作(11-18)把

寫作考慮到h(N－1－n)=h(n)，可以把寫作于是，從式（11-18）右邊兩項中提取e－jω(N/2－1）得到：(11-19)（3）奇數(shù)長度奇對稱。有（11-20）把寫作考慮到h(N－1－n)=h(n)，可以把(11-21)于是，從上式右邊兩項中提取e－jω(N/2－1)得到：考慮到ejωn－e－jωn=j2sin(ωn)，上式可以寫作(11-22)明顯地，H(jω)的滯后相位為ω(N－1)/2+π/2，它也是線性的。（4）偶數(shù)長度奇對稱。有寫作把(11-23)考慮到h(N－1－n)=－h(huán)(n)，可以把于是，從式(11-23)右邊兩項中提取e－jω(N/2－1)得到：11.1.3線性相位FIR濾波器的窗口設(shè)計

1.線性相位FIR濾波器窗口設(shè)計的基本思路

線性相位FIR濾波器窗口設(shè)計的基本思路是：瞄準(zhǔn)目標(biāo)是理想頻率特性，可以是理想低通、高通、帶通或帶阻，根據(jù)需要而定。

第一步，由傅里葉反變換求其頻率特性Hd(ω)的離散時間序列hd(n)。

第二步，選擇窗口的大小，也就是FIR濾波器系數(shù)的個數(shù)N，并截取hd(n)的以n=0為中心的中間N項，將時間序列右移延時floor(N/2)。函數(shù)floor(x)表示不大于x的整數(shù)。第三步，選擇窗函數(shù)，并求取長度為N的窗函數(shù)系數(shù)。將hd(n)的中間N項與窗函數(shù)系數(shù)對應(yīng)項相乘，即得到FIR濾波器的系數(shù)h(n)。

第四步，求h(n)的頻譜H(ω)，該頻譜與期望的理想頻率特性有差異，差異與截取的長度和所選擇的窗函數(shù)有關(guān)。第五步，根據(jù)采用窗函數(shù)的特點分析獲得的頻譜過渡帶斜率及截止角頻率ωc，若不滿足設(shè)計要求，則調(diào)整窗函數(shù)重復(fù)第二步至第四步，直至滿足要求。

2.線性相位低通濾波器的設(shè)計

期望的理想低通濾波器頻率特性Hd(ω)為(1) 對Hd(ω)進(jìn)行傅里葉反變換，以獲取時間序列hd(n)。離散傅里葉反變換的形式為(11-25)(11-26)將頻域擴(kuò)展到－ωc，式(11-25)的值代入式(11-26)可得：(11-27)式中sinc(nωc)=sin(nωc)/nωc，稱之為sinc函數(shù)，它是以2π為周期的。當(dāng)0≤n<π/ωc時，sinc函數(shù)大于0；當(dāng)n=π/ωc時，sinc函數(shù)第一次經(jīng)過0點；當(dāng)2π/ωc>n>π/ωc時，sinc函數(shù)小于0。由于隨著n的增大，sinc函數(shù)的分母成比例增大，于是sinc函數(shù)隨著n的增加而作衰減振蕩。另外，對于給定的截止頻率ωc，Tc是一個常數(shù)，因此，hd(n)是以n=0為中心對稱，其幅值呈sinc函數(shù)形式衰減的無限長時間序列，其形狀如圖11-2(b)所示。圖11-2理想矩形幅頻特性及無限長序列的插入(a)理想矩形幅頻特性|Hd(ω)|；

(b)|Hd(ω)|對應(yīng)放入無限長時間序列hd(n)（2)將無限長序列hd(n)有限截短，并向右移延時。如圖11-2(b)所示，已知hd(n)是以n=0為中心對稱的無限長序列，因此物理實現(xiàn)是不可能的。為了讓FIR濾波成為可能，需要從hd(n)中提取一段有限個系數(shù)來近似hd(n)，這個過程稱之為截短。從圖11-2(b)還可以看出，n=0時，hd(n)的絕對值最大，然后向n的正負(fù)方向振蕩衰減，直至趨于0，

因此，截短一般是截取以n=0為中心的一段長度為N的系數(shù)hN(n)。由于hd(n)是關(guān)于hd(0)偶對稱的，N一般取為正的奇數(shù)。另外，負(fù)時間(－n)的序列h(－n)與正時間(n)的序列h(n)以n=0為中心對稱，但負(fù)時間序列是真實物理世界不存在的，且又不可見，這稱之為非因果的，也是不可實現(xiàn)的，因此需要對截短后的序列hN(n)向右移動(N－1)/2，也就是將原來的序列hN［(N－1)/2］,…,hN(－1),hN(0),hN(1),…,

hN［(N－1)/2］變成hN(0),hN(1),…,hN(N－1)。

（3)加窗。hd(n)的有限截短實際上是對hd(n)加了一個矩形窗，在這個窗之內(nèi)的系數(shù)可以全部看到，在這個窗之外的系數(shù)全部被置為0，這等效于給FIR濾波器系數(shù)另外再乘了一個系數(shù)，只是在窗內(nèi)的系數(shù)為1，窗外的系數(shù)為0。在數(shù)學(xué)上，這個過程可以寫作(11-28)式中RN(n)為矩形窗函數(shù)，其函數(shù)式為(11-29)除了矩形窗之外，為了獲得不同指標(biāo)的濾波器，有關(guān)研究工作者對矩形窗做了改進(jìn)，窗內(nèi)的系數(shù)并非全為1，而是從中心點開始，濾波器系數(shù)向左、向右逐步衰減，從而獲得了其它形式的窗函數(shù)。對于選用矩形窗來說，加窗的過程實質(zhì)上在截短的過程中已經(jīng)完成。若選用的是其它窗函數(shù)，則首先要求取窗函數(shù)的系數(shù)，然后代入到式(11-28)中，才能求得濾波器系數(shù)h(n)。（4)時間序列h(n)的頻率H(ω)。

根據(jù)傅里葉變換中頻域卷積特性可知，時域兩個函數(shù)相乘，其積的頻譜是各自頻譜函數(shù)的卷積，因此，經(jīng)過將無限長時間序列hd［(N－1)/2］與窗函數(shù)相乘后，被截短為有限長時間序列h(n)，其頻譜H(ω)可由頻域卷積定理來求。下面分別求時域兩函數(shù)各自的頻譜函數(shù)，再求它們頻譜函數(shù)的卷積。①無限長時間序列的頻譜函數(shù)。

根據(jù)傅里葉變換中的時移特性，將信號在時域中沿時間軸平移一個常值而成為時延信號，在頻域其幅值不變，而相頻譜中相位角的改變與頻率成正比。已知未延時時間

序列的頻譜為Hd(ω)時，時延序列

的頻譜為(11-30)其中j′(ω)=j(ω)+ω(N－1)/2，相位角的改變量為Δj=j′(ω)－j(ω)=ω(N－1)/2，當(dāng)N為常數(shù)時，相位角的改變與頻率ω成正比。當(dāng)初始相位角j(ω)=0時有

特別地，理想低通濾波器的頻譜函數(shù)具有式(11-27)所示的形式。(11-31a)(11-31b)②窗函數(shù)的頻譜。已知未時延窗函數(shù)時間序列即j(ω)=0。經(jīng)時延(N－1)/2后有(11-32)按傅里葉正變換定義式可求WN(ω)。特別地，對于矩形窗，時延(N－1)/2的頻譜為(11-33)矩形窗函數(shù)及其幅值譜如圖11-3所示。圖11-3矩形窗函數(shù)及幅值譜(a)矩形窗；(b)矩形窗頻譜；(c)矩形窗幅值譜③求兩時間序列hd(n)與wN(n)之乘積h(n)的頻譜。

根據(jù)傅里葉變換中的頻域卷積特性，將該兩時間序列的頻譜進(jìn)行卷積。已知式(11-31)和式(11-32)成立，則乘積(11-34)將式(11-31)和式(11-32)代入式(11-34)，由卷積定義有(11-35a)特別地，根據(jù)式(11-25)可知，在－ωc至ωc區(qū)間，Hd(ω)=1，將式(11-31a)和式(11-33)代入式(11-35a)可得窗函數(shù)為矩形窗的低通濾波器系數(shù)h(n)的頻譜為(11-35b)式(11-35b)卷積的結(jié)果如圖11-4所示,圖中曲線表示離散時刻幅值的包絡(luò)。圖11-4卷積結(jié)果圖(a)Hd(ω)——無限長時間序列hd(n)的頻譜；(b)WR(ω)——矩形窗的頻譜；(c)H(ω)——有限長時間序列h(n)，n=0,1,…,N－1的頻譜（5)有限長時間序列h(n)頻譜H(ω)的相關(guān)參數(shù)。

由圖11-4可看出，加矩形窗后使實際頻率響應(yīng)偏離理想頻率響應(yīng)，主要影響有三個方面：首先是理想幅頻特性陡直邊緣處形成過渡帶，過渡帶寬取決于矩形窗函數(shù)頻率響應(yīng)的主瓣寬度；其次是過渡帶兩側(cè)形成肩峰和波紋，這是矩形窗函數(shù)頻率響應(yīng)的旁瓣引起的，旁瓣相對值越大，旁瓣越多，波紋越多；第三是隨窗函數(shù)寬度N的增大，矩形窗函數(shù)頻率響應(yīng)的主瓣寬度減小，但不改變旁瓣的相對值。以矩形窗為例，從圖11-4(c)所示的H(ω)可以看出，當(dāng)窗函數(shù)寬度為N時，加窗截取后的頻譜圖存在過渡帶及泄漏(即截止頻率之外的幅頻特性不為0)。

①過渡帶寬定義為

Δω=ωs－ωp

式中ωp為通帶截止頻率，ωs為阻帶允許最小衰減頻率。②相對衰耗定義為

A(ω)=20lg［|H(ejω)|/|H(ej0)|］=20lg[H(ω)/H(0)](11-36)濾波器的幾乎所有重要指標(biāo)都是由窗函數(shù)決定的，改進(jìn)濾波器的關(guān)鍵在于改進(jìn)窗函數(shù)。

為了便于研究，我們把窗譜W(ω)也化為相對衰耗

20lg［W(ejω)/W(ej0)］的形式。

窗譜的兩個最重要的指標(biāo)是主瓣寬度和旁瓣峰值衰耗。旁瓣峰值衰耗定義為旁瓣峰值衰耗＝20lg(第一旁瓣峰值/主瓣峰值)(11-37)

對于矩形窗而言，旁瓣峰值衰耗為－13dB。為了改善FIR濾波器性能，必須修改窗函數(shù)，使其具有更好的窗譜。一個好的窗譜，應(yīng)滿足以下兩方面的條件：主瓣盡可能窄，以使設(shè)計出來的濾波器有較陡的過渡帶；

第一副瓣面積相對主瓣面積盡可能小，即能量盡可能集中在主瓣，外泄少，這樣設(shè)計出來的濾波器才能尖峰和余振小。

對任一個具體的窗函數(shù)而言，以上兩個條件互相矛盾，不能同時滿足。我們所能做的是根據(jù)具體設(shè)計指標(biāo)，選擇一種能兼顧各項指標(biāo)的相對最佳的窗口。

3.幾種常用的窗口函數(shù)簡介

(1)矩形窗:

ω(n)=RN(n)

(11-38)

矩形窗的特點是旁瓣較大，尤其是第一個旁瓣峰太高，達(dá)到主瓣高度的21%，所以泄漏很大。矩形窗的優(yōu)點是容易獲得，而且主瓣寬度小，其等效帶寬為1/T。(11-39)它是由兩個長度為N/2的矩形窗進(jìn)行線性卷積而得到的。

(3)漢寧（Hanning）窗，又稱為余弦窗:漢寧窗的主要思路是：通過矩形窗譜的合理疊加減小旁瓣面積。因此可寫出其窗譜(11-4１)(11-40)(11-42)式中WR(ejω)為矩形窗譜。

(4)海明（Hamming）窗:式（11-43）可轉(zhuǎn)化為(11-43)由此式可知，海明窗是在矩形窗上拼結(jié)一個漢寧窗而形成的，它包括一個高為0.08的矩形窗和一個最大高度為0.92的漢寧窗。（5)布萊克曼（Blackman）窗：

(11-44)它是5個矩形窗譜的疊加。前兩項與漢寧窗類似，是兩個矩形窗譜的疊加；后一項構(gòu)成兩個調(diào)整副譜，使旁瓣進(jìn)一步壓低，但主瓣寬度增加到12π/N。（6)凱塞(Kaiser)窗的表達(dá)式是:式中，I0［·］是第一類修正的零階貝塞爾函數(shù)，可以用級數(shù)展開來計算它的值：(11-45)(11-46)一般取15～25次計算，便可達(dá)到足夠精度。I0［x］在［0,∞］上是單調(diào)增函數(shù)，I0［0］=1。β是決定窗函數(shù)形狀從而決定窗函數(shù)特性的參數(shù)。不管β是多少，窗口的中心點n=(N－1)/2處總是1：ω(N－1)/2=I0［β］/I0［β］=1。而在窗的兩邊n=0和n=N－1處，ω(0)=ω(N－1)=1/I0［β］。β越大，I0［β］就越大，窗包絡(luò)兩邊ω(0)=ω(N－0)就壓得越低。改變β值，就能得到不同特性的窗：β=0時

凱塞窗相當(dāng)于矩形窗，β=5.44時相當(dāng)于海明窗，β=8.5時接近于布萊克曼窗。從頻譜看，β越大，主瓣越寬，旁瓣峰值衰耗越大，實踐中常取4<β<9。幾種常用窗口的包絡(luò)形狀見圖11-5，性能比較見表11-1。其中，旁瓣衰減速度（第一旁瓣高度/主瓣高度）是指旁瓣

峰值下降的漸進(jìn)速度，單位是dB/oct（倍頻程）。從表11-1可明顯地看出：窗函數(shù)的主瓣寬度和旁瓣峰值衰耗是矛盾的，一項指標(biāo)的提高總是以另一項指標(biāo)的下降為代價，窗口選擇實際上是對兩項指標(biāo)作權(quán)衡。以上介紹的前五種窗中，這兩項指標(biāo)是跳變的。于是有人提出可調(diào)整窗，適當(dāng)修改參數(shù)，可在這兩項指標(biāo)間作連續(xù)的選擇。常用的可調(diào)整窗是凱塞窗。圖11-5五種窗函數(shù)波形及其窗譜11.1.4［示例11-1］用窗口法設(shè)計一個線性相位低通FIR濾波器

要求：截止頻率為fc，采樣頻率是8fc；通帶范圍內(nèi)，衰減度不超過5.8dB。

(1)寫出理想的頻響：(11-47)式中ωc=2πfcT=2πfc/fs=0.25π;a=(N－1)/2。

(2)選窗的形狀與大?。毫頝=11，此時，

①如用矩形窗

wR(n)=RN(n)

②如用海明窗,把N＝11和n＝0,1,…,10代入式(11-43)得

w(n)={0.08,0.0168,0.399,0.682,0.912,1,0.912,

0.682,0.399,0.168,0.008}(11-48)

(3)加窗：

h(n)=hd(n)·ω(n)

序列相乘等于對應(yīng)項相乘。

如用矩形窗，由于矩形窗的系數(shù)均為1，因此直接取式(11-27)的中間11項即可:

h(n)={－0.045,0,0.075,0.159,0.225,0.25,0.225,

0.159,0.075,0,－0.045}(11-49)如用海明窗，則用式(11-48)和式(11-47)的中間11項的對應(yīng)項相乘，例如－0.045×0.08=－0.0036，于是得到：

h(n)={－0.0036,0,0.03,0.108,0.205,0.25,0.205,

0.108,0.03,0,－0.0036}(11-50)

h(n)就是所要求的FIR濾波器系數(shù)。從式(11-49)和(11-50)中可以看出，這兩個濾波器的系數(shù)的確是奇數(shù)偶對稱的。

(4)檢驗。h(n)是否合乎要求需要通過檢驗來判別，

特別地，還要檢驗濾波器在截止頻率處的幅頻特性是否滿

足要求。于是，令ω=0.25π，分別把式(11-49)和式(11-50)代入并計算－20lg10(|H(ejω)|)，分別得到5.58dB和5.98dB。因此，矩形窗設(shè)計的濾波器滿足要求，海明窗不滿足要求，若選用海明窗，需要增大濾波器長度。11.1.5［示例11-2］線性相位FIR濾波器在Matlab中的實現(xiàn)

1.函數(shù)介紹

1)firceqrip函數(shù)

調(diào)用格式有多種，其調(diào)用格式與參數(shù)意義分別如下：(1)h=firceqrip(n,wo,del):用來設(shè)計n階具有線性相位的有限沖擊響應(yīng)低通濾波器，該濾波器在阻帶具有等紋波。參數(shù)wo指定截止頻率;del=［d1,d2］指定峰值或者通帶與阻帶的最大允許誤差，其中d1設(shè)置通帶誤差，d2指定阻帶誤差。由于firceqrip采用的是正則化頻率，也就是采樣頻率為fs，低于采樣頻率的信號頻率f正則化為w=2f/fs。顯然，根據(jù)采樣定律，信號頻率不能高于采樣頻率的1/2，即w≤1，因此wo的值必須設(shè)置在(0，1)范圍內(nèi)。

(2)h=firceqrip(…,′slope′,r)：

用輸入關(guān)鍵詞′slope′和r設(shè)計阻帶不具有等紋波特性的濾波器。r以dB為單位，確定阻帶的傾斜度，r>0。

（3）h=firceqrip(…,′passedge′)：指定通帶起始頻率。

（4）h=firceqrip(…,′stopedge′)：指定阻帶起始頻率。

（5）h=firceqrip(…,′high′)：設(shè)計高通濾波器。

（6）h=firceqrip(…,′min′)：設(shè)計最小相位FIR濾波器。（7）h=firceqrip(…,′invsinc′,c))：設(shè)計具有sinc函數(shù)形狀的低通濾波器，關(guān)鍵詞invsinc采用

逆sinc函數(shù)，它由c是標(biāo)量還是二元素向量來確定。

當(dāng)c為標(biāo)量時，通帶采用函數(shù)1/sin(c*w)，其中w是正則化頻率。

當(dāng)c為二元素向量[cp]時，通帶采用函數(shù)1/sin(c*w)p,其中w是正則化頻率。

2)filter函數(shù)

filter函數(shù)根據(jù)濾波器參數(shù)對信號進(jìn)行濾波，并返回濾波后得到的信號。filter函數(shù)也有多種調(diào)用格式，分別介紹如下。

（1）y=filter(b,a,x)：參數(shù)a和b分別指定濾波器分母和分子中的系數(shù)向量。該函數(shù)的作用實際上就是實現(xiàn)如

下運算：

a(1)y(n)=b(1)x(n)+b(2)x(n－1)+…+b(nb+1)x(n－nb)

－a(2)y(n－1)－…－a(na+1)y(n-na)如果a(1)不等于1，濾波器對濾波器系數(shù)正則化；如果a(1)=0，返回錯誤。如果x是矩陣，那么該函數(shù)對矩陣的每一列進(jìn)行濾波；如果x是一個多維陣列，那么函數(shù)沿著x的第一個非單維陣列濾波。

（2）[y,zf]=filter(b,a,x)和[y,zf]=filter(b,a,x,zi)：

濾波后產(chǎn)生一個附加輸出zf，包含從0初始狀態(tài)計算得到的最后狀態(tài)向量。（3）y=filter(Hq,x)：

用濾波器Hq對輸入數(shù)據(jù)x進(jìn)行濾波，得到濾波后的輸出數(shù)據(jù)y，向量x和y具有相同的長度。如果x是矩陣，那么該函數(shù)對矩陣的每一列進(jìn)行濾波；如果x是一個多維陣列，

那么函數(shù)沿著x的第一個非單維陣列濾波。

（4）［y,zf］=filter(Hq,x)：

濾波后產(chǎn)生一個附加輸出zf，包含從0初始狀態(tài)計算得到的最后狀態(tài)向量。

3)fvtool函數(shù)

fvtool函數(shù)可用來分析數(shù)字濾波器的圖形化用戶界面，它也有多種調(diào)用格式。

（1）fvtool(B,A)：裝載濾波器可視化工具，并計算濾波器的幅頻響應(yīng)，公式如下：分子和分母的系數(shù)分別在向量B和A中。（2）fvtool(B,A,B1,A1,…)：實現(xiàn)多濾波器分析，如果用來實現(xiàn)有限沖擊響應(yīng)濾波器的分析，把Ai(i=1,2,…)設(shè)為1即可。如果只對一個有限沖擊響應(yīng)濾波器進(jìn)行分析，可直接寫作fvtool(h)，h為FIR濾波器的系數(shù)向量。數(shù)字濾波器的圖形化用戶界面如圖11-6所示，點擊

圖標(biāo)，只顯示濾波器的幅頻特性曲線，點擊只顯示相頻特性，點擊則同時顯示濾波器的幅頻特性曲線和相頻特性曲線。點擊顯示濾波器的沖擊響應(yīng)，點擊則顯示濾波器的階越響應(yīng)，點擊顯示濾波器的零點和極點，點擊則可查看濾波器的系數(shù)。圖11-6數(shù)字濾波器的圖形化用戶界面從圖11-6中還可以看出，相位曲線成鋸齒波狀，這是將相位歸算到0~－360°的結(jié)果。從左到右，每多一個鋸齒，其實際相位是顯示的相位多增加一次360°，即從左邊開始，第i(i=1,2,3,…)個鋸齒中顯示的相位為θ，其實際相位為－(i－1)×360°+θ。例如，若第二個鋸齒中顯示

相位為－180°，其實際滯后的相位為－(2－1)×360°－180°=－540°。另外，除了第一鋸齒以外，其它鋸齒的最大值往往不是0，最小值也不是－360°，這是由于濾波器長度有限所致。上一個鋸齒的結(jié)束點和下一個鋸齒的起始點，這兩個相鄰離散化相位之間有一定的間隔，這個間隔分散在上一個鋸齒的結(jié)束處和下一個鋸齒的起始處。濾波器長度越大，則間隔越小，鋸齒的起始點越接近0，結(jié)束點越接近－360°。

2.應(yīng)用實例

程序如下：

%生成方波

t=0:0.01:20;

forj=1:length(t)

fori=1:10

ift(j)<=2*i&t(j)>2*i-1

a(j)=1+0.2*(rand-0.5);%producerandom

noisewhoseamplitudeis0.2

end

ift(j)<=2*i－1&t(j)>2*i-2

a(j)=－1+0.2*(rand－0.5);

end

end

end

%設(shè)定濾波器參數(shù)，獲得線性相位濾波器系數(shù)

n=19;%設(shè)定濾波器階次為19

wo=0.1;%設(shè)定截止頻率為0.1

del=［0.001,0.01］;%指定通帶與阻帶的最大允許誤差分別為0.001,0.01

h=firceqrip(n,wo,del);

%得到的濾波器系數(shù)為：

［－0.00150.00910.01580.02690.04070.0561

0.07160.08520.09550.10100.10100.0955

0.08520.07160.05610.04070.02690.0158

0.0091－0.0015］;

%卷積實現(xiàn)濾波

fori=1:length(t)

c=0;

forj=1:n+1

ifi>j

c=c+h(j)*a(i-j+1);

end

end

b(i)=c;%b為濾波后得到的信號

end

%直接調(diào)用函數(shù)filter實現(xiàn)濾波

y=filter(h,1,a);

%顯示信號

figure(1)

fvtool(h)%顯示濾波器的幅頻特性和相頻特性

figure(2)

subplot(4,1,1)%在同一張圖中設(shè)定4行1列的4個子圖，并指定第一行的子圖

plot(t,a,′r′)

%在第一行的子圖中用紅色線條′r′指定紅色)繪出方波a的波形，橫坐標(biāo)為時間軸t

subplot(4,1,2)

plot(t,y,′b′)%在第二行的子圖中用藍(lán)色線條(′b′指定藍(lán)色)繪出用filter函數(shù)濾波后的信號y的波形，橫坐標(biāo)為時間軸t

subplot(4,1,3)

plot(t,b,′c′)%在第三行的子圖中用青色線條(′c′指定青色)繪出自行編程濾波后的信號b的波形，橫坐標(biāo)為時間軸t

subplot(4,1,4)

plot(t,b-y,′k′)%在第四行的子圖中用黑色線條(′k′指定黑色)繪出兩種濾波方法得到的信號的差異，橫坐標(biāo)為時間軸t

最后得到的濾波器幅頻特性(藍(lán)線)與相頻特性(綠線)結(jié)果如圖11-6所示。從圖中可以看出，在頻率高于0.27以后，即在阻帶區(qū)，濾波器具有等紋波，而在通帶內(nèi)，濾波器具有線性相位。初始方波信號、濾波后的信號以及兩種實現(xiàn)濾波計算得到的濾波信號的差異如圖11-7所示。對比圖11-7(a)與(b)和(c)可以看出，濾波后的信號具有一定程度的信號滯后。

而由圖11-7(d)可以看出，兩種實現(xiàn)濾波計算的方法得到的信號的差異非常小，其幅值不及初始信號幅值的1/1015，這說明filter函數(shù)所完成的工作實質(zhì)上就是濾波器系數(shù)與初始信號的卷積過程。兩種實現(xiàn)濾波計算方法得到的結(jié)果之間的差異是由于計算的截斷誤差所帶來的。圖11-7濾波結(jié)果(a)初始方波信號；(b)用filter函數(shù)計算得到的濾波信號；(c)自行編程計算得到的濾波信號；(d)兩種濾波信號的差異11.2自適應(yīng)濾波器

11.2.1自適應(yīng)濾波器的結(jié)構(gòu)

圖11-8示出了FIR自適應(yīng)濾波器的一般結(jié)構(gòu)。圖中的x(n)為濾波器輸入信號，也就是初始的測試信號，y(n)為濾波器輸出信號，d(n)為參考信號或期望信號，e(n)則是d(n)和y(n)的誤差信號。自適應(yīng)濾波器的濾波器系數(shù)受誤差信號e(n)控制，根據(jù)e(n)的值和自適應(yīng)算法自動調(diào)整。圖11-8FIR自適應(yīng)濾波器結(jié)構(gòu)11.2.2自適應(yīng)濾波理論與算法

對于一個含有噪聲的觀測信號，它可以表示為

x(n)=s(n)+v(n)(11-51)

式中s(n)為信號的真值，v(n)為噪聲信號的值，x(n)為含有噪聲的觀測值。即

與經(jīng)典濾波的區(qū)別在于：經(jīng)典濾波的前提是把v(n)看做頻率可預(yù)測的，而且是一般不會變換的，其幅值往往是有一定范圍的，而自適應(yīng)濾波則把v(n)看做未知的，其頻率、幅值都是時變的。在自適應(yīng)濾波中，最廣泛采用的目標(biāo)函數(shù)之一是均方誤差(MSE)，其定義為

F(e(n))=ξ(n)=E(e2(n))

=E(d2(n)－2d(n)y(n)+y2(n))(11-52)

假設(shè)自適應(yīng)濾波器是由線性組合器構(gòu)成的，即輸出信號是由來自于陣列的信號的線性組合構(gòu)成，在許多應(yīng)用中，輸入信號向量的每個元素是由同一個信號的時延形式構(gòu)成的，此時，信號y(n)是輸入信號經(jīng)過FIR濾波的結(jié)果：(11-53)式中:x(n)=［x(n),x(n－1),…,x(n－N)］T，w(n)=［w0(n),w1(n),…,wN(n)］T。此時，E(e2(n))=ξ(n)=E［d2(n)－2d(n)wT(n)x(n)+wT(n)x(n)xT(n)w(n)］=E［d2(n)］－2E［d(n)wT(n)x(n)］+E［wT(n)x(n)xT(n)w(n)］(11-54)對于固定系數(shù)的濾波器而言，MSE函數(shù)為

ξ(n)=E［d2(n)］－2wTE[d(n)x(n)］+wTE[x(n)xT(n)]w

=E［d2(n)］－2wTP+wTRw(11-55)

式中P=E[d(n)x(n)]，是輸入信號與期望信號之間的互相關(guān)向量；R=E［x(n)xT(n)］，是輸入信號的自相關(guān)矩陣。MSE函數(shù)對向量w的梯度向量為(11-56)如果P向量和矩陣R已知，令該梯度向量為0，可以直接求解w：

w0=R－1P(11-57)

上式中的解稱為維納解。遺憾的是，P向量和矩陣R很難

精確估計，當(dāng)輸入和期望信號都是遍歷性過程時，能利用時間平均估計R和P，大多數(shù)自適應(yīng)算法都隱含地利用了這一點。

1.最小均方算法

最小均方算法(LeastMeanSquare,LMS)是由Widrow和Hoff(1960)開發(fā)并命名的，它是隨機(jī)梯度算法族中的一員。LMS算法的一個顯著特點是它的簡單性，此外，它不需要計算有關(guān)的相關(guān)函數(shù)，也不需要矩陣求逆運算，這使得它成為其它線性自適應(yīng)濾波算法的參照標(biāo)準(zhǔn)。

LMS算法包含兩個過程：

(1)濾波過程,包括計算線性濾波器輸出對輸入信號的響應(yīng)以及通過比較輸出結(jié)果與期望響應(yīng)產(chǎn)生的估計誤差。（2)自適應(yīng)過程，即根據(jù)估計誤差自動調(diào)整濾波器

參數(shù)。

這兩個過程一起工作組成一個反饋環(huán)。首先，我們有一個用來完成濾波過程的橫向濾波器；其次，我們有一個調(diào)節(jié)濾波器系數(shù)的自適應(yīng)控制算法。

對于式(11-55)，如果可以得到P向量和矩陣R的較好的估計，分別記為和，則可以利用如下最陡下降算法搜索維納解：(11-58)式中表示目標(biāo)函數(shù)對濾波器系數(shù)的梯度向量估計值。一種可能的解是利用R和P的瞬態(tài)估計值來估計梯度向量，即(11-59)(11-60)于是得到的梯度估計值為(11-61)如果用瞬態(tài)平方誤差e2(n)代替MSE作為目標(biāo)函數(shù)，梯度估計值則代表了真實的梯度向量，因為(11-62)由于得到的梯度算法使得平方誤差的均值最小化，因此被稱為LMS算法，其更新方程為

w(n+1)=w(n)+2μe(n)x(n)(11-63)式中，收斂因子μ應(yīng)該在一個小的范圍內(nèi)取值，以保證算法的收斂性。

2.RLS算法

最小二乘(LeastSquare,LS)算法旨在使期望信號與模型濾波器輸出之差的平方和達(dá)到最小。當(dāng)在每次迭代中接收到輸入信號的新采樣值時，可采用遞歸形式求解最小二乘問題，得到遞歸最小二乘(RecursiveLeastSquare,RLS)算法。RLS算法能實現(xiàn)快速收斂，即使是在輸入信號相關(guān)矩陣的特征值擴(kuò)展比較大的情況下。對于最小二乘算法，目標(biāo)函數(shù)是確定性的，并且由下式給出：(11-64)式中λ為指數(shù)加權(quán)因子或遺忘因子，其值應(yīng)選擇在［0,1］范圍內(nèi)。應(yīng)該注意的是，在推導(dǎo)LMS算法和基于RLS準(zhǔn)則的算法時，利用了先念誤差。在RLS算法中，用e(n)表示后念誤差，e′(n)表示先念誤差，在推導(dǎo)基于RLS準(zhǔn)則算法的過程中，將首先選擇后念誤差。

每一個誤差都是由期望信號和采用最近的系數(shù)w(n)得到的濾波器輸出之差所組成的。將ξd(n)相對于w(n)求導(dǎo)，可以得到：(11-65)令式(11-65)等于0，則可以通過如下關(guān)系式找到使得最小二乘誤差最小的最優(yōu)向量w(n)：(11-66)從而得到最優(yōu)系數(shù)向量w(n)的表達(dá)式為(11-67)式中:稱為輸入信號與期望信號的確定性互相關(guān)向量；(11-68)稱為輸入信號的確定性相關(guān)矩陣。為了避免求矩陣RD(n)的逆，我們可以利用如下的矩陣求逆引理：

[A+BCD]－1=A－1－A－1B[DA－1B+C－1]－1DA－1(11-69)

式中A、

B、C和D是具有合適維數(shù)的矩陣，并且矩陣A和C是非奇異矩陣。于是由矩陣求逆引理可得：(11-70)利用式(11-70)來計算，與每次迭代過程中直接計算的逆矩陣(乘法次數(shù)在N3數(shù)量級)相比，更新方程的計算復(fù)雜度更低(乘法次數(shù)在N2數(shù)量級)。當(dāng)n=0時，只要x(0)≠0，則有：當(dāng)n=0時，有：，當(dāng)n=0時，有：，利用歸納法，可以證明在任意時刻n，有：(11-71)這樣，在計算初始的N+1個濾波器系數(shù)時，采用式(11-71)即可，不需要計算式(11-70)。11.2.3Matlab中的自適應(yīng)濾波函數(shù)

1.s=initlms(w0,mu,zi,lf)

LMS算法的矩陣形式是：

w(k+1)=w(k)+μe(k)x(k)

為了確保良好的收斂速度和穩(wěn)定性，μ(mu)一般選擇在如下范圍內(nèi)：式中N為信號中樣本的數(shù)量。輸入?yún)?shù)zi指定濾波器初始條件，若忽略zi或者指定為空，initkalman默認(rèn)zi為0向量，其長度等于length(w0)－1。對于限定的處理條件，例如for循環(huán)中使用adaptlms，指定初始條件非常重要，LMS算法的每一次迭代都使用上一次迭代的權(quán)值，由于提供了初始條件，因此第一次迭代具有一組先驗濾波器權(quán)值，所以開始進(jìn)行迭代運算。輸入?yún)?shù)lf為泄漏因子。指定泄漏因子可以改善算法的特性。泄漏權(quán)值w(k)強制算法進(jìn)一步適應(yīng)運算，就算已經(jīng)達(dá)到了最小值。這可能意味著泄漏LMS沒有獲得非常精確的最小均方誤差的測度時，使用泄漏因子算法可以減小誤差的敏感性，或者對輸入的小數(shù)值的敏感性。典型地，lf設(shè)置為0.9至1.0之間，表示沒有泄漏。如果指定lf為空，它默認(rèn)為1。

在使用initlms函數(shù)后，如果核對s的內(nèi)容，Matlab會顯示結(jié)構(gòu)元素，而不是輸入?yún)?shù)的名稱。為了幫助記住各元素對應(yīng)的initlms的輸入?yún)?shù)，表11-2給出了該映射關(guān)系。

2.［y,e,s］=adaptlms(x,d,s)

該函數(shù)的功能是：用FIR最小均方自適應(yīng)濾波器作用于x和期望信號d，濾波后的信號返回給y。s是一個結(jié)構(gòu)體，它包含定義所采用的LMS的初始化設(shè)置，這同一些濾波器適應(yīng)過程的輸出一樣。表11-2詳細(xì)列出了