調(diào)頻初相位補償?shù)膹V義Warblet 變換時頻分析

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

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

0 引言

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

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

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

1 廣義Warblet 變換

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

GWT 算法的典型應(yīng)用場景之一是處理正弦調(diào)頻信號，信號模型為

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

定義旋轉(zhuǎn)算子φR和平移算子φS分別為

GWT 算法可以表示為[8]

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

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

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

2 調(diào)頻初相位補償?shù)腉WT 方法

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

2.1 調(diào)頻參數(shù)估計與初相位補償

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

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

傅里葉變換系數(shù)與傅里葉變換間關(guān)系為

其中：n0為各調(diào)制頻率所在頻域離散采樣點序數(shù)。

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

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

2.2 迭代提升參數(shù)估計精度

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

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

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

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

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

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

3 仿真驗證

3.1 算法仿真實現(xiàn)

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

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

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

表1 GWT 前3 次迭代調(diào)頻參數(shù)估計結(jié)果 Table 1 Frequency modulation parameter estimation results for the first 3 iterations of the GWT method

圖1 GWT 前3 次迭代瞬時頻率恢復(fù)值、真實值及提取值對比 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 次迭代結(jié)果 Fig.2 The simulation results for the first 3 iterations of the GWT method

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

表2 FMIPC-GWT 前3 次迭代調(diào)頻參數(shù)估計結(jié)果 Table 2 Frequency modulation parameter estimation results for the first 3 iterations of the FMIPC-GWT method

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

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

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

保持信號參數(shù)不變，通過調(diào)整噪聲強度考察FMIPC-GWT 方法性能。利用蒙特卡洛（Monte Carlo）方法研究FMIPC-GWT 方法在不同信噪比條件下瞬時頻率提取性能，譜級信噪比為?10～5 dB，按照1 dB 間隔分別進行100 次試驗，Monte Carlo試驗結(jié)果如圖5 所示。圖5 中，縱坐標表示可以準確提取瞬時頻率的試驗次數(shù)占比。由圖5 可知，當譜級信噪比為?3 dB 時，F(xiàn)MIPC-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 實驗數(shù)據(jù)處理

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

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

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

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

5 結(jié)論

本文針對GWT 算法在長時間窗處理時算法失效問題，提出一種調(diào)頻初相位補償?shù)膹V義Warblet變換時頻分析方法。通過在調(diào)頻參數(shù)估計時將調(diào)頻成分在1/2 時間窗長內(nèi)經(jīng)過的相位補償?shù)秸{(diào)頻初相位估計值中，提高調(diào)頻參數(shù)估計的準確性，避免因長時間窗處理導致旋轉(zhuǎn)算子和平移算子與信號失配問題的出現(xiàn)，提高了瞬時頻率時變信號線譜時頻能量聚集度和頻率分辨率。仿真和海上實驗數(shù)據(jù)處理結(jié)果表明，F(xiàn)MIPC-GWT 方法能夠避免GWT 方法在長時間窗分析時出現(xiàn)參數(shù)失配問題，在瞬時頻率變化信號線譜檢測中具有比STFT 和GWT 方法更優(yōu)的性能。

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

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