基于相位差法和FRFT的Chirp信號參數估計算法

(重慶郵電大學 通信與信息工程學院，重慶 400065)

0 引 言

啁啾(Chirp)信號也稱為線性調頻信號，在信號處理領域有著廣泛的應用，如在超聲波成像[1]、雷達[2]、水下探測[3]以及電力線載波通信[4]等領域。Chirp信號的檢測與估計在各個應用領域中具有重要的意義和實用價值，近年來越來越受重視。最大似然估計[5]對含有噪聲的Chirp信號檢測具有較好的性能，但其復雜度較高，不適合工程實現；基于短時傅里葉變換(Short Time Fourier Transform, STFT)[6]的估計方法，時頻域的分辨率與選擇的窗長密切相關；分數階傅里葉變換(Fractional Fourier Transform, FRFT)隨著Ozaktas H M提出的采樣型算法[7]出現，可利用快速傅里葉變換(Fast Fourier Transform, FFT)實現，其運算復雜度與FFT相當且計算精度高，被廣泛應用于Chirp信號的檢測與估計。通過FRFT進行信號檢測其實質是一個二維搜索的過程。文獻[8]提出通過FFT對Chirp信號的調頻率進行粗估計，再通過FRFT進行精估計，該方法降低了信號檢測的計算量，但其抗噪聲性能較差；文獻[9]提出通過STFT對Chirp信號參數進行粗估計并聯合FRFT的方法，其計算量依舊很大。

基于上述問題，本文提出首先通過改變窗長相位差法來預估待測信號的調頻斜率，通過調頻斜率粗估計值將FRFT搜索階次限定在一個極小的區間內，然后再通過二分法在這個區間搜索待檢測信號的最佳階次FRFT，從而求得Chirp信號的初始頻率和調頻斜率，為Chirp信號參數估計提供了一種新的思路。

1 Chirp信號特征

Chirp信號的時域表達式為

式中：α(t)為Chirp信號的包絡，通常選取矩形脈沖； j為虛數單位；f0為Chirp的初始頻率，單位為Hz；t為時間軸，單位為s；k為Chirp信號的調頻斜率，單位為Hz/s；T為脈沖寬度，又稱掃頻時間，單位為s。

Chirp信號的瞬時頻率表達式為

由式(2)可知，在一個脈沖寬度內其瞬時頻率與時間成線性關系。B=kt為Chirp信號的調頻帶寬。

圖1所示為包絡為1、帶寬為3.6 MHz、調頻斜率為1.08×1012Hz/s、采樣頻率為76.8 MHz和時寬為20 μs的Chirp信號時域表示圖。

圖1 Chirp信號時域表示圖

2 改變窗長相位差估計調頻斜率

2.1 改變窗長相位差估計瞬時頻率

改變窗長相位差法[10]估計頻率的核心思想是，將長度為N的待測信號x1(n)先進行FFT得到X1(n)，通過X1(n)中的最大值得到實際峰值對應的譜線號k1max，并計算出該點的相位φ1=arg(X1(k1max))；再利用該段長度為前M=b·N(0

式中，n=0,1,…,M。

由φ2=arg(X2(k2max))計算出第2段數據中x2(k2max)的相位，則φ1與φ2的相位差Δφ為

由上式可知，歸一化的頻率估計值為

則最終得到的該段FFT峰值所對應實際估計頻率為

式中，fs為信號采樣頻率。

2.2 最小二乘法進行數據擬合得到調頻斜率

通過改變窗長相位差得到待測信號不同時段的瞬時頻率，再利用最小二乘法[11]對不同時段的瞬時頻率進行擬合，從而求得檢測信號的調頻斜率。圖2所示為采用最小二乘法進行數據擬合求得的調頻斜率及未經最小二乘直接通過前后兩段數據瞬時頻率擬合求得的調頻斜率，再通過1 000次Monte Carlo后取平均得到的歸一化均方誤差(Normalized Mean Square Error，NMSE)與信噪比(Signal to Noise Ratio，SNR)的關系圖。

圖2 最小二乘擬合與直接擬合的NMSE與SNR的關系曲線圖

NMSE計算公式為

3 二分法FRFT

3.1 FRFT估計Chirp信號參數理論

FRFT處理Chirp信號的原理是通過一組新的Chirp正交基函數對檢測信號進行重構。如果Chirp信號的調頻斜率與其中一個正交基相同，就會在分數階傅里葉域出現沖激特性，從而可以對Chirp信號進行參數估計和檢測[12]。

Chirp信號經過p階的FRFT為

式中：Kp(u,t)為FRFT的核函數；u為分數階傅里葉域的坐標軸；α=pπ/2為旋轉角度，p為FRFT的階次。當Chirp信號旋轉到α=-arccotk時，Chirp信號在FRFT域里表現出Sinc函數形式的沖激，故可以利用FRFT對Chirp信號進行參數估計：

(1) 通過設置旋轉角度α，連續地對待檢測的信號進行FRFT，從而獲得在(α,u)平面上以信號能量為縱軸的二維分布。

(2) 在(α,u)平面上搜索峰值點，即當k=-cotα時，Chirp信號的能量高度聚集在u=f0/cscα，此時為峰值點，此時在二維平面內搜索：

Chirp信號在各階次的FRFT域三維波形和最佳階次下的FRFT域波形分別如圖3和4所示。

圖3 Chirp信號的FRFT三維圖

圖4 Chirp信號在最佳階次下的FRFT波形

3.2 FRFT階次搜索區間設置

通過改變窗長相位差法將待檢測的Chirp信號進行粗估計可以得到調頻斜率k1,利用粗估計調頻斜率k1計算所對應的FRFT階次p1：

因粗估計所得調頻斜率有一定的誤差，為了進一步精確估計出Chirp信號的參數，需要設置一個階次搜索區間以及搜索階次步長。區間設置越大且階次步長設置越小，則進行FRFT的次數就越多，使得計算量增大，并且會使FRFT的結果受到頻譜偽峰的干擾；區間設置太小且搜索階次步長太大，會導致估計精度下降，甚至會檢測不到峰值而導致估計參數錯誤。故需要設置一個合理的區間及步長。

3.3 二分法搜索

二分法也稱折半查找法，是對一組有序數據進行搜索的算法，其是在一個區間內通過不斷取中間值使得區間范圍逐漸縮小，從而使得區間的中間值逐漸逼近實際值。二分法體現了無限逼近的極限思想，也可以說二分法是一種區間迭代的數值算法[13]。

用粗估計的調頻斜率設定一個調頻斜率區間范圍，將二分法運用到FRFT中，然后在對應的FRFT階次范圍[pbegin,pend]內用二分法搜索，每經過一次搜索，階次范圍就縮小一半，若滿足(pend-pbegin)

步驟1:令pmid=(pbegin+pend)/2，分別計算出階次pbegin、pmid和pend所對應的FRFT的最大值Fbegin、Fmid和Fend。

步驟2：比較Fbegin、Fmid和Fend的大小，若Fbegin最大，則令pend=pmid，Fend=Fmid，pmid=(pbegin+pend)/2，重新計算Fmid；若Fend最大，則令pbegin=pmid，Fbegin=Fmid，同時令pmid=(pbegin+pend)/2，重新計算Fmid；若Fmid最大，為了減小誤差影響，則需要令pbegin=(pbegin+pmid)/2、pend=(pend+pmid)/2，重新計算Fbegin和Fend。

步驟3：若(pend-pbegin)>r，則重復執行Step2。

4 算法設計流程

對于一段待檢測的Chirp信號，其頻率是呈線性變換的。利用改變窗長相位差法和FRFT聯合檢測的方法設計步驟如下：

(1)將待檢測長度為N的信號分成i=N/M段長度為M的數據，分別利用改變窗長相位差法求出各段待檢測信號所對應的瞬時頻率fi；

(2)利用最小二乘法將得到的瞬時頻率進行線性擬合，計算出該段Chirp信號的調頻斜率k′；

(3)根據粗估計的調頻斜率k′，將FRFT搜索階次限定在一個極小的區間[pbegin,pend]；

設計算法流程圖如圖5所示。

圖5 本文算法流程圖

5 仿真與分析

在本節中，我們主要通過NMSE來評估Chirp信號參數估計的性能，Chirp信號的仿真參數如表1所示。

表1 Chirp信號參數

圖6 調頻斜率粗估計的NMSE與SNR的關系

為了更好地對本文算法進行性能分析，將使用本文提出的改變窗長相位差法和文獻[8]提出的將信號進行FFT后進行5點平滑再取模平方處理方法、文獻[9]使用的STFT算法、文獻[14]提出的Rife算法以及文獻[15]提出的修正的M-Rife算法對Chirp信號調頻斜率進行粗估計后再進行比較。圖6所示為SNR為[-20,5] dB、采用1 000次Monte Carlo仿真取平均得到的Chirp信號調頻斜率粗估計的NMSE仿真圖。由圖可知，FFT 5點平滑處理方法的性能最差，其抗噪聲性能只能達到-7 dB，且估計的準確度受噪聲影響較大；Rife算法是通過FFT之后在最大值和次大值之間進行插值，從而求得瞬時頻率，該算法比FFT 5點平滑處理方法具有更高的估計精度；M-Rife算法是修正的Rife算法，通過對Rife算法估計之后的頻率進行頻譜搬移，再進行一次Rife算法，以提高Rife算法估計頻率的準確度，其估計精度要比Rife算法好，但其運算量要比Rife算法大；STFT算法的估計性能在SNR>-6 dB時比M-Rife算法要好。相較于上述4種算法，改變窗長相位差法的估計準確度是最好的。改變窗長相位差法通過一次FFT和一次單點DFT運算，可以在SNR>-10 dB時較準確地對Chirp信號參數進行估計，且其在SNR=-8 dB時相較于M-Rife算法，其NMSE提升了9 dB，故利用改變窗長相位差法來對Chirp信號進行粗估計是一種很好的選擇。

圖7和8分別是本文提出的通過改變窗長相位差法進行粗估計再通過二分法FRFT進行精確估計的算法、采用FFT 5點平滑處理方法進行粗估計聯合FRFT進行精估計的算法以及通過STFT聯合FRFT算法在SNR為[-20,5] dB、采用1 000次Monte Carlo仿真取平均得到的Chirp信號調頻斜率估計值和初始頻率估計值的NMSE仿真圖。其中改變窗長相位差的點數及STFT窗長設置為256，搜索階次步長為0.000 1。由圖可知，本文算法與通過STFT聯合FRFT算法的性能接近，在SNR>-9 dB時，估計Chirp信號調頻斜率以及初始頻率可以接近克拉美羅下界(Cramer-Rao Lower Bound，CRLB),并且本文算法的抗噪聲性能相比于采用FFT 5點平滑處理方法后再采用FRFT的算法提升了2 dB。

圖7 調頻斜率估計值的NMSE與SNR的關系

圖8 初始頻率估計值的NMSE與SNR的關系

為了進一步評估本文通過改變窗長相位差法聯合二分法FRFT算法對Chirp信號參數進行估計的計算復雜度和實用性，將本文算法和FFT聯合FRFT算法以及STFT聯合FRFT算法的時間復雜度進行對比與分析。設待檢測信號的長度為N,則本文算法中在對信號調頻斜率進行粗估計時，涉及到N/M次FFT和N/M次單點DFT運算，故改變窗長相位差法的復雜度為N/M(O(Mlog2M)+O(M))=O(Nlog2M)。因進行一次FRFT復雜度為O(Nlog2N),當階次搜索區間寬度為Δp，搜索精度為r時,通過二分法搜索的復雜度為O(log2(Δp/r))，故改變窗長相位差聯合二分法FRFT算法的復雜度為O((Nlog2N)log2(Δp/r)+Nlog2M)=O((Nlog2N)log2(Δp/r))；文獻[8]通過FFT 5點平滑處理方法進行粗估計包含一次N點FFT運算及5(N-1)次加法和N-1次除法，則其時間復雜度為O(Nlog2N+5(N-1)+(N-1))=O(Nlog2N)，再通過FRFT以搜索精度r進行局部等步長搜索的復雜度為O((Δp/r)Nlog2N)，故FFT 5點平滑和FRFT聯合算法的時間復雜度為

O((Δp/r)Nlog2N+Nlog2N)=O((Δp/r)Nlog2N)；

通過STFT和FRFT對Chirp信號進行參數估計時涉及到N-M次M點的FFT以及Δp/r次FRFT，故其復雜度為O((N-M)Mlog2M+(Δp/r)Nlog2N)，可化簡為O((Δp/r)Nlog2N)。由此可知，本文算法的計算復雜度要比上述兩種方法小很多，當SNR=-7 dB、Δp=0.2、r=0.000 1時，本文算法只需要進行log22 000≈11次FRFT就可以估計出Chirp信號的參數。而通過FFT 5點平滑處理方法進行粗估計，再通過局部等步長窮舉搜索的話，需要進行2 000次FRFT。同理，通過STFT聯合FRFT的方法也需要進行2 000次FRFT。故可知本文通過改變窗長相位差法進行粗估計和二分法FRFT對Chirp信號進行精估計的算法極大地降低了計算復雜度。

6 結束語

本文針對FRFT對Chirp信號進行盲檢測時存在運算量大的問題,提出了一種通過改變窗長相位差法聯合二分法FRFT對Chirp信號參數進行估計的方法。在該方法中，先通過改變窗長相位差對待檢測的Chirp信號進行粗估計，再通過二分法進行搜索和FRFT對Chirp信號參數進行精估計，可以在低SNR下對Chirp信號進行參數估計，同時也降低了計算復雜度。綜合估計精度和時間復雜度來說，本文所提算法對Chirp信號參數估計具有很好的性能。