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

1、通信系统的计算机模拟第九讲1随机信号的产生与处理 信道噪声、干扰和衰落等随机波形级精确地仿真这些系统，随机影响建立准确的模型其基本构建模块是随机数发生器。如何产生采样后的随机波形（信号、干扰和噪声等），以用于仿真。在仿真环境下，所有的随机过程必须用随机变量序列来表示。随机数发生器产生的是“伪随机序列”，因为尽管它们是确定的，但在具体应用中会呈现随机性。精度随应用不同：产生一个波形来表示锁相环鉴相器输入端的噪声，对一个输入SNR为50dB的系统建立噪声波形模型所要求的精度要比输入端SNR为8dB的系统高得多。2问题在（0，1）上产生均匀分布且不相关的随机数将不相关且均匀分布的随机数映射成不相关的

2、，具有任意（需要的）概率密度函数（pdf）的随机产生不相关的、具有高斯型pdf的随机数产生相关的、具有高斯型pdf的随机数产生相关的、具有任意（需要的）pdf的随机数伪噪声（PN）序列的产生以及几种应用于随机数序列的计算方法。3平稳与遍历性过程信号、噪声和干扰而产生的样本函数通常假设为各态历经的。原因：依次处理系统中波形的时域样本，系统中的每个点上只有一个波形（样本函数）。假设仿真所处理的波形是其内在统计模型定义的总体（Emsemble）中的一个典型成员，这时各种统量如各阶矩、信噪比和误比特率，就可以当作时间平均来计算。仿真计算得到的时间平均等于总体均值。隐含的假设条件是对应的随机过程为遍历性

3、的。4例7-1 假设随机过程的样本函数表达式为 是对应随机试验的样本空间中的一个输出，每一个映射为一个相位。再假定对应的随机试验就是从均匀随机数发生器输出端抽取一个数，抽取得到的结果 ，这里ui 在（0，1）区间上均匀分布。然后映射成相位 。当幅度A和频率f固定， 的值便决定了波形。K-5f = 1; % frequency of sinusoidfs = 100; % sampling frequencyt = (0:200)/fs; % time vectorfor i=1:20 x(:,i) = cos(2*pi*f*t+rand(1)*2*pi); y(:,i) = cos(2*pi*

4、f*t+rand(1)*pi/2); z(:,i) = (1+rand(1)*cos(2*pi*f*t);endsubplot(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)67由x(t)、y(t)和z(t)组成的所有样本函数的时间平均都等于零。大量满足0ti2的，计算出来的的总体均值近似0；Y(t） 0.875 1.875 0.375 1.375 0.125 0.625 1.125 1.625 -周期平稳Z（t） 8

5、例7-2 在后面的工作中，我们经常要用到数字调制器的模型。这些调制器的基本构建模块是函数random-binary，它产生电平值为+1或-1的二进制波形，产生的比特数以及每比特的采样数是该函数的参数。 9function 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,n

6、bits);for m=1:nbits for n=1:nsamples index = (m-1)*nsamples + n; x(1,index) = (-1)bits(m); endend10函数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 = r

7、andom_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)1213均匀随机数发生器具有均匀概率密度函数的随机变量很容易转换成具有其他所需pdf的随

8、机变量，因此，要产生一个具有特定pdf的随机变量，首先产生一个在（0，1）区间均匀分布的随机变量。通常，先产生一个列介于0和M之间的数（整数），然后将序列中每个元素除以M。实现随机数发生器最常用的方法是线性同余。（linear congruence）。14线性同余线性同余发生器定义为如下运算a和c分别称作乘子和增量，m叫做模数。这是一个确定性的序列算法，能依次产生连续的值。初始值记为x0，称作线性同余发生器的种子数（seed nuber）如果a、c、m和都是整数，则LCG产生的所有数也都是整数。由于对进行运算，式（7-3）至多可产生m个不同的整数。发生器输出的一个理想特征是它应具备很长的周期，

9、从而在序列重复前，输出序列能产生最多数目的整数。对于给定值，当周期最大时，我们称发生器是全周期（full period）的。此外，具体仿真程序的应用对LCG会提出其他的要求，我们通常要求样本和互不相关。可能还要求LCG的输出能通过其他统计测试。LCG可以采用多种不同的形式15方法A：混合同余算法最通用的同余算法就是c0的“混合”同余算法。之所以称之为混合算法，是因为在求解的过程中要同时用到乘法与加法。混合线性算法具有式（7-4）的形式c 0时，发生器的最大周期为m。当且仅当满足以下特性时才能达到这个周期：*增量c与m互质。换句话说，与没有素公因子（prime factor）。*a-1是p的倍数

10、，则p是m的任一素因子。*如果m是4的倍数，则a1是4的倍数。（7-4）16例7-3我们要设计一个周期m5000的混合同余发生器。因为令c等于除2和5之外的素数的乘积，可保证c和m互质，这满足了第一条特性。其中可能的一个设定是：现在必须选定a的值。若满足则满足第二特性，p1=2，p2=5,k1,k2任意17cont因为4是m=5000的因数，设定其中k是整数。当 ，则 ，于是可以满足第三特性18例7-4这个例子将说明前一例子中设计的LCG确实具有周期。在以下MATLAB程序中，输入种子数，并运行程序直到种子数再次出现。设产生n个整数，如果nm而种子数没有再次出现，则认为发生器进入了一个重复产生

11、短序列的循环中。MATLAB程序如下：19a = 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)end2021方法B：具有素模数的乘性算法乘性发生器的定义式为它是增量等于零的混合算

12、法。注意，c=0时，xi不能为零，因此，全周期是而不是前面那种情况下的。若满足以下特性，乘性算法能产生全周期序列；*m是素数（通常要求m取较大值）*m为的mod(m)本原元素如果除i=m-1了外没有更小的i值，使ai-1是m的倍数，则a为的mod(m)本原元素。换句话说，a为mod(m)的本原元素，如果满足22方法C：具有非素模数的乘性算法模数不是素数的同余算法中最重要的情况是m等于2的幂，即 这里n为整数。对式（7-16）定义的算法，最大周期为2n/4=2n-2。如果满足以下条件，可取得这个周期：*乘子 结果为3或5*种子数x0是奇数（7-16）23由于两个奇数的乘积是奇数，可以推出，若x0

13、是奇数，则由式（7-16）产生的所有值都是奇数，也就是说，产生的值都不是偶数，这使周期减为原来的一半。式（7-16）产生的奇数分成两组，其中只有一组由给定的种子数产生，这样又使周期减小了一倍。产生的这组奇数一般跟所选的种子数有关采用m=2n的好处是整数溢出能用于运算，这样缩短了计算时间，这样的结果确实比较理想，但所得的程序不易移植。24随机数发生器的测试 前面一节为我们提供了工具，来产生在0与1之间均匀分布的伪随机数。到目前为止，我们只考虑了由LCG产生的序列的周期。于特定的应用来说，还要满足其他理想的特性，至少，我们希望序列要满足相关（白噪声）。具体的应用可能还必须满足其他的要求。现已开发了

14、多个程序来测试某一给定序列的随机性（randomness），其中最常用的Chi-方（Chi-square）测试、Kolomogorov-Smirnov测试和谱测试。 25谱测试似乎是功能最强大的。在后面很多的例子中，要满足的一个最重要的特性是给定序列的元素相互独立或至少互不相关。为此，我们来考虑两个非常简单的测试：散点图（scatterplots）和Durbin-Watson测试。需要指出的是给定序列的性质适用于完整的（全周期）序列。如果只用到序列的一部分，则性质不再存在。 26散点图 所谓散点图就是x(i+1)与x(i)的函数关系图，它表示了随机数发生器经验式的质量指标。本例中所考虑的两个随

15、机数发生器定义如下：27可以验证是全周期的m = 2048; c = 1; seed = 1; % default values of m and ca1 = 65; a2 = 1229; % multiplier valuesix1 = seed; ix2 = seed; % initialize algorithmx1 = zeros(1,m); x2 = zeros(1,m); % initialize arraysfor i=1:m ix1 = rem(ix1*a1+c),m); x1(i) = ix1/m; ix2 = rem(ix2*a2+c),m); x2(i) = ix2/m;

16、endsubplot(1,2,1)y1 = x1(1,2:m),x1(1,1);plot(x1,y1,.) % plot results for a1subplot(1,2,2)y2 = x2(1,2:m),x2(1,1);plot(x2,y2,.) % plot results for a22829Durbin-Watson测试Durbin-Watson独立性测试可以通过计算如下Durbin参数来完成：其中Xn是一个零均值（zero-mean）随机变量4。下面会讲到当D值在2附近时，Xn和Xn-1相关性很小。30Durbin-Watson测试为了说明Durbin-Watson测试的性质，假设

17、和相关并且是一个遍历性过程。为了简化记号，假设N足够大，使得N=N-1，式（7-19）可写成： （7-20）由假定Xn和Xn-1相关，可令其中X和Z不相关，是X和Y的相关系数。注意，X、Y和Z有相同的方差，用表示。将式（7-21）代入式（7-20）得：31X和Z不相关且均值为零，所以中间项等于零。由于和方差相同，所以有因为，Durbin参数在0与4之间，且有=0时D=2。D0表示值为负。以下的MATLAB函数计算Durbin参数值：32function D = durbin(x)N = length(x); % length of input vectory = x-mean(x); % re

18、move dcydiff = y(2:N)-y(1:(N-1); % numerator summandNum = sum(ydiff.*ydiff); % numerator factor of D Den = sum(y.*y); % denominator factor of DD = Num/Den; % Durbin factor33例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);fo

19、r i=1:m ix1 = rem(ix1*a1+c),m); x1(i) = ix1; ix2 = rem(ix2*a2+c),m); x2(i) = ix2;endD1 = c7_Durbin(x1); D2 = c7_Durbin(x2); % calculate Durbin parametersrho1 = 1-D1/2; rho2 = 1- D2/2; % calculate correlationtext1 = The value of D1 is ,num2str(D1), and rho1 is ,num2str(rho1),.;text2 = The value of D2

20、 is ,num2str(D2), and rho2 is ,num2str(rho2),.; disp(text1)disp(text2)34 c7_LCDemo2The 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）算法。最低标准算法具有以下性质：*全周

21、期*能通过所有适用的随机性统计测试*易于从一台计算机移植到另一台计算机如果使用一个最低标准算法，不用担心算法本身的正确性，但必须保证在给定计算环境下算法得到正确的实现。编程时要特别注意的是，由算法产生的所有的数都要能唯一地表示出来。36Lewis, Goodman和Miller最低标准其中m是Mersenne素数231-1，这个值最初由Lehmer推荐使用，很容易在32位的计算机以上以整型运算实现，如果尾数超过31位，则以浮点运算实现。37Wichmann-Hill算法希望获得一些具有长周期的随机数发生器。要构造具有长周期的波型，一个行之有效的办法是将几个周期相差很小的周期波型相加。Cos2(

22、1)t的周期是1s， Cos2(1.0001)t的周期是1000/10001s，略小于1s 。两者的复合波形可写成如下形式：其周期为10000s或约合2.78个小时。复合波形的一个周期内，第一个分量经过10000个周期，第二个分量经过10001个周期。还可再添加更多的分量。38可以将同样的方法用于LCG，即把几个周期相差不大的随机数发生器合成在一起。Wichmann-Hill算法可能是最有名的一个合成随机数发生器的例子。该算法用到的三个分量发生器，定义如下： （7-26） （7-28）39这三个分量发生器都是全周期发生器，它们合成得到输出为40Wichmann-Hill算法等于价于乘性LCG，它的倍数为模数为显然，m不是素数，周期小于m-