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

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

（1.中國(guó)科學(xué)院聲學(xué)研究所，聲場(chǎng)聲信息國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100190；2.中國(guó)科學(xué)院大學(xué)，北京 100049）

0 引言

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

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

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

1 廣義Warblet 變換

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

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

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

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

GWT 算法可以表示為[8]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

（1）對(duì)信號(hào)進(jìn)行STFT，高斯時(shí)間窗長(zhǎng)L、滑動(dòng)長(zhǎng)度為0.02L。在一定帶寬內(nèi)利用式（5）得到瞬時(shí)頻率的估計(jì)。

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

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

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

將STFT 作為算法第1 次迭代，兩次瞬時(shí)頻率估計(jì)相對(duì)差值為

3 仿真驗(yàn)證

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

設(shè)頻率時(shí)變信號(hào)滿(mǎn)足式（1）并包含2個(gè)頻率調(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，信號(hào)時(shí)長(zhǎng)1 000 s。對(duì)信號(hào)添加譜級(jí)信噪比為?2 dB 的高斯白噪聲，譜級(jí)信噪比定義為全帶寬信噪比與工作帶寬B的乘積，對(duì)數(shù)表示為[14]

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

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

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

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

表2 FMIPC-GWT 前3 次迭代調(diào)頻參數(shù)估計(jì)結(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 各次迭代瞬時(shí)頻率恢復(fù)值、真實(shí)值及 提取值對(duì)比 Fig.4 Comparison of the recovery,real and extracted values of the instantaneous frequency for each iteration of the FMIPC-GWT method

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

保持信號(hào)參數(shù)不變，通過(guò)調(diào)整噪聲強(qiáng)度考察FMIPC-GWT 方法性能。利用蒙特卡洛（Monte Carlo）方法研究FMIPC-GWT 方法在不同信噪比條件下瞬時(shí)頻率提取性能，譜級(jí)信噪比為?10～5 dB，按照1 dB 間隔分別進(jìn)行100 次試驗(yàn)，Monte Carlo試驗(yàn)結(jié)果如圖5 所示。圖5 中，縱坐標(biāo)表示可以準(zhǔn)確提取瞬時(shí)頻率的試驗(yàn)次數(shù)占比。由圖5 可知，當(dāng)譜級(jí)信噪比為?3 dB 時(shí)，F(xiàn)MIPC-GWT 方法仍有71%的可能正確輸出，表明算法在較低信噪比下仍有較好性能。

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

4 實(shí)驗(yàn)數(shù)據(jù)處理

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

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

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

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

5 結(jié)論

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

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

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