通信系统的计算机模拟第九讲

,*,TechniTemplate,Test,Test,Test,Test,通信系统的计算机模拟,第九讲,1,随机信号的产生与处理,信道噪声、干扰和衰落等随机,波形级精确地仿真这些系统，随机影响建立准确的模型,其基本构建模块是随机数发生器。如何产生采样后的随机波形（信号、干扰和噪声等），以用于仿真。在仿真环境下，所有的随机过程必须用随机变量序列来表示。,随机数发生器产生的是“伪随机序列”，因为尽管它们是确定的，但在具体应用中会呈现随机性。,精度随应用不同：产生一个波形来表示锁相环鉴相器输入端的噪声，对一个输入,SNR,为,50dB,的系统建立噪声波形模型所要求的精度要比输入端,SNR,为,8dB,的系统高得多。,2,问题,在（,0,，,1,）上产生均匀分布且不相关的随机数,将不相关且均匀分布的随机数映射成不相关的，具有任意（需要的）概率密度函数（,pdf,）的随机,产生不相关的、具有高斯型,pdf,的随机数,产生相关的、具有高斯型,pdf,的随机数,产生相关的、具有任意（需要的）,pdf,的随机数,伪噪声（,PN,）序列的产生以及几种应用于随机数序列的计算方法。,3,平稳与遍历性过程,信号、噪声和干扰而产生的样本函数通常假设为各态历经的。,原因：依次处理系统中波形的时域样本，系统中的每个点上只有一个波形（样本函数）。,假设仿真所处理的波形是其内在统计模型定义的总体（,Emsemble,）中的一个典型成员，这时各种统量如各阶矩、信噪比和误比特率，就可以当作时间平均来计算。仿真计算得到的时间平均等于总体均值。隐含的假设条件是对应的随机过程为遍历性的。,4,例,7-1,假设随机过程的样本函数表达式为,是对应随机试验的样本空间中的一个输出，每一个映射为一个相位。,再假定对应的随机试验就是从均匀随机数发生器输出端抽取一个数，抽取得到的结果 ，这里,u,i,在（,0,，,1,）区间上均匀分布。然后映射成相位 。当幅度,A,和频率,f,固定， 的值便决定了波形。,K-,5,f = 1; % frequency of sinusoid,fs = 100; % sampling frequency,t = (0:200)/fs; % time vector,for i=1:20,x(:,i) = cos(2*pi*f*t+rand(1)*2*pi);,y(:,i) = cos(2*pi*f*t+rand(1)*pi/2);,z(:,i) = (1+rand(1)*cos(2*pi*f*t);,end,subplot(3,1,1); plot(t,x,k); ylabel(x(t),subplot(3,1,2); plot(t,y,k); ylabel(y(t),subplot(3,1,3); plot(t,z,k); ylabel(z(t),6,7,由,x(t),、,y(t),和,z(t),组成的所有样本函数的时间平均都等于零。,大量满足,0t,i,2,的，计算出来的的总体均值近似,0,；,Y(t,）,0.875 1.875 0.375 1.375 0.125 0.625 1.125 1.625 -,周期平稳,Z,（,t,）,8,例,7-2,在后面的工作中，我们经常要用到数字调制器的模型。这些调制器的基本构建模块是函数,random-binary,，它产生电平值为,+1,或,-1,的二进制波形，产生的比特数以及每比特的采样数是该函数的参数。,9,function x, bits = random_binary(nbits,nsamples),%,This function genrates a random binary waveform of length nbits,% sampled at a rate of nsamples/bit.,x = zeros(1,nbits*nsamples);,bits = round(rand(1,nbits);,for m=1:nbits,for n=1:nsamples,index = (m-1)*nsamples + n;,x(1,index) = (-1)bits(m);,end,end,10,函数,random-binary,的应用,函数,random-binary,可以仿真多个数字调制器，例如，可用如下,MATLAB,语句仿真一个,QPSK,调制器：,x=random_binary(nbits, nsamples)+i*random_binary (nbits, nsaples);,11,产生一个长,10,比特的,QPSK,信号，采样频率为每比特,8,个采样点,nbits = 10; nsamples = 8;,x = random_binary(nbits,nsamples)+i*random_binary(nbits,nsamples);,xd = real(x); xq = imag(x);,subplot(2,1,1),stem(xd,.); grid; axis(0 80 -1.5 1.5);,xlabel(Sample Index); ylabel(xd),subplot(2,1,2),stem(xq,.); grid; axis(0 80 -1.5 1.5);,xlabel(Sample Index); ylabel(xq),12,13,均匀随机数发生器,具有均匀概率密度函数的随机变量很容易转换成具有其他所需,pdf,的随机变量，因此，要产生一个具有特定,pdf,的随机变量，首先产生一个在（,0,，,1,）区间均匀分布的随机变量。通常，先产生一个列介于,0,和,M,之间的数（整数），然后将序列中每个元素除以,M,。实现随机数发生器最常用的方法是线性同余。（,linear congruence,）。,14,线性同余,线性同余发生器定义为如下运算,a,和,c,分别称作乘子和增量，,m,叫做模数。,这是一个确定性的序列算法，能依次产生连续的值。,初始值记为,x,0,，称作线性同余发生器的种子数（,seed nuber,）,如果,a,、,c,、,m,和都是整数，则,LCG,产生的所有数也都是整数。由于对进行运算，式（,7-3,）至多可产生,m,个不同的整数。发生器输出的一个理想特征是它应具备很长的周期，从而在序列重复前，输出序列能产生最多数目的整数。,对于给定值，当周期最大时，我们称发生器是全周期（,full period,）的。,此外，具体仿真程序的应用对,LCG,会提出其他的要求，,我们通常要求样本和互不相关。,可能还要求,LCG,的输出能通过其他统计测试。,LCG,可以采用多种不同的形式,15,方法,A,：混合同余算法,最通用的同余算法就是,c0,的“混合”同余算法。之所以称之为混合算法，是因为在求解的过程中要同时用到乘法与加法。混合线性算法具有式（,7-4,）的形式,c ,0,时，发生器的最大周期为,m,。当且仅当满足以下特性时才能达到这个周期：,*增量,c,与,m,互质。换句话说，与没有素公因子（,prime factor,）。,*,a-1,是,p,的倍数，则,p,是,m,的任一素因子。,*如果,m,是,4,的倍数，则,a,1,是,4,的倍数。,（,7-4,）,16,例,7-3,我们要设计一个周期,m,5000,的混合同余发生器。因为,令,c,等于除,2,和,5,之外的素数的乘积，可保证,c,和,m,互质，这满足了第一条特性。其中可能的一个设定是：,现在必须选定,a,的值。若满足,则满足第二特性，,p1=2,，,p2=5,k1,k2,任意,17,cont,因为,4,是,m=5000,的因数，设定,其中,k,是整数。当 ，则 ，于是,可以满足第三特性,18,例,7-4,这个例子将说明前一例子中设计的,LCG,确实具有周期。在以下,MATLAB,程序中，输入种子数，并运行程序直到种子数再次出现。设产生,n,个整数，如果,nm,而种子数没有再次出现，则认为发生器进入了一个重复产生短序列的循环中。,MATLAB,程序如下：,19,a = input(Enter multiplier a );,c = input(Enter offset c );,m = input(Enter modulus m );,seed = input(Enter seed );,n=1; ix = rem(seed*a+c),m);,while (ix=seed)&(nm,disp(Caught in a loop.),else,text = The period is ,num2str(n,15),.;,disp(text),end,20,21,方法,B,：具有素模数的乘性算法,乘性发生器的定义式为,它是增量等于零的混合算法。注意，,c=0,时，,x,i,不能为零，因此，全周期是而不是前面那种情况下的。若满足以下特性，乘性算法能产生全周期序列；,*,m,是素数（通常要求,m,取较大值）,*,m,为的,mod(m),本原元素,如果除,i=m-1,了外没有更小的,i,值，使,a,i,-1,是,m,的倍数，则,a,为的,mod(m),本原元素。换句话说，,a,为,mod(m),的本原元素，如果满足,22,方法,C,：具有非素模数的乘性算法,模数不是素数的同余算法中最重要的情况是,m,等于,2,的幂，即,这里,n,为整数。对式（,7-16,）定义的算法，最大周期为,2,n,/4=2,n-2,。如果满足以下条件，可取得这个周期：,*乘子 结果为,3,或,5,*,种子数,x,0,是奇数,（,7-16,）,23,由于两个奇数的乘积是奇数，可以推出，若,x,0,是奇数，则由式（,7-16,）产生的所有值都是奇数，也就是说，产生的值都不是偶数，这使周期减为原来的一半。,式（,7-16,）产生的奇数分成两组，其中只有一组由给定的种子数产生，这样又使周期减小了一倍。产生的这组奇数一般跟所选的种子数有关,采用,m=2,n,的好处是整数溢出能用于运算，这样缩短了计算时间，这样的结果确实比较理想，但所得的程序不易移植。,24,随机数发生器的测试,前面一节为我们提供了工具，来产生在,0,与,1,之间均匀分布的伪随机数。,到目前为止，我们只考虑了由,LCG,产生的序列的周期。,于特定的应用来说，还要满足其他理想的特性，至少，我们希望序列要满足,相关（白噪声）。,具体的应用可能还必须满足其他的要求。,现已开发了多个程序来测试某一给定序列的随机性（,randomness,），其中最常用的,Chi-,方（,Chi-square,）测试、,Kolomogorov-Smirnov,测试和谱测试。,25,谱测试似乎是功能最强大的。,在后面很多的例子中，要满足的一个最重要的特性是给定序列的元素相互独立或至少互不相关。,为此，我们来考虑两个非常简单的测试：散点图（,scatterplots,）和,Durbin-Watson,测试。,需要指出的是给定序列的性质适用于完整的（全周期）序列。如果只用到序列的一部分，则性质不再存在。,26,散点图,所谓散点图就是,x(i+1),与,x(i),的函数关系图，它表示了随机数发生器经验式的质量指标。本例中所考虑的两个随机数发生器定义如下：,27,可以验证是全周期的,m = 2048; c = 1; seed = 1; % default values of m and c,a1 = 65; a2 = 1229; % multiplier values,ix1 = seed; ix2 = seed; % initialize algorithm,x1 = zeros(1,m); x2 = zeros(1,m); % initialize arrays,for i=1:m,ix1 = rem(ix1*a1+c),m);,x1(i) = ix1/m;,ix2 = rem(ix2*a2+c),m);,x2(i) = ix2/m;,end,subplot(1,2,1),y1 = x1(1,2:m),x1(1,1);,plot(x1,y1,.) % plot results for a1,subplot(1,2,2),y2 = x2(1,2:m),x2(1,1);,plot(x2,y2,.) % plot results for a2,28,29,Durbin-Watson,测试,Durbin-Watson,独立性测试可以通过计算如下,Durbin,参数来完成：,其中,Xn,是一个零均值（,zero-mean,）随机变量,4,。,下面会讲到当,D,值在,2,附近时，,Xn,和,Xn-1,相关性很小。,30,Durbin-Watson,测试,为了说明,Durbin-Watson,测试的性质，假设和相关并且是一个遍历性过程。为了简化记号，假设,N,足够大，使得,N=N-1,，式（,7-19,）可写成：,（,7-20,）,由假定,Xn,和,Xn-1,相关，可令,其中,X,和,Z,不相关，,是,X,和,Y,的相关系数。注意，,X,、,Y,和,Z,有相同的方差，用,表示。将式（,7-21,）代入式（,7-20,）得：,31,X,和,Z,不相关且均值为零，所以中间项等于零。由于和方差相同，所以有,因为，,Durbin,参数在,0,与,4,之间，且有,=0,时,D=2,。,D0,表示,值为负。以下的,MATLAB,函数计算,Durbin,参数值：,32,function D = durbin(x),N = length(x); % length of input vector,y = x-mean(x); % remove dc,ydiff = y(2:N)-y(1:(N-1); % numerator summand,Num = sum(ydiff.*ydiff); % numerator factor of D,Den = sum(y.*y); % denominator factor of D,D = Num/Den; % Durbin factor,33,例,7-6,本例将计算例,7-5,中讨论到的两个噪声发生器的值。,m = 2048; c = 1; seed = 1;,a1 = 65; a2 = 1229;,ix1 = 1; ix2 = 1;,x1 = zeros(1,m); x2 = zeros(1,m);,for i=1:m,ix1 = rem(ix1*a1+c),m);,x1(i) = ix1;,ix2 = rem(ix2*a2+c),m);,x2(i) = ix2;,end,D1 = c7_Durbin(x1); D2 = c7_Durbin(x2); % calculate Durbin parameters,rho1 = 1-D1/2; rho2 = 1- D2/2; % calculate correlation,text1 = The value of D1 is ,num2str(D1), and rho1 is ,num2str(rho1),.;,text2 = The value of D2 is ,num2str(D2), and rho2 is ,num2str(rho2),.;,disp(text1),disp(text2),34, c7_LCDemo2,The value of D1 is 1.9925 and fho1 is 0.0037273,。,The value of D2 is 1.6037 and rho2 is 0.19814.,35,最低标准,表明给定,LCG,能通过各种随机性的统计测试，以便彻底测试它的质量，这是一项首要任务。对于产生长的序列来说更是如此。为了部分地解决这个问题，多种算法被确定为最低标准（,Minimum standard,）算法。最低标准算法具有以下性质：,*全周期,*能通过所有适用的随机性统计测试,*易于从一台计算机移植到另一台计算机,如果使用一个最低标准算法，不用担心算法本身的正确性，但必须保证在给定计算环境下算法得到正确的实现。编程时要特别注意的是，由算法产生的所有的数都要能唯一地表示出来。,36,Lewis, Goodman,和,Miller,最低标准,其中,m,是,Mersenne,素数,2,31,-1,，这个值最初由,Lehmer,推荐使用，很容易在,32,位的计算机以上以整型运算实现，如果尾数超过,31,位，则以浮点运算实现。,37,Wichmann-Hill,算法,希望获得一些具有长周期的随机数发生器。,要构造具有长周期的波型，一个行之有效的办法是将几个周期相差很小的周期波型相加。,Cos2,(1)t,的周期是,1s,，,Cos2,(1.0001)t,的周期是,1000/10001s,，略小于,1s,。两者的复合波形可写成如下形式：,其,周期为,10000s,或约合,2.78,个小时。复合波形的一个周期内，第一个分量经过,10000,个周期，第二个分量经过,10001,个周期。还可再添加更多的分量。,38,可以将同样的方法用于,LCG,，即把几个周期相差不大的随机数发生器合成在一起。,Wichmann-Hill,算法可能是最有名的一个合成随机数发生器的例子。该算法用到的三个分量发生器，定义如下：,（,7-26,）,（,7-28,）,39,这三个分量发生器都是全周期发生器，它们合成得到输出为,40,Wichmann-Hill,算法等于价于乘性,LCG,，它的倍数为,模数为,显然，,m,不是素数，周期小于,m-1,。周期近似为,710,12,，尽管小于,m,，但仍然很长。,Wichmann-Hill,被证明通过了所有标准的统计测试，仍然认为它是最低标准均匀随机数发生器。,41,MATLAB,实现,MATLAB,版本,5,以前，,rand,是由式（,7-24,）定义的最低标准随机数发生器要产生的是浮点数而不是缩放的整数，,Math Works,声称这个随机数发生器的周期超过,2,1492,，并且“相当确信”可以产生,eps,与,1-eps/2,之间的所有浮点数，其中,MATLAB,常数,eps,，如第,3,章中所介绍，等于,2,52,。新的随机数发生器只用到加法和减法运算，没有用到乘法或除法运算。,42,种子数与种子向量,MATLAB,运用种子数的方式。,“旧版本”,MATLAB,随机数发生器，使用单个种子数,“新版本”随机数发生器使用种子向量，称为随机数发生器的状态。,该向量由,35,个表示随机数发生器状态的元素组成，它们分别是,32,个浮点数、两个整数和一个标志位。,上面两种随机数发生器都可以使用，其中新的随机数发生器设定为默认的发生器。,旧版的随机数发生器可通过命令,RAND,（,seed 0,）,or RAND(seed, J),来调用。整数种子和状态向量：,43,注意事项,*使用者可以使用默认的种子数，也可以自己给定一个种子数。,*关闭并重新打开,MATLAB,可以将种子数复位到默认值。如果调用一个随机数发生顺次后，关闭,MATLAB,，接着重新打开，再调用随机数发生器次，那么两种情况下产生的个数是相同的。这一特性可以用来发挥,MATLAB,的优势，它允许我们重新产生结果完全相同的序列。这对测试很有用。,*系统时钟可用来随机化初始的种子数（具体细节见,MATLAB help,）。,*种子数存在一个缓冲器中，而不是,MATLAB,工作区。命令,clear all,对种子数不起作用。,44,45,