零相位數(shù)字濾波器在非平穩(wěn)信號處理中的應(yīng)用

零相位數(shù)字濾波器在非平穩(wěn)信號處理中的應(yīng)用

數(shù)字濾波是數(shù)字數(shù)據(jù)處理的一般方法。普通數(shù)字濾波在濾波過程中會發(fā)生一定的相位位移。為了解決這個問題，引入了零相位數(shù)字濾波。兩個主要的零相位數(shù)字濾波算法（之前的序列、欺詐和結(jié)果序列）和rr（不安全、不安全、結(jié)果序列和結(jié)果序列）。首先，我們研究了采用合適的濾波初期模型、信號延伸度、frr和rr的方法，改善了零相位數(shù)字濾波的邊界問題，并對穩(wěn)定信號濾波產(chǎn)生了良好的效果。在對非平穩(wěn)信號進行濾波時,與小波分解重構(gòu)和經(jīng)驗?zāi)B(tài)分解等無相移的非平穩(wěn)信號處理方法相比,零相位數(shù)字濾波器具有截止頻率明確、計算量小的優(yōu)點.尤其在計算能力較弱的嵌入式應(yīng)用系統(tǒng)中,零相位數(shù)字濾波器的這一優(yōu)點更加突出.調(diào)幅信號、調(diào)頻信號是兩類典型的非平穩(wěn)信號,很多工程實際信號可以借助這兩種信號進行模擬.例如,電力系統(tǒng)的次同步振蕩、電壓暫態(tài)振蕩、開關(guān)設(shè)備動作時的振動、發(fā)電機組軸系的扭振等就可用調(diào)幅信號模擬;旋轉(zhuǎn)機械調(diào)速過程中的轉(zhuǎn)速則多為調(diào)頻信號.文獻利用線性調(diào)頻信號對零相位數(shù)字濾波器進行了分析.因此,本文作者將主要針對調(diào)幅信號來討論零相位數(shù)字濾波器在非平穩(wěn)信號處理中的特點.利用調(diào)幅信號對一種典型的帶通零相位數(shù)字濾波器進行了仿真分析;對小波包分解重構(gòu)、經(jīng)驗?zāi)B(tài)分解及零相位數(shù)字濾波器方法的濾波效果進行了比較;探討了在調(diào)幅信號處理中存在的過渡過程,及其對零相位數(shù)字濾波的幅值特性以及頻率特性的影響;論述了影響濾波誤差的因素;提出了利用分段零相位濾波的方法改善零相位數(shù)字濾波器過渡過程的方法;用實測的次同步振蕩信號對兩種方法進行了驗證.1零相位濾波、ir數(shù)字濾波器對于一維時間序列,零相位數(shù)字濾波器的基本原理是,利用正向時間序列和翻轉(zhuǎn)時間序列通過濾波器時的相移相互抵消,從而實現(xiàn)濾波結(jié)果的零相移.利用傅里葉變換的時域位移性質(zhì)和時間翻轉(zhuǎn)性質(zhì),即可從理論上論證其正確性.以FRR方法為例,將長度為N+1的單邊時間序列x[n]通過沖激響應(yīng)為h[n]的濾波器,得到一次濾波序列y1[n],再將一次濾波序列進行時序反轉(zhuǎn)、平移得到數(shù)據(jù)序列y2[n],y2[n]通過濾波器h[n]得到二次濾波序列y3[n],最后將二次濾波序列進行時序反轉(zhuǎn)、平移得到最終濾波輸出y[n].從而,有y1[n]=x[n]*h[n](1)y2[n]=y1[Ν-n](2)y3[n]=y2[n]*h[n](3)y[n]=y3[Ν-n](4)設(shè)x[n]?X(ω)、y1[n]?Y1(ω)、y2[n]?Y2(ω)、y3[n]?Y3(ω)、y[n]?Y(ω)、h[n]?H(ω)構(gòu)成傅里葉變換對,且考慮到對實序列h[n]有H(-ω)=H*(ω).依式(1～4),FRR濾波輸出的頻域表達如下式推導Y(ω)=e-jωΝY3(-ω)=e-jωΝY2(-ω)Η(-ω)=e-jωΝejωΝY1(ω)Η(-ω)=Y1(ω)Η(-ω)=X(ω)Η(ω)Η(-ω)=X(ω)|Η(ω)|2(5)因而,FRR方法實現(xiàn)了零相移濾波,其傳遞函數(shù)頻域表示如下ΗFRR(ω)=Η(ω)Η(-ω)(6)同理可得,RRF方法實現(xiàn)零相移濾波時,其傳遞函數(shù)頻域表示如下ΗRRF(ω)=Η(-ω)Η(ω)(7)選定沖激響應(yīng)序列為h[n]的普通數(shù)字濾波器,實際使用時可以通過兩個分立的濾波過程來實現(xiàn)零相位濾波,也可以直接構(gòu)造零相位數(shù)字濾波器,其沖激響應(yīng)序列分別為hFRR[n]和hRRF[n].hFRR[n]=h[n]*h[-n](8)hRRF[n]=h[-n]*h[n](9)其中h[-n]為h[n]的時間翻轉(zhuǎn)序列.h[n]的選擇可有FIR濾波器和IIR濾波器兩種類型.在零相位濾波器的應(yīng)用中,IIR濾波器在相頻特性上的劣勢可以忽略,且在相同階數(shù)下IIR濾波器比FIR濾波器具有更好的截止特性,因而h[n]一般選用IIR數(shù)字濾波器.23khs.自適應(yīng)變化為比較小波包分解重構(gòu)、經(jīng)驗?zāi)B(tài)分解和零相位數(shù)字濾波器3種方法的差異.選取含有3個衰減振蕩分量的輸入信號f(t),信號采樣率為1kHz.設(shè)定3個分量的振蕩頻率分別為10Hz、23Hz、39Hz,其計算公式如下f1(t)=e-3(t-1)u(t-1)sin(20π(t-1))(10)f2(t)=e-2(t-1)u(t-1)sin(46π(t-1))(11)f3(t)=e-(t-1)u(t-1)sin(78π(t-1))(12)輸入信號f(t)的計算公式如下f(t)=f1(t)+f2(t)+f3(t)(13)1小波包的濾波、并存在重要差異利用小波包可將信號分解到不同節(jié)點,各節(jié)點對應(yīng)了不同頻帶的子波,對目標節(jié)點系數(shù)進行重構(gòu)即可得到所需的濾波輸出.考慮到重構(gòu)信號的不失真和光滑,采用了bior、sym和db小波基對多分量信號進行處理,最終選定db8小波基對f(t)進行5層分解重構(gòu),由于信號采樣率為1kHz,小波包5層分解的各子波頻帶寬度為15.625Hz,使10Hz、23Hz、39Hz3個頻率成分位于、、節(jié)點對應(yīng)的頻帶內(nèi),濾波結(jié)果如圖1所示(圖中虛線為濾波輸出,實線為輸入分量).小波包分解重構(gòu)的濾波輸出,突變點附近存在過渡過程,其余部分信號不光滑,在幅值和相位上都存在抖動,其中23Hz和39Hz分量抖動比較嚴重.2mt經(jīng)驗?zāi)B(tài)分解mt模式,m經(jīng)驗?zāi)B(tài)分解方法是自適應(yīng)的將信號分解為不同時間尺度的本征模態(tài)函數(shù)(IMF)分量和趨勢項,本征模態(tài)函數(shù)是一類零均值的調(diào)幅調(diào)頻信號.對含隨機寬帶噪聲的信號,經(jīng)驗?zāi)B(tài)分解是類小波的二進濾波器組,第一個本征模態(tài)函數(shù)的頻帶最高.由于m(t)由3個調(diào)幅正弦信號組成,經(jīng)驗?zāi)B(tài)分解得到的前3個本征模態(tài)函數(shù)分量即對應(yīng)其濾波輸出.其中,第一個本征模態(tài)函數(shù)對應(yīng)39Hz分量、第二個本征模態(tài)函數(shù)對應(yīng)23Hz分量、第三個本征模態(tài)函數(shù)對應(yīng)10Hz分量,濾波結(jié)果如圖2所示(圖中虛線為濾波輸出,實線為輸入分量).經(jīng)驗?zāi)B(tài)分解方法的濾波輸出,突變點附近存在過渡過程,但其余部分信號很光滑,幅值和相位的穩(wěn)定性強于小波包分解重構(gòu)方法,但是仍然存在抖動.我們又采用標準差為0.005、聚合次數(shù)為50次的聚合經(jīng)驗?zāi)B(tài)分解對f(t)進行濾波,發(fā)現(xiàn)能改善過渡過程.但是,同時造成了39Hz分量混疊到23Hz分量中,為減少混疊效應(yīng)需要增加聚合次數(shù),極大的增加了計算量.3過渡過程中的誤差檢驗數(shù)字濾波器在對有限長信號進行濾波時,相當于給信號加了一個矩形窗,在濾波輸出信號上造成衰減振蕩.零相位數(shù)字濾波器進行了兩次濾波,在輸出信號的首尾兩端形成失真.FRR方法第二次濾波的首端是輸入信號的尾端,濾波輸出的首端失真較小;RRF方法第二次濾波的首端是輸入信號的首端,濾波輸出的尾端失真較小.對于本文考慮的調(diào)幅正弦信號,通過在輸入信號首尾兩端進行信號延拓的方法,來消除端點效應(yīng).依據(jù)對端點效應(yīng)誤差的要求,延拓的數(shù)據(jù)長度應(yīng)比擬于過渡過程時間,本文取1s長度的數(shù)據(jù)點.消除端點效應(yīng)的FRR方法和RRF,在使用中的差異可忽略不計.本文的后續(xù)討論中,零相位數(shù)字濾波器均采用FRR方法實現(xiàn),且采用了首尾端的端點延拓,不再特別說明.綜合考慮濾波器階躍響應(yīng)的過渡過程時間和振蕩次數(shù)等因素,選定用4階巴特沃斯帶通濾波器(直接型系數(shù)見表1)構(gòu)造零相位數(shù)字濾波器,中心頻率分別為10Hz、23Hz、39Hz,帶寬為3Hz,記為hbp10、hbp23、hbp39.濾波結(jié)果見圖3(圖中虛線為濾波輸出,實線為輸入分量).濾波輸出很光滑,除了突變點附近,幅值和相位與原始信號重合,效果好于小波包分解重構(gòu)和經(jīng)驗?zāi)B(tài)分解方法.不考慮信號的突變,則零相位數(shù)字濾波器實現(xiàn)了較嚴格的零相移濾波.3種方法都存在過渡過程.以10Hz分量的誤差為例來考察過渡過程,將3種方法所得結(jié)果前3s的誤差示于圖4中.從上至下依次為小波包分解重構(gòu)、經(jīng)驗?zāi)B(tài)分解、零相位濾波器誤差,實線為原始誤差,虛線為幅值放大100倍后的誤差.其中,經(jīng)驗?zāi)B(tài)分解誤差極值最大,零相位濾波器與小波包分解重構(gòu)的誤差極值較小;零相位濾波與經(jīng)驗?zāi)B(tài)分解誤差穩(wěn)定時間較短,小波包分解重構(gòu)的誤差穩(wěn)定時間較長;但零相位濾波器反時序的特性造成突變點前有信號的提前輸出.下面我們將對零相位濾波器的過渡過程詳細探討,并研究其幅值和頻率特性,從而對其過渡過程性能進行改善.3過渡過程3.1同上轉(zhuǎn)臺濾波過濾算法零相位數(shù)字濾波器解決了濾波中存在的相移問題,但是在對非平穩(wěn)信號進行濾波時,仍需注意濾波中可能存在的過渡過程問題.尤其是在信號發(fā)生突變的位置影響較大,這是因為突變點的存在類似于將信號在該點截斷,引入了新的端點.考察一個在第2秒時刻幅值由0躍變?yōu)?在第8秒再跳變?yōu)?的10Hz正弦信號x(t),其計算公式如下x(t)=sin(20π(t-2))[u(t-2)-u(t-8)](14)其中u(t)為單位階躍信號.將x(t)通過hbp10進行濾波,濾波結(jié)果如圖5所示,存在過渡過程.為了更清楚地表達過渡過程,將其幅值響應(yīng)曲線示于圖6.輸入信號的突變在濾波輸出的對應(yīng)位置上存在過渡過程,穩(wěn)定到±2%的時間約0.48s,最大超調(diào)量為3.4%;進入平穩(wěn)階段后,濾波輸出與輸入信號能很好的吻合,誤差均值在10-5量級,而誤差的方差在10-7量級,平穩(wěn)階段時間越長,誤差越小.對這類幅值發(fā)生階躍的正弦信號,可以采用濾波輸出信號最大幅值的中值時刻來近似定位突變點,從而去除突變點前的虛假信號.采用該方法估計x(t)的突變點的估計誤差只有0.005s,實際應(yīng)用中受波形形態(tài)及噪聲的影響,誤差變大.事實上,不同特點的信號通過同一零相位數(shù)字濾波器,濾波信號幅值和頻率表現(xiàn)出的特性也會不同.為進行具體的仿真分析,選取衰減振蕩、發(fā)散振蕩、發(fā)散振蕩經(jīng)采取措施后又衰減振蕩這3種情況進行研究,利用幅值衰減、幅值發(fā)散、幅值先發(fā)散后衰減這3種不同調(diào)幅正弦信號進行模擬.3.23濾波器hbp10濾波特性先構(gòu)造調(diào)幅正弦信號的幅值包絡(luò)a1(t)、a2(t)、a3(t),其中a1(t)為衰減幅值、a2(t)為發(fā)散幅值、a3(t)為先發(fā)散后衰減幅值,其計算公式如下a1(t)=e-(t-1)u(t-1)(15)a2(t)=e(t-9)[1-u(t-9)](16)a3(t)=e(t-5)[1-u(t-5)]+e-(t-5)u(t-5)(17)對應(yīng)的調(diào)幅正弦信號分別為m1(t)、m2(t)、m3(t),其計算公式如下m1(t)=a1(t)sin(20πt)(18)m2(t)=a2(t)sin(20πt)(19)m3(t)=a3(t)sin(20πt)(20)采用一個高斯白噪聲信號wn(t)作為噪聲干擾信號,wn(t)的統(tǒng)計分布為N(0,0.01).對應(yīng)的含噪聲的輸入信號分別為s1(t)、s2(t)、s3(t),其表達式分別為下列公式所示s1(t)=m1(t)+wn(t)(21)s2(t)=m2(t)+wn(t)(22)s3(t)=m3(t)+wn(t)(23)s1(t)、s2(t)、s3(t)的波形如圖7所示,經(jīng)過零相位帶通數(shù)字濾波器hbp10濾波后,濾波輸出信號如圖8所示.可明顯看出,這3種信號的濾波輸出的最大幅值小于的x(t)的濾波輸出的最大幅值;過渡過程時間也小于x(t)的過渡過程時間;其中s3(t)的濾波結(jié)果最好,誤差均值和方差都在10-5量級.可見輸入信號的形態(tài)特點對濾波結(jié)果有影響.4振幅值和頻率特性4.1濾波結(jié)果的頻率方差濾波輸出的幅值和頻率可利用希爾伯特變換求得.其基本原理如下.對信號ν(t),其希爾伯特變換由下式定義ν?(t)=1π∫-∞+∞ν(τ)t-τdτ(24)對應(yīng)的幅值包絡(luò)a(t)可按下式計算a(t)=ν2(t)+ν?2(t)(25)對應(yīng)的相位θ(t)可按下式計算θ(t)=arctanν?(t)ν(t)(26)對應(yīng)的瞬時頻率定義為f(t)=12πdθ(t)d(t)(27)下面利用希爾伯特變換獲取濾波結(jié)果的幅值包絡(luò)和頻率,對零相位數(shù)字濾波器的幅值響應(yīng)和頻率響應(yīng)特點進行詳細分析.4.2放大電路一般特性s1(t)、s2(t)、s3(t)經(jīng)過零相位帶通數(shù)字濾波器hbp10濾波后,理想的幅值應(yīng)為a1(t)、a2(t)、a3(t).實際按式(25)計算出的濾波輸出信號幅值分別如圖9～圖11所示.可明顯看出,s1(t)、s2(t)、s3(t)的濾波輸出的最大幅值分別只有0.81、0.80和0.93,明顯小于的x(t)的濾波輸出的最大幅值1.034;s3(t)的最大值沒有明顯延遲,s1(t)、s2(t)過渡過程時間為0.22s、0.22s,也明顯小于對x(t)濾波時的過渡過程時間0.48s;信號幅值較小的部分噪聲影響比較大.由于不像x(t)有穩(wěn)定的信號階段可以使濾波輸出完成過渡達到穩(wěn)定,上述3類信號衰減系數(shù)越大幅值損失越大.將不同衰減系數(shù)的10Hz分量濾波輸出的幅值示于圖12中,其結(jié)果印證了該結(jié)論.對于零相位帶通數(shù)字濾波器而言,與普通帶通數(shù)字濾波器相似,其幅值響應(yīng)具有如下特點:1)濾波器類型一定時,較寬的頻帶可以獲得較大的幅值和較小的時延;提高濾波器階數(shù)可加速阻頻帶衰減,減少臨近頻帶的干擾,但也會帶來通帶幅值的振蕩.2)輸入信號幅值包絡(luò)變化時,濾波輸出的時延從幅值包絡(luò)變化率大的部分偏向變化率小的部分;幅值包絡(luò)變化越平緩,濾波后的時延也越小;濾波器階數(shù)大于等于4時,信號幅值包絡(luò)越接近階躍,超調(diào)量越大.3)因為白噪聲無法完全濾除,幅值較小的部分噪聲影響比較明顯.4.3濾波、濾波效率及信號背景分析按式(27)式定義對m1(t)、m2(t)、m3(t)的濾波輸出信號進行計算,得到的瞬時頻率分別示于圖13～圖15中.其中,首尾兩端的振蕩是由于希爾伯特變換的邊界效應(yīng)造成的.另外需要注意的地方是,在信號突變點附近(約為一個過渡過程時間),濾波輸出存在較大的頻率突變,說明濾波輸出信號在該位置存在畸變,在圖13～圖14中分別位于0.755s、9.246s處.取中間頻率穩(wěn)定的3～8s時間段進行分析,可得頻率均值為10.000Hz,方差小于10-5量級.其中,m3(t)在中間突變點處的頻率處于9.99Hz與10.01Hz之間,明顯好于m1(t)、m2(t)在突變點處的濾波結(jié)果.為考察噪聲對頻率響應(yīng)的影響,取s1(t)信號為例.圖16給出了s1(t)經(jīng)濾波后的瞬時頻率曲線,可見濾波結(jié)果變差,信號幅值較小的部位上噪聲造成相位的波動,使波形相對發(fā)生了較大的畸變.仍取3～8s時間段進行分析,可得頻率均值為9.962Hz,方差在10-1量級.零相位帶通數(shù)字濾波器的頻率響應(yīng)具有如下特點:1)中心頻率和頻帶寬度變化對穩(wěn)定的頻率段影響很小,主要對信號突變點附近的頻率有影響.2)當輸入信號的幅值包絡(luò)變化時,濾波輸出的頻率畸變點位置與幅值時延的方向相反,位于突變點偏幅值包絡(luò)變化率較大的一側(cè).3)噪聲的存在,造成信號的畸變,濾波輸出頻率存在波動,幅值較小的部分噪聲影響比較明顯.5分段零相位濾波由上述仿真結(jié)果可知,m3(t)信號的濾波結(jié)果明顯好于m1(t)與m2(t).m3(t)與m1(t)、m2(t)的不同之處在于,其突變點兩邊的幅值是對稱的.基于此,作者提出分段零相位濾波方法,用以改善零相位濾波器在處理存在突變的調(diào)幅信號時的過渡過程性能,其具體步驟如下:1)對信號進行零相位濾波.2)在濾波幅值包絡(luò)的中值點附近,搜索濾波輸出中絕對值最小的點,進行突變點的近似定位.3)將輸入信號按定位結(jié)果分段,對分段后的信號進行端點延拓,再分別進行零相位濾波.4)將分段濾波結(jié)果組合,形成濾波輸出.零相位濾波器采用hbp10,利用零相位濾波和分段零相位濾波方法分別對s1(t)進行濾波,濾波結(jié)果示于圖17中.從圖中可以明顯觀察到,分段零相位濾波的結(jié)果不但比原始的零相位濾波結(jié)果過渡過程縮短,而且也解決了反時序提前輸出的問題.可見,分段零相位濾波能縮短過渡過程時間,消除反時序提前輸出,極大改善零相位濾波器的濾波效果.其缺點是,增加了濾波次數(shù)和端點延拓次數(shù),增加了計算量.6零相位數(shù)字濾波器和分段零相位數(shù)字濾波器為從次同步振蕩中分離出主要的振蕩分量,將該方法應(yīng)用到次同步振蕩信號的濾波中.在某次試驗中測得振蕩信號,采樣率為1kHz,數(shù)據(jù)長度為95s,如圖18所示.通過時頻分析知該振蕩信號中主要包含有13Hz、50Hz、25Hz等幾個較大的頻率成分,其中在次同步頻域內(nèi)13Hz振幅最大.以13Hz分量為例,仍采用4階巴特沃斯帶通濾波器(直接型的中心頻率為13Hz,帶寬3Hz),構(gòu)造零相位