三角波調頻信號參數快速精確估計方法

，

(機電動態控制重點實驗室，陜西 西安 710065)

0 引言

近年來，國內外學者針對三角波調頻信號的參數估計問題進行了大量的研究。文獻[1]和文獻[2]提出一種基于 Radon-Ambiguity (RAT)變換與分數階傅里葉變換相結合的參數估計方法，該方法相比于采用二維搜索的Radon-Wigner變換，以較小的運算量實現了對信號初始頻率的估計，但利用RAT變換估計信號調頻斜率的運算量非常大。文獻[3]和文獻[4]采用高階累積量與多相濾波器組相結合的方法實現對信號的參數估計，該方法無需先驗知識且能較好地抑制高斯噪聲，然而該方法的估計精度會隨著多相濾波器組數的增加而下降。文獻[5]提出一種分數階傅里葉變換與聚類分析相結合的參數估計方法，該方法克服了信號尖峰必須高于噪聲幅度的限制，可以實現低信噪比條件下信號的檢測與參數估計，但該方法要求截獲信號的起點位于調制周期的起點，因此不適合無先驗知識的信號參數估計。文獻[6]提出一種基于短時傅里葉變換，解線性調頻和最大似然估計相結合的參數估計方法，估計精度高，但引入最大似然估計校正精度，運算量非常大。文獻[7]和文獻[8]提出了基于最佳頻率分辨率窗寬的高斯窗短時傅里葉變換與相關解調相結合的參數估計方法，該方法實時性強，但并未使用最佳高斯窗進行短時傅里葉變換，且在進行相關解調時粗估了信號線性調頻段的起點和終點，導致該方法參數估計精度不足。本文針對Radon-Ambiguity變換、解線性調頻、最大似然估計及短時傅里葉變換等三角波調頻信號參數估計方法均無法兼顧運算量小和估計精度高的問題，提出了基于信號特征的最佳高斯窗短時傅里葉變換與分段相關解調相結合的參數估計方法。

1 三角波調頻信號參數估計方法

1.1 基于短時傅里葉變換與最大似然估計的三角波調頻信號參數估計方法

對于三角波調頻信號的參數估計問題，文獻[6]提出了基于短時傅里葉變換和最大似然估計的參數估計方法，該方法首先對信號解線性調頻并做短時傅里葉變換，求出中心頻率，調制帶寬和調制周期的粗略估計值，最后在時頻平面進行最大似然搜索，求出上述參數的精確估計值。方法流程如圖1所示。

1.2 基于短時傅里葉變換與相關解調的三角波調頻信號參數估計方法

對于三角波調頻信號的參數估計問題，文獻[7]提出了基于最佳頻率分辨率窗寬的高斯窗短時傅里葉變換與相關解調相結合的參數估計方法，該方法首先根據時頻分析的目的，選擇高斯窗作為短時傅里葉變換的窗函數，其次使用基于瞬時頻率的窗寬倍增搜索方法求出最佳頻率分辨率窗寬，然后對三角波調頻信號進行短時傅里葉變換，求得信號的中心頻率和調制帶寬，再使用相關解調求得信號的調制周期。方法流程如圖2所示。

高斯窗gs(t)的定義如下：

(1)

式(1)中，σ為高斯窗的標準差，σ的取值將對高斯窗的時頻分辨率產生影響。由Heisenberg不確定性原理可知，提高高斯窗的時頻分辨率，可以減小信號的參數估計誤差。1.2節方法雖然求得了最佳高斯窗寬，但在對算法進行Matlab仿真時選擇了Matlab自帶的高斯窗函數，該函數使用固定的標準差σ1=2.5進行短時傅里葉變換，并未考慮信號波形特性與高斯窗標準差的關系，因此選擇固定標準差σ1=2.5的高斯窗短時傅里葉變換并不能使窗內信號的時頻聚集性達到最佳，這將導致參數估計精度不足。所以選擇固定標準差σ1=2.5的高斯窗并非最佳高斯窗，要想使用最佳高斯窗，不僅需要求得最佳高斯窗寬，還需要根據信號波形特性求出高斯窗最佳時頻聚集性標準差。

2 基于信號特征的三角波調頻信號參數快速精確估計方法

本文提出了基于信號特征的最佳高斯窗短時傅里葉變換與分段相關解調相結合的三角波調頻信號參數估計方法。該方法首先選擇固定的標準差σ1=2.5，使用1.2節提出的方法求得三角波調頻信號的調頻斜率k，根據k求出高斯窗最佳時頻聚集性標準差σbest，進而選擇σbest并使用1.2節提出的基于最佳頻率分辨率窗寬的高斯窗短時傅里葉變換求出信號的中心頻率fc和調制帶寬ΔF，其次利用分段相關解調確定三角波調頻信號嚴格線性調頻段的起點和終點，最后在嚴格線性調頻段內使用相關解調估計出信號的調制周期T。方法流程如圖3所示。

2.1 基于信號特征的最佳高斯窗短時傅里葉變換

基于信號特征的高斯窗最佳時頻聚集性標準差推導步驟如下：

給定線性調頻信號xLFM(t)：

xLFM(t)=ej2π(f0t+0.5kt2)

(2)

引入變量β，對式(2)配平方，整理得式(3)，

xLFM(t)=ej2π[f0t+0.5k(t-β)2-0.5kβ2+ktβ]

(3)

任選窗函數L(t)，對式(3)進行短時傅里葉變換，得式(4)，

(4)

對式(4)進行整理，得式(5)

(5)

定義函數γ(t)，

(6)

定義函數ρ(f)，

(7)

定義加窗的線性調頻信號λ(t)：

λ(t)=ejπkt2L(t)

(8)

由ρ(f)的定義可知，函數ρ(f)是λ(t)的傅里葉變換，將式(8)代入式(6)可得，

(9)

根據傅里葉變換的時移和頻移性質，整理式(9)可得，

γ(t)=e-j2πf′t·FT{ej2π[(f0+kt)t′]λ(t′)}=
{e-j2πf′t·FT[λ(t′)]}|f′=f-f0-kt

(10)

將式(7)代入式(10)可得，

γ(t)=e-j2π(f-f0-kt)t·ρ(f-f0-kt)

(11)

將式(11)代入式(5)可得，

STFTx(t,f)=e-jπkt2γ(t)=
e-jπkt2e-j2π(f-f0-kt)t·ρ(f-f0-kt)

(12)

則線性調頻信號xLFM(t)的譜圖SPECx(t,f)為：

SPECx(t,f)=|e-jπkt2e-j2π(f-f0-kt)t·
ρ(f-f0-kt)|2=|ρ(f-f0-kt)|2

(13)

由式(13)可知，譜圖SPECx(t,f)相當于λ(t)頻譜的模平方ρ(f)2沿直線f=f0+kt的排列。因此ρ(f)2的寬度決定了線性調頻信號xLFM(t)譜圖的時頻聚集性。

將式(1)定義的高斯窗函數gs(t)代入式(7)的ρ(f)中，可得，

(4)

由于信號x(t)的α轉角分數階傅里葉變換公式為：

(15)

其中，

(16)

因此，高斯窗函數gs(t)的α轉角分數階傅里葉變換為：

(17)

比較式(17)與式(14)可知，當cotα=2πk且ucscα=2πf時，

(18)

對式(18)進行整理得，

(19)

取cotα=2πk且ucscα=2πf，進行三角函數代換，整理式(19)并代入式(18)可得，

(20)

整理式(20)可得，

(21)

整理式(21)可得，

(22)

由式(22)可知，ρ(f)2為高斯型函數，因此其寬度由z(σ)決定：

(23)

z(σ)min=k/2π

(24)

2.2 分段相關解調

分段相關解調具體步驟如下：

1)將區間[tmin,tmax]內的信號均分為r段，為計算方便，取r為4的整數倍，各段信號編號依次為1,2,…，r，為對每段信號進行相關解調，求得每段信號的調頻斜率依次為ki,1≤i≤r。

2)計算各段信號調頻斜率的均值：因區間[tmin,tmax]端點處的信號為非嚴格線性調頻信號，為減小誤差，使用區間中間編號為第0.25r+1小段到編號為第0.75r小段的信號計算每段信號調頻斜率的均值ka：

(25)

3)計算各段信號調頻斜率ki與調頻斜率均值ka的絕對值差Δki：

Δki=ki-ka,1≤i≤20

(26)

綜上，本文提出了基于信號特征的最佳高斯窗短時傅里葉變換與分段相關解調相結合的參數估計方法，該方法根據三角波調頻信號的波形特性，通過理論推導給出了高斯窗短時傅里葉變換最佳時頻聚集性標準差，選擇該標準差進行短時傅里葉變換對信號的中心頻率和調制帶寬進行了參數估計；使用分段相關解調確定了信號線性調頻段的起點和終點，求出了信號的嚴格線性調頻區間，在嚴格線性調頻區間內進行相關解調，估計出了信號的調制周期。

3 運算量分析與仿真驗證

3.1 運算量分析

假定信號線性調頻段內信號點數為M點，采樣信號共包含c個線性調頻段，以最佳頻率分辨率窗寬qbest進行短時傅里葉變換，設步長為a，則初次估計時，短時傅里葉變換的運算量V1為：

(27)

初次估計時相關解調的運算量P1為：

(28)

二次估計時，依然使用最佳頻率分辨率窗寬qbest進行短時傅里葉變換，則二次估計時短時傅里葉變換的運算量V2為：

V2=V1

(29)

假定將近似線性調頻區間分為r段，令U=M/(2r)，則分段相關解調的運算量P2為：

(30)

因為窗寬倍增搜索方法和相關解調的運算量遠小于短時傅里葉變換本身的運算量，而分段相關解調的運算量又小于相關解調的運算量，因此若采樣信號包含c個線性調頻段，則本文方法的總運算量V為：

(31)

對比1.1節方法，其運算量主要為解線性調頻，短時傅里葉變換和最大似然估計的運算量，依然假定線性調頻段內信號點數為M點，采樣信號包含c個線性調頻段，則解線性調頻的運算量X為：

X=M(4+2lbM)

(32)

以窗寬q進行短時傅里葉變換，設步長為b，則短時傅里葉變換的運算量Y為：

(33)

最大似然估計的運算量Z為：

Z=M(M+4+2lbM)

(34)

1.1節方法總的運算量R為：

(35)

由以上分析可知，1.2節方法和本文方法的運算量均由M決定，而1.1節方法的運算量由M2決定，因此本文方法運算量與1.2節方法的運算量在同一數量級，遠小于1.1節方法的運算量，適用于信號的實時處理。

3.2 實驗仿真

為驗證本文方法的參數估計精度，選擇三組三角波調頻信號，在相同條件下，分別使用1.1節方法，1.2節方法和本文方法對信號的中心頻率，調制帶寬，調制周期進行參數估計。

給定三組三角波調頻信號：

信號1：中心頻率fc=270 MHz，調制帶寬ΔF=40 MHz，調制周期T=20 μs。

信號2：中心頻率fc=290 MHz，調制帶寬ΔF=30 MHz，調制周期T=30 μs。

信號3：中心頻率fc=310 MHz，調制帶寬ΔF=20 MHz，調制周期T=40 μs。

仿真條件設定：三組信號長度均為5個周期，采樣率均為4倍采樣，仿真時均加入高斯白噪聲，令信噪比從0～20 dB，每個信噪比進行1 000次蒙特卡洛實驗，均采用歸一化均方誤差(NRMSE)作為參數估計精度的誤差衡量標準。

圖4分別給出了不同信噪比下，采用1.1節方法，1.2節方法以及本文方法時，信號1、信號2、信號3中心頻率的估計誤差；圖5分別給出了不同信噪比下，采用1.1節方法，1.2節方法以及本文方法時，信號1、信號2、信號3調制帶寬的估計誤差。圖6分別給出了不同信噪比下，采用1.1節方法，1.2節方法以及本文方法時，信號1、信號2、信號3調制周期的估計誤差。

以上仿真結果表明，在相同信噪比下，本文方法相比于1.1節方法，中心頻率的估計誤差下降了一個數量級，調制帶寬的估計誤差下降了50%左右，調制周期的估計誤差下降了40%左右；本文方法相比于1.2節方法，中心頻率的估計誤差下降了一個數量級，調制帶寬的估計誤差下降了70%左右，調制周期的估計誤差下降了50%左右。

4 結論

本文針對三角波調頻信號的參數估計問題，提出了基于信號特征的最佳高斯窗短時傅里葉變換與分段相關解調相結合的參數估計方法。該方法根據三角波調頻信號的波形特征確定了高斯窗最佳窗寬和最佳時頻聚集性標準差，進而使用短時傅里葉變換估計出了信號的中心頻率和調制帶寬，其次使用分段相關解調確定了信號的嚴格線性調頻段，并估計出了信號的調制周期。理論分析和仿真表明本文方法以較小的運算量大幅降低了三角波調頻信號中心頻率，調制帶寬和調制周期估計的均方誤差，實現了對三角波調頻信號參數的快速精確估計。

參考文獻:

[1]袁偉明，王敏，吳順君.三角波調頻信號的檢測與參數估計[J].電波科學學報，2005，20(5)：594-597.

[2]WANG MSH, CHAIN AK. Linear frequency-modulated signal detection using time-frequency analysis [J]. IEEE Transactions on Signal Processing, 2011, 46(3)： 471-486.

[3]馮志紅，賴濤，趙擁軍.低信噪比下對稱三角線性調頻信號參數估計[J].電波科學學報，2012，27(3)：520-525.

[4]Stank L. Quadratic and higher order time-frequency analysis based on the short time Fourier transform [J]. Signal Processing and its Applications, 2014, 10(14), 581-588.

[5]劉鋒，徐會法，陶然.基于FRFT的對稱三角線性LFMCW信號檢測與參數估計[J].電子與信息學報，2011，33(8)： 364-370.

[6]華旭東，臧小剛，唐斌.一種三角波調頻信號參數估計實時系統的設計[J].信息技術，2007，25(7)： 1009-1018.

[7]燕天，鄒金龍，馬捷.基于短時傅里葉變換與相關解調的參數估計方法[J].探測與控制學報，2016，38(2)： 18-22.

[8]Jones D L. Signal dependent time-frequency analysis using Gaussian window [J]. Applied Signal Processing, 2014,12(5)：379-388.