鄭 鵬 ， 張 鑫 ， 劉 鋒 ，， 陶 然

（1.海軍航空工程學(xué)院 電子信息工程系，山東 煙臺(tái) 264001；2.北京理工大學(xué) 信息科學(xué)技術(shù)學(xué)院，北京 100081）

譜相關(guān)理論由Gardner W A等人深入研究并發(fā)展[1-2]，其優(yōu)點(diǎn)主要表現(xiàn)在分辨率高、抗干擾能力強(qiáng)、雙頻率域信息豐富且易于提取等。在通信、雷達(dá)、遙測(cè)、水聲、電子對(duì)抗等領(lǐng)域中，由于某些人為的或無意的調(diào)制，而使大量信號(hào)呈現(xiàn)出循環(huán)平穩(wěn)性，因此譜相關(guān)理論在信號(hào)檢測(cè)、參數(shù)估計(jì)、系統(tǒng)辨識(shí)、調(diào)制識(shí)別、故障診斷等現(xiàn)代信號(hào)處理諸多領(lǐng)域中得到了越來越深入的研究。循環(huán)譜估計(jì)是譜相關(guān)理論的主要研究對(duì)象，對(duì)循環(huán)譜的估計(jì)是研究循環(huán)平穩(wěn)過程的重要前提，也是其應(yīng)用基礎(chǔ)，循環(huán)譜估計(jì)方法中工程應(yīng)用最廣的主要是基于FFT的時(shí)域平滑法中的FAM和SSCA、頻域平滑法中的FSM[3-4]。其中SSCA算法因其計(jì)算量最小，計(jì)算效率最高，因此在工程中得到了廣泛應(yīng)用，其采用方法[5]都是文獻(xiàn)[3-4]中提到的經(jīng)典算法，但是這種算法存在循環(huán)譜泄露的問題，而且會(huì)產(chǎn)生錯(cuò)誤的循環(huán)特征[6-7]，其循環(huán)譜估計(jì)效果并不理想，本文在已有文獻(xiàn)基礎(chǔ)上提出了一種改進(jìn)的SSCA算法，該算法能夠改善譜估計(jì)，便于實(shí)現(xiàn)檢測(cè)、參數(shù)估計(jì)、調(diào)制識(shí)別等后續(xù)信號(hào)處理。

1 循環(huán)譜基本理論

設(shè)x（t）為循環(huán)平穩(wěn)信號(hào)，其均值和自相關(guān)函數(shù)呈現(xiàn)以T0為周期的周期性，即

那么可以將（2）式展開成Fourier級(jí)數(shù)的形式，即

其中α為循環(huán)頻率，循環(huán)自相關(guān)函數(shù)Rαx（τ）是Fourier級(jí)數(shù)的系數(shù)，

Rαx（τ）的 Fourier變換 Sαx（ f）稱為循環(huán)譜密度，簡(jiǎn)稱循環(huán)譜：

由上述分析可以看出，Rαx（τ）可以看成 x（t）的兩個(gè)復(fù)數(shù)頻移分量 u（t）和 v（t）的互相關(guān)

由式（5）可以看出，Sαx（f）=Sαuv（f）。由互譜分析可以得到

其中，式（7）是循環(huán)譜密度估計(jì)式，式（8）是循環(huán)周期圖，式（9）是接收數(shù)據(jù)的短時(shí)傅里葉變換。Δt是接收數(shù)據(jù)的長(zhǎng)度，T是短時(shí)傅里葉變換的窗長(zhǎng)，上式*表示復(fù)共軛。

循環(huán)周期圖不能直接作為譜相關(guān)函數(shù)的估計(jì)值，它的方差很大，而且隨著采集時(shí)間的增大并不趨為零，同時(shí)，它也不是無偏估計(jì)。為了減小估值誤差，減小譜泄露，得到循環(huán)譜可靠估計(jì)，對(duì)于有限時(shí)間長(zhǎng)度的接收數(shù)據(jù)，必須對(duì)循環(huán)周期圖進(jìn)行平滑處理。平滑分為時(shí)域平滑和頻域平滑兩種，其中時(shí)域平滑因具有較高的計(jì)算效率和較小的計(jì)算量而在工程中被廣泛采用。在這里采用時(shí)域平滑的循環(huán)周期圖來估計(jì)信號(hào)的循環(huán)譜密度。循環(huán)譜密度估計(jì)的時(shí)域平滑SSCA算法可以表示為

其中 XT（r，f）是復(fù)解調(diào)，又被稱為循環(huán)周期圖，a（n）為數(shù)據(jù)衰減窗，gc（n）為平滑窗，q=-P/2，…，P/2-1，L 為抽取因子，滿足L＜N′，當(dāng)L=N′/4時(shí)效果最好。N為總的采樣數(shù)，LP=N。SSCA算法的實(shí)現(xiàn)原理圖如圖1所示。

圖1 SSCA算法實(shí)現(xiàn)Fig.1 Realization of SSCA

SSCA算法的頻率分辨率為Δf=fs/N′，循環(huán)頻率分辨率為Δα=1/Δt=fs/N。要得到循環(huán)譜可靠估計(jì)，使循環(huán)譜密度函數(shù)的估計(jì)與時(shí)間幾乎不相關(guān)，需滿足可靠性條件：

由上式可見要實(shí)現(xiàn)可靠估計(jì)，循環(huán)頻率分辨率Δα必須比頻率分辨率Δf高得多。

（1）從產(chǎn)業(yè)轉(zhuǎn)移角度分析，建筑業(yè)、工業(yè)與交通運(yùn)輸、倉儲(chǔ)和郵政業(yè)在區(qū)域間的轉(zhuǎn)移趨勢(shì)明顯。東北、華北和華東地區(qū)表現(xiàn)為相對(duì)穩(wěn)定的轉(zhuǎn)出，而中南、西南和西北地區(qū)表現(xiàn)為穩(wěn)定的轉(zhuǎn)入，產(chǎn)業(yè)結(jié)構(gòu)的空間調(diào)整同時(shí)也對(duì)能源強(qiáng)度的空間分布產(chǎn)生影響。

2 算法改進(jìn)及數(shù)字實(shí)現(xiàn)

算法改進(jìn)之處在于線性內(nèi)插和數(shù)據(jù)加窗：

1）復(fù)解調(diào)實(shí)現(xiàn)的內(nèi)插方法 由于信道化處理的采樣處理，復(fù)解調(diào)的采樣速率為fs/L。在SSCA算法的實(shí)現(xiàn)過程中要求復(fù)解調(diào)與未進(jìn)行濾波處理的原始信號(hào)共軛進(jìn)行相乘，而原始信號(hào)的采樣速率fs。為了完成乘法運(yùn)算需要對(duì)復(fù)解調(diào)進(jìn)行內(nèi)插處理以使其采樣速率與原始信號(hào)共軛相匹配。傳統(tǒng)方法采用的是“保持”的處理方法，即將每一個(gè)復(fù)解調(diào)復(fù)制L次，但是這種方法會(huì)產(chǎn)生錯(cuò)誤的循環(huán)特征。因此，本文采取了線性內(nèi)插方法來代替“保持”的方法。

2）數(shù)據(jù)加窗 算法的第2個(gè)改進(jìn)之處在于窗函數(shù)的選擇。對(duì)窗函數(shù)的選擇應(yīng)滿足以下要求：a（n）應(yīng)具有大的衰減，gc（n）應(yīng)具有小的帶寬。在Gardner的經(jīng)典文獻(xiàn)中提到兩者分別采用漢明窗和矩形窗，這種方法被廣泛采用。實(shí)際上這種選擇并不十分合適，對(duì)循環(huán)譜泄露和虛假循環(huán)特征的產(chǎn)生的抑制效果并不是很好。基于窗函數(shù)的選擇要求，本文采取了適應(yīng)性更強(qiáng)的凱澤窗，通過調(diào)整窗函數(shù)的參數(shù)來達(dá)到最好別的效果。

改進(jìn)算法的實(shí)現(xiàn)原理圖如圖2所示，數(shù)字實(shí)現(xiàn)過程如下：

圖2 改進(jìn)的SSCA算法實(shí)現(xiàn)Fig.2 Realization of innovatory SSCA

1）n時(shí)刻N(yùn)′的點(diǎn)數(shù)據(jù)加窗，然后進(jìn)行N′點(diǎn)FFT；

2）折疊運(yùn)算。求出一個(gè)周期的頻譜，得到信號(hào)的復(fù)解調(diào)表達(dá)式，即 fk=k（fs/N′），k=-N′/2，…，N′/2-1；

3）求出 n 時(shí)刻復(fù)解調(diào)的相關(guān)，即 XT（n， fk）x*n，得到 N′個(gè)數(shù)據(jù)；

4）采用線性內(nèi)插對(duì)采樣后的復(fù)解調(diào)進(jìn)行賦值，重復(fù)1，2，3直到 n=N，得到 N×N′矩陣；

5）對(duì)N×N′矩陣中的每一列求N點(diǎn)FFT，得到頻率和循環(huán)頻率分別為

的N個(gè)數(shù)據(jù)。

通過上面步驟，可以看出SSCA算法主要有3個(gè)階段：

第1階段：利用N′點(diǎn)FFT計(jì)算數(shù)據(jù)的復(fù)解調(diào)。

第2階段：計(jì)算復(fù)解調(diào)與信號(hào)的共軛相關(guān)。

第3階段：計(jì)算N點(diǎn)FFT，得到循環(huán)譜估計(jì)值。

需要指出的是SSCA算法由于利用原始信號(hào)代替復(fù)解調(diào)信號(hào)的復(fù)共軛，采用了1個(gè)復(fù)解調(diào)濾波器，它的輸出信噪比相對(duì)于FAM算法有所減小，但換來了較小的計(jì)算量。

3 仿真實(shí)驗(yàn)及結(jié)果分析

若噪聲為加性平穩(wěn)噪聲，那么MPSK信號(hào)的模型表示為：

其中：xl（t）為 x（t）的復(fù)包絡(luò)；f0為信號(hào)載頻；θk為調(diào)制相位信息，θk∈{2π（k-1）/M}，k=1，…，M；q（t）為成型脈沖；Tc為碼元寬度，其倒數(shù)fc為碼元速率。

MPSK信號(hào)循環(huán)譜的一般形式[8-9]如下，當(dāng)M=2時(shí)，即BPSK信號(hào)的循環(huán)譜為：其中

可以看出，對(duì)所有MPSK信號(hào)，在α=n/Tc處均取得較大非零值，而在其它地方取值為零或接近于零。BPSK信號(hào)的循環(huán)譜除具有上述特征外，在f=0，α=2f0n/Tc處也取得較大非零值。本文采用BPSK、QPSK信號(hào)作為仿真對(duì)象，噪聲為高斯白噪聲。仿真中參數(shù)設(shè)置如下：采樣頻率fs=8 192 Hz，載頻f0=2 048 Hz，碼片速率為fc=512 Hz，信號(hào)初始時(shí)間和初始相位均為0。

圖3是算法改進(jìn)前后BPSK信號(hào)循環(huán)譜密度。圖4是取f=0其中切面循環(huán)譜包絡(luò)。圖5是算法改進(jìn)前后QPSK信號(hào)的循環(huán)譜密度。圖6是取f=fc其中切面的循環(huán)譜包絡(luò)。由圖可以看出算法改進(jìn)前信號(hào)的循環(huán)譜在錯(cuò)誤的循環(huán)頻率處出現(xiàn)了峰值，循環(huán)譜泄露現(xiàn)象仍舊比較明顯，而改進(jìn)算法則對(duì)上述兩種現(xiàn)象抑制效果明顯，這對(duì)于信號(hào)循環(huán)特征的提取時(shí)十分有利的，這說明改進(jìn)算法對(duì)循環(huán)譜估計(jì)效果要好于原算法。

圖3 BPSK信號(hào)的循環(huán)譜密度Fig.3 Cyclic spectrum of BPSK signal

圖4 BPSK信號(hào)f=0處循環(huán)譜包絡(luò)Fig.4 The amplitude of BPSK signal’s cyclic spectrum

4 結(jié) 論

隨著對(duì)信號(hào)處理研究的深入，譜相關(guān)理論的應(yīng)用會(huì)越來越廣。本文針對(duì)SSCA算法存在譜泄露以及產(chǎn)生錯(cuò)誤的循環(huán)頻率信息等問題，提出了一種改進(jìn)算法，這種算法實(shí)現(xiàn)簡(jiǎn)單，能夠有效地改善循環(huán)譜估計(jì)，這對(duì)于信號(hào)檢測(cè)、參數(shù)估計(jì)、調(diào)制識(shí)別等后續(xù)信號(hào)處理是十分有益的。

[1]Gardner W A.Statistical spectral analysis：a non-probabilistic theory[M].NJ：Prentice-Hall， Englewood Cliffs，1988.

圖5 Cyclic spectrum of QPSK signal Fig.5 QPSK信號(hào)的循環(huán)譜密度

圖6 QPSK信號(hào)f=0處循環(huán)譜包絡(luò)Fig.6 The amplitude of QPSK signal’s cyclic spectrum

[2]Gardner W A.Exploitation of spectral redundancy in cyclostationary signals[J].IEEE SP MAG，1991，（04）：14-36.

[3]Roberts R S，Brown W A，Loomis H H， et al.Computionally efficient algorithms for cyclic spectral analysis[J].IEEE Signal Processing， 1991，8（2）：38-49.

[4]Brown W A，Loomis H，Jr.Digital implementations of spectral correlation analyzers[J].IEEE Trans.Signal Processing，1993，4（2）： 703-720.

[5]高玉龍，張中兆.SSCA算法改進(jìn)及實(shí)現(xiàn)[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2008，40（9）：1374-1377.

GAO Yu-long，ZHANG Zhong-zhao.Modification and realization ofSSCA algorithm[J].JournalofHarbin Institute of Technology，2008，40（9）：374-1377.

[6]Gillman A M.Non co-operative detection of LPI/LPD signals via cyclic spectral analysis[D].ADA 361720/XAB，1999.

[7]Eric April.On the implementation of the strip spectral correlation algorithm forcyclic spectrum estimation[M].Ottawa：Defence Research Establishment，1995.