李景慧，李國軍，葉昌榮，王尊立

(1.重慶郵電大學(xué) 通信與信息工程學(xué)院，重慶 400065；2.重慶郵電大學(xué) 超視距可信信息傳輸研究所，重慶 400065；3.重慶郵電大學(xué) 光電工程學(xué)院，重慶 400065)

0 引 言

在現(xiàn)代數(shù)字通信系統(tǒng)中，同步性能是評估系統(tǒng)性能的關(guān)鍵指標(biāo)，是可靠數(shù)字通信的重要保證。傳統(tǒng)同步方法一般通過檢測位于數(shù)據(jù)幀頭部的同步碼實(shí)現(xiàn)，其性能與同步碼的長度相關(guān)。同步碼太短，在低信噪比條件下將無法被識別；而同步碼太長，又將占用過多的信道資源，影響通信效率。因此，為惡劣環(huán)境下的短波通信系統(tǒng)設(shè)計(jì)更好的同步碼和更穩(wěn)定的檢測算法已經(jīng)成為研究同步的核心問題。實(shí)現(xiàn)同步的主要方法有相關(guān)法[1]、最大似然法(maximum likelihood, ML)[2-3]和似然比檢驗(yàn)(likelihood ratio test, LRT)[4]等3種方法。其中，相關(guān)法具有最低的復(fù)雜度，且最容易實(shí)現(xiàn)[5]，但是在低信噪比條件下的性能最差；最大似然法和似然比檢驗(yàn)效果好于相關(guān)法，但計(jì)算比相關(guān)法復(fù)雜[6]。提高相關(guān)法在惡劣環(huán)境下的性能，降低算法實(shí)現(xiàn)難度，對于提高惡劣環(huán)境下的通信質(zhì)量具有重要意義。

同步碼中使用較多的是偽隨機(jī)碼(PN碼)[7-8]，延遲相關(guān)法和本地相關(guān)法是兩種常用的同步算法。延遲相關(guān)法具有較大的頻率偏移容限，但是在低信噪比和多徑信道下，其自相關(guān)函數(shù)曲線在峰值處變化平緩，無法精確指示同步起始點(diǎn)。本地相關(guān)法具有尖銳的自相關(guān)峰值，準(zhǔn)確度高于延遲相關(guān)法，能精確地指示同步起始點(diǎn)，抗噪聲和多徑能力強(qiáng)，但是對頻偏敏感。

本文提出了一種基于Costas序列的時(shí)頻聯(lián)合同步方法，使用Costas序列作為同步碼，將傳統(tǒng)的滑動(dòng)窗口法擴(kuò)展到時(shí)頻域中，可以同時(shí)估計(jì)時(shí)延和頻偏，具有更好的抗噪聲和抗頻偏性能。

1 Costas序列特性分析

1.1 Costas序列定義及構(gòu)造方法

任意每行每列有且僅有1個(gè)元素等于1、其余元素都為0的n×n階矩陣稱為置換矩陣。Costas序列是一類特殊的置換矩陣，它與自身任意方向的平移副本之間都至多有1個(gè)元素“1”重合，如矩陣A所示(其序列表示為[4,1,6,7,5,8,3,2])。

(1)

(2)

由于置換矩陣A是有限維的，元素ai,j構(gòu)成的序列是有限長的，因此在非循環(huán)差異函數(shù)C(r,s)中，當(dāng)i+r或j+s超出區(qū)間[1,N]時(shí)，取ai+r,j+s為0。非循環(huán)差異函數(shù)滿足(3)式的n×n階置換矩陣稱為n階Costas序列[9]。

(3)

Costas序列構(gòu)造方法[9-11]分為計(jì)算機(jī)窮舉搜索法和基于有限域理論的構(gòu)造法兩種方法。計(jì)算機(jī)窮舉搜索法會窮舉出所有的置換矩陣，并檢測矩陣是否滿足Costas序列的定義，滿足則輸出該Costas序列。計(jì)算機(jī)窮舉搜索法可以輸出滿足所需階數(shù)的所有Costas序列，但由于窮舉法的復(fù)雜度較高，計(jì)算機(jī)窮舉法所能搜索的階數(shù)是有限制的[9]，在構(gòu)造階數(shù)較大的Costas序列時(shí)需要使用基于有限域理論的構(gòu)造方法。基于有限域理論的構(gòu)造法有Welch構(gòu)造法、Golomb構(gòu)造法和Lempel構(gòu)造法[9-10]等多種方法。

本文使用8階Costas序列，計(jì)算機(jī)窮舉搜索法可以滿足使用。如無特別說明，本文中出現(xiàn)的Costas序列均由計(jì)算機(jī)窮舉搜索法生成。

1.2 Costas序列相關(guān)特性分析

對于一個(gè)Costas序列([4,1,6,7,5,8,3,2])，稱[4,1,6,7,5,8,3,2]為其序列表示，(1)式為其矩陣表示。Costas序列的自相關(guān)性體現(xiàn)在其矩陣表示上。為了分析矩陣的相關(guān)特性，定義函數(shù)RAB(a,b)為n×n階矩陣A和B的相關(guān)函數(shù)。計(jì)算相關(guān)函數(shù)的方法如圖1所示，矩陣A不動(dòng)，矩陣B沿x軸和y軸進(jìn)行移動(dòng)，其中a為矩陣B沿x軸移動(dòng)的單位數(shù)(-n+1≤a≤n-1)，b為矩陣B沿y軸移動(dòng)的單位數(shù)(-n+1≤b≤n-1)，則函數(shù)RAB(a,b)的值為矩陣B移動(dòng)后，兩矩陣中“1”單元格重合的個(gè)數(shù)。當(dāng)矩陣A和矩陣B相等時(shí)，稱RAB(a,b)為自相關(guān)函數(shù)[12]。

圖1 矩陣相關(guān)計(jì)算示意圖

圖2是一個(gè)8階Costas序列的自相關(guān)函數(shù)計(jì)算結(jié)果，峰值幅度為其階數(shù)n，非峰值處的最大幅值為1。不同Costas序列做同步碼會表現(xiàn)出不同的性能。圖3和圖4分別為Costas序列1和序列2自相關(guān)函數(shù)示意圖。圖3中，序列1存在“相鄰1”，導(dǎo)致自相關(guān)函數(shù)存在臨峰(序列中相鄰兩個(gè)元素的差值為1的情況稱為“相鄰1”，序列自相關(guān)函數(shù)中與峰值相鄰的不為0的值稱為臨峰)，在同步過程中，臨峰會對同步造成干擾。而圖4 中Costas序列2由于沒有“相鄰1”的存在，自相關(guān)函數(shù)中不存在臨峰，因此其在同步時(shí)受到的干擾會小于Costas序列1。

圖2 8階Costas序列自相關(guān)函數(shù)

圖3 Costas序列1自相關(guān)函數(shù)示意圖

圖4 Costas序列2自相關(guān)函數(shù)示意圖

2 基于Costas序列的MFSK信號時(shí)頻聯(lián)合同步方法

2.1 粗同步方法

Costas序列經(jīng)過FSK調(diào)制后轉(zhuǎn)換到時(shí)頻域，信號在時(shí)頻域的能量分布符合Costas序列自相關(guān)特性。本文通過離散短時(shí)傅里葉變換(discrete short time Fourier transform, DSTFT)將信號映射到時(shí)頻域中，如圖5所示。

圖5 時(shí)頻譜轉(zhuǎn)換框圖

在沒有信道干擾的情況下，x(t)的表達(dá)式為

(4)

(4)式中：A為幅度；g(t)為窗函數(shù)；M為FSK調(diào)制階數(shù)；M與Costas序列的階數(shù)n相同；Ts為碼元時(shí)間；fi為第i個(gè)碼元符號的調(diào)制頻率；θi為第i個(gè)碼元的初始相位。對x(t)采樣后得到x(n)，采樣率為fs。將時(shí)域信號按照碼元時(shí)間Ts進(jìn)行分割，分割后的信號集合為{x1(k),x2(k),…,xM(k)}，其中k為一個(gè)碼元的采樣點(diǎn)數(shù)，k=Ts·fs。分別對{x1(k),x2(k),…,xM(k)}進(jìn)行傅里葉變換可得

Xi(W)=F[xi(k)]

(5)

(5)式中，F(xiàn)[·]表示傅里葉變換。將|Xi(W)|進(jìn)行組合得到信號x(t)的時(shí)頻譜，表示為

[|X1(W)|,|X2(W)|,…,|XM(W)|]=

TF[x(t)]

(6)

(6)式中，TF[·]表示時(shí)頻譜變換，其結(jié)果為矩陣形式，示意圖如圖6所示。

圖6 時(shí)頻譜示意圖

假設(shè)FTi,j為時(shí)頻譜在碼元時(shí)間i、頻率j處的能量值，定義在碼元時(shí)間t、頻率f處的判決量為

(7)

(7)式中：n為Costas序列階數(shù)；Ci為Costas序列第i位的值。如果在時(shí)頻譜范圍內(nèi)判決量大于判決門限，則認(rèn)為信號存在(在一般情況下取最大相關(guān)檢測值的60%—90%作為判決門限)。圖7所示為iturHFMD信道下的判決量仿真結(jié)果，iturHFMD信道參數(shù)如表1所示。在信號嚴(yán)重衰落的情況下，通過(7)式定義的判決量依然可以準(zhǔn)確指示信號起止位置。

圖7 判決量示意圖

表1 信道參數(shù)

定義第一個(gè)大于判決門限的自相關(guān)值坐標(biāo)為[tes,fes]，則tes·Ts為時(shí)間粗估計(jì)值，fes·fs/N為頻率粗估計(jì)值，其中N為傅里葉變換點(diǎn)數(shù)。

滑動(dòng)相關(guān)后的時(shí)間估計(jì)精度為Ts，頻率估計(jì)精度為fs/N=1/NΔT，其中ΔT為采樣間隔，NΔT為碼元時(shí)間。提高時(shí)延估計(jì)精度，需要減小碼元時(shí)間Ts，而Ts的減小會使得頻率估計(jì)精度減小。為了同時(shí)提高時(shí)間估計(jì)精度和頻偏估計(jì)精度，使短時(shí)傅里葉變換中窗函數(shù)的步進(jìn)長度ΔL為Ts/K，其中K為正整數(shù)。同時(shí)使用過采樣(補(bǔ)零)提高頻譜分辨率，過采樣系數(shù)定義為E。K=4時(shí)的窗函數(shù)步進(jìn)示意圖如圖8所示。

圖8 DSTFT步進(jìn)示意圖

圖9為E=2、K=4的參數(shù)條件下，信號x(t)的時(shí)頻譜。此時(shí)的時(shí)間估計(jì)精度為Ts/K，頻率估計(jì)精度為1/ENΔT。通過這種方式提高估計(jì)精度會導(dǎo)致計(jì)算量成倍增加。因此，在適當(dāng)調(diào)整參數(shù)K和E后，使用離散頻譜分析方法進(jìn)行精估計(jì)。

圖9 K=4, E=2的時(shí)頻譜示意圖

2.2 精同步方法

由2.1節(jié)的tes可以獲得一個(gè)碼元時(shí)間的頻譜|Xi(W)|，fes為|Xi(W)|的頻譜峰值坐標(biāo)，示意圖如圖10所示。

圖10 基于二倍過采樣的非線性頻偏估計(jì)示意圖

圖10中，x為經(jīng)過傅里葉變換的頻譜峰值位置，x0為頻偏Δf影響后實(shí)際的頻譜峰值位置。假設(shè)窗函數(shù)為矩形窗，定義其功率譜函數(shù)表達(dá)式為[13]

(8)

根據(jù)窗函數(shù)的重心法則，在n趨近于無窮時(shí)有

(9)

根據(jù)(8)—(9)式可以計(jì)算出x0的位置坐標(biāo)為[13]

(10)

定義k為峰值譜線坐標(biāo)，由(10)式可以推導(dǎo)出頻偏估計(jì)值為

(11)

最終的頻率估計(jì)為

(12)

在估計(jì)時(shí)延之前，使用Δf對信號頻譜進(jìn)行糾正。取一個(gè)碼元時(shí)間的信息進(jìn)行分析，對信號重新進(jìn)行時(shí)頻譜變換。由實(shí)驗(yàn)數(shù)據(jù)可得，在時(shí)頻譜的時(shí)域取8個(gè)數(shù)據(jù)點(diǎn)，數(shù)據(jù)點(diǎn)的值大致服從線性分布，如圖11所示。

圖11 基于1/4步進(jìn)滑動(dòng)窗口的線性時(shí)延估計(jì)示意圖

假設(shè)峰值坐標(biāo)為x，次高值坐標(biāo)為x-1。定義x到x+3表示的直線函數(shù)為Y1(x)，x-4到x-1表示的直線函數(shù)為Y2(x)。則Y1(x)與Y2(x)定義如下。

(13)

數(shù)據(jù)幀實(shí)際開始位置x0可以通過計(jì)算Y1(x)與Y2(x)的交點(diǎn)所得。

即求解

Y1(x)=Y2(x)

(14)

f(x)=2kx+b1-b2

(15)

所得唯一可行解即為x0，因此，時(shí)延精估計(jì)值為

Δt=|x-x0|Ts/K

(16)

(16)式中：Ts為碼元時(shí)間；K為窗函數(shù)步進(jìn)系數(shù)。

2.3 窗函數(shù)影響分析

DSTFT中窗函數(shù)的不同會對同步性能產(chǎn)生影響。令窗函數(shù)步進(jìn)長度ΔL為Ts，假設(shè)f1為碼元1所代表的頻率，f2為碼元2所代表的頻率。當(dāng)存在時(shí)延時(shí)，第一個(gè)碼元符號的部分信息會延后到x2(k)中。

x2(n)=A1sin[2πf1(n-τ)]+A2sin[2πf2(n-τ)]

(17)

(17)式中：A1表示第一個(gè)碼元時(shí)間的幅度；A2表示第二個(gè)碼元時(shí)間的幅度；τ為時(shí)間延遲。此時(shí)對x2(k)進(jìn)行傅里葉變換后，頻率f1所造成的頻譜泄漏現(xiàn)象會極大地干擾頻率f2的主瓣譜線，如圖12所示。這種會嚴(yán)重影響頻偏估計(jì)的干擾，可以通過選擇不同的窗函數(shù)或者全相位傅里葉變換來減小。

圖12 時(shí)間延遲分析

3 仿真驗(yàn)證

仿真過程使用8階Costas序列[1,5,7,4,2,8,3,6]作為同步碼，窗函數(shù)為矩形窗函數(shù)，使用MFSK(M=8)方式對傳輸序列進(jìn)行調(diào)制，采樣率設(shè)置為12 kHz，為2倍過采樣，窗函數(shù)步進(jìn)長度為1/4碼元時(shí)間。根據(jù)國際電信聯(lián)盟給出的Watterson模型建議參數(shù)[14]，選擇中緯度地區(qū)的3種短波信道仿真。3種信道參數(shù)如表1所示。

圖13—圖15給出了本文方法在不同判決門限下虛警率和檢測率的ROC曲線。在AWGN信道條件下，隨著判決門限的逐漸增大，虛警率不斷降低，但檢測率變化很小。在短波信道中，不添加噪聲時(shí)的檢測性能與添加-20 dB噪聲干擾的檢測性能相近。綜上所述，本文同步方法對閾值不敏感，且具有優(yōu)秀的抗噪聲性能，雖在短波信道中的性能有所下降，但依然可以保證其在短波信道的可用性。

圖13 ROC曲線(AWGN信道)

圖14 ROC曲線(iturHF信道)

圖15 ROC曲線(iturHF信道, -20 dB)

以最大相關(guān)檢測值的70%作為判決門限，本文同步方法與傳統(tǒng)PN碼同步方法在iturHFMQ信道下的檢測性能如圖16所示。傳統(tǒng)的PN碼同步方法，在信噪比逐漸降低時(shí)性能出現(xiàn)陡降；本文同步方法在信噪比為-20 dB時(shí)，依然具有可觀的檢測率，可以有效提高惡劣環(huán)境下短波通信的質(zhì)量。

圖16 基于Costas序列同步方法的漏檢率

圖17—圖18為本文同步方法在頻偏影響下的自相關(guān)峰值變化圖，仿真信道為iturHFMQ。由于自相關(guān)函數(shù)是關(guān)于時(shí)間和頻率的參數(shù)，需要從兩個(gè)角度來分析其峰值變化，圖17為頻域角度的自相關(guān)峰值變化圖，圖18為時(shí)域角度的自相關(guān)峰值變化圖。本文同步方法的自相關(guān)峰值幅度主要受頻譜泄漏的影響，當(dāng)頻偏越接近1/2頻譜分辨率的奇數(shù)倍時(shí)，影響越大。從時(shí)域分析，自相關(guān)峰值附近的值相對較低，出現(xiàn)了一個(gè)凹陷，驗(yàn)證了本文1.2節(jié)中所描述的“相鄰1”的影響。在小頻偏影響下，相關(guān)結(jié)果峰值幾乎不受頻偏的影響，在頻偏值造成的影響較大時(shí)(接近1/2頻譜分辨率的奇數(shù)倍)，其相關(guān)結(jié)果峰值依然可以近似達(dá)到無頻偏時(shí)的70%左右。對比文獻(xiàn)[15]中的方法，本文同步方法可以有效減小大頻偏對同步造成的影響。

圖17 同步相關(guān)峰示意圖(頻域)

圖18 同步相關(guān)峰示意圖(時(shí)域)

圖19為本文同步方法與PN碼相關(guān)法的定時(shí)精度對比。圖19中，PN碼長度與本文使用的同步碼長度相同(1.28 s)。由圖19可以看出，本文所提出的同步方法在信噪比較低時(shí)具有明顯的優(yōu)勢。

圖19 定時(shí)誤差分析

由于使用了離散頻譜分析方法進(jìn)行精估計(jì)，本文同步方法的精度在低信噪比情況下呈現(xiàn)出明顯的周期性。設(shè)置信噪比為-17 dB，圖20和圖21為本文同步方法在iturHFMQ信道中的時(shí)間估計(jì)誤差和頻率估計(jì)誤差仿真圖。在低信噪比短波信道條件下，本文同步方法依然具有很高的估計(jì)精度，但估計(jì)誤差出現(xiàn)周期性變化的現(xiàn)象。這是由于在精同步過程中，依賴頻譜或自相關(guān)函數(shù)中最高值坐標(biāo)和次高值坐標(biāo)的判定，當(dāng)時(shí)間延遲或者頻率偏移接近其歸一化值-0.5、0、+0.5時(shí)會出現(xiàn)最高值坐標(biāo)或次高值坐標(biāo)判定錯(cuò)誤的情況。以歸一化頻率偏差接近±0.5為例，其峰值譜線的兩側(cè)譜線在實(shí)際中趨近于相等，由于受到噪聲的干擾，其值相對于真實(shí)值來說出現(xiàn)波動(dòng)，會造成次高值判定錯(cuò)誤的情況，因此，估計(jì)誤差會出現(xiàn)一個(gè)周期性的變化(在本文的仿真條件中，頻譜分辨率為3.125 Hz，時(shí)間分辨率為0.04 s，歸一化是指一個(gè)分辨率)。

圖20 時(shí)間估計(jì)誤差分析

圖21 頻率估計(jì)誤差分析

4 結(jié)束語

本文提出了一種基于Costas序列的時(shí)頻聯(lián)合同步方法，經(jīng)仿真驗(yàn)證，該方法與基于PN碼的同步方法相比，具有更好的抗噪聲和頻偏能力，可以有效提高短波通信中的捕獲率和同步精度，且相同長度下同步性能更好，更加適用于短波通信，可以有效提高惡劣信道條件下的同步性能，為提高短波通信的可靠性奠定基礎(chǔ)。