蘇紹璟，伍文君，黃芝平，劉純武

(國防科技大學 機電工程與自動化學院，湖南 長沙410073)

0 引言

m 序列在通信領(lǐng)域中有著廣泛應用。如在SDH 通信系統(tǒng)中，m 序列用來對數(shù)據(jù)流進行擾碼處理，以增強數(shù)據(jù)的隨機性，便于時鐘恢復以及增強通信的保密性;在擴頻通信中，m 序列也是應用最廣泛的擴頻碼序列。因此，有必要研究如何檢測與恢復m 序列。

國內(nèi)外一些學者對m 序列的檢測進行了探討。Adams，Peter 等［1-4］提出高階統(tǒng)計的方法對m 序列進行檢測，并提到根據(jù)高階統(tǒng)計函數(shù)峰值的位置可以測定m 序列的本原多項式。文獻［5］也研究了基于高階統(tǒng)計的m 序列檢測，提出了基于偏三階相關(guān)函數(shù)峰值特性的m 序列的檢測方法和識別標準，證明偏TCF 具有與三階相關(guān)函數(shù)所對應截取區(qū)域相同的峰值特性，根據(jù)這一特性可以對m 序列進行檢測及識別。

現(xiàn)有文獻的研究重點是m 序列的檢測，而對m序列的本原多項式的測定研究的并不深入。Adams等提到可以通過計算峰值位置的重復性來確定m序列的本原多項式。但該方法在本原多項式的階數(shù)較高，而截斷的序列長度有限的情況下，并不適用，且計算量很大。實際中常用的測定本原多項式的方法是Walsh 變換法［6］。但該方法的性能受制于本原多項式的抽頭數(shù)，且當截斷序列長度一定，而誤碼率增加時，其能準確測定的概率將顯著降低。為了克服Walsh 變換法的這一不足，本文在Adams 等的研究成果的基礎(chǔ)上，利用概率分析，深入研究了基于高階統(tǒng)計原理的本原多項式測定算法。新算法克服Walsh 變換法性能受制于本原多項式抽頭數(shù)的影響，但計算復雜度相對更高。

1 m 序列及高階統(tǒng)計原理

1.1 m 序列的產(chǎn)生與性質(zhì)

m 序列是最大周期線性反饋移位寄存器序列。圖1為生成線性反饋移位寄存器序列的原理圖。從圖中可以看出，線性反饋移位寄存器序列由寄存器的抽頭系數(shù)及寄存器的初始狀態(tài)決定。設(shè)序列的生成多項式為

f(x)的多項式系數(shù)對應于反饋移位寄存器的抽頭系數(shù)。因此f(x)就確定了線性反饋移位寄存器的結(jié)構(gòu)。當f(x)為二元域上本原多項式時，其生成的序列具有最大周期，也就是m 序列。

圖1 線性反饋移位寄存器原理圖Fig.1 Schematic of liner feedback shift register

m 序列具有一些優(yōu)良特性:

1)平衡性:m 序列中一個周期內(nèi)“1”的數(shù)目比“0”的數(shù)目多1 個。

2)平移相加性:m 序列與其移位后的序列逐位模2 相加后，所得序列仍然是一m 序列，只是相移不同而已。

3)二值相關(guān)性:m 序列具有優(yōu)良的二值自相關(guān)特性，碼位數(shù)越長越接近于隨機噪聲的自相關(guān)特性。但是，2 個相同長度的m 序列的互相關(guān)特性并不是很好。

1.2 高階統(tǒng)計原理

以3 階統(tǒng)計為例進行介紹高階統(tǒng)計原理。m 序列的3 階相關(guān)函數(shù)(triple correlation function，TCF)定義為

式中:R(t1，t2)為3 階相關(guān)函數(shù);z(t)為m 序列;t1，t2=1，2，…，L -1，L 為序列的周期。實際分析中，z(t)由±1 組成，而不是0，1，對應關(guān)系如下:(0，1)→(1，-1).故式(2)中的乘法，其實就是模2加。3 階相關(guān)值在實際中由下式近似計算

由m 序列的平移相加性，對任意的相移n1，z(n)z(n+n1)仍是m 序列。假設(shè)得到的m 序列相對原序列的相移為n2，則?n，z(n)z(n +n1)z(n +n2)=1，此時，R(n1，n2)=1.而對于不滿足上述的條件的n1，n2，z(n)z(n +n1)z(n +n2)也仍然是m序列。由m 序列的平衡性可知，在m 序列的一個周期內(nèi)，”1”的數(shù)目比”0”的數(shù)目多1.因此，此時R(n1，n2)=-1/L.所以，m 序列的TCF 為

因此，m 序列的平移相加性導致m 序列的TCF有峰值出現(xiàn)。根據(jù)這一特點，就可以通過計算3 階相關(guān)值來檢測m 序列。更高階的相關(guān)統(tǒng)計同樣存在這樣的特征。

實際上，大多數(shù)情況下我們并不知道m(xù) 序列的周期長度，而且，通常截取的m 序列長度小于m 序列的周期長度，那么，截取的這一部分m 序列的TCF 是否能體現(xiàn)m 序列的TCF 的峰值特性呢?文獻［5］說明這是可能的，并提出了偏3 階相關(guān)函數(shù)(偏TCF)的概念，定義如下

式中，N 為截取的序列長度。按照式(5)計算得到的1～N 區(qū)間的TCF 峰值位置與按照整個序列周期長度計算的TCF 相對應的局部區(qū)域的峰值位置是一致的。

2 基于高階統(tǒng)計的本原多項式測定

2.1 基本思路

令z(n)(0≤n＜N)為m 序列，f(x)為其本原多項式，其階數(shù)為l，截斷序列長度為N.計算該序列的偏d 階相關(guān)函數(shù)

若R 在某位置(s1，s2，…，sd-1)出現(xiàn)峰值，則有

觀察上式可知，多項式g(x)=1 + xs1+ xs2+…+xsd-1的反轉(zhuǎn)多項式g'(x)符合m 序列z(n).由m 序列的基本理論可知，f(x)整除g'(x)，也即f(x)是g'(x)的一個因式。這樣，若能找到R 的幾個峰值，通過求這幾個峰值位置對應多項式的最大公約式，就能測定該序列的本原多項式。

在實際應用中，得到的截斷m 序列的長度通常要遠小于該序列的周期。此時，序列高階相關(guān)函數(shù)并不一定存在峰值。例如，對于本原多項式f(x)=x17+x15+x14+x13+x11+x10+x9+x8+x6+x5+x4+x2+1，截斷序列長度N=400，其3 階相關(guān)函數(shù)并不存在峰值，而只有更高階的相關(guān)函數(shù)才存在峰值。因此，我們有必要對m 序列的高階相關(guān)函數(shù)的峰值分布進行研究。

文獻［7］提到，由于m 序列具有很好的偽隨機特性，假定其高階相關(guān)函數(shù)峰值也平均分布，則

式中，m(d)為d 階相關(guān)函數(shù)峰值數(shù)的估計值。文獻［7］還以本原多項式f1(x)=x17+x3+1 及f2(x)=x17+x15+x14+x13+x11+x10+x9+x8+x6+x5+x4+x2+1 做了仿真實驗，實驗結(jié)果見表1［7］.實驗表明這個公式是合理的。該公式也是本文算法的基本假設(shè)之一。

表1 式(8)估計值與實際值對比Tab.1 The comparison between approximate value with Function (8)and actual value

有了上述假設(shè)，我們就可以設(shè)計算法了。設(shè)有含錯m 序列z(n)(0≤n＜N)，N 為序列長度，誤碼率為p.?(n1，n2，…，nd-1)，(1＜n1＜n2＜…＜nd-1＜N1)，計算其d 階相關(guān)函數(shù)

式中，N1為向量(n1，n2，…，nd-1)對應多項式的最高階數(shù)，N2為統(tǒng)計次數(shù)，顯然N1+N2≤N.搜索R(n1，n2，…，nd-1)的峰值，通過求解這些峰值位置所對應的多項式的最大公約式，就可以導出該序列的本原多項式了。

2.2 參數(shù)選取

式(9)中有3 個參數(shù)需要確定:相關(guān)函數(shù)階數(shù)d，多項式最高階數(shù)N1以及統(tǒng)計次數(shù)N2.下面分別討論它們的選取。首先對于特定相關(guān)函數(shù)階數(shù)d，N1的選取至少應滿足m(d)≥2，當m(d)=2 時，得到的最大公約式仍然可能為本原多項式的倍式。實踐中發(fā)現(xiàn)，當m(d)≥4 時，得到的最大公約式通常為所求的本原多項式。另外考慮到誤碼率較大時，峰值位置可能出錯或不能得到，本算法暫取m(d)≥10.故:

N2的選取應滿足在高誤碼率情況下，也能將部分d 階相關(guān)函數(shù)R 的峰值區(qū)分出來。對于含錯m序列，由puling-up 定理［8］可知，式z(n)z(n+n1)z(n+n2)…z(n+nd-1)成立概率s 為

式中，ε=1/2 -p.這樣，當向量(n1，n2，…，nd-1)處于R 的峰值位置時，Pr(z(n)z(n +n1)z(n +n2)…z(n+nd-1)=1)=s.而對于非峰值位置的向量(n'1，n'2，…，n'd-1)，由m 序列的偽隨機特性，有Pr(z(n)z(n+n'1)z(n +n'2)…z(n +n'd-1)=1)=0.5.假設(shè)?n，z(n)z(n +n1)z(n +n2)…z(n +nd-1)相互獨立，則R(n1，n2，…，nd-1)可以看成N2次伯努利試驗。由De Moivre-Laplace 中心極限定理，當N2較大時，R(n1，n2，…，nd-1)漸進于正態(tài)分布。

對處于峰值位置的向量(n1，n2，…，nd-1)，有E(z(n)z(n+n1)z(n +n2)…z(n +nd-1))=2s -1，D(z(n)z(n+n1)z(n+n2)…z(n +nd-1))=4s(1 -s)，則

而對于其他向量(n'1，n'2，…，n'd-1)，有E(z(n)z(n+n'1)z(n+n'2)…z(n+n'd-1))=0，D(z(n)z(n +n'1)z(n+n'2)…z(n+n'd-1))=1.則

考察式(12)及式(13)可以發(fā)現(xiàn)，當誤碼率p 越大時，s 趨近于0.5，此時峰值位置與非峰值位置的期望差越小，也就是兩者的概率空間重疊部分越多，峰值就越難區(qū)分出來。另一方面，N2越大，兩者的方差越小，此時概率空間重疊越小，峰值就更容易區(qū)分。圖2，3 反映了這一過程。圖2，3 為誤碼率p =0.25，N2分別為300 和800 時，峰值位置的了概率分布與非峰值的R(n1，n2，…，nd-1)概率分布情況。

圖2 p=0.25，N2 =300 時，R(n1，n2，…，nd-1)概率分布Fig.2 The probability distribution of R(n1，n2，…，nd-1)when p=0.25 and N2 =300

由以上分析可知，對于特定的p 與d，N2的取值決定了是否能區(qū)分峰值。設(shè)定閾值h，當R(n1，n2，…，nd-1)＞h 時，我們認為這個值為d 階相關(guān)函數(shù)的峰值。h 應使非峰值被誤認為峰值的平均數(shù)目小于1，這樣我們就最多可能得到一個錯誤的峰值，有

另外，h 還應滿足使大多數(shù)的峰值能區(qū)分出來。由于計算Pr(R(n1，n2，…，nd-1)＞h)的解析式較為復雜，為了簡化計算，我們?nèi) 小于峰值位置R(n1，n2，…，nd-1)的期望值，即

圖3 p=0.25，N2 =800 時，R(n1，n2，…，nd-1)概率分布Fig.3 The probability disribution of R(n1，n2，…，nd-1)when p=0.25 and N2 =800

綜合式(14)和式(15)，我們就能確定參數(shù)N2.

現(xiàn)在，我們可以說明之前為什么選擇m(d)≥10.由式(15)可知，峰值被找到的概率大于0.5，這樣，很容易算得:從10 個峰值中找出3 個峰值的概率就大于0.95.因此，選擇m(d)≥10，使得我們至少能找到3 個峰值，以推導出本原多項式。

相關(guān)函數(shù)的階數(shù)d 是需要探討的另一個重要的參數(shù)。式(10)表明，峰值的數(shù)量m 隨著d 的增長而增長。當截斷m 序列的長度遠低于序列周期時，可以通過增加d 以找到算法成功所需的峰值數(shù)。然而，另一方面，式(11)也表明，d 的增加將導致s 的降低。而由式(15)，s 的降低將導致N2的急劇增加，并可能使得N2＞N，算法失敗。由以上分析可知，d，N1和N2這3 個參數(shù)是相互聯(lián)系，相互制約的，需要綜合權(quán)衡才能確定它們的最佳取值。在本算法中，初始化d=3，當算得N1＞N 時，再增加d 重新再算，直到N1＜N.

2.3 算法步驟

有了以上的這些探討，下面給出新算法的完整步驟:

1)通過觀察m 序列的游程，估計本原多項式的階數(shù)l.初始化d=3，由式(10)計算N1.如果N1＞N，則更新d =d +1，并重新計算N1.重復這一步驟直到N1＜N.

2)估計序列的誤碼率為p，由式(14)、式(15)并查正態(tài)分布表，確定閾值h 和N2.如果N1+N2≥N，則轉(zhuǎn)步驟5.

3)?(n1，n2，…，nd-1)，0≤n1≤n2≤…≤nd-1＜N1，計算其d 階相關(guān)函數(shù)R(n1，n2，…，nd-1)，當R(n1，n2，…，nd-1)＞h，儲存該向量(n1，n2，…，nd-1).如果找到了4 個峰值，則轉(zhuǎn)到步驟4.a，如果只找到了3 個峰值，則轉(zhuǎn)步驟4.b，否則轉(zhuǎn)步驟5.

4)求解儲存向量(n1，n2，…，nd-1)所對應的反轉(zhuǎn)多項式(1 +xn1+xn2+…+xnd-1)'的最大公約式。

a.如果算得的最大公約式不為1，則認為該最大公約式就是序列的本原多項式。如果其中一個多項式與其他多項式互素，則認為該多項式是錯誤的，重新計算其他多項式的最大公約式，如果這個最大公約式不為1，則同樣認為該最大公約式就是序列的最大公約式。否則，轉(zhuǎn)步驟5.

b.如果算得的最大公約式不為1，則該最大公約式即為序列的本原多項式，否則轉(zhuǎn)步驟5.

5)未能測得序列的本原多項式，算法失敗。

3 仿真計算及性能分析

3.1 仿真計算

使用Matlab 軟件，在Pentium4 2.8 G 機器上進行了大量的仿真計算。例如，以f(x)=x17+ x15+x14+x13+x11+x10+x9+x8+x6+x5+x4+x2+1 為本原多項式，誤碼率0.25，生成截斷序列長度10 000的含錯m 序列。首先通過觀察含錯序列的游程，估計其本原多項式階數(shù)l 為20.令d =3，由式(10)算得N1=4 096.由式(14)、式(15)及查正態(tài)分布表，算得h =0.125 和N2=1 600.窮舉所有向量(n1，n2)，0≤n1≤n2＜N1，計算其3 階相關(guān)函數(shù)R(n1，n2)，儲存R(n1，n2)≥h =0.125 的向量。這樣，我們得到(130，1258)，(210，813)，(278，1543)和(216，483)，其對應的反轉(zhuǎn)多項式分別是x1258+x1127+1，x813+x602+1，x1543+x1264+1 及x483 +x266 +1，它們的最大公約式便是其本原多項式f(x).上述過程耗時約35 s.由于每次仿真找到的峰值位置都不相同，計算時間也就不同。在上面這個例子中，計算時間大約在15～60 s 之間。

3.2 計算復雜度分析

計算新算法的計算復雜度。在最壞的情況下，新算法需要窮舉所有可能的向量(n1，n2，…，nd-1)，0≤n1≤n2≤…≤nd-1＜N1，總共有個向量。由式(10)，有Nd-11/(d -1)!=10 ×2l，所以總共有10 ×2l個向量。對于每個向量都需要做N2次運算，因此新算法的計算復雜度為o(2lN2).相比Walsh 變換法的計算復雜度o(2ll))［6］，新算法的計算復雜度要更高。在實際應用中，通常并不需要窮舉所有可能的向量，因此計算時間通常要遠小于最壞情況下的計算時間。

3.3 算法成功概率分析

Walsh 變換法成功測定本原多項式的概率與本原多項式的抽頭數(shù)相關(guān)，本原多項式的抽頭數(shù)越多，其成功測定的概率越低。在上一小節(jié)的例子中，Walsh 變換法就未能成功找到本原多項式，但如果將f(x)=x17+x3+1 替代上例使用的本原多項式，則Walsh 變換法也能成功求出本原多項式。由此可見，新算法的性能與本原多項式的抽頭數(shù)無關(guān)。

為進一步考察新算法的性能，分別以C1(x)=x17+x3+1 和C2(x)=x17+ x15+ x14+ x13+ x11+x10+x9+x8+x6+x5+x4+x2+1 為生成多項式，在不同的誤碼率及不同的截斷序列長度情況下做了大量的仿真試驗。圖4和圖5分別給出了生成多項式為C1(x)和C2(x)時，在不同的誤碼率下，Walsh 變換法和新算法能成功測定生成多項式的概率。圖中截斷序列長度分別為20 000 和80 000，算法成功概率值為隨機試驗500 次得到的估計值。從圖中可以看出，當生成多項式的抽頭數(shù)為2 時，Walsh 變換法和新算法的容錯性能相當，但隨著生成多項式抽頭數(shù)的增加，Walsh 變換法的容錯能力急劇降低，而新算法的容錯性能不隨生成多項式抽頭數(shù)的增加而降低。此外，截斷序列長度越長，相同誤碼率下，算法成功的概率越高。

圖4 新算法與Walsh 算法容錯能力對比(生成多項式抽頭數(shù)為2)Fig.4 The error tolerance ability comparison between the new algorithm and the Walsh-Hadamard transformation when the number of tap is 2

由以上對新算法的性能分析可知:相比實際中常用的Walsh 變換法，新算法克服了其性能受制本原多項式抽頭數(shù)的不足，并且在本原多項式階數(shù)為17，截斷序列長度為50 000 及誤碼率達到0.36 的情況下，依然可以算得正確的本原多項式，體現(xiàn)出了很好的容錯能力。新算法不足之處在于，算法的計算復雜度相對較高。

圖5 新算法與Walsh 算法容錯能力對比(生成多項式抽頭數(shù)為12)Fig.5 The error tolerance ability comparison between the new algorithm and the Walsh-Hadamard transformation when the number of tap is 12

4 小結(jié)

實際中常用Walsh 變換法［6］來測定m 序列的本原多項式，但是該方法的性能受制于本原多項式的抽頭數(shù)。當本原多項式的階數(shù)較高時，其抽頭數(shù)也往往較多，這時Walsh 變換法常常無法得到本原多項式。本文基于這點考慮，對基于高階統(tǒng)計的m序列本原多項式測定做了深入的研究。本文給出了基于高階統(tǒng)計的m 序列本原多項式測定的基本思路，詳細探討了各個參數(shù)的確定方法，并給出算法的完整步驟。該算法與常用的Walsh 變換法相比，性能不再受制于本原多項式的抽頭數(shù)，且具有很好的容錯能力。新算法不足之處在于，計算復雜度相對較高。

References)

［1］Adams E R，Gouda M，Hill P C J.Detection ＆characterisation of DS/SS signals using higher-order correlation［C］∥Proc IEEE ISSSTA’96.Mainz，1996:7 -31.

［2］Batty K E，Adams E R.Detection and blind identification of msequence codes using higher order statistics［C］∥Proceedings of IEEE on Signal Processing Workshop，Caesarea.1999:16 -20.

［3］Adams E R，Gouda M，Hill P C J.Statistical techniques for blind detection ＆ discrimination of m-sequence codes in DS/SS systems［C］∥Proc IEEE 5th ISSSTA symposium’98.Sun City，1998:853 -857.

［4］Adams E R，Hill P C J.Detection of direct sequence spread spectrum signals using higher-order statistical processing［C］∥Proc Int.Conference on Acoustics Speech and Signal Processing，ICASSP’97.Munich，1997:3849 -3852.

［5］俎云霄.基于高階統(tǒng)計處理技術(shù)的m-序列檢測及識別［J］.電子與信息學報，2007，29(7):1576 -1579.ZU Yun-xiao.The detection and recognition of m-sequence using higher-order statistical processing［J］.Journal of Electronics and Information Technology，2007，29(7):1576 -1579.(in Chinese)

［6］游凌，朱中梁.Walsh 函數(shù)在解二元域方程組上的應用［J］.信號處理，2000，16:27 -30.YOU Ling，ZHU Zhong-liang.The application of walsh function in resolving of F(2)equations［J］.Signal Processing(chinese)，2000，16:27 -30.(in Chinese)

［7］Canteaut A，Trabbia M.Improved fast correlation attacks using parity-check equations of weight 4 and 5［G］∥Bart Preneel，Advances in Cryptology-EUROCRYPT 2000 LNCS.Germany:Springer，2000，1807:573 -588.

［8］Matsui M.Linear cryptanalysis method for DES cipher［G］∥Tor Helleseth，Advances in Cryptology-EUROCRYPT 1993.LNCS.Germany:Springer，1994，765:386 -397.