多相濾波器組譜分析方法的性能討論與天文觀測應用*

朱凱，甘恒謙，朱巖，金乘進

(中國科學院國家天文臺，北京 100012)

對天體輻射信號的功率譜分析可以獲得射電源輻射機制、運動性質等諸多信息。在射電天文數字終端設備中，對信號進行譜分析主要采用自相關或周期圖方法［1］;另一方面，以多相結構實現的數字濾波器組能夠實現信號的均勻子帶分解，在數據流與計算量上相對于濾波器組的直接實現有很大的降低，因此應用多相濾波器組［2］進行譜分析，也是一種可選擇的方案，如GALFA中性氫巡天項目所用的頻譜儀即基于多相濾波器組進行研制。本文擬對這種方法的原理進行分析，并從仿真和實測上與通常的周期圖方法性能進行比較，以得到有參考意義的結果。

本文內容組織如下，第1節簡要介紹多相濾波器組的原理，第2節對PFB譜估計器的統計性能和計算量進行討論，并與具有相同通道數及相同單位樣值響應長度(即原型FIR濾波器長度)的周期圖方法分別進行比較，第3、4節介紹仿真與實測的結果。

1 PFB原理簡述

多相濾波器組是數字濾波器組的一種高效實現形式［3］，用于對信號進行均勻子帶分解。所謂均勻子帶分解，是指濾波器組中的每一個子帶濾波器具有相同的通帶寬度，依次占據信號整個頻帶，從而將信號分解為若干個相同帶寬的子帶信號，其頻譜分解過程如圖1。

圖1 信號的均勻子帶分解示意Fig.1 Uniform sub band decomposition of a signal

其中h0(n)為低通濾波器的單位樣值響應，hk(n)為第k個子帶濾波器的單位樣值響應，由于各子帶濾波器均通過低通濾波器得到，因而此低通濾波器也被稱為原型低通濾波器。

通常，對每一個子帶信號，經帶通濾波之后都需要進行混頻、抽取處理。可以證明，對均勻劃分的多通道濾波器組，上述子帶處理過程可以用圖2中更為高效的多相結構形式實現［4］，也即多相濾波器組。

上圖所示的均勻子帶分解過程中，低通濾波器通帶為［－f1/2D，f1/2D］，第k個子帶濾波器的通帶為［(2k－1)f1/2D，(2k+1)f1/2D］，并且上述數字濾波器均為有限沖激響應(FIR)濾波器。由于濾波器組的通帶劃分具有均勻性，因此各子帶濾波器的頻率響應可以由其中的低通濾波器經頻譜搬移得到，此頻譜搬移過程可表示為:

圖2 多相濾波器組結構Fig.2 Structure of the polyphase filter bank

圖中抽樣時間T2=DT1，第ρ個子帶的濾波器單位樣值響應為:

第ρ個子帶輸入的信號變為x(nD－ρ)，是對原信號進行ρ點延遲和D倍抽取后的結果。

相對于直接實現的濾波、混頻、抽取處理，多相濾波器組將信號的抽取提前，同時對濾波器單位樣值響應也進行了抽取和重排列，子帶濾波器的單位樣值響應長度從KD下降為K，從而大大降低了資源量。

2 譜估計性能

隨機信號x(n)的功率譜定義為:

其中，r(k)是輸入信號x(n)的自相關函數，理論上信號的自相關函數需要通過無限長信號的計算得到，但是由于實際采集到的信號樣本長度總是有限的，因此需要利用有限樣本估計無限長信號的功率譜，即構成一個信號參數估計問題。

實際中應用最普遍的譜估計算法是基于FFT的周期圖方法，周期圖方法在大樣本量場合下是一個良好的譜估計器，將多相濾波器組用于功率譜估計，也構成一個譜估計器，下面對其性能進行討論，并與周期圖方法進行比較。

下面的討論中，如無特別說明，均設定多相濾波器組通道數為D，多相結構下的子濾波器長度(又稱多相因子)為K，對應的原型濾波器長度為N=KD。

對于譜估計方法性能的討論總是在有限樣本長度的條件下進行的，需要區分大樣本量與小樣本量情況。對于PFB方法，顯然在小樣本量(如L＜KD)的情況下是不適用的。此處僅考慮樣本數L大于原型濾波器長度N=KD的情況，設樣本數為N的整數倍，L=MN。輸入信號記為:

2.1 PFB譜估計原理

應用濾波器組進行信號功率譜估計的原理如圖3。

圖3 應用濾波器組進行功率譜估計Fig.3 Estimation of power spectrum with filter banks

圖中上部為信號處理的流程，下部為信號頻譜的變化。通過對帶通濾波后的信號進行功率計算，功率計算是一個積分平均的過程，由此得到該通帶內信號功率的估計值。設計一組帶通濾波器，使其通帶之和覆蓋整個信號帶寬，就可以估計出輸入信號的功率譜。

對信號進行帶通濾波和功率計算，實際上是對信號功率譜進行了平滑。而對于不同的帶通濾波器形式、以及不同的濾波器組實現方案，會具有不同的估計性能。

D通道PFB譜估計器首先將所有樣本數據有重疊的分成MK個子段，每一子段含KD個樣本數據，每一子段相對于上一子段新增D個樣本，分段方法如圖4。

圖4 D通道PFB譜估計器的樣本分段方式Fig.4 The segmentation of the D-channel PFB estimator

然后對上述的每一子段樣本依次進行D通道的子帶濾波，由PFB處理過程與圖2的等效性可知，PFB是采用單位樣值響應長度為KD的原型濾波器進行子帶濾波，因此對于每個子段KD個樣本輸入，每個頻率通道僅得到一個穩定的濾波輸出，在頻率ωk位置，第j子段的輸出值可表示為:

其中xj(m)=x［(j－1)N+m］，為第j子段的第m個樣本;h0(n)為原型低通濾波器的單位樣值響應。利用所有子段共MK個穩定輸出值進行功率計算，頻率ωk處功率計算的結果為所有子段輸出值的均值，可表示為:

在濾波器組通道數D確定的情況下，對信號功率譜的估計總是有偏差的。在樣本數趨于無限時，頻率ωk處功率譜的估計結果收斂到此通道處帶通濾波器的頻率響應與輸入信號真實功率譜乘積的積分，從而估計偏差也可由此得到。

對于PFB濾波器組，所有子帶濾波器的頻率響應都取決于原型低通濾波器的設計，而通常通過選擇合適的多相因子K和濾波器設計方法，可以獲得很好的濾波器頻響特性，因此相對于同通道數的周期圖方法其估計偏差更低，下面對此進行比較。

2.2 偏差比較

作為比較對象，考慮常見的周期圖(即FFT)譜估計方法，在相同的樣本量L=MN的條件下，分相同通道數和相同樣值響應長度兩種情況，與PFB譜估計法進行比較。

(1)相同通道數D

D通道Bartlett周期圖方法將數據無重疊的分成MK個子段，每個子段含D個樣本，首先計算FFT得到每個子段上的功率譜觀測值，然后，利用所有子段共MK個觀測值進行功率計算。可將周期圖等效為濾波器組，則對應的周期圖濾波器組第k通道濾波器的單位樣值響應為:

此濾波器的樣值響應長度為D，通帶寬度為1/D，對應的傳遞函數為:

例如對于D=128，周期圖低通濾波器的頻響曲線如圖5中的細實線所示，圖中橫坐標為歸一化圓頻率，抽樣頻率歸一化為2π，縱坐標為濾波器幅頻響應。圖中粗實線同時也給出了D=128，多相因子K=8，采用Parks-McClellan等紋波方法設計的PFB低通濾波器的幅頻響應，點劃線為通道數D的理想通帶形狀。D通道PFB與周期圖濾波器性能參數如表1。

圖5 相同通道數的周期圖與PFB濾波器頻率響應(D=128;K=8)Fig.5 Frequency responses of the D-channel periodogram and PFB filters(—both have D=128;K=8)

圖6 相同通道數的加Hamming時間窗周期圖與PFB濾波器頻率響應(D=128;K=8)Fig.6 Frequency responses of the D-channel periodogram and PFB filters with Hammming time windows(—both have D=128;K=8)

表1 PFB與周期圖濾波器參數比較Table 1 Comprison between PFB and periodogram filters

表中矩形系數定義為30 dB帶寬與3 dB帶寬比值，同時列出的包括KD通道周期圖的濾波器性能參數。

在相同的譜分辨率即通道數D條件下，譜估計器的性能取決于濾波器的通帶特性(通帶矩形系數與旁瓣抑制水平)，由此來看，PFB方法能夠獲得比周期圖濾波器更好的通帶矩形系數和旁瓣抑制性能，因而從理論上來說應該能夠獲得更好的譜估計質量。

圖6和表1中給出了加Hamming時間窗(直接加在樣本數據而非自相關函數上)的周期圖與同通道數PFB的比較結果。加窗周期圖在提高旁瓣抑制比性能的同時，降低了通帶選擇性能。可以看到，在加Hamming時間窗的情況下，對應的旁瓣抑制比提高到43 dB，但濾波器主瓣寬度增加了約1倍。周期圖方法的旁瓣抑制比和主瓣寬度是一對矛盾，在保持譜通道數不變的情況下，加窗可以改善旁瓣抑制水平，但要以損失一定的主瓣選擇性能為代價。而PFB可以較好地解決這個矛盾，在相同譜通道數情況下，在增加旁瓣抑制比的同時，仍然保持原有的主瓣3 dB帶寬，并且具有更好的通帶矩形系數。

不過，與周期圖方法相比，PFB在提高通帶性能的同時，也付出了額外的代價，因為在相同通道數的限制下，PFB使用了具有更長單位樣值響應的濾波器;而一般來說，單位樣值響應更長的濾波器總是對應著更多的計算量和資源，因此下面對相同單位樣值響應長度條件下的PFB與周期圖方法性能進行比較。

(2)相同單位樣值響應長度KD

由(8)式知，周期圖濾波器的單位樣值響應長度與其譜通道數總是相等的，因此在濾波器樣值響應長度KD的條件下，周期圖方法可達到的通道數為KD，這是PFB方法的K倍。對于KD通道周期濾波器，其性能與D通道周期圖相似，只需將參數D換成KD即可。

對應D=128、K=8，具有相同單位樣值響應長度KD=1024的PFB與周期圖濾波器頻率響應如圖7，而對應濾波器性能參數列于表1中。

圖7 相同單位樣值響應長度KD=1024的周期圖與PFB濾波器頻率響應(D=128;K=8)Fig.7 Frequency responses of periodogram and PFB filters with the same unit of sampling response time(—both have D=128;K=8)

其中粗實線為128通道的PFB原型濾波器的幅頻響應，細實線為相同單位樣值響應長度周期圖濾波器的幅頻響應，點劃線和點線分別128通道和1024通道的矩形通帶示意。KD通道周期圖的濾波器性能參數見表1。

從圖7和表1可見，在相同的單位樣值響應長度內，周期圖方法可以取得此情況下能達到的最高譜分辨率，PFB方法則是通過放松對分辨率的要求，以換取更好的通帶矩形系數和旁瓣抑制性能，在這一點上PFB與加窗的效果是相似的。

2.3 估計方差

在樣本量足夠大的情況下，上述幾種方法的估計方差可以降低到足夠小，但是方差下降速度有所不同。PFB與D通道周期圖都是在每D個濾波輸出中抽取了一個輸出用于功率計算，而KD通道周期圖方法則是在KD個輸出中抽取一個輸出用于功率計算，因此對應于相同的樣本量，KD通道周期圖方差約為PFB和D通道周期圖方法的K倍。或者說，PFB方法的方差降速是相同單位樣值響應長度周期圖方法的K倍。

從上面的分析可以看出，PFB的估計質量與直接實現的濾波器組是一樣的，它相對于直接實現結構的優點在于計算量和存儲資源上的大幅降低。下面對PFB方法實現所需的計算量、存儲資源進行考察。

2.4 計算量與資源

下表列出了在KD個樣本數據輸入的情況下，上述3種方法的計算量、存儲資源以及主瓣通道寬度、方差降低的水平。根據FFT實現方法不同，FFT暫存可能需要不同的存儲量，此處以全并行方式計算。

表2 PFB與周期圖的計算量、資源需求、方差水平Table2 Calculation amounts，resource requirements，and variance levels of the PFB and periokogram filters

從上表中可以看到，在KD點樣本輸入的情況下，KD通道的周期圖方法有更高的譜分辨率，而D通道的PFB與周期圖方法的方差降低到其1/K的水平，有更快的方差下降降速，這意味實用中達到相同的信噪比需要更低的采集時間。依據上述表達式繪出3種方法計算量的曲線，如圖8。

圖8 3種譜估計方法計算量的比較Fig.8 Computation amounts of three estimation methods

可以看到，對于相同的輸入樣本量，3種方法中，D通道周期圖的乘法與加法計算量均為最小，PFB的乘法量稍大一些。例如在16 384通道的情況下，PFB乘法次數為D通道周期圖的1.38倍，為KD通道周期圖計算量的1.14倍，幾種方法的計算量大體在同一數量級上。

2.5 小結

從上述分析中得到的幾個結論如下:

(1)與同通道數的周期圖譜分析相比，PFB能夠得到更接近理想帶通響應的幅頻響應，從而降低估計偏差，理論上可以獲得更好的譜估計質量，但實現時需要增加一定的資源和計算量;

(2)與同單位樣值響應長度的周期圖譜分析相比，PFB仍然具有幅頻響應和旁瓣抑制性能的優勢，但譜分辨率較低，降低的倍數取決于多相因子K;

(3)與直接實現的濾波器組譜分析相比，PFB通過抽取大幅度降低了資源與計算量，這也是多相結構的主要優點之一。

需要說明的是，上面的結論都基于大數據樣本譜分析的情形，對于小樣本情況，無論是PFB還是周期圖方法都需要在方差性能和分辨率之間進行折中，有時可采用參數化方法來提高分辨率［5］。

3 仿真

為檢驗PFB譜分析方法的效果，應用MATLAB進行仿真，構造一個4階ARMA連續譜加單頻正弦構成的輸入信號，其中正弦信號的歸一化頻率為0.54π(采樣頻率為2π)，此信號的頻譜如圖9。

圖9 輸入信號的功率譜Fig.9 Power spectrum of the input signal

在仿真中模擬FAST密云模型接收機基帶系統的處理過程(見下節)，對輸入信號采用正交混頻，將信號中心頻率混至零頻，進行雙通道采集，獲得原信號的IQ分量。由得到雙通道IQ分量組成復信號只含有原信號的單邊頻譜。

分別應用K=8，D=2048通道PFB、2048通道周期圖以及16 384通道周期圖方法對此單邊譜復信號進行譜估計，得到的功率譜如圖10。

圖10 PFB與周期圖方法得到的功率譜Fig.10 Estimated power spectra with the PFB and periodogram filters

圖10從上至下依次為2048通道PFB、2048通道周期圖、16 384通道周期圖方法得到的功率譜。在等量數據樣本輸入的情況下，2048通道數的PFB與周期圖方法得到的功率譜估計結果沒有大的差別，與相同單位樣值響應長度(16 384通道)的周期圖方法相比，方差更小。

在比較仿真結果中譜線強度與連續譜基線的差值上，下圖要比上兩圖高8 dB，這說明16 384通道的周期圖方法由于具有更高的分辨率，對于噪聲中微弱的單頻譜線的檢測而言是有好處的。因此，在相同的資源消耗下，PFB有更好的頻率選擇特性，而周期圖能夠獲得更高的譜分辨率，從這個意義上說，兩種方法各有所長，實際中可以根據天文觀測的科學目標進行選用。例如需要研究譜線的精細結構或譜線的漂移時，即可選用可獲得更高分辨率的周期圖方法。

4 實測

FAST接收機組利用FAST密云模型L波段接收機以及新研制的基帶處理和數據采集系統進行了銀河系中性氫譜線的觀測實驗，下面對此進行介紹。

FAST密云模型L波段接收機為單偏振接收機，中頻帶寬100 MHz，基帶系統采用正交混頻，如圖11，模擬基帶低通濾波器LPF帶寬為5 MHz，之后進行雙通道復采樣信號，采樣頻率為10 MHz，得到有效帶寬10 MHz的復信號。

圖11 FAST模型接收機基帶系統信號處理流程Fig.11 Processing procedure of the baseband system of the MyFAST Receiver

基帶處理系統實物如圖12。

圖12 FAST模型接收機基帶系統實物Fig.12 Baseband system of the MyFAST Receiver

利用上述基帶處理與數據采集系統采集得到的中性氫數據，取出nKD=8000×8×2048個樣本，分別用2048通道PFB和周期圖方法對采樣數據進行處理。圖13給出了樣本信號分析得到的相對功率譜密度圖，其中上部為PFB得到的功率譜，下部為2048通道周期圖的結果。

圖13 PFB與周期圖得到的HI信號功率譜Fig.13 HI line profiles observed with the PFB(upper)and Periodogram(lower)filters

從圖中可見，HI譜線峰值位于頻率1 420.3 MHz處。在相同樣本數輸入的情況下，兩種方法得到的功率譜結果方差均為0.1 dB，可見兩者的方差性能是相同的。而在圖中所顯示的頻率范圍內，上述兩種方法得到的功率譜估計值之間的差值序列(PFB－FFT)如圖14。

圖14 PFB與周期圖方法估計結果的差值Fig.14 Difference between the power spectra obtained with the PFB and periodogram filters

在圖中所示的頻率范圍內，此差值序列的均值為－0.05 dB，兩種方法估計結果之間存在很小的系統差值，這個差值實際上顯示了頻譜泄漏的大小。也即，在相同的通道數情況下，由于頻譜泄漏，周期圖方法的估計結果中多數通道的估計值較PFB方法更大，估計出的功率譜相對于PFB更加平滑，因此兩者差值序列的均值小于0，這驗證了PFB相對于周期圖方法，獲得了估計偏差性能上的改善。

5 結論

本文對多相濾波器組應用于譜分析的性能進行了討論，利用多相濾波器組進行譜估計，由于濾波器相對于同通道數的周期圖具有更好的頻率響應、更少的頻譜泄漏，理論上能夠獲得更好的估計質量。

在銀河系21 cm中性氫譜線實測中，PFB估計效果與同通道數的周期圖方法相比，方差性能相同，并降低了功率譜估計結果的偏差。由于快速增長的數字芯片性能與集成度使得應用PFB所需的資源變得越來越容易獲得，因而PFB是一種很好的數字式頻譜儀可選方案。另一方面，多相濾波器組作為一種信號子帶分解和多抽樣率處理的方法，在信號處理中的應用并不限于功率譜分析，如脈沖星消色散、VLBI基帶轉換［6］等場合也有很好的應用。

［1］胡廣書.數字信號處理——理論、算法與實現［M］.北京:清華大學出版社，2005:496－526.

［2］Ronald E Crochiere，Lawrence R Rabiner.Multirate Digital Signal Processing［M］.Prentice-Hall，Inc，1983:290 －325.

［3］Fredric J Harris，Chris Dick，Michael Rice.Digital Receivers and Transmitters Using Polyphase Filter Banks for Wireless Communications［J］.Microwave Theory and Techniques，2003，51(4):1395－1412.

［4］董暉，姜秋喜，畢大平.多相濾波寬帶信道化數字接收機［J］.雷達科學與技術.2007，5(1):73－77.Dong Hui， Jiang Qiuxi， Bi Daping. Polyphase Filtering Wide-Band Channelized Digital Receiver［J］.Radar Science and Technology，2007，5(1):73 －77.

［5］Petre Stoica，Randolph Moses.現代信號譜分析［M］.吳仁彪，韓萍，馮青，譯.北京:電子工業出版社，2007:69－87.

［6］陳嵐，張秀忠.用于VLBI數字基帶轉換的多相濾波技術研究［J］.天文學進展，2008，26(1):87－94.Chen Lan，Zhang Xiuzhong.The Study of DBBC Based on Poly-phase Filter Banks and FFT in VLBI［J］.Progress in Astronomy，2008，26(1):87 －94.