李 偉, 楊 斌
(西南交通大學(xué)信息科學(xué)與技術(shù)學(xué)院,四川成都610031)
隨機(jī)噪聲在通信等過程中往往是有害并需設(shè)法消除的無用信號,然而在一些情況下,也需要人為產(chǎn)生噪聲,如自適應(yīng)控制系統(tǒng)、系統(tǒng)辨識等領(lǐng)域。正態(tài)分布的高斯噪聲是最常用的一種噪聲。隨機(jī)噪聲由隨機(jī)數(shù)表達(dá),因此隨機(jī)噪聲的產(chǎn)生就歸結(jié)于隨機(jī)數(shù)的產(chǎn)生問題。由隨機(jī)信號理論可知,在(0,1)上服從均勻分布的隨機(jī)數(shù)經(jīng)過一定的變換,可產(chǎn)生服從N(0,1)正態(tài)分布的高斯隨機(jī)數(shù)[1]。
隨機(jī)數(shù)的產(chǎn)生大致可分為查表法、物理法(如真空管射槍噪聲)、數(shù)學(xué)遞推法[2]。計(jì)算機(jī)產(chǎn)生隨機(jī)數(shù)用到的是數(shù)學(xué)遞推法,字長和精度有限的計(jì)算機(jī)或DSP產(chǎn)生真正的隨機(jī)數(shù)是不可能的,因而現(xiàn)有的隨機(jī)數(shù)產(chǎn)生算法產(chǎn)生的隨機(jī)數(shù)均是偽隨機(jī)數(shù),即所得隨機(jī)數(shù)都有一定的周期性。實(shí)際應(yīng)用中也不必采用真正隨機(jī)的數(shù),達(dá)到一定程度隨機(jī)性的偽隨機(jī)數(shù)就可以滿足一般要求了。
文中介紹了3種產(chǎn)生(0,1)均勻分布隨機(jī)數(shù)的算法,隨后介紹了2種用(0,1)均勻分布隨機(jī)數(shù)生成正態(tài)分布隨機(jī)數(shù)的快速生成方法,最后在ADSP-TS201S上實(shí)現(xiàn)了實(shí)時高斯噪聲生成算法,能以10kHz的頻率連續(xù)產(chǎn)生兩路2048對不重復(fù)的正交高斯噪聲數(shù)據(jù)。
均勻分布隨機(jī)數(shù)產(chǎn)生的算法很多,如乘同余法、線性同余法、移位寄存器法、迭代取中法、加同余法、二次剩余法、進(jìn)位加/借位減/進(jìn)位乘(AWC/SWB/MWC)發(fā)生器法[3]等。根據(jù)實(shí)時性要求,主要介紹3種相對產(chǎn)生隨機(jī)數(shù)速度較快的算法,3種算法產(chǎn)生的隨機(jī)數(shù)均在(0,2L)區(qū)間上服從均勻分布,結(jié)果除以常數(shù)2L,即可得到(0,1)均勻分布的隨機(jī)數(shù)。
乘同余法用如下的遞推同余式產(chǎn)生隨機(jī)序列:
其中,xi為當(dāng)前隨機(jī)數(shù),xi-1為前一時刻隨機(jī)數(shù);M為模值,取2的冪,即 M=2L,L>2的整數(shù);A為乘因子,一般取小于M 且模8等于3或5的奇數(shù);初值 x0取0~2L之間任意奇數(shù)。可以證明,該算法產(chǎn)生的隨機(jī)數(shù)序列的周期為2L-2[2]。
線性同余法又稱為混合同余法,遞推公式如下:
這是一種類似乘同余的算法,只多了個增量c。根據(jù)隨機(jī)數(shù)理論,A的后3位一般為“偶數(shù)-2-1”形式,如421、821等;初值 x0取0~2L之間任意奇數(shù);Coveyou給出增量c滿足的一個等式可以在 A取值不大的條件下,有效地減小序列的相關(guān)系數(shù),而且不占用過多的機(jī)器時間。用該方法產(chǎn)生的隨機(jī)序列的周期為2L[1]。
這是一種基于m序列理論的算法,可以利用線性反饋移位寄存器(LFSR)產(chǎn)生隨機(jī)序列。根據(jù)m序列本原特征多項(xiàng)式,采取合適的回饋,M級LFSR產(chǎn)生的隨機(jī)數(shù)的最大周期為2M-1。M階m序列本原特征多項(xiàng)式可以表示為:
⊕表示模2或者異或運(yùn)算,系數(shù)ci取0或1,轉(zhuǎn)化為最大周期LFSR就是根據(jù)系數(shù)在LFSR相應(yīng)位抽頭。例如(16,12,3,1,0)表示16級LFSR的本原特征多項(xiàng)式為:
則應(yīng)在LFSR16、12、3、1位抽頭,如圖1所示。每次移位運(yùn)算后,把整個的16位數(shù)作為隨機(jī)數(shù)輸出。同樣,也可以根據(jù)(32,7,5,3,2,1,0)表示的32級LFSR的本原特征多項(xiàng)式得到32位隨機(jī)數(shù)。
正態(tài)分布隨機(jī)數(shù)可由(0,1)均勻分布隨機(jī)數(shù)產(chǎn)生,下面介紹兩種算法,一種基于中心極限定理,疊加均分分布隨機(jī)數(shù)產(chǎn)生;另一種基于公式變換,其產(chǎn)生的隨機(jī)數(shù)統(tǒng)計(jì)分布特性更好,但生成速度相對較慢。
統(tǒng)計(jì)近似抽樣[2]法的理論依據(jù)是著名的中心極限定理,它通過疊加均勻分布數(shù)據(jù)產(chǎn)生高斯噪聲,計(jì)算速度快,能夠滿足實(shí)時性要求。
當(dāng)隨機(jī)變量相互獨(dú)立且服從同一分布時,中心極限定理可表述為:設(shè)隨機(jī)變量X1,X2,…,Xn,…相互獨(dú)立,服從同一分布,且數(shù)學(xué)期望和方差均存在且方差大于零,即E(Xk)=μ、D(Xk)=σ2≠0(k=1,2,3…),則隨機(jī)變量近似服從正態(tài)分布N(nμ,nσ2),即近似服從標(biāo)準(zhǔn)正態(tài)分布N(0,1)[3]。
圖1 16級LFSR實(shí)現(xiàn)原理圖
服從(0,1)均勻分布的隨機(jī)序列X1,X2,…,Xn,…的數(shù)學(xué)期望和方差分別為則根據(jù)上述定理,服從N(0,1)正態(tài)分布。取合適的n值,就能產(chǎn)生N(0,1)正態(tài)分布的隨機(jī)數(shù),試驗(yàn)表明,n取6時就基本滿足要求,點(diǎn)數(shù)越大,高斯分布效果越好,實(shí)際實(shí)現(xiàn)方法要考慮硬件資源的限制,在效率和結(jié)果之間作出折中。
要得到兩路相互獨(dú)立、服從正態(tài)分布的高斯噪聲,可將上面方法得到的隨機(jī)數(shù)乘以正交的隨機(jī)相位即可。即若x是上述方法得到的一個隨機(jī)數(shù),y是服從(0,1)均勻分布的隨機(jī)數(shù),那么
是相互獨(dú)立、服從正態(tài)分布的隨機(jī)變量,這樣就得到了兩路正交高斯噪聲。
變換抽樣法[2]基于Box-Muller變換[4],該方法可以產(chǎn)生分布特性非常好的高斯噪聲,但與統(tǒng)計(jì)近似抽樣法相比,其計(jì)算相對復(fù)雜,可用于低實(shí)時要求的系統(tǒng)中。
Box-Muller變換可表述為:設(shè) x,y是兩個相互獨(dú)立的服從(0,1)均勻分布的隨機(jī)數(shù),作如下變換:
則可得到兩個相互獨(dú)立、服從N(0,1)正態(tài)分布的隨機(jī)數(shù) m,n,即得到了兩路正交高斯噪聲。
ADSP-TS201S是ADI公司一款性能極高、600MHz時鐘頻率、靜態(tài)超標(biāo)量處理器[5]。處理器集成了24Mbits大容量存儲器,兼有ASIC和FPGA的信號處理性能,具有高度的可編程性與靈活性。處理器將非常寬的存儲器和雙運(yùn)算模塊組合在一起,具有優(yōu)異的并行處理能力,廣泛應(yīng)用于無線通信、圖像處理等領(lǐng)域。
TS201S內(nèi)部結(jié)構(gòu)可分成內(nèi)核和I/O接口兩部分,如圖2。這兩部分通過4條總線(J、K、I、S)傳送數(shù)據(jù)、地址和控制信號。內(nèi)核部分包括程序控制器、數(shù)據(jù)地址產(chǎn)生器和XY雙運(yùn)算模塊。程序控制器提供完全可中斷的編程模式,支持匯編語言、C/C++編程,支持10指令周期流水線,IAB可預(yù)存5條指令,BTB可減小分支跳轉(zhuǎn)延遲;數(shù)據(jù)地址產(chǎn)生器包含兩個ALU,支持立即尋址和間接尋址,支持位反序和環(huán)形緩沖尋址;雙運(yùn)算模塊能夠獨(dú)立或同時工作以實(shí)現(xiàn)SIMD,每個周期每個運(yùn)算模塊可執(zhí)行2條運(yùn)算指令。I/O接口包括內(nèi)部存儲器、JTAG口、外部接口、DMA控制器和鏈路口。內(nèi)部存儲器為大小是24Mbits的DRAM;JTAG接口可用于片上仿真;外部接口包括主機(jī)接口、SRAM接口、多處理器接口、EPROM接口;DMA通道共14個,可無需處理器干預(yù)下完成設(shè)備間的數(shù)據(jù)傳輸;雙向鏈路口采用低壓差分信號(LVDS)鏈路口技術(shù),可達(dá)4Gbps的數(shù)據(jù)吞吐率[5]。
TS201S支持32位和40位浮點(diǎn)運(yùn)算以及8、16、32、64位定點(diǎn)運(yùn)算。每周期最多可執(zhí)行 4條指令,在600MHz的時鐘頻率下,可達(dá)到每秒48億次乘加運(yùn)算(GMACS)和每秒36億次浮點(diǎn)運(yùn)算(GFLOPS)的速度[5]。
圖2 ADSP-TS201S內(nèi)部結(jié)構(gòu)框圖
根據(jù)實(shí)時性要求,選擇線性同余法和統(tǒng)計(jì)近似抽樣法組合來產(chǎn)生實(shí)時高斯噪聲。線性同余法產(chǎn)生(0,65536)均勻分布的16位隨機(jī)整數(shù),結(jié)果除以65536即得(0,1)均勻分布隨機(jī)數(shù),此步驟合并到統(tǒng)計(jì)近似抽樣階段一次完成,以節(jié)約時間。試驗(yàn)表明,統(tǒng)計(jì)近似抽樣法疊加點(diǎn)數(shù)n取6已基本達(dá)到要求。
為了進(jìn)一步提高噪聲生成速度,此算法還采取了如下優(yōu)化措施:
(1)正交隨機(jī)相位正余弦值通過查表獲得,事先建立了4096間隔[0,2π]正余弦表文件“cos-sin-4096.dat”,用到時直接隨機(jī)讀取;
(2)統(tǒng)計(jì)近似抽樣法6點(diǎn)均勻分布隨機(jī)數(shù)累加分兩組并行計(jì)算,節(jié)約處理器時間;
(3)利用ADSP-TS201S一指令行周期最多可以執(zhí)行4條指令的特點(diǎn)[5],縮減程序長度,將可以并發(fā)執(zhí)行的指令放在同一指令行,節(jié)約指令執(zhí)行時間;
(4)利用ADSP-TS201S兩個運(yùn)算塊可同時獨(dú)立進(jìn)行計(jì)算的特點(diǎn)[5],減少循環(huán)次數(shù),要得到正交兩路各2048點(diǎn),共4096點(diǎn)的噪聲數(shù)據(jù),只要循環(huán)1024次,每次生成兩對正交共4點(diǎn)噪聲。
具體程序代碼如下:
圖3 正交兩路高斯噪聲統(tǒng)計(jì)分布圖和時域波形圖
經(jīng)試驗(yàn),程序在600MHz ADSP-TS201S上執(zhí)行時間約100μ s,即能以10kHz的頻率連續(xù)產(chǎn)生兩路各2048點(diǎn)不重復(fù)的正交高斯噪聲。通過調(diào)試器dump出-noiseout[4096]的噪聲數(shù)據(jù),用Matlab畫出的高斯噪聲的統(tǒng)計(jì)分布圖和時域波形如圖3所示,可見算法滿足實(shí)時性要求和分布特性要求。
討論了幾種均勻分布隨機(jī)數(shù)和正態(tài)分布隨機(jī)數(shù)的快速生成算法,并在ADSP-TS201S上實(shí)現(xiàn)了其中一種算法,實(shí)時產(chǎn)生了滿足分布要求的高斯噪聲。該方法已成功應(yīng)用到一款基于ADSP-TS201S的雷達(dá)模擬源中,并滿足參數(shù)要求。今后的工作可在程序優(yōu)化部分作出改進(jìn),以進(jìn)一步提高實(shí)時性。
[1]王衛(wèi)江.實(shí)時高斯噪聲產(chǎn)生方法[J].Chinese Journal of Scientific Instrument,2006,27(11):1523-1526.
[2]方崇智,蕭德云.過程識別[M].北京:清華大學(xué)出版社,1994:45-60.
[3]Marsaglia G.A Current View of Random Number Generator[J].Ann.Appl.Prob.,1991,(3):461-481.
[4]李裕奇,何平.概率論與數(shù)理統(tǒng)計(jì)[M].北京:國防工業(yè)出版社,2004.
[5]劉書明,羅勇江.ADSP TS20XS系列DSP原理與應(yīng)用設(shè)計(jì)[M].北京:電子工業(yè)出版社,2007.