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

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

0 引 言

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

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

1 Chirp信號特征

Chirp信號的時域表達(dá)式為

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

Chirp信號的瞬時頻率表達(dá)式為

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

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

圖1 Chirp信號時域表示圖

2 改變窗長相位差估計調(diào)頻斜率

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

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

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

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

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

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

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

2.2 最小二乘法進(jìn)行數(shù)據(jù)擬合得到調(diào)頻斜率

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

圖2 最小二乘擬合與直接擬合的NMSE與SNR的關(guān)系曲線圖

NMSE計算公式為

3 二分法FRFT

3.1 FRFT估計Chirp信號參數(shù)理論

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

Chirp信號經(jīng)過p階的FRFT為

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

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

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

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

圖3 Chirp信號的FRFT三維圖

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

3.2 FRFT階次搜索區(qū)間設(shè)置

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

因粗估計所得調(diào)頻斜率有一定的誤差，為了進(jìn)一步精確估計出Chirp信號的參數(shù)，需要設(shè)置一個階次搜索區(qū)間以及搜索階次步長。區(qū)間設(shè)置越大且階次步長設(shè)置越小，則進(jìn)行FRFT的次數(shù)就越多，使得計算量增大，并且會使FRFT的結(jié)果受到頻譜偽峰的干擾；區(qū)間設(shè)置太小且搜索階次步長太大，會導(dǎo)致估計精度下降，甚至?xí)z測不到峰值而導(dǎo)致估計參數(shù)錯誤。故需要設(shè)置一個合理的區(qū)間及步長。

3.3 二分法搜索

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

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

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

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

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

4 算法設(shè)計流程

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

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

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

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

設(shè)計算法流程圖如圖5所示。

圖5 本文算法流程圖

5 仿真與分析

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

表1 Chirp信號參數(shù)

圖6 調(diào)頻斜率粗估計的NMSE與SNR的關(guān)系

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

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

圖7 調(diào)頻斜率估計值的NMSE與SNR的關(guān)系

圖8 初始頻率估計值的NMSE與SNR的關(guān)系

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

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

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

6 結(jié)束語

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