基于插值短時(shí)分?jǐn)?shù)階傅里葉變換- 變權(quán)擬合的線性調(diào)頻信號(hào)參數(shù)估計(jì)

曹偉浩，姚直象，夏文杰，閆肅

(1.海軍工程大學(xué) 電子工程學(xué)院, 湖北 武漢 430033;2.國(guó)家海洋技術(shù)中心漳州基地籌建辦公室，福建 廈門 361000; 3.91388部隊(duì)，廣東 湛江 524000)

0 引言

線性調(diào)頻(LFM)信號(hào)是主動(dòng)聲納應(yīng)用較廣泛的一種信號(hào)形式，由于水聲環(huán)境的復(fù)雜性及水中可利用的信號(hào)頻帶有限，對(duì)噪聲污染的LFM信號(hào)進(jìn)行瞬時(shí)頻率估計(jì)變得非常重要。LFM信號(hào)高精度的瞬時(shí)頻率估計(jì)，不僅在水下目標(biāo)探測(cè)、水聲通信多普勒補(bǔ)償?shù)确矫姘缪萘酥匾巧?，而且在雷達(dá)、通信等領(lǐng)域得到了廣泛應(yīng)用。此外，提高LFM信號(hào)參數(shù)估計(jì)的精度，對(duì)水下高精度目標(biāo)識(shí)別也具有重要意義。

短時(shí)傅里葉變換(STFT)[1-2]和維格納- 威利分布(WVD)[3-4]是兩種常用的對(duì)LFM信號(hào)進(jìn)行時(shí)頻分析的方法，但是STFT對(duì)于寬帶信號(hào)頻率估計(jì)效果并不理想，WVD是雙線性變換，時(shí)頻分析存在交叉項(xiàng)問(wèn)題。分?jǐn)?shù)階傅里葉變換(FRFT)對(duì)LFM信號(hào)具有較好的能量聚集性，利用此變換對(duì)LFM信號(hào)進(jìn)行時(shí)頻分析的應(yīng)用前景廣闊[5-10]。

本文針對(duì)參數(shù)估計(jì)問(wèn)題，利用聲納、雷達(dá)等領(lǐng)域常用的一種頻譜修正算法——Rife插值法，提出一種基于Rife插值短時(shí)FRFT(STFRFT)方法來(lái)估計(jì)LFM信號(hào)瞬時(shí)頻率，通過(guò)變權(quán)擬合(VWSF)估計(jì)信號(hào)的初始頻率和調(diào)頻率，進(jìn)一步推導(dǎo)了初始頻率和調(diào)頻率估計(jì)的克拉美- 羅限(CRLB)，分析了CRLB隨參數(shù)變化的規(guī)律。

1 時(shí)頻分析和參數(shù)估計(jì)方法原理

1.1 短時(shí)分?jǐn)?shù)階傅里葉變換

LFM信號(hào)作為水聲信號(hào)處理領(lǐng)域中一種常用的信號(hào)，其表達(dá)式如下：

x(t)=Aexp {-j(2πf0t+πμt2)}，

式中：A為信號(hào)幅度；f0和μ分別為信號(hào)初始頻率和調(diào)頻率；t為信號(hào)持續(xù)時(shí)間。實(shí)際處理中，信號(hào)采用離散形式，故給出離散FRFT的表達(dá)式[10]為

(1)

式中：Aα=exp {-jπsgn(sinα)/4+jα/2}/|sinα|1/2，α=(p-1)π/2，p為FRFT最優(yōu)階數(shù)，0<|p|<2；x為量綱歸一化后新的橫坐標(biāo)；Δx為量綱歸一化寬度；X1(·)為信號(hào)的傅里葉變換；N為采樣點(diǎn)數(shù)，|m0|≤N，|n0|≤N.

若對(duì)信號(hào)進(jìn)行加窗處理，則由(1)式可得第i個(gè)短時(shí)窗內(nèi)離散FRFT的表達(dá)式為

(2)

式中：i=1,2,…,K，K為窗移動(dòng)的總次數(shù)，ΔK為移動(dòng)步進(jìn)；m=0,1,2,…,M-1；h(·)為長(zhǎng)度為M的短時(shí)窗；αi=(pi-1)π/2，pi為第i個(gè)窗內(nèi)FRFT的最優(yōu)階數(shù)。通?？紤]到估計(jì)精度與計(jì)算量等因素，移動(dòng)步進(jìn)取窗長(zhǎng)的0.1～0.3倍之間比較合適。

圖1 LFM信號(hào)歸一化功率譜密度

FRFT作為傅里葉變換的一種推廣，直觀上可以認(rèn)為是Chirp基分解，非常適用于處理Chirp信號(hào)。對(duì)一個(gè)窗內(nèi)的LFM信號(hào)分別作傅里葉變換和FRFT，得到歸一化功率譜密度如圖1所示。從圖1中可以明顯看出，經(jīng)過(guò)傅里葉變換后頻譜幅度在一定帶寬內(nèi)起伏變化，而經(jīng)過(guò)FRFT后頻譜幅度位于某一頻率，表明FRFT對(duì)LFM信號(hào)具有更好的能量聚集性。

1.2 基于STFRFT的LFM瞬時(shí)頻率估計(jì)

對(duì)信號(hào)進(jìn)行瞬時(shí)頻率估計(jì)時(shí)，在一定頻率范圍內(nèi)，短時(shí)窗內(nèi)的信號(hào)可近似認(rèn)為是單頻信號(hào)。STFT中通常取信號(hào)頻譜幅度的均值作為一個(gè)窗內(nèi)的頻譜幅度，而實(shí)際頻譜幅度高于此均值。因此，可以令第i個(gè)短時(shí)窗內(nèi)FRFT幅度最大值對(duì)應(yīng)的頻率為該窗內(nèi)LFM信號(hào)瞬時(shí)頻率fi的估計(jì)值i，

i=kiΔf，

由于FRFT滿足Parseval關(guān)系，可得

(3)

式中：Efi、Eti分別為信號(hào)頻域和時(shí)域能量；Ni為第i個(gè)短時(shí)窗內(nèi)時(shí)域采樣點(diǎn)數(shù)；Nαi為第i個(gè)短時(shí)窗內(nèi)分?jǐn)?shù)階域采樣點(diǎn)數(shù)，Nαi=Nfisinαi=μT2sinαi，Nfi為第i個(gè)短時(shí)窗內(nèi)頻域采樣點(diǎn)數(shù)，αi=piπ/2，pi為第i個(gè)短時(shí)窗內(nèi)FRFT的最優(yōu)階數(shù)，μ為調(diào)頻率，T為短時(shí)窗時(shí)長(zhǎng)；Ai為第i個(gè)短時(shí)窗內(nèi)信號(hào)的頻譜幅度。由(3)式可以得到第i個(gè)短時(shí)窗內(nèi)信號(hào)經(jīng)過(guò)FRFT后的頻譜幅度為

式中：fs為采樣率。

對(duì)于離散頻譜，考慮柵欄效應(yīng)后的頻譜幅度為

(4)

當(dāng)信號(hào)頻率不是頻率分辨率的整數(shù)倍時(shí)，不可避免地會(huì)存在頻譜泄露情況，例如1個(gè)主瓣內(nèi)會(huì)存在2個(gè)峰值。采用較常用的Rife插值法進(jìn)行頻譜修正，該方法的頻差表示[11]為

δi=

分別采用STFT、插值STFT(ISTFT)、STFRFT和插值STFRFT(ISTFRFT)4種方法估計(jì)LFM信號(hào)的瞬時(shí)頻率，結(jié)果如圖2所示。由圖2可見(jiàn)，整體而言，STFRFT的方法較STFT方法更優(yōu)，對(duì)于STFRFT方法而言，ISTFRFT方法比STFRFT方法估計(jì)的頻率更加接近實(shí)際頻率。

圖2 LFM信號(hào)瞬時(shí)頻率估計(jì)

在得到LFM信號(hào)每一時(shí)刻的瞬時(shí)頻率后，利用曲線擬合方法估計(jì)其初始頻率和調(diào)頻率。最小二乘擬合(CLSF)是將觀測(cè)值賦予相同的權(quán)值，使觀測(cè)值與擬合值的誤差平方和達(dá)到最小。當(dāng)Rife插值估計(jì)的頻率存在異常點(diǎn)時(shí)，將每次觀測(cè)值的影響同等看待顯然不合理，應(yīng)該對(duì)誤差平方和較大的項(xiàng)賦予較小的權(quán)值，誤差平方和較小的項(xiàng)賦予較大的權(quán)值。因此采用VWSF方法估計(jì)信號(hào)的初始頻率和調(diào)頻率。

1.3 基于VWSF的LFM初始頻率和調(diào)頻率估計(jì)

式中：f′i=a1+a2t，a1、a2為2個(gè)待估計(jì)參量；Wi為第i個(gè)時(shí)刻估計(jì)瞬時(shí)頻率的加權(quán)值。分別求RW(f)關(guān)于a1、a2的偏導(dǎo)并令其偏導(dǎo)為0，可解得2個(gè)參量的值[12]如下：

以上兩個(gè)公式中的權(quán)值可由預(yù)先設(shè)定的精度控制實(shí)現(xiàn)。

分別利用CLSF和VWSF兩種方法估計(jì)的瞬時(shí)頻率結(jié)果如圖3所示。從圖3中可以看出，當(dāng)估計(jì)頻率出現(xiàn)異常點(diǎn)時(shí)，利用常規(guī)CLSF方法得到的曲線偏離實(shí)際頻率，而VWSF方法能夠很好地?cái)M合實(shí)際頻率。

圖3 LFM信號(hào)瞬時(shí)頻率擬合曲線

2 初始頻率和調(diào)頻率估計(jì)的CRLB研究

2.1 插值STFRFT頻率估計(jì)的CRLB

在一定條件下，任何估計(jì)量都存在一個(gè)均方誤差下限，即CRLB.對(duì)于未知確定參量f0和μ的均方誤差等于估計(jì)量的方差，該估計(jì)為無(wú)偏估計(jì)。因此，求得估計(jì)量的方差限即為其CRLB[13].

本文采用Rife插值法，參照文獻(xiàn)[14]中的方法，給出一個(gè)窗內(nèi)插值STFRFT方法估計(jì)頻率的均方誤差為

(5)

假設(shè)估計(jì)的LFM信號(hào)瞬時(shí)頻率為

f′i=f′0+μ′i+wi，

式中：f′0和μ′分別為估計(jì)的初始頻率和調(diào)頻率；wi為高斯白噪聲，其方差為σ2.則瞬時(shí)頻率的似然函數(shù)為

P(x|f′0,μ′)=

關(guān)于未知參量的Fisher信息矩陣為

式中：Pf′0f′0為函數(shù)P對(duì)f′0的2階偏導(dǎo)；Pf′0μ′、Pμ′f′0為函數(shù)P對(duì)f′0和μ′的2階偏導(dǎo)；Pμ′μ′為函數(shù)P對(duì)μ′的2階偏導(dǎo)。由此可得

(6)

2.2 初始頻率估計(jì)的CRLB分析

根據(jù)(6)式可以得到初始頻率估計(jì)的CRLB為

(7)

在輸入信噪比一定的情況下，取δ=0.5，初始頻率的CRLB與滑窗次數(shù)、FRFT的最優(yōu)階數(shù)有關(guān)。對(duì)(7)式的K求偏導(dǎo)，有

(8)

考慮實(shí)際情況，當(dāng)K大于1時(shí)(8)式恒小于0，表明初始頻率的估計(jì)精度隨滑窗次數(shù)增加而提高，即CRLB隨滑窗次數(shù)的增加而減小。同理，對(duì)α求偏導(dǎo)，有

(9)

由于階數(shù)p是以4為周期的參數(shù)，只需考察區(qū)間p∈[-2,2]即可，根據(jù)對(duì)稱性，取階數(shù)范圍為p∈[0,2]，則α∈[0,π]。當(dāng)α∈[0，π/2]時(shí)cosα>0，(9)式為正值；當(dāng)α∈[π/2,π]時(shí)cosα<0，(9)式為負(fù)值。因此初始頻率估計(jì)精度隨階數(shù)的增加先降低、后提高，即CRLB隨階數(shù)的增加先增大、后減小。當(dāng)階數(shù)為1時(shí)，初始頻率估計(jì)性能最低，此時(shí)CRLB達(dá)到最大。當(dāng)階數(shù)不變時(shí)，初始頻率估計(jì)性能與滑窗次數(shù)呈正比，但是滑窗次數(shù)的確定還需兼顧頻率分辨率、計(jì)算量等因素。因此，分析CRLB的變化規(guī)律對(duì)窗長(zhǎng)等參數(shù)選取具有一定的指導(dǎo)意義。

初始頻率估計(jì)的CRLB仿真結(jié)果歸一化后如圖4所示。其中，圖4(a)的階數(shù)取值為[1，2]。由圖4(a)可以看出CRLB隨滑窗次數(shù)增加呈遞減趨勢(shì)。圖4(b)的階數(shù)取值為[0，2]，由圖4(b)可以看出CRLB隨階數(shù)的增加先增大、后減小，階數(shù)為1時(shí)估計(jì)性能最低。這一結(jié)論與理論分析結(jié)果一致。

2.3 調(diào)頻率估計(jì)的CRLB分析

根據(jù)(6)式可以得到調(diào)頻率估計(jì)的CRLB為

(10)

式中：τ=ΔK/fs，即短時(shí)窗移動(dòng)步進(jìn)對(duì)應(yīng)時(shí)長(zhǎng)。

在輸入信噪比一定的情況下，取δ=0.5，移動(dòng)步進(jìn)取窗長(zhǎng)的一半時(shí)，調(diào)頻率估計(jì)的CRLB與滑窗次數(shù)、信號(hào)時(shí)長(zhǎng)和FRFT的最優(yōu)階數(shù)有關(guān)。調(diào)頻率的CRLB表達(dá)式為

(11)

式中：Tx為信號(hào)時(shí)長(zhǎng)，為便于分析，取Tx=1 s，此時(shí)調(diào)頻率估計(jì)的CRLB與滑窗次數(shù)、FRFT的最優(yōu)階數(shù)有關(guān)，分析方法與初始頻率估計(jì)的CRLB分析相同。分別求(11)式關(guān)于K和α的偏導(dǎo)，有

(12)

(13)

當(dāng)K>1時(shí)，(12)式的值恒小于0，故調(diào)頻率的估計(jì)精度隨滑窗次數(shù)的增加而提高，即CRLB隨滑窗次數(shù)的增加而減小。(13)式的正負(fù)由cosα決定，變化規(guī)律與初始頻率估計(jì)的CRLB情況一致。

對(duì)調(diào)頻率估計(jì)的CRLB仿真結(jié)果歸一化后如圖5所示。同樣地，圖5(a)的階數(shù)取值為[1，2]。由圖5(a)可以明顯看出，調(diào)頻率估計(jì)的CRLB隨滑窗次數(shù)增加而減小。圖5(b)的階數(shù)取值為[0,2]。由圖5(b)可以看出，調(diào)頻率估計(jì)的CRLB隨階數(shù)的增加先增大、后減小，階數(shù)為1時(shí)估計(jì)性能最低。這一結(jié)論與理論分析結(jié)果一致。

圖5 歸一化調(diào)頻率估計(jì)的CRLB

3 實(shí)驗(yàn)驗(yàn)證

3.1 仿真實(shí)驗(yàn)

仿真信號(hào)采用LFM信號(hào)，初始頻率為600 Hz，調(diào)頻率為800 Hz/s，采樣率為4 096 Hz，脈寬為1 s，噪聲為高斯白噪聲，信噪比變化區(qū)間為-12～0 dB.經(jīng)過(guò)5 000次蒙特卡洛仿真，得到估計(jì)的調(diào)頻率均值及其均方根誤差如圖6所示；估計(jì)的初始頻率均值及其均方根誤差如圖7所示。圖6和圖7中，VWSF-ISTFT表示利用插值STFT估計(jì)信號(hào)瞬時(shí)頻率，然后通過(guò)VWSF估計(jì)參數(shù)；VWSF-ISTFRFT表示利用插值STFRFT估計(jì)信號(hào)瞬時(shí)頻率，然后通過(guò)VWSF估計(jì)參數(shù)。

圖6 調(diào)頻率估計(jì)

從圖6(a)可以看出，輸入信噪比為-11 dB時(shí)，VWSF-ISTFRFT方法估計(jì)的調(diào)頻率已經(jīng)接近真值，與VWSF-ISTFT方法相比，所要求的輸入信噪比降低了1 dB；相比于CLSF-ISTFRFT方法，VWSF-ISTFRFT方法所要求的信噪比降低了約3 dB.從圖6(b)可以看出，輸入信噪比為-11 dB時(shí)，VWSF-ISTFRFT方法估計(jì)的調(diào)頻率精度已經(jīng)接近CRLB;與VWSF-ISTFT方法相比，相同估計(jì)精度下，VWSF-ISTFRFT方法要求的輸入信噪比降低了1 dB；達(dá)到相同的估計(jì)精度時(shí)，VWSF方法同CLSF方法相比，要求的輸入信噪比降低了約3 dB.

圖7 初始頻率估計(jì)

從圖7(a)可以看出，輸入信噪比為-11 dB時(shí)，VWSF-ISTFRFT方法估計(jì)的初始頻率已經(jīng)接近真值，與VWSF-ISTFT方法相比信噪比降低了1 dB；相比于CLSF-ISTFRFT方法，VWSF-ISTFRFT方法所要求的信噪比降低了約3 dB.從圖7(b)可以看出，相同估計(jì)精度下，VWSF-ISTFRFT方法估計(jì)初始頻率所要求的輸入信噪比較VWSF-ISTFT方法降低了1 dB，VWSF方法所要求的輸入信噪比較CLSF方法降低了約3 dB.

針對(duì)被動(dòng)聲納檢測(cè)到的脈沖信號(hào)參數(shù)估計(jì)問(wèn)題，為了證明本文方法的普適性，在被動(dòng)聲納工作頻率范圍內(nèi)隨機(jī)選取參數(shù)進(jìn)行仿真。根據(jù)文獻(xiàn)[15]可知，現(xiàn)有的被動(dòng)聲納工作的大致頻率范圍為0.1～10 kHz，調(diào)頻信號(hào)通常作為窄帶信號(hào)應(yīng)用，脈寬由幾十毫秒到幾秒，在合理的參數(shù)取值范圍內(nèi)取3組參數(shù)進(jìn)行仿真。信號(hào)均采用LFM信號(hào)，噪聲為高斯白噪聲，信噪比變化范圍為-12～0 dB，其他仿真條件及得到的結(jié)果如下。

1)仿真條件1.采樣率為2 kHz，初始頻率為50 Hz，帶寬為80 Hz，脈寬為4 s，經(jīng)過(guò)5 000次蒙特卡洛仿真，得到調(diào)頻率和初始頻率的估計(jì)誤差如圖8和圖9所示。

圖8 調(diào)頻率估計(jì)誤差(仿真條件1)

圖9 初始頻率估計(jì)誤差(仿真條件1)

2)仿真條件2.采樣率為8 kHz，初始頻率為2 kHz，帶寬為100 Hz，脈寬為1 s，經(jīng)過(guò)5 000次蒙特卡洛仿真，得到調(diào)頻率和初始頻率的估計(jì)誤差如圖10和圖11所示。

圖10 調(diào)頻率估計(jì)誤差(仿真條件2)

圖11 初始頻率估計(jì)誤差(仿真條件2)

3)仿真條件3.采樣率為20 kHz，初始頻率為8 kHz，帶寬為400 Hz，脈寬為0.4 s，經(jīng)過(guò)5 000次蒙特卡洛仿真，得到調(diào)頻率和初始頻率的估計(jì)誤差如圖12和圖13所示。

圖12 調(diào)頻率估計(jì)誤差(仿真條件3)

圖13 初始頻率估計(jì)誤差(仿真條件3)

從上述3個(gè)不同條件的仿真結(jié)果中均可以看出，輸入信噪比為-12 dB時(shí)，VWSF-ISTFRFT方法仍有較高的估計(jì)精度，與基于STFT方法相比，VWSF-ISTFRFT方法的初始頻率和調(diào)頻率估計(jì)精度更高，估計(jì)性能更優(yōu)；達(dá)到相同的估計(jì)精度時(shí)，VWSF-ISTFRFT方法所要求的輸入信噪比更低。

3.2 海試驗(yàn)證

2017年12月在某地海試，發(fā)射信號(hào)為L(zhǎng)FM信號(hào)，采樣率為64 kHz，脈寬為128 ms，歸一化后初始頻率為0.96，信噪比為-8 dB.分別采用VWSF-ISTFT、CLSF-ISTFRFT和VWSF-ISTFRFT方法3種方法估計(jì)的瞬時(shí)頻率如圖14所示。從圖14中可以看出，VWSF-ISTFRFT方法估計(jì)的瞬時(shí)頻率最接近實(shí)際頻率。將初始頻率和調(diào)頻率歸一化，估計(jì)的參數(shù)如表1所示。從表1中可以看出，VWSF-ISTFRFT方法估計(jì)的值最接近真實(shí)值，CLSF-ISTFRFT方法估計(jì)參數(shù)誤差最大。

圖14 LFM信號(hào)瞬時(shí)頻率估計(jì)

從海試數(shù)據(jù)的處理結(jié)果中可以看出，VWSF-ISTFRFT方法估計(jì)的瞬時(shí)頻率最接近實(shí)際頻率；從表1中數(shù)據(jù)可以看出，VWSF-ISTFRFT方法估計(jì)的參數(shù)精度最高，與CLSF方法相比，VWSF方法的參數(shù)估計(jì)精度由百分之一量級(jí)提高到千分之一量級(jí)，提高了一個(gè)數(shù)量級(jí)。

表1 估計(jì)參數(shù)及相對(duì)誤差

注：括號(hào)內(nèi)的值為參數(shù)估計(jì)相對(duì)誤差。

4 結(jié)論

本文針對(duì)低信噪比下LFM信號(hào)參數(shù)估計(jì)問(wèn)題，利用Rife插值STFRFT方法估計(jì)LFM信號(hào)的瞬時(shí)頻率，通過(guò)VWSF估計(jì)初始頻率和調(diào)頻率，并推導(dǎo)分析了初始頻率和調(diào)頻率估計(jì)的CRLB影響因素。通過(guò)仿真分析和海試數(shù)據(jù)驗(yàn)證，得到主要結(jié)論如下：

1)插值STFRFT方法估計(jì)的LFM信號(hào)瞬時(shí)頻率較插值STFT精度更高，有著更優(yōu)的估計(jì)性能，估計(jì)精度相同時(shí)，對(duì)輸入信噪比的要求前者較后者更低。

2)利用瞬時(shí)頻率進(jìn)行參數(shù)估計(jì)時(shí)，VWSF方法可以消除插值異常點(diǎn)的影響，與CLSF方法相比，得到的參數(shù)估計(jì)精度更高。

3)初始頻率和調(diào)頻率估計(jì)的CRLB受滑窗次數(shù)和FRFT階數(shù)影響，隨滑窗次數(shù)的增加而減小，隨FRFT階數(shù)的增加先增大、后減小，階數(shù)為1時(shí)，CRLB達(dá)到最大值，此時(shí)參數(shù)的估計(jì)性能最低。