高斯噪声实时产生算法及在ADSP-TS201S上的实现

2010-06-29李伟,杨斌

成都信息工程大学学报 2010年2期

李伟, 杨斌

(西南交通大学信息科学与技术学院,四川成都610031)

1 引言

随机噪声在通信等过程中往往是有害并需设法消除的无用信号,然而在一些情况下,也需要人为产生噪声,如自适应控制系统、系统辨识等领域。正态分布的高斯噪声是最常用的一种噪声。随机噪声由随机数表达,因此随机噪声的产生就归结于随机数的产生问题。由随机信号理论可知,在(0,1)上服从均匀分布的随机数经过一定的变换,可产生服从N(0,1)正态分布的高斯随机数[1]。

随机数的产生大致可分为查表法、物理法(如真空管射枪噪声)、数学递推法[2]。计算机产生随机数用到的是数学递推法,字长和精度有限的计算机或DSP产生真正的随机数是不可能的,因而现有的随机数产生算法产生的随机数均是伪随机数,即所得随机数都有一定的周期性。实际应用中也不必采用真正随机的数,达到一定程度随机性的伪随机数就可以满足一般要求了。

文中介绍了3种产生(0,1)均匀分布随机数的算法,随后介绍了2种用(0,1)均匀分布随机数生成正态分布随机数的快速生成方法,最后在ADSP-TS201S上实现了实时高斯噪声生成算法,能以10kHz的频率连续产生两路2048对不重复的正交高斯噪声数据。

2 均匀分布随机数的产生

均匀分布随机数产生的算法很多,如乘同余法、线性同余法、移位寄存器法、迭代取中法、加同余法、二次剩余法、进位加/借位减/进位乘(AWC/SWB/MWC)发生器法[3]等。根据实时性要求,主要介绍3种相对产生随机数速度较快的算法,3种算法产生的随机数均在(0,2L)区间上服从均匀分布,结果除以常数2L,即可得到(0,1)均匀分布的随机数。

2.1 乘同余法

乘同余法用如下的递推同余式产生随机序列：

其中,xi为当前随机数,xi-1为前一时刻随机数;M为模值,取2的幂,即 M=2L,L＞2的整数;A为乘因子,一般取小于M 且模8等于3或5的奇数;初值 x0取0～2L之间任意奇数。可以证明,该算法产生的随机数序列的周期为2L-2[2]。

2.2 线性同余法

线性同余法又称为混合同余法,递推公式如下：

这是一种类似乘同余的算法,只多了个增量c。根据随机数理论,A的后3位一般为“偶数-2-1”形式,如421、821等;初值 x0取0～2L之间任意奇数;Coveyou给出增量c满足的一个等式可以在 A取值不大的条件下,有效地减小序列的相关系数,而且不占用过多的机器时间。用该方法产生的随机序列的周期为2L[1]。

2.3 移位寄存器法

这是一种基于m序列理论的算法,可以利用线性反馈移位寄存器(LFSR)产生随机序列。根据m序列本原特征多项式,采取合适的回馈,M级LFSR产生的随机数的最大周期为2M-1。M阶m序列本原特征多项式可以表示为：

⊕表示模2或者异或运算,系数ci取0或1,转化为最大周期LFSR就是根据系数在LFSR相应位抽头。例如(16,12,3,1,0)表示16级LFSR的本原特征多项式为：

则应在LFSR16、12、3、1位抽头,如图1所示。每次移位运算后,把整个的16位数作为随机数输出。同样,也可以根据(32,7,5,3,2,1,0)表示的32级LFSR的本原特征多项式得到32位随机数。

3 正态分布随机数的产生

正态分布随机数可由(0,1)均匀分布随机数产生,下面介绍两种算法,一种基于中心极限定理,叠加均分分布随机数产生;另一种基于公式变换,其产生的随机数统计分布特性更好,但生成速度相对较慢。

3.1 统计近似抽样法

统计近似抽样[2]法的理论依据是著名的中心极限定理,它通过叠加均匀分布数据产生高斯噪声,计算速度快,能够满足实时性要求。

当随机变量相互独立且服从同一分布时,中心极限定理可表述为：设随机变量X1,X2,…,Xn,…相互独立,服从同一分布,且数学期望和方差均存在且方差大于零,即E(Xk)=μ、D(Xk)=σ2≠0(k=1,2,3…),则随机变量近似服从正态分布N(nμ,nσ2),即近似服从标准正态分布N(0,1)[3]。

图1 16级LFSR实现原理图

服从(0,1)均匀分布的随机序列X1,X2,…,Xn,…的数学期望和方差分别为则根据上述定理,服从N(0,1)正态分布。取合适的n值,就能产生N(0,1)正态分布的随机数,试验表明,n取6时就基本满足要求,点数越大,高斯分布效果越好,实际实现方法要考虑硬件资源的限制,在效率和结果之间作出折中。

要得到两路相互独立、服从正态分布的高斯噪声,可将上面方法得到的随机数乘以正交的随机相位即可。即若x是上述方法得到的一个随机数,y是服从(0,1)均匀分布的随机数,那么

是相互独立、服从正态分布的随机变量,这样就得到了两路正交高斯噪声。

3.2 变换抽样法

变换抽样法[2]基于Box-Muller变换[4],该方法可以产生分布特性非常好的高斯噪声,但与统计近似抽样法相比,其计算相对复杂,可用于低实时要求的系统中。

Box-Muller变换可表述为：设 x,y是两个相互独立的服从(0,1)均匀分布的随机数,作如下变换：

则可得到两个相互独立、服从N(0,1)正态分布的随机数 m,n,即得到了两路正交高斯噪声。

4 高斯噪声实时产生算法在ADSP-TS201S上的实现

4.1 ADSP-TS201S介绍

ADSP-TS201S是ADI公司一款性能极高、600MHz时钟频率、静态超标量处理器[5]。处理器集成了24Mbits大容量存储器,兼有ASIC和FPGA的信号处理性能,具有高度的可编程性与灵活性。处理器将非常宽的存储器和双运算模块组合在一起,具有优异的并行处理能力,广泛应用于无线通信、图像处理等领域。

TS201S内部结构可分成内核和I/O接口两部分,如图2。这两部分通过4条总线(J、K、I、S)传送数据、地址和控制信号。内核部分包括程序控制器、数据地址产生器和XY双运算模块。程序控制器提供完全可中断的编程模式,支持汇编语言、C/C++编程,支持10指令周期流水线,IAB可预存5条指令,BTB可减小分支跳转延迟;数据地址产生器包含两个ALU,支持立即寻址和间接寻址,支持位反序和环形缓冲寻址;双运算模块能够独立或同时工作以实现SIMD,每个周期每个运算模块可执行2条运算指令。I/O接口包括内部存储器、JTAG口、外部接口、DMA控制器和链路口。内部存储器为大小是24Mbits的DRAM;JTAG接口可用于片上仿真;外部接口包括主机接口、SRAM接口、多处理器接口、EPROM接口;DMA通道共14个,可无需处理器干预下完成设备间的数据传输;双向链路口采用低压差分信号(LVDS)链路口技术,可达4Gbps的数据吞吐率[5]。

TS201S支持32位和40位浮点运算以及8、16、32、64位定点运算。每周期最多可执行 4条指令,在600MHz的时钟频率下,可达到每秒48亿次乘加运算(GMACS)和每秒36亿次浮点运算(GFLOPS)的速度[5]。

图2 ADSP-TS201S内部结构框图

4.2 算法实现

根据实时性要求,选择线性同余法和统计近似抽样法组合来产生实时高斯噪声。线性同余法产生(0,65536)均匀分布的16位随机整数,结果除以65536即得(0,1)均匀分布随机数,此步骤合并到统计近似抽样阶段一次完成,以节约时间。试验表明,统计近似抽样法叠加点数n取6已基本达到要求。

为了进一步提高噪声生成速度,此算法还采取了如下优化措施：

(1)正交随机相位正余弦值通过查表获得,事先建立了4096间隔[0,2π]正余弦表文件“cos-sin-4096.dat”,用到时直接随机读取;

(2)统计近似抽样法6点均匀分布随机数累加分两组并行计算,节约处理器时间;

(3)利用ADSP-TS201S一指令行周期最多可以执行4条指令的特点[5],缩减程序长度,将可以并发执行的指令放在同一指令行,节约指令执行时间;

(4)利用ADSP-TS201S两个运算块可同时独立进行计算的特点[5],减少循环次数,要得到正交两路各2048点,共4096点的噪声数据,只要循环1024次,每次生成两对正交共4点噪声。

具体程序代码如下：

图3 正交两路高斯噪声统计分布图和时域波形图

经试验,程序在600MHz ADSP-TS201S上执行时间约100μ s,即能以10kHz的频率连续产生两路各2048点不重复的正交高斯噪声。通过调试器dump出-noiseout[4096]的噪声数据,用Matlab画出的高斯噪声的统计分布图和时域波形如图3所示,可见算法满足实时性要求和分布特性要求。