調頻初相位補償的廣義Warblet 變換時頻分析

張 地，周士弘，戚聿波，杜淑媛

（1.中國科學院聲學研究所，聲場聲信息國家重點實驗室，北京 100190；2.中國科學院大學，北京 100049）

0 引言

時頻分析能夠體現被分析信號的時間-頻率二維能量分布，是常用的準穩態信號處理方法。時頻分析方法主要有兩類：不需要信號先驗知識的非參數化方法和需要已知信號先驗知識的參數化方法。非參數化方法包含短時傅里葉變換（Short Time Fourier Transform,STFT）、Wigner-Vile 分布（Wigner-Vile Distribution,WVD）和連續小波變換（Continuous Wavelet Transform,CWT）及其改進算法等。STFT方法將信號劃分為若干重疊或不重疊的小段并對每一小段信號加窗后做傅里葉變換，性能由所選取的時間窗決定，若時間窗較長則頻率分辨率高但時間分辨率低，反之則時間分辨率高但頻率分辨率低。Kwok 等[1]提出了自適應STFT（Adaptive STFT,ASTFT）方法，利用最大似然準則在不同時間選用不同的時間窗，進而得到最高的時頻輸出峰值。WVD方法對各段信號自相關函數做傅里葉變換，調頻信號自相關能夠得到比信號與傅里葉基函數相關輸出更大的模值，但自相關操作會引入交叉項，對瞬時頻率估計帶來很大影響。呂軍等[2]提出了一種改進的STFT-WVD 方法，利用在STFT 輸出中不存在WVD 輸出交叉項的特性，將兩種時頻分析輸出矩陣對應位置處元素相乘并設定合適的閾值進行后續分析。CWT 方法通過改變時頻分辨尺度提高分析性能，在低頻段具有更高的頻率分辨率，在高頻段則具有更高的時間分辨率[3]。與STFT 類似，CWT性能與被分析信號和選取的母小波有關。

參數化方法可利用廣義時頻變換（General Parameterized Time-Frequency Transform,GPTF）[4]框架統一描述，包含多項式Chirplet 變換（Polynomial Chirplet Transform,PCT）[5]、多普勒Chirplet 變換（Doppler Chirplet Transform,DCT）[6]、樣條基Chirplet 變換（Spline-Kernelled Chirplet Transform,SCT）[7]、廣 義 Warblet 變換（Generalized Warblet Transform,GWT）[8-9]等。這些方法從不同角度描述信號時頻走勢，PCT、SCT 和GWT 分別采用多項式、三次樣條插值函數和三角函數描述信號的時頻走勢；DCT 則直接解算信號源與接收點間徑向運動引起的多普勒頻移關系得到瞬時頻率表達式。GPTF 的核心思想是利用旋轉算子抵消信號相位中除頻率一次項之外的部分，在時頻圖上體現為將時變線譜旋轉為時不變線譜，然后利用平移算子修正旋轉算子帶來的頻率偏移，再利用單頻信號作為核函數對算子作用后的信號做傅里葉變換，此時每個時間分析窗內信號頻率均不隨時間變化，能夠與傅里葉基函數完全匹配，使能量集中在某一頻點，從而獲得更高的時頻能量聚集度，更有利于線譜檢測和瞬時頻率估計。PCT 和DCT 方法適用于頻率非周期變化信號，在處理頻率周期變化信號方面存在局限性；SCT 方法雖能處理頻率周期變化信號，但算法涉及到的樣條基函數較復雜；GWT 方法利用三角函數描述信號瞬時頻率，算法較易實現且對頻率周期變化和非周期變化類信號均有較好的處理效果。傳統的Warblet 變換（Warblet Transform,WT）[9]方法在每一個時間窗內確定調頻相位，算法較難實現。GWT 在此基礎上，利用旋轉算子和平移算子對信號整體操作后再進行STFT，算法更易實現。GWT 方法采用三角函數描述信號瞬時頻率，適用于瞬時頻率隨時間近似呈三角函數變化的信號，例如微多普勒雷達信號[10]以及頻率周期變化的水聲信號[11]。Zhou 等[12]提出修正的廣義參數化時頻分析（Modified Generalized Parameterized Time-Frequency,MGPTF）方法和相干單程多普勒干涉 MGPTF（Coherent Single Range Doppler Interferometry-MGPTF,CSRDI-MGPTF）方法，MGPTF 利用三角函數描述信號瞬時頻率并在GWT 基礎上增加了一層關于調制成分個數的循環，使得瞬時頻率的估計更為準確；CSRDI-MGPTF 方法能夠適用時頻軌跡部分中斷的情況。王曉[13]利用時域平移校正法對瞬時頻率傅里葉變換結果進行校正，提高時頻軌跡參數的估計精度，進而提高瞬時頻率估計精度。

現有算法在調制成分估計時均沒有考慮時頻輸出起始時刻與信號起始時刻之間存在由1/2 個時間窗長度引起的相位差。對于頻率時變信號，當時頻分析所需時間窗較長時，如果不對調頻初相位進行補償則無法對調頻成分進行準確估計，GWT 算法性能將下降甚至失效。本文提出一種調頻初相位補償的GWT（Frequency Modulation Initial Phase Compensation-GWT,FMIPC-GWT）方法，在長時間窗時頻分析時依舊能夠準確估計瞬時頻率，對時變線譜檢測和瞬時頻率估計具有重要意義。

1 廣義Warblet 變換

當信號頻率隨時間變化時，對整段信號做頻譜分析無法體現信號頻率隨時間變化情況。通常假設信號在各小段時間內頻率不發生變化，對各小段時間做頻譜分析，再將所有結果表示為能量在時間-頻率二維平面上的分布。然而頻率時變信號的數學形式決定了在時頻分析時間窗內，信號頻率依舊是變化的，傅里葉變換輸出無法將信號全部能量集中到某一頻點上，造成時頻能量聚集度下降、輸出信噪比降低以及瞬時頻率估計準確性下降。

GWT 算法的典型應用場景之一是處理正弦調頻信號，信號模型為

其中：A0為信號幅度；f0為中心頻率；m為調制成分序數；Fm和ψm分別為調制頻率和調頻初相位；?0為載波初相位；Bm描述頻率調制成分的強度，為方便后續推導，令，am為頻率調頻幅度。信號瞬時頻率可以表示為

定義旋轉算子φR和平移算子φS分別為

GWT 算法可以表示為[8]

其中：t0為各時間窗中心時刻；h(t)為時間窗函數，常用時間窗函數有高斯窗、漢寧窗等。

由式（4）可以看出，旋轉算子將信號相位項中時變部分抵消掉，使得信號頻率不隨時間變化，之后通過平移算子將各時間窗內的信號頻率修正到原信號頻率。與STFT 相比所不同的是經過旋轉算子作用后，信號頻率在每一個分析時間窗內都是恒定不變的，能夠與傅里葉基函數exp (i 2πft)的共軛完全匹配，從而使得時頻分析算法輸出信噪比最大，提高了時變線譜的檢測能力。

對于不是正弦調頻信號的其他非線性調頻信號，雖然其瞬時頻率不再具有式（2）的形式，但依舊可以通過傅里葉變換將瞬時頻率隨時間變化的情況用若干三角函數的線性組合來表示。因此，GWT類方法依舊適用[4]。

2 調頻初相位補償的GWT 方法

GWT 法需利用瞬時頻率軌跡估計調頻參數。由于時頻分析結果對應時刻為所選取的各時間窗對應的中心時刻，當時間窗較長時，調頻初相位估計結果將出現較大偏差。FMIPC-GWT 方法對調頻初相位估計結果進行修正，避免了調頻參數估計結果與信號原有調頻成分間的失配。同時FMIPC-GWT 瞬時頻率估計結果可再次作為參數輸入FMIPC-GWT 算法進行進一步的時頻分析，經過多次迭代可進一步提高線譜瞬時頻率估計準確性。

2.1 調頻參數估計與初相位補償

GWT 算法的核心思想是用一系列三角函數刻畫信號瞬時頻率，這需要對調頻參數進行估計，確定各調制成分幅度am、調制頻率Fm和調頻初相位ψm。首先對信號進行φR=1、φS=1 的GWT，此時GWT 退化為STFT。瞬時頻率的估計可表示為

其中：n為整數，?f=1 (?t?N)為離散頻率間隔，?t為各時間窗中心時刻的間隔，N為離散傅里葉變換點數。

傅里葉變換系數與傅里葉變換間關系為

其中：n0為各調制頻率所在頻域離散采樣點序數。

需要注意的是，上述對信號調頻參數的估計是建立在對信號時頻分析結果的基礎上，而時頻分析輸出時刻對應各時間窗的中心時刻，因此調頻初相位估計實際上是對信號瞬時頻率軌跡在第一個時間窗中心時刻處相位的估計。由于旋轉算子和平移算子是作用到起始時刻為0 的整段信號上，因此式（9）中對調頻初相位的估計與實際調頻初相位之間存在由1/2 個時間窗長度所引入的相位差。因此，相位補償后調頻初相位估計可以表示為

需要指出的是，GWT 算法中將信號瞬時頻率寫作一系列余弦函數和正弦函數的線性組合[8]，可改寫為式（2）帶有初相位的余弦函數形式，但此處初相位依然是沒有考慮1/2 個時間窗長所引入的相位。

2.2 迭代提升參數估計精度

利用STFT 結果估計調頻參數后進行FMIPC-GWT 時頻分析，可提高瞬時頻率估計精度，并可利用FMIPC-GWT 結果進行調頻參數估計，將估計結果作為參數再次輸入FMIPC-GWT 算法做進一步的時頻分析。多次迭代的FMIPC-GWT 的步驟為：

（1）對信號進行STFT，高斯時間窗長L、滑動長度為0.02L。在一定帶寬內利用式（5）得到瞬時頻率的估計。

（2）利用式（6）對瞬時頻率估計結果進行傅里葉變換得到Fins，可以得到各調制成分的調制頻率、調制幅度和調制頻率初相位。

（3）將步驟（2）中得到的調頻參數估計結果代入式（3）得到旋轉算子φR和平移算子φS，并利用得到的算子對原信號進行FMIPC-GWT 時頻分析。

（4）利用式（5）得到新的瞬時頻率估計，并重復步驟（2）～（4）直到滿足終止條件。

將STFT 作為算法第1 次迭代，兩次瞬時頻率估計相對差值為

3 仿真驗證

3.1 算法仿真實現

設頻率時變信號滿足式（1）并包含2個頻率調制成分，調制頻率Fm分別為3 mHz 和9 mHz，調頻幅度am均為0.3 Hz，調頻初相位ψm分別為π/6和π/3，中心頻率f0=30 Hz，A0=1，?0=0，信號時長1 000 s。對信號添加譜級信噪比為?2 dB 的高斯白噪聲，譜級信噪比定義為全帶寬信噪比與工作帶寬B的乘積，對數表示為[14]

其中：和RSN分別為譜級信噪比和全帶寬信噪比的對數形式，，fs=250 Hz 為信號采樣頻率。

對信號進行STFT，高斯時間窗長 80 sL=，滑動長度1.6 s。如果不對調頻初相位進行補償，直接在STFT 基礎上利用式（5）～（9）對信號進行GWT，迭代終止門限δ設為0.1%，最大迭代次數為5，算法因觸發最大迭代次數約束條件而終止，其中STFT 作為第1 次迭代。每次迭代后利用瞬時頻率軌跡進行調頻參數估計，前3 次參數估計結果如表1 所示。表1 中，調頻幅度和調制頻率誤差均為估計值與真值差值的絕對值與真實值的百分比，調頻初相位誤差為估計值與真值差值的絕對值與2π的百分比。誤差保留小數點后1 位有效數字。調頻參數估計結果，特別是調頻初相位估計結果誤差較大。利用表1 中調頻參數估計的各次迭代結果恢復的瞬時頻率與真實值以及瞬時頻率提取值對比如圖1 所示。圖1 中黑色虛線為時頻分析后利用式（5）提取到的瞬時頻率軌跡（簡稱提取值），藍色點劃線為利用對瞬時頻率軌跡進行參數估計到的調頻參數恢復出的瞬時頻率軌跡（簡稱恢復值），紅色實線為瞬時頻率真實值。由圖1 可以看到紅色實線與藍色點劃線之間存在較大誤差，這主要由相位差引起，利用藍色點劃線所表示的參數進行下一步迭代便會造成算法性能下降甚至失效。GWT 算法前3 次迭代輸出結果如圖2（a）～2（c）所示，其中圖2（a）為STFT 結果，第4～5 次迭代時頻輸出結果與第2～3次結果（圖2（b）～2（c））類似（圖略）。由于調頻參數估計誤差較大，圖2（b）～2（c）中是許多雜亂的亮線，無法分辨出線譜，時頻分析效果甚至不如STFT 方法。

表1 GWT 前3 次迭代調頻參數估計結果 Table 1 Frequency modulation parameter estimation results for the first 3 iterations of the GWT method

圖1 GWT 前3 次迭代瞬時頻率恢復值、真實值及提取值對比 Fig.1 Comparison of the recovery,real and extracted values of the instantaneous frequency for the first 3 iterations of the GWT method

圖2 仿真GWT 前3 次迭代結果 Fig.2 The simulation results for the first 3 iterations of the GWT method

利用第2 節中介紹的FMIPC-GWT 方法對信號進行時頻分析，算法經過3 次迭代后終止，各次迭代的輸出結果如圖3（a）～3（c）所示。圖3（a）為第1次迭代即STFT 方法輸出結果，由于存在較強噪聲且傅里葉分析基函數不能很好地匹配頻率時變信號，導致信號時頻軌跡存在許多斷點，很難分辨出線譜所在位置。表2 為FMIPC-GWT 算法各次迭代調頻參數估計值結果，可以看出對調頻初相位進行補償能夠有效提高調頻參數特別是調頻初相位估計的準確性。利用相對準確的調頻參數估計結果得到的旋轉算子和平移算子對信號做 GWT 即FMIPC-GWT，第2～3 次迭代的輸出結果如圖3（b）～3（c）所示，第2 次迭代輸出相較于STFT 結果有很大改善，能夠比較清晰地看到信號瞬時頻率軌跡，第3 次迭代輸出結果進一步提升了瞬時頻率軌跡的清晰程度。圖4 為利用FMIPC-GWT 時頻分析算法各次迭代調頻參數估計值恢復出的瞬時頻率軌跡與真實值和在時頻圖上提取的瞬時頻率軌跡的對比，可以看到經過對調頻初相位補償后的調頻參數能夠很好地描述信號瞬時頻率走勢，經過迭代后瞬時頻率提取值與真實值基本一致。

表2 FMIPC-GWT 前3 次迭代調頻參數估計結果 Table 2 Frequency modulation parameter estimation results for the first 3 iterations of the FMIPC-GWT method

圖3 仿真FMIPC-GWT 各次迭代輸出結果 Fig.3 The simulation results for each interation of the FMIPC-GWT method

圖4 FMIPC-GWT 各次迭代瞬時頻率恢復值、真實值及 提取值對比 Fig.4 Comparison of the recovery,real and extracted values of the instantaneous frequency for each iteration of the FMIPC-GWT method

3.2 不同信噪比下算法性能分析

保持信號參數不變，通過調整噪聲強度考察FMIPC-GWT 方法性能。利用蒙特卡洛（Monte Carlo）方法研究FMIPC-GWT 方法在不同信噪比條件下瞬時頻率提取性能，譜級信噪比為?10～5 dB，按照1 dB 間隔分別進行100 次試驗，Monte Carlo試驗結果如圖5 所示。圖5 中，縱坐標表示可以準確提取瞬時頻率的試驗次數占比。由圖5 可知，當譜級信噪比為?3 dB 時，FMIPC-GWT 方法仍有71%的可能正確輸出，表明算法在較低信噪比下仍有較好性能。

圖5 不同信噪比下FMIPC-GWT 算法準確提取瞬時頻率的占比 Fig.5 The proportion of accurately extracting instantaneous frequency by FMIPC-GWT method under different signal to noise ratios

4 實驗數據處理

2018 年夏季，我們在某海域開展了一次水聲實驗觀測過往船只輻射噪聲，某時段接收噪聲信號線譜具有瞬時頻率時變特征。高斯時間窗長為16 000個數據點，每次向前滑動320 個數據點，時頻分析算法共輸出476 個時間采樣點。利用式（5）在STFT結果中提取瞬時頻率，對瞬時頻率減去均值后做頻譜分析如圖6 所示，由高到低依次選取7 個峰值點所在頻率作為調制頻率。

圖6 STFT 提取瞬時頻率的頻譜 Fig.6 The frequency spectrum of the instantaneous frequency extracted by STFT method

如不補償調頻初相位，直接對信號使用GWT算法即在調頻參數估計時利用式（8）～（9）求調頻幅度和調頻初相位，迭代終止門限δ設為0.1%，最大迭代次數設為5，算法經過2 次迭代終止，分析結果如圖7（a）～7（b）所示。圖7（a）為STFT 結果，由于海洋環境噪聲及信號瞬時頻率隨時間變化使線譜時頻能量聚集度下降等因素，STFT 結果中線譜存在許多不連續的地方。由于未對1/2 個時間窗長引起的相位差進行補償，使得調頻參數估計結果與信號真實調頻參數存在較大誤差，旋轉算子和平移算子與信號失配，GWT 法后續迭代結果并未改善STFT 的輸出效果，甚至加劇了線譜的不連續性和時頻能量分散程度。GWT 第2 次迭代的結果如圖7（b）所示。對信號做相同時間窗長、相同滑動時間長度的FMIPC-GWT 處理即在調頻參數估計時利用式（8）、（10）求調頻幅度和調頻初相位，迭代終止門限δ設為0.1%，最大迭代次數設為5，算法經2次迭代后終止，結果如圖7（c）所示。相比STFT 和GWT 結果，FMIPC-GWT 結果線譜連續性更好，強度更高，FMIPC-GWT 算法可有效提高瞬時頻率變化線譜時頻能量聚集度，有利于提高時變線譜時間分辨率和頻率分辨率。

圖7 實驗數據GWT/FMIPC-GWT 時頻分析結果 Fig.7 Time-frequency analysis results of the experimental data obtained by GWT/FMIPC-GWT method

5 結論

本文針對GWT 算法在長時間窗處理時算法失效問題，提出一種調頻初相位補償的廣義Warblet變換時頻分析方法。通過在調頻參數估計時將調頻成分在1/2 時間窗長內經過的相位補償到調頻初相位估計值中，提高調頻參數估計的準確性，避免因長時間窗處理導致旋轉算子和平移算子與信號失配問題的出現，提高了瞬時頻率時變信號線譜時頻能量聚集度和頻率分辨率。仿真和海上實驗數據處理結果表明，FMIPC-GWT 方法能夠避免GWT 方法在長時間窗分析時出現參數失配問題，在瞬時頻率變化信號線譜檢測中具有比STFT 和GWT 方法更優的性能。

FMIPC-GWT 方法需要的調頻成分個數需要根據其瞬時頻率軌跡的傅里葉變換結果給出。若選取的調制成分個數較少，則不足以描述瞬時頻率的變化，算法輸出信噪比下降。而如果選取的調制成分個數較多，則可能出現過擬合，降低瞬時頻率估計的準確性。此外，當分析頻帶內存在較強的干擾線譜時，調頻參數估計將出現錯誤導致算法失效，提高算法抗干擾能力將是未來的研究方向。

致謝感謝2018年夏季參加海上實驗的全體同志，他們的辛勤工作為本文提供了寶貴的實驗數據。