通信系統(tǒng)的計(jì)算機(jī)模擬第九講

1、通信系統(tǒng)的計(jì)算機(jī)模擬第九講1隨機(jī)信號(hào)的產(chǎn)生與處理 信道噪聲、干擾和衰落等隨機(jī)波形級(jí)精確地仿真這些系統(tǒng)，隨機(jī)影響建立準(zhǔn)確的模型其基本構(gòu)建模塊是隨機(jī)數(shù)發(fā)生器。如何產(chǎn)生采樣后的隨機(jī)波形（信號(hào)、干擾和噪聲等），以用于仿真。在仿真環(huán)境下，所有的隨機(jī)過(guò)程必須用隨機(jī)變量序列來(lái)表示。隨機(jī)數(shù)發(fā)生器產(chǎn)生的是“偽隨機(jī)序列”，因?yàn)楸M管它們是確定的，但在具體應(yīng)用中會(huì)呈現(xiàn)隨機(jī)性。精度隨應(yīng)用不同：產(chǎn)生一個(gè)波形來(lái)表示鎖相環(huán)鑒相器輸入端的噪聲，對(duì)一個(gè)輸入SNR為50dB的系統(tǒng)建立噪聲波形模型所要求的精度要比輸入端SNR為8dB的系統(tǒng)高得多。2問(wèn)題在（0，1）上產(chǎn)生均勻分布且不相關(guān)的隨機(jī)數(shù)將不相關(guān)且均勻分布的隨機(jī)數(shù)映射成不相關(guān)的

2、，具有任意（需要的）概率密度函數(shù)（pdf）的隨機(jī)產(chǎn)生不相關(guān)的、具有高斯型pdf的隨機(jī)數(shù)產(chǎn)生相關(guān)的、具有高斯型pdf的隨機(jī)數(shù)產(chǎn)生相關(guān)的、具有任意（需要的）pdf的隨機(jī)數(shù)偽噪聲（PN）序列的產(chǎn)生以及幾種應(yīng)用于隨機(jī)數(shù)序列的計(jì)算方法。3平穩(wěn)與遍歷性過(guò)程信號(hào)、噪聲和干擾而產(chǎn)生的樣本函數(shù)通常假設(shè)為各態(tài)歷經(jīng)的。原因：依次處理系統(tǒng)中波形的時(shí)域樣本，系統(tǒng)中的每個(gè)點(diǎn)上只有一個(gè)波形（樣本函數(shù)）。假設(shè)仿真所處理的波形是其內(nèi)在統(tǒng)計(jì)模型定義的總體（Emsemble）中的一個(gè)典型成員，這時(shí)各種統(tǒng)量如各階矩、信噪比和誤比特率，就可以當(dāng)作時(shí)間平均來(lái)計(jì)算。仿真計(jì)算得到的時(shí)間平均等于總體均值。隱含的假設(shè)條件是對(duì)應(yīng)的隨機(jī)過(guò)程為遍歷性

3、的。4例7-1 假設(shè)隨機(jī)過(guò)程的樣本函數(shù)表達(dá)式為 是對(duì)應(yīng)隨機(jī)試驗(yàn)的樣本空間中的一個(gè)輸出，每一個(gè)映射為一個(gè)相位。再假定對(duì)應(yīng)的隨機(jī)試驗(yàn)就是從均勻隨機(jī)數(shù)發(fā)生器輸出端抽取一個(gè)數(shù)，抽取得到的結(jié)果 ，這里ui 在（0，1）區(qū)間上均勻分布。然后映射成相位 。當(dāng)幅度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)組成的所有樣本函數(shù)的時(shí)間平均都等于零。大量滿足0ti2的，計(jì)算出來(lái)的的總體均值近似0；Y(t） 0.875 1.875 0.375 1.375 0.125 0.625 1.125 1.625 -周期平穩(wěn)Z（t） 8

5、例7-2 在后面的工作中，我們經(jīng)常要用到數(shù)字調(diào)制器的模型。這些調(diào)制器的基本構(gòu)建模塊是函數(shù)random-binary，它產(chǎn)生電平值為+1或-1的二進(jìn)制波形，產(chǎn)生的比特?cái)?shù)以及每比特的采樣數(shù)是該函數(shù)的參數(shù)。 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函數(shù)random-binary的應(yīng)用函數(shù)random-binary可以仿真多個(gè)數(shù)字調(diào)制器，例如，可用如下MATLAB語(yǔ)句仿真一個(gè)QPSK調(diào)制器：x=random_binary(nbits, nsamples)+i*random_binary (nbits, nsaples);11產(chǎn)生一個(gè)長(zhǎng)10比特的QPSK信號(hào)，采樣頻率為每比特8個(gè)采樣點(diǎn) 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均勻隨機(jī)數(shù)發(fā)生器具有均勻概率密度函數(shù)的隨機(jī)變量很容易轉(zhuǎn)換成具有其他所需pdf的隨

8、機(jī)變量，因此，要產(chǎn)生一個(gè)具有特定pdf的隨機(jī)變量，首先產(chǎn)生一個(gè)在（0，1）區(qū)間均勻分布的隨機(jī)變量。通常，先產(chǎn)生一個(gè)列介于0和M之間的數(shù)（整數(shù)），然后將序列中每個(gè)元素除以M。實(shí)現(xiàn)隨機(jī)數(shù)發(fā)生器最常用的方法是線性同余。（linear congruence）。14線性同余線性同余發(fā)生器定義為如下運(yùn)算a和c分別稱作乘子和增量，m叫做模數(shù)。這是一個(gè)確定性的序列算法，能依次產(chǎn)生連續(xù)的值。初始值記為x0，稱作線性同余發(fā)生器的種子數(shù)（seed nuber）如果a、c、m和都是整數(shù)，則LCG產(chǎn)生的所有數(shù)也都是整數(shù)。由于對(duì)進(jìn)行運(yùn)算，式（7-3）至多可產(chǎn)生m個(gè)不同的整數(shù)。發(fā)生器輸出的一個(gè)理想特征是它應(yīng)具備很長(zhǎng)的周期，

9、從而在序列重復(fù)前，輸出序列能產(chǎn)生最多數(shù)目的整數(shù)。對(duì)于給定值，當(dāng)周期最大時(shí)，我們稱發(fā)生器是全周期（full period）的。此外，具體仿真程序的應(yīng)用對(duì)LCG會(huì)提出其他的要求，我們通常要求樣本和互不相關(guān)。可能還要求LCG的輸出能通過(guò)其他統(tǒng)計(jì)測(cè)試。LCG可以采用多種不同的形式15方法A：混合同余算法最通用的同余算法就是c0的“混合”同余算法。之所以稱之為混合算法，是因?yàn)樵谇蠼獾倪^(guò)程中要同時(shí)用到乘法與加法?；旌暇€性算法具有式（7-4）的形式c 0時(shí)，發(fā)生器的最大周期為m。當(dāng)且僅當(dāng)滿足以下特性時(shí)才能達(dá)到這個(gè)周期：*增量c與m互質(zhì)。換句話說(shuō)，與沒(méi)有素公因子（prime factor）。*a-1是p的倍數(shù)

10、，則p是m的任一素因子。*如果m是4的倍數(shù)，則a1是4的倍數(shù)。（7-4）16例7-3我們要設(shè)計(jì)一個(gè)周期m5000的混合同余發(fā)生器。因?yàn)榱頲等于除2和5之外的素?cái)?shù)的乘積，可保證c和m互質(zhì)，這滿足了第一條特性。其中可能的一個(gè)設(shè)定是：現(xiàn)在必須選定a的值。若滿足則滿足第二特性，p1=2，p2=5,k1,k2任意17cont因?yàn)?是m=5000的因數(shù)，設(shè)定其中k是整數(shù)。當(dāng) ，則 ，于是可以滿足第三特性18例7-4這個(gè)例子將說(shuō)明前一例子中設(shè)計(jì)的LCG確實(shí)具有周期。在以下MATLAB程序中，輸入種子數(shù)，并運(yùn)行程序直到種子數(shù)再次出現(xiàn)。設(shè)產(chǎn)生n個(gè)整數(shù)，如果nm而種子數(shù)沒(méi)有再次出現(xiàn)，則認(rèn)為發(fā)生器進(jìn)入了一個(gè)重復(fù)產(chǎn)生

11、短序列的循環(huán)中。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：具有素模數(shù)的乘性算法乘性發(fā)生器的定義式為它是增量等于零的混合算

12、法。注意，c=0時(shí)，xi不能為零，因此，全周期是而不是前面那種情況下的。若滿足以下特性，乘性算法能產(chǎn)生全周期序列；*m是素?cái)?shù)（通常要求m取較大值）*m為的mod(m)本原元素如果除i=m-1了外沒(méi)有更小的i值，使ai-1是m的倍數(shù)，則a為的mod(m)本原元素。換句話說(shuō)，a為mod(m)的本原元素，如果滿足22方法C：具有非素模數(shù)的乘性算法模數(shù)不是素?cái)?shù)的同余算法中最重要的情況是m等于2的冪，即 這里n為整數(shù)。對(duì)式（7-16）定義的算法，最大周期為2n/4=2n-2。如果滿足以下條件，可取得這個(gè)周期：*乘子 結(jié)果為3或5*種子數(shù)x0是奇數(shù)（7-16）23由于兩個(gè)奇數(shù)的乘積是奇數(shù)，可以推出，若x0

13、是奇數(shù)，則由式（7-16）產(chǎn)生的所有值都是奇數(shù)，也就是說(shuō)，產(chǎn)生的值都不是偶數(shù)，這使周期減為原來(lái)的一半。式（7-16）產(chǎn)生的奇數(shù)分成兩組，其中只有一組由給定的種子數(shù)產(chǎn)生，這樣又使周期減小了一倍。產(chǎn)生的這組奇數(shù)一般跟所選的種子數(shù)有關(guān)采用m=2n的好處是整數(shù)溢出能用于運(yùn)算，這樣縮短了計(jì)算時(shí)間，這樣的結(jié)果確實(shí)比較理想，但所得的程序不易移植。24隨機(jī)數(shù)發(fā)生器的測(cè)試 前面一節(jié)為我們提供了工具，來(lái)產(chǎn)生在0與1之間均勻分布的偽隨機(jī)數(shù)。到目前為止，我們只考慮了由LCG產(chǎn)生的序列的周期。于特定的應(yīng)用來(lái)說(shuō)，還要滿足其他理想的特性，至少，我們希望序列要滿足相關(guān)（白噪聲）。具體的應(yīng)用可能還必須滿足其他的要求。現(xiàn)已開(kāi)發(fā)了

14、多個(gè)程序來(lái)測(cè)試某一給定序列的隨機(jī)性（randomness），其中最常用的Chi-方（Chi-square）測(cè)試、Kolomogorov-Smirnov測(cè)試和譜測(cè)試。 25譜測(cè)試似乎是功能最強(qiáng)大的。在后面很多的例子中，要滿足的一個(gè)最重要的特性是給定序列的元素相互獨(dú)立或至少互不相關(guān)。為此，我們來(lái)考慮兩個(gè)非常簡(jiǎn)單的測(cè)試：散點(diǎn)圖（scatterplots）和Durbin-Watson測(cè)試。需要指出的是給定序列的性質(zhì)適用于完整的（全周期）序列。如果只用到序列的一部分，則性質(zhì)不再存在。 26散點(diǎn)圖 所謂散點(diǎn)圖就是x(i+1)與x(i)的函數(shù)關(guān)系圖，它表示了隨機(jī)數(shù)發(fā)生器經(jīng)驗(yàn)式的質(zhì)量指標(biāo)。本例中所考慮的兩個(gè)隨

15、機(jī)數(shù)發(fā)生器定義如下：27可以驗(yàn)證是全周期的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測(cè)試Durbin-Watson獨(dú)立性測(cè)試可以通過(guò)計(jì)算如下Durbin參數(shù)來(lái)完成：其中Xn是一個(gè)零均值（zero-mean）隨機(jī)變量4。下面會(huì)講到當(dāng)D值在2附近時(shí)，Xn和Xn-1相關(guān)性很小。30Durbin-Watson測(cè)試為了說(shuō)明Durbin-Watson測(cè)試的性質(zhì)，假設(shè)

17、和相關(guān)并且是一個(gè)遍歷性過(guò)程。為了簡(jiǎn)化記號(hào)，假設(shè)N足夠大，使得N=N-1，式（7-19）可寫(xiě)成： （7-20）由假定Xn和Xn-1相關(guān)，可令其中X和Z不相關(guān)，是X和Y的相關(guān)系數(shù)。注意，X、Y和Z有相同的方差，用表示。將式（7-21）代入式（7-20）得：31X和Z不相關(guān)且均值為零，所以中間項(xiàng)等于零。由于和方差相同，所以有因?yàn)?，Durbin參數(shù)在0與4之間，且有=0時(shí)D=2。D0表示值為負(fù)。以下的MATLAB函數(shù)計(jì)算Durbin參數(shù)值：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 本例將計(jì)算例7-5中討論到的兩個(gè)噪聲發(fā)生器的值。 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最低標(biāo)準(zhǔn)表明給定LCG能通過(guò)各種隨機(jī)性的統(tǒng)計(jì)測(cè)試，以便徹底測(cè)試它的質(zhì)量，這是一項(xiàng)首要任務(wù)。對(duì)于產(chǎn)生長(zhǎng)的序列來(lái)說(shuō)更是如此。為了部分地解決這個(gè)問(wèn)題，多種算法被確定為最低標(biāo)準(zhǔn)（Minimum standard）算法。最低標(biāo)準(zhǔn)算法具有以下性質(zhì)：*全周

21、期*能通過(guò)所有適用的隨機(jī)性統(tǒng)計(jì)測(cè)試*易于從一臺(tái)計(jì)算機(jī)移植到另一臺(tái)計(jì)算機(jī)如果使用一個(gè)最低標(biāo)準(zhǔn)算法，不用擔(dān)心算法本身的正確性，但必須保證在給定計(jì)算環(huán)境下算法得到正確的實(shí)現(xiàn)。編程時(shí)要特別注意的是，由算法產(chǎn)生的所有的數(shù)都要能唯一地表示出來(lái)。36Lewis, Goodman和Miller最低標(biāo)準(zhǔn)其中m是Mersenne素?cái)?shù)231-1，這個(gè)值最初由Lehmer推薦使用，很容易在32位的計(jì)算機(jī)以上以整型運(yùn)算實(shí)現(xiàn)，如果尾數(shù)超過(guò)31位，則以浮點(diǎn)運(yùn)算實(shí)現(xiàn)。37Wichmann-Hill算法希望獲得一些具有長(zhǎng)周期的隨機(jī)數(shù)發(fā)生器。要構(gòu)造具有長(zhǎng)周期的波型，一個(gè)行之有效的辦法是將幾個(gè)周期相差很小的周期波型相加。Cos2(

22、1)t的周期是1s， Cos2(1.0001)t的周期是1000/10001s，略小于1s 。兩者的復(fù)合波形可寫(xiě)成如下形式：其周期為10000s或約合2.78個(gè)小時(shí)。復(fù)合波形的一個(gè)周期內(nèi)，第一個(gè)分量經(jīng)過(guò)10000個(gè)周期，第二個(gè)分量經(jīng)過(guò)10001個(gè)周期。還可再添加更多的分量。38可以將同樣的方法用于LCG，即把幾個(gè)周期相差不大的隨機(jī)數(shù)發(fā)生器合成在一起。Wichmann-Hill算法可能是最有名的一個(gè)合成隨機(jī)數(shù)發(fā)生器的例子。該算法用到的三個(gè)分量發(fā)生器，定義如下： （7-26） （7-28）39這三個(gè)分量發(fā)生器都是全周期發(fā)生器，它們合成得到輸出為40Wichmann-Hill算法等于價(jià)于乘性LCG，它的倍數(shù)為模數(shù)為顯然，m不是素?cái)?shù)，周期小于m-