ch7 隨機(jī)信號的產(chǎn)生與處理

第七講隨機(jī)信號的產(chǎn)生與處理1-平穩(wěn)與遍歷性過程2-均勻隨機(jī)數(shù)發(fā)生器3-將均勻分布隨機(jī)變量映射成任意pdf4-產(chǎn)生不相關(guān)的高斯隨機(jī)變量5-產(chǎn)生相關(guān)的高斯隨機(jī)變量6-同時(shí)確定pdf和PSD7-PN序列發(fā)生器8-信號處理1引言以前所討論的是仿真中的確定性信號。在所有具有實(shí)際意義的通信系統(tǒng)中，在信息承載信號從信源端到終端用戶的傳送過程中，信道噪聲、干擾和衰落等隨機(jī)影響會使它造成損失。要想在波形級精確地仿真這些系統(tǒng)，首先要對這些隨機(jī)影響建立準(zhǔn)確的模型，因此，需要產(chǎn)生這些隨機(jī)影響的算法，其基本構(gòu)建模塊是隨機(jī)數(shù)發(fā)生器。雖然隨機(jī)數(shù)發(fā)生器方面的內(nèi)容很豐富，但本章強(qiáng)調(diào)的是隨機(jī)數(shù)發(fā)生器在通信系統(tǒng)仿真中的運(yùn)用。因此，我們僅限于討論如何產(chǎn)生采樣后的隨機(jī)波形（信號、干擾和噪聲等），以用于仿真。2在仿真環(huán)境下，所有的隨機(jī)過程必須用隨機(jī)變量序列來表示，本章主要內(nèi)容便是產(chǎn)生和測試這些隨機(jī)序列，許多適用于仿真程序開發(fā)的程序設(shè)計(jì)語言，如MATLAB，都將隨機(jī)數(shù)發(fā)生器含于其“內(nèi)置”函數(shù)庫中。了解這些隨機(jī)數(shù)發(fā)生器的原理，有助于深入理解整個(gè)仿真程序。確保隨機(jī)數(shù)發(fā)生器設(shè)計(jì)合理并適用于給定的應(yīng)用，方不失為明智之舉。我們將看到，嚴(yán)格來講，隨機(jī)數(shù)發(fā)生器所產(chǎn)生的并不是真正的隨機(jī)序列，而是在觀察（仿真）區(qū)間上“呈現(xiàn)隨機(jī)性”的序列，以此來近似在給定仿真程序中的隨機(jī)過程的樣本函數(shù)。所謂“呈現(xiàn)隨機(jī)性”指的是產(chǎn)生的序列在仿真區(qū)間上具有某些特性，能對具體應(yīng)用中的隨機(jī)過程建立準(zhǔn)確模型并達(dá)到要求的精度，我們稱這種序列為“偽隨機(jī)序列”。這樣命名是因?yàn)楸M管它們是確定的，但在具體應(yīng)用中會呈現(xiàn)隨機(jī)性。3所要求的模型精度取決于具體的應(yīng)用。例如，如果我們要產(chǎn)生一個(gè)波形來表示鎖相環(huán)鑒相器輸入端的噪聲，對一個(gè)輸入端SNR為50dB的系統(tǒng)建立噪聲波形模型所要求的精度要比輸入端SNR為8dB的系統(tǒng)高得多。要建立數(shù)字通信系統(tǒng)中噪聲分量的模型，誤比特率為10-7時(shí)所需的精度會比誤比特率為10-3時(shí)所需的精度要高。

本章首先介紹隨機(jī)過程樣本函數(shù)的產(chǎn)生，在仿真環(huán)境下考察平穩(wěn)的概念，然后簡要介紹一下數(shù)字調(diào)制器的仿真模型。有了這些預(yù)備知識后，我們再轉(zhuǎn)向本章的重點(diǎn)，來討論如下問題：4※在（0，1)上產(chǎn)生均勻分布且不相關(guān)的隨機(jī)數(shù)※將不相關(guān)且均勻分布的隨機(jī)數(shù)映射成不相關(guān)的、具有任意（需要的）概率密度函數(shù)（pdf）的隨機(jī)數(shù)※產(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ì)算方法。5一、平穩(wěn)與遍歷性過程

在仿真通信系統(tǒng)時(shí)，用于表示信號、噪聲和干擾而產(chǎn)生的樣本函數(shù)通常假設(shè)為各態(tài)歷經(jīng)的。作這一假設(shè)，是因?yàn)槲覀兺ǔＲ来翁幚硐到y(tǒng)中波形的時(shí)域樣本，而且系統(tǒng)中的每個(gè)點(diǎn)上只有一個(gè)波形（樣本函數(shù)）。我們還假設(shè)仿真所處理的波形是由其內(nèi)在的統(tǒng)計(jì)模型定義的總體(Ensemble）中的一個(gè)典型成員，這時(shí)各種統(tǒng)計(jì)量如各階矩、信噪比和誤比特率，就可以當(dāng)作時(shí)間平均來計(jì)算。將仿真結(jié)果與對應(yīng)的理論結(jié)果進(jìn)行比較時(shí)，常有一個(gè)潛在的假設(shè)：仿真計(jì)算得到的時(shí)間平均等于總體均值。因此，隱含的假設(shè)條件是對應(yīng)的隨機(jī)過程為遍歷性的。6遍歷性過程總是平穩(wěn)的，所以，總是可以假設(shè)仿真中產(chǎn)生的樣本函數(shù)屬于平穩(wěn)隨機(jī)過程。由基本的隨機(jī)過程理論可知，平穩(wěn)性是這樣定義的：所有統(tǒng)計(jì)量與時(shí)間原點(diǎn)無關(guān)。為了說明其中的幾點(diǎn)思想，先來看一個(gè)簡單的例子。例7-1假設(shè)隨機(jī)過程的樣本函數(shù)表達(dá)式為式中，i是對應(yīng)隨機(jī)試驗(yàn)的樣本空間中的一個(gè)輸出，每一個(gè)i映射為一個(gè)相位i。再假定對應(yīng)的隨機(jī)試驗(yàn)就是從均勻隨機(jī)數(shù)發(fā)生器輸出端抽取一個(gè)數(shù)，抽取得到的結(jié)果i=ui，這里ui在(0,1)區(qū)間上均勻分布。然后ui映射成相位i＝kui。當(dāng)幅度A和頻率f固定，i的值便決定了波形。7這個(gè)例子中我們感興趣的是k的兩個(gè)值：k=2，對應(yīng)的相位均勻分布在（0，2）區(qū)間，和k＝/2，對應(yīng)的相位均勻分布在（0，/2）區(qū)間。作為第二個(gè)例子，假設(shè)隨機(jī)過程由以下表達(dá)式描述：在這種情況下，幅度在（A，2A）上均勻分布。下面MATLAB程序產(chǎn)生一個(gè)隨機(jī)過程的三組樣本函數(shù)。1）波形表示為x(t)，對應(yīng)式（7-1)中k＝2的情況。2）波形表示為y(t)，對應(yīng)式（7-1)中k＝/2的情況。3）波形表示為z(t)，由式（7-2）所定義。對于所有的波形有A＝1，f＝1。每個(gè)仿真過程產(chǎn)生兩秒鐘長度的數(shù)據(jù)以及20個(gè)樣本函數(shù)。8%File:c7_sinewave.mf=1; %frequencyofsinusoidfs=100; %samplingfrequencyt=(0:200)/fs; %timevectorfori=1:20x(:,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)';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)')%Endofscriptfile.9

由x(t)、y(t)和z(t)組成的所有樣本函數(shù)的時(shí)間平均都等于零。很容易便可證明：對于大量滿足0≤ti≤2的ti,計(jì)算出來的x(t)的總體均值近似為零；若樣本函數(shù)的個(gè)數(shù)趨近∞，時(shí)間平均收斂于零。然而，對于y(t)，當(dāng)t在0.875和1.875附近時(shí)總體均值近似為1，當(dāng)t在0.375和1.375附近時(shí)總體均值近似為-1，而當(dāng)t在0.125、0.625、1.125和1.625附近時(shí)總體均值近似為零。這是周期平穩(wěn)過程的一個(gè)例子，這種過程的矩是周期的。

z(t)表示的樣本函數(shù)也是周期平穩(wěn)過程的樣本函數(shù)。需要注意的是，若以t＝0.5k對該過程進(jìn)行采樣，得到的隨機(jī)變量在k為偶數(shù)時(shí)其均值近似為+1.5，而在k為奇數(shù)時(shí)其均．值近似為-1.5。■10

前面的例子使用了MATLAB中均勻隨機(jī)數(shù)發(fā)生器rand,接下來的例子將說明如何用隨機(jī)數(shù)發(fā)生器來建立數(shù)字調(diào)制器的模型。下一節(jié)是實(shí)現(xiàn)均勻隨機(jī)數(shù)發(fā)生器的算法。圖7-1三個(gè)不同隨機(jī)過程的樣本函數(shù)20個(gè)實(shí)現(xiàn)11圖7-1-三個(gè)不同隨機(jī)過程的樣本函數(shù)500個(gè)實(shí)現(xiàn)12

例7-2在后面的工作中，我們經(jīng)常要用到數(shù)字調(diào)制器的模型。這些調(diào)制器基本構(gòu)建模塊是函數(shù)randombinary，它產(chǎn)生電平值為+1或-1的二進(jìn)制波形，產(chǎn)生的比特?cái)?shù)以及每比特的采樣數(shù)是該函數(shù)的參數(shù)。程序如下：%File:random_binary.mfunction[x,bits]=random_binary(nbits,nsamples)%Thisfunctiongenratesarandombinarywaveformoflengthnbits%sampledatarateofnsamples/bit.x=zeros(1,nbits*nsamples);bits=round(rand(1,nbits));form=1:nbitsforn=1:nsamplesindex=(m-1)*nsamples+n;x(1,index)=(-1)^bits(m);endend%Endoffunctionfile.13函數(shù)randombinary可以仿真多個(gè)數(shù)字調(diào)制器，例如，可用如下MATLAB語句仿真一個(gè)QPSK調(diào)制器：x=randombinary(nbits,nsamples)+i*randombinary(nbits,nsamples);下列MATLAB程序產(chǎn)生一個(gè)長10比特的QPSK信號，采樣頻率為每比特8個(gè)采樣點(diǎn)：%File:c7_example2.mnbits=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([080-1.51.5]);xlabel('SampleIndex');ylabel('xd')subplot(2,1,2)stem(xq,'.');grid;axis([080-1.51.5]);xlabel('SampleIndex');ylabel('xq')%Endofscriptfile.14運(yùn)行該程序，得到QPSK信號的同相分量和正交分量，如圖7-2所示。圖7-2QPSK信號的同相分量和正交分量（10比特）15圖7-2-QPSK信號的同相分量和正交分量（100比特）16具有均勻概率密度函數(shù)的隨機(jī)變量很容易轉(zhuǎn)換成具有其他所需pdf的隨杭變量，因此，要產(chǎn)生一個(gè)具有特定pdf的隨機(jī)變量，首先需要產(chǎn)生一個(gè)在（0，1)區(qū)間均勻分布的隨機(jī)變量。通常，先產(chǎn)生一列介于0和M之間的數(shù)（整數(shù)），然后將序列中每個(gè)元素除以M。實(shí)現(xiàn)隨機(jī)數(shù)發(fā)生器最常用的方法是線性同余（linearcongruence).（一）線性同余

線性同余發(fā)生器（linearcongruencegenerator，LCG）定義為如下運(yùn)算：二、均勻隨機(jī)數(shù)發(fā)生器17其中，a和c分別稱作乘子（multiplier）和增量，m叫做模數(shù)。當(dāng)然，這是一個(gè)確定性的序列算法，能依次產(chǎn)生連續(xù)的x值。x的初始值記為x0，稱作線性同余發(fā)生器的種子數(shù)（seednumber）。如果x0、a、c和m都是整數(shù)，則LCG產(chǎn)生的所有數(shù)也都是整數(shù)。由于對[axi＋c］進(jìn)行mod(m)運(yùn)算，式（7-3）至多可產(chǎn)生m個(gè)不同的整數(shù)。發(fā)生器輸出的一個(gè)理想特性是它應(yīng)具備很長的周期，從而在序列重復(fù)前，輸出序列能產(chǎn)生最多數(shù)目的整數(shù)。對于給定m值，當(dāng)周期最大時(shí)，我們稱發(fā)生器是全周期（fullperiod）的。此外，具體仿真程序的應(yīng)用對LCG會提出其他的要求，比如，我們通常要求樣本xi和xi+1，互不相關(guān)。另外，根據(jù)具體的應(yīng)用，可能還要求LCG的輸出能通過其他統(tǒng)計(jì)測試。LCG可以采用多種不同的形式，本節(jié)只考慮最常用的算法。18方法A：混合同余算法最通用的同余算法就是c≠0的“混合”同余算法。之所以稱之為混合算法，是因?yàn)樵谇蠼鈞i+1的過程中要同時(shí)用到乘法與加法?；旌暇€性算法具有式（7-4）形式c≠0時(shí)，發(fā)生器的最大周期為m。當(dāng)且僅當(dāng)滿足以下特性時(shí)才能達(dá)到這個(gè)周期：

·增量c與m互質(zhì)。換句話說，c與m沒有素公因子（primefactor）。

·a-1是p的倍數(shù)，這里p表示模數(shù)m的任意一個(gè)素因子。

·如果m是4的倍數(shù)，則a-1是4的倍數(shù)。以上命題的證明見Knuth[1]。19

例7-4前一例子中設(shè)計(jì)的LCG確實(shí)具有周期m＝50000在以下MATLAB程序中，輸入種子數(shù)，并運(yùn)行程序直到種子數(shù)再次出現(xiàn)。設(shè)產(chǎn)生n個(gè)整數(shù)，如果n＞m而種子數(shù)沒有再次出現(xiàn)，則認(rèn)為發(fā)生器進(jìn)入了一個(gè)重復(fù)產(chǎn)生短序列的循環(huán)中。MATLAB程序如下：%Filec7_LCGperiod.ma=input('Entermultipliera>');c=input('Enteroffsetc>');m=input('Entermodulusm>');seed=input('Enterseed>');n=1;ix=rem((seed*a+c),m);while(ix~=seed)&(n<m+2)n=n+1;ix=rem((ix*a+c),m);endifn>mdisp('Caughtinaloop.')elsetext=['Theperiodis',num2str(n,15),'.'];disp(text)end%Endofscriptfile.執(zhí)行程序得以下對話：>>c7l_LCGperiodEntermultipliera>241Enteroffsetc>1323Entermodulusm>5000Enterseed>1Theperiodis5000.可以看到，跟期望的一樣，周期確實(shí)是5000?！?0方法B：具有素模數(shù)的乘性算法·乘性發(fā)生器的定義式為

xi+1＝[axi]mod(m)（7-13）它是增量c等于零的混合算法。注意，c＝0時(shí)，xi不能為零，因此，全周期是m-1而不是前面那種情況下的m。若滿足以下特性，乘性算法能產(chǎn)生全周期序列：

．m是素?cái)?shù)（通常要求m取較大值）

．m為mod(m）的本原元素我們知道，素?cái)?shù)只能被1和它本身整除，這里可能需要對第二條特性作一下說明。如果除了i＝m-1外沒有更小的i值，使ai-1是m的倍數(shù)，則a為mod(m)的本原元素。換句話說，a為mod(m)的本原元素，如果滿足21其中k是一個(gè)任意整數(shù)。

方法C：具有非素模數(shù)的乘性算法模數(shù)m不是素?cái)?shù)的同余算法中最重要的情況是m等于2的冪，即：

xi+l＝＝[axi]mod(2n)（7-16)（這里n為整數(shù)。對式（7-16）定義的算法，最大周期為2n/4＝2n-2。如果滿足以下條件，可取得這個(gè)周期：

．乘子amod（8）結(jié)果為3或5

．種子數(shù)x0是奇數(shù)22由于兩個(gè)奇數(shù)的乘積是奇數(shù)，可以推出,若x0是奇數(shù)，則由式（7-16）產(chǎn)生的所有值都是奇數(shù)，也就是說，產(chǎn)生的xi值都不是偶數(shù)，這使周期減為原來的一半。式（7-16）產(chǎn)生的奇數(shù)分成兩組，其中只有一組由給定的種子數(shù)產(chǎn)生，這樣又使周期減小了一倍。產(chǎn)生的這組奇數(shù)一般跟所選的種子數(shù)有關(guān)。采用m＝2k的好處是整數(shù)溢出能用于mod(m)運(yùn)算，這樣縮短了計(jì)算時(shí)間。這樣的結(jié)果確實(shí)比較理想，但所得程序不易移植。23(二)隨機(jī)數(shù)發(fā)生器的測試前面一節(jié)為我們提供了工具，來產(chǎn)生在0與1之間均勻分布的偽隨機(jī)數(shù)。到目前為止，我們只考慮了由LCG產(chǎn)生的序列的周期。雖然我們很希望序列的周期越長越好，但對于特定的應(yīng)用來說，還要滿足其他理想的特性，至少，我們希望序列要滿足相關(guān)（白噪聲）。具體的應(yīng)用可能還必須滿足其他的要求。現(xiàn)已開發(fā)了多個(gè)程序來測試某一給定序列的隨機(jī)性（randomness），其中最常用的Chi-方(Chi-square）測試、Kolomogorov-Smirnov測試和譜測試。在這些測試中，譜測試似乎是功能最強(qiáng)大的。在后面很多的例子中，要滿足的一個(gè)最重要的特性是給定序列的元素相互獨(dú)立或至少互不相關(guān)。24為此，我們來考慮兩個(gè)非常簡單的測試：散點(diǎn)圖（scatterplots）和Durbin-Watson測試。需要指出的是給定序列的性質(zhì)適用于完整的（全周期）序列。如果只用到序列的一部分，則性質(zhì)不再存在。散點(diǎn)圖下面通過一個(gè)例子來說明散點(diǎn)圖。例7-5所謂散點(diǎn)圖就是xi+1與xi的函數(shù)關(guān)系圖，它表示了隨機(jī)數(shù)發(fā)生器經(jīng)驗(yàn)式的質(zhì)量指標(biāo)。本例中所考慮的兩個(gè)隨機(jī)數(shù)發(fā)生器定義如下：

xi十1＝(65xi＋1)mod(2048)（7-17）和

xi＋1＝(1229xi＋1)mod(2048）（7-18）用例7-4中的程序c7_LCGperiod.m驗(yàn)證，這兩個(gè)發(fā)生器都是全周期的。每個(gè)發(fā)生器產(chǎn)生散點(diǎn)圖的MATLAB代碼如下：25%File:c7_LCDemo1.mm=2048;c=1;seed=1; %defaultvaluesofmandca1=65;a2=1229; %multipliervaluesix1=seed;ix2=seed; %initializealgorithmx1=zeros(1,m);x2=zeros(1,m);%initializearraysfori=1:mix1=rem((ix1*a1+c),m);x1(i)=ix1/m;ix2=rem((ix2*a2+c),m);x2(i)=ix2/m;endsubplot(1,2,1)y1=[x1(1,2:m),x1(1,1)];plot(x1,y1,'.') %plotresultsfora1subplot(1,2,2)y2=[x2(1,2:m),x2(1,1)];plot(x2,y2,'.') %plotresultsfora2%Endofscriptfile.執(zhí)行程序產(chǎn)生的散點(diǎn)圖見圖7-3。人們希望得到的散點(diǎn)圖中縱坐標(biāo)xi+1和橫坐標(biāo)xi的各種組合都26會出現(xiàn)，這種情形下的散點(diǎn)圖缺乏結(jié)構(gòu)性，從圖7-3看來，a＝65的發(fā)生器比a＝1229的發(fā)生器具有更小的序列相關(guān)性。在例7-6中我們將看到情況確實(shí)如此?！鰣D7-3a1=65（左圖）和a2=1229（右圖）時(shí)的散點(diǎn)圖27Durbin-Watson測試Durbin-Watson獨(dú)立性測試可以通過計(jì)算Durbin參數(shù)來完成：其中X[n]是一個(gè)零均值(zero-mean)隨機(jī)變量。下面會講到當(dāng)D值在2附近時(shí)，X[n]和X[n-1］相關(guān)性很小。為了說明Durbin-Watson測試的性質(zhì)，假設(shè)X[n]和X[n-1]相關(guān)并且X[n]是一個(gè)遍歷性過程。為了簡化記號，假設(shè)N足夠大，使得N-1≈N,式（7-19)可寫成：28式中X表示X[n]，Y表示X[n-1]，E{·｝表示數(shù)學(xué)期望。由于假定X[n]和X[n-1]相關(guān)，可令其中X和Z不相關(guān)，是X和Y的相關(guān)系數(shù)。注意，X、Y和Z有相同的方差，用2表示。將式（7-21）代入式（7-20）得：X和Z不相關(guān)且均值為零，所以中間項(xiàng)等于零。由于X和Z方差相同，所以有29因?yàn)?1≤≤1，Durbin參數(shù)D在O與4之間，且有＝O時(shí)D＝2。D＜2表示正相關(guān)，而D＞2表示值為負(fù)。以下的MATLAB函數(shù)計(jì)算Durbin參數(shù)值：%File:c7_DurbinfunctionD=durbin(x)N=length(x); %lengthofinputvectory=x-mean(x); %removedcydiff=y(2:N)-y(1:(N-1)); %numeratorsummandNum=sum(ydiff.*ydiff); %numeratorfactorofDDen=sum(y.*y); %denominatorfactorofDD=Num/Den; %Durbinfactor%Endoffunctionfile.30例7-6本例將計(jì)算例7-5中討論到的兩個(gè)噪聲發(fā)生器的D值。MATLAB代碼如下：%File:c7_LCDemo2m=2048;c=1;seed=1;a1=65;a2=1229;ix1=1;ix2=1;x1=zeros(1,m);x2=zeros(1,m);fori=1:mix1=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);%calculateDurbinparametersrho1=1-D1/2;rho2=1-D2/2; %calculatecorrelationtext1=['ThevalueofD1is',num2str(D1),'andrho1is',num2str(rho1),'.'];text2=['ThevalueofD2is',num2str(D2),'andrho2is',num2str(rho2),'.'];disp(text1)disp(text2)%Endofscriptfile.31執(zhí)行程序得到：>>c71.CDemo2ThevalueofD1is1.9925andrholis0.0037273.ThevalueofD2is1.6037andrho2is0.19814.a1＝65對應(yīng)的相關(guān)值約為0;a2＝1229對應(yīng)的相關(guān)值約為0.2。因此，從Durbin測試結(jié)果得出a1=65給出的結(jié)果優(yōu)于a2＝1229，這個(gè)結(jié)論與圖7-3所示的散點(diǎn)圖一致?！觯ㄈ┳畹蜆?biāo)準(zhǔn)表明給定LCG能通過各種隨機(jī)性的統(tǒng)計(jì)測試，以便徹底測試它的質(zhì)量，這是一項(xiàng)首要任務(wù)。對于產(chǎn)生長的序列來說更是如此。為了部分地解決這個(gè)問題，多種算法被確定為最低標(biāo)準(zhǔn)（minimumstandard）算法。最低標(biāo)準(zhǔn)算法具有以下性質(zhì)：·全周期/·能通過所有適用的隨機(jī)性統(tǒng)計(jì)測試、·易于從一臺計(jì)算機(jī)移植到另一臺計(jì)算機(jī)32一旦算法經(jīng)確定并適當(dāng)?shù)赜涗浵聛?，它就成了一個(gè)最低標(biāo)準(zhǔn)算法，于是這個(gè)算法可以在別處放心地使用而無須額外的測試。如果使用一個(gè)最低標(biāo)準(zhǔn)算法，不用擔(dān)心算法本身的正確性，但必須保證在給定計(jì)算環(huán)境下算法得到正確的實(shí)現(xiàn)編程時(shí)要特別注意的是，由算法產(chǎn)生的所有的數(shù)都要能唯一地表示出來。33Lewis，Goodman和Miller最低標(biāo)準(zhǔn)

Lewis，Goodman和Miller最低標(biāo)準(zhǔn)定義為：

xi+1＝(16807xi)mod(2147483647)（7-24）其中m是Mersenne素?cái)?shù)231-1，這個(gè)值最初由Lehmer推薦使用，他在半個(gè)多世紀(jì)前完成了有關(guān)LCG的大量基礎(chǔ)研究工作。這種算法應(yīng)用廣泛，很容易在32位的計(jì)算機(jī)上以整型運(yùn)算實(shí)現(xiàn)，如果尾數(shù)超過31位，則以浮點(diǎn)運(yùn)算實(shí)現(xiàn)畫。

Wichmann-Hill算法前面已說明我們希望獲得一些具有長周期的隨機(jī)數(shù)發(fā)生器。要構(gòu)造具有長周期的波形，一個(gè)行之有效的辦法是將幾個(gè)周期相差很小的周期波形相加。例如，cos2(1)t的周期是1s，cos2(1.0001)t的周期是10000/10OOls，略小于1s。兩者的復(fù)合波形如下：34其周期為10OOOs或約合2.78個(gè)小時(shí)。復(fù)合波形的一個(gè)周期內(nèi)，第一個(gè)分量經(jīng)過10000個(gè)周期，第二個(gè)分量經(jīng)過10001個(gè)周期。如果有需要，還可再添加更多分量?？梢詫⑼瑯拥姆椒ㄓ糜贚CG，即把幾個(gè)周期相差不大的隨機(jī)數(shù)發(fā)生器合成在一起，Wichmann-Hill算法可能是最有名的一個(gè)合成隨機(jī)數(shù)發(fā)生器的例子。Wichmann-Hill算法可以有多種變形。該算法用到的三個(gè)分量發(fā)生器，定義如下：35這三個(gè)分量發(fā)生器都是全周期發(fā)生器。它們合成得到的輸出為Wichmann-Hill算法等價(jià)于乘性LCG，它的倍數(shù)為模數(shù)為顯然，m不是素?cái)?shù)，周期小于m-1。[3]中有講到周期近似為7.0×1012,盡管小于m，但仍然很長。36

Wichmann-Hill算法雖然在結(jié)構(gòu)上與前面講到的最低標(biāo)準(zhǔn)有些不同，但由于它已經(jīng)得到廣泛的測試，被證明通過了所有標(biāo)準(zhǔn)的統(tǒng)計(jì)測試，并且易于從一臺機(jī)子上移植到另一臺機(jī)子上，所以仍然認(rèn)為它是最低標(biāo)準(zhǔn)均勻隨機(jī)數(shù)發(fā)生器。（四）MATLAB實(shí)現(xiàn)

MATLAB版本5以前，函數(shù)庫內(nèi)的均勻隨機(jī)數(shù)發(fā)生器rand是由式（7-24）定義的最低標(biāo)準(zhǔn)隨機(jī)數(shù)發(fā)生器。MATLAB版本5和版本6中的rand則是基于由Marsaglia開發(fā)的一種方法。這個(gè)隨機(jī)數(shù)發(fā)生器要產(chǎn)生的是浮點(diǎn)數(shù)而不是縮放的整數(shù)。MathWorks聲稱這個(gè)隨機(jī)數(shù)發(fā)生器的周期超過21492，并且“相當(dāng)確信”可以產(chǎn)生eps與1-eps／2之間的所有浮點(diǎn)數(shù)，其中MATLAB常數(shù)eps等于2-52。新的隨機(jī)數(shù)發(fā)生器只用到加法和減法運(yùn)算，沒有用到乘法或除法運(yùn)算，因此算法執(zhí)行起來比LCG快得多。37（五）種子數(shù)與種子向量本書的仿真實(shí)例都是基于MATLAB，因此有必要簡單了解一下MATLAB運(yùn)用種子數(shù)的方式?！芭f版本”MATLAB隨機(jī)數(shù)發(fā)生器（MATLAB版本5以前由式（7-24）定義的）使用單個(gè)種子數(shù)，“新版本”隨機(jī)數(shù)發(fā)生器使用種子向量，稱為隨機(jī)數(shù)發(fā)生器的狀態(tài)。該向量由35個(gè)表示隨機(jī)數(shù)發(fā)生器狀態(tài)的元素組成，它們分別是32個(gè)浮點(diǎn)數(shù)、兩個(gè)整數(shù)和一個(gè)標(biāo)志位。到了版本5以及后來的新版本，上面兩種隨機(jī)數(shù)發(fā)生器都可以使用，其中新的隨機(jī)數(shù)發(fā)生器設(shè)定為默認(rèn)的發(fā)生器。由式(7-24)定義的舊版的隨機(jī)數(shù)發(fā)生器可通過命令RAND（’seed’，0)或RAND（’seed’，J）來調(diào)用。跟其他所有MATLAB命令一樣，使用者應(yīng)該仔細(xì)研究由help命令提供的信息。此外，使用者還須注意以下幾點(diǎn)（種子數(shù)這一術(shù)語可用來表示整數(shù)種子和狀態(tài)向量）：38·使用者可以使用默認(rèn)的種子數(shù)，也可以自己給定一個(gè)種子數(shù)。·關(guān)閉并重新打開MATLAB可以將種子數(shù)復(fù)位到默認(rèn)值。因此，如果調(diào)用一個(gè)隨機(jī)數(shù)發(fā)生器N次后，關(guān)閉MATLAB，接著重新打開，再調(diào)用隨機(jī)數(shù)發(fā)生器N次，那么兩種情況下產(chǎn)生的N個(gè)數(shù)是相同的。這一特性可以用來發(fā)揮MATLAB的優(yōu)勢，它允許我們重新產(chǎn)生結(jié)果完全相同的序列。這對測試很有用。·系統(tǒng)時(shí)鐘可用來隨機(jī)化初始的種子數(shù)（具體細(xì)節(jié)見MATLABhelp).·種子數(shù)存在一個(gè)緩沖器中，而不是MATLAB工作區(qū)。因此，命令clearall對種子數(shù)不起作用。39

RAND('state',0)resetsthegeneratortoitsinitialstate.RAND('state',J),forintegerJ,resetsthegeneratortoitsJ-thstate.RAND('state',sum(100*clock))resetsittoadifferentstateeachtime.==============================================

RANDN('state',0)resetsthegeneratortoitsinitialstate.RANDN('state',J),forintegerJ,resetsthegeneratortoitsJ-thstate.RANDN('state',sum(100*clock))resetsittoadifferentstateeachtime.40有多種方法可將一個(gè)均勻分布的隨機(jī)變量映射成一個(gè)具有非均勻pdf的目標(biāo)隨機(jī)變量。基本有三種情況：

1.已知目標(biāo)隨機(jī)變量的概率分布函數(shù)（CDF）具有閉合形式。可以看到，如果目標(biāo)隨機(jī)變量的CDF具有閉合形式，就可以采用一種叫做逆變換法的簡單方法。

2.已知目標(biāo)隨機(jī)變量的pdf具有閉合形式，但CDF不具有閉合形式。常用的高斯隨機(jī)變量就屬于這種類型。對于這種情況，可以用許多專門方法來解決，此外，還可以用舍棄法。

3.已知pdf和CDF都不具有閉合形式。在尋求一個(gè)隨機(jī)數(shù)發(fā)生器，使其pdf與實(shí)驗(yàn)數(shù)據(jù)的pdf相一致時(shí)，經(jīng)常會出現(xiàn)這種情況?，F(xiàn)在來討論可用于以上三種情況的各種方法。三、將均勻分布隨機(jī)變量映射成任意pdf41逆變換方法可以將一個(gè)不相關(guān)的均勻分布的隨機(jī)序列U變換為一個(gè)具有概率分布函數(shù)Fx(x)的不相關(guān)的（獨(dú)立的樣本）序列X。如圖7-4所示，這個(gè)變換使用了一個(gè)無記憶的非線性器件。由于非線性器件是無記憶的，在輸入序列不相關(guān)時(shí)，能保證輸出序列也是不相關(guān)的。當(dāng)然，根據(jù)維納--辛欽定理，不相關(guān)隨機(jī)數(shù)序列的功率譜密度是常數(shù)（白噪聲）。逆變換方法簡單地設(shè)定：（一）逆變換法并求解X，得42使用逆變換方法時(shí)要求已知分布函數(shù)Fx(x)為閉合形式。容易看出，逆變換方法能產(chǎn)生具有所需分布函數(shù)的隨機(jī)變量。我們知道，分布函數(shù)Fx(x）是自變量x的一個(gè)非減函數(shù)，如圖7-5所示。根據(jù)分布函數(shù)的定義43令X＝F-1（U），考慮到FX(x)單調(diào)，所以有即得要求的結(jié)果?，F(xiàn)在先通過一個(gè)例子來說明逆變換方法。44例7-7本例中，將一個(gè)均勻分布的隨機(jī)變量變換成具有單邊指數(shù)分布的隨機(jī)變量其中u(x)是單位階躍函數(shù)，定義如下我們首先是要求出CDF，也就是令分布函數(shù)與均勻隨機(jī)變量U相等，有45求解X，上式可化為由于隨機(jī)變量1-U與U等價(jià)(Z＝U與Z＝1-U具有相同的pdf)，所以X可寫成實(shí)現(xiàn)均勻分布一指數(shù)分布變換的MATLAB代碼如下：46%File:c7_uni2expclearall %besafen=input('Enternumberofpoints>');b=3; %setpdfparameteru=rand(1,n); %generateUy_exp=-log(u)/b; %transformation[N_samp,x]=hist(y_exp,20); %gethistogramparameterssubplot(2,1,1)bar(x,N_samp,1) %plothistogramylabel('NumberofSamples')xlabel('IndependentVariable-x')subplot(2,1,2)y=b*exp(-3*x); %calculatepdfdel_x=x(3)-x(2); %determinebinwidthp_hist=N_samp/n/del_x; %probabilityfromhistogramplot(x,y,'k',x,p_hist,'ok') %compareylabel('ProbabilityDensity')xlabel('IndependentVariable-x')legend('truepdf','samplesfromhistogram',1)%Endofscriptfile.47圖7-6均勻分布轉(zhuǎn)換為指數(shù)分布（N=100）48圖7-7均勻分布轉(zhuǎn)換為指數(shù)分布（N=2000）49＝3，N＝100對應(yīng)的結(jié)果如圖7-6所示。圖7-6中上面部分是直方圖，第二部分同時(shí)給出理論上的pdf和使用1000個(gè)樣本所得的試驗(yàn)值。N=100時(shí)所得的效果相對比較差，這促使我們?nèi)「嗟臉颖军c(diǎn)重新進(jìn)行試驗(yàn)。N＝2000對應(yīng)的結(jié)果如圖7-7所示，可以看到明顯的改善?！隼?-8作為第二個(gè)例子，考慮瑞利隨機(jī)變量，其pdf表示為跟前面一樣，單位階躍函數(shù)確定了pdf是單邊的。瑞利隨機(jī)變量的CDF為50令FR（R）＝U,則有上式等價(jià)于這里也考慮到了1-U與U等價(jià)。求解R得到這個(gè)變換是Box-Muller算法中的第一步。而Box-Muller算法是產(chǎn)生高斯隨機(jī)變量的一個(gè)基本算法。51從前面的例子看出，我們感興趣的是實(shí)現(xiàn)這個(gè)變換，并且根據(jù)變換的點(diǎn)數(shù)評估其性能。MATLAB程序如下，N＝10000時(shí)的執(zhí)行結(jié)果如圖7-8所示。可以選用其他N值，并將執(zhí)行結(jié)果與圖7-8進(jìn)行比較。■以上兩個(gè)例子闡明了如何將逆變換方法應(yīng)用于連續(xù)隨機(jī)變量，然而，逆變換方法同樣也適用于離散隨機(jī)變量。接下來要介紹的直方圖法是逆變換法的數(shù)字（離散數(shù)據(jù)）版本。52圖7-8均勻分布轉(zhuǎn)換為瑞利分布（N=10000）53%File:c7_uni2rayclearall %besafen=input('Enternumberofpoints>');varR=3; %setpdfparameteru=rand(1,n); %generateUy_exp=sqrt(-2*varR*log(u)); %transformation[N_samp,r]=hist(y_exp,20); %gethistogramparameterssubplot(2,1,1)bar(r,N_samp,1) %plothistogramylabel('NumberofSamples')xlabel('IndependentVariable-x')subplot(2,1,2)term1=r.*r/2/varR; %exponentray=(r/varR).*exp(-term1); %Rayleighpdfdel_r=r(3)-r(2); %determinebinwidthp_hist=N_samp/n/del_r; %probabilityfromhistogramplot(r,ray,'k',r,p_hist,'ok') %compareresultsylabel('ProbabilityDensity')xlabel('IndependentVariable-x')legend('truepdf','samplesfromhistogram',1)%Endofscriptfile.54假設(shè)有一組通過實(shí)驗(yàn)測得的數(shù)據(jù)。在這種情況下，盡管pdf可由數(shù)據(jù)的直方圖近似，但pdf和CDF都未知，我們的問題就是如何開發(fā)出一種算法，來產(chǎn)生一組樣本點(diǎn)使之與實(shí)驗(yàn)數(shù)據(jù)具有相似的pdf。第一步是產(chǎn)生實(shí)驗(yàn)數(shù)據(jù)的直方圖，假設(shè)結(jié)果如圖7-9中所示。一旦得出直方圖，就大體知道了pdf和CDF,因此可以運(yùn)用逆變換方法。這里介紹的方法是前一節(jié)中討論的逆變換方法的簡單擴(kuò)展。（二）直方圖法55樣本值x落在直方圖第i個(gè)直方（bin）上的概率為在x點(diǎn)處估算的CDF值可表示為這里接下來，令Fx(X)＝U,其中U是均勻隨機(jī)變量。于是有求解X得到56求取所需隨機(jī)數(shù)發(fā)生器的算法可分為以下三步實(shí)現(xiàn)：1.從一個(gè)能產(chǎn)生(0,1)區(qū)間均勻分布的隨機(jī)數(shù)發(fā)生器中取出樣本，產(chǎn)生隨機(jī)變量U。2.確定滿足如下條件的i值：其中Fi由式（7-49）定義。3.根據(jù)式（7-51）產(chǎn)生X,并將X返回給主調(diào)程序。

隨機(jī)數(shù)發(fā)生器的重現(xiàn)精度顯然取決于對應(yīng)直方圖的精度，而直方圖的精度又取決于可獲得的樣本數(shù)目。57(三)舍棄法用于產(chǎn)生具有期望或“目標(biāo)”pdffX(x)的隨機(jī)變量的舍棄（或接受）方法，基本都涉及到用函數(shù)Mgx(x)界定目標(biāo)pdf,這里gX(x)表示一個(gè)易于產(chǎn)生的隨機(jī)變量的pdf,M是一個(gè)滿足下式的適當(dāng)大的常數(shù)

在最簡單的形式下，gx(x)均勻分布在(0,a)區(qū)間。如果目標(biāo)pdffX(x),在(0,a)以外為零，則有由于Mgx(x)界定fX(x)，在上式中有58這一點(diǎn)如圖7-10中所示。用于產(chǎn)生pdf為fX(x)的隨機(jī)變量X的算法包括四步：1.產(chǎn)生在（0，1)上均勻分布的U1和U2。2.產(chǎn)生在（0，a)上均勻分布的V1，其中a是X的最大值。3.產(chǎn)生在（0，b)上均勻分布的V2，其中b不小于fX(x)的最大值。4.如果V2≤fx(V1)，令X＝V1;如果不等式不滿足，則丟棄V1和V2，從第一步開始重復(fù)以上過程。下面證明，以上算法產(chǎn)生的隨機(jī)變量X具有目標(biāo)pdffx(x)。5960因?yàn)閂1和V2都是均勻分布，所以樣本對（V1，V2）產(chǎn)生的點(diǎn)等概率地分布于ab區(qū)域內(nèi)的任何位置。V1的接受概率等于fX(x)所覆蓋的部分在區(qū)域ab所占比重，即圖7-10陰影部分面積跟總面積ab的比值，因此有由概率密度函數(shù)的定義，上式中分子等于1，于是有：因此有61這就定義了目標(biāo)pdf。例7-9在本例中，我們將剛討論的方法應(yīng)用于以下pdf:令a＝R，b＝M/R，就可獲得產(chǎn)生pdf為fx(x)的隨機(jī)變量X的算法。對這種特殊情況，圖7-10變?yōu)閳D7-11。MATLAB代碼如下：62取N＝3000點(diǎn)的情況執(zhí)行該程序，所得結(jié)果如圖7-12所示。3000點(diǎn)中，接受2301點(diǎn)，舍棄699點(diǎn)，效率達(dá)到76.70%，接近理論值（N->∞時(shí)）78.54%?！?3%File:c7_rejex1.mR=7; %defaultvalueofRM=4/pi; %valueofMN=input('InputnumberofpointsN>'); %setNfx=zeros(1,N); %arrayofoutputsamplesu1=rand(1,N);u2=rand(1,N); %generateu1andu2v1=R*u1; %generatev1v2=(M/R)*rand(1,N); %generatev2(g(x))kpts=0; %initializecounterfork=1:Nifv2(k)<(M/(R*R))*sqrt(R*R-v1(k)*v1(k));kpts=kpts+1; %incrementcounterfx(kpts)=v1(k); %saveoutputsampleendendfx=fx(1:kpts);[N_samp,x]=hist(fx,20); %gethistogramparameterssubplot(2,1,1)bar(x,N_samp,1) %plothistogramylabel('NumberofSamples')xlabel('IndependentVariable-x')subplot(2,1,2)yt=(M/R/R)*sqrt(R*R-x.*x); %calculatepdfdel_x=x(3)-x(2); %determinebinwidthp_hist=N_samp/kpts/del_x; %probabilityfromhistogramplot(x,yt,'k',x,p_hist,'ok') %compareylabel('ProbabilityDensity')xlabel('IndependentVariable-x')legend('truepdf','samplesfromhistogram',3)text=['Thenumberofpointsacceptedis',...num2str(kpts,15),'andNis',num2str(N,15),'.'];disp(text)%Endofscriptfile.64圖7-12-1運(yùn)用舍棄法得到的結(jié)果N=300065圖7-12-2運(yùn)用舍棄法得到的結(jié)果N=3000066盡管我們前面介紹的舍棄法是針對界定pdf為均勻分布這種情況，這里所講的方法只要稍加修改，便可適用于界定pdf非均勻的情況。理想情況下，應(yīng)該緊密界定目標(biāo)pdf,從而最小化舍棄概率。如果目標(biāo)pdf具有有限支撐（只在有限范圍內(nèi)有非零值），舍棄法的應(yīng)用效果良好。然而，有限支撐不是必需的，這里闡述的舍棄法可以很容易地?cái)U(kuò)展到無限支撐的情況。67從學(xué)習(xí)中可知，在通信系統(tǒng)中經(jīng)常會碰到高斯隨機(jī)變量，它恰當(dāng)?shù)亟o出了熱噪聲和許多其他現(xiàn)象的模型。在很多仿真中，高斯噪聲發(fā)生器都是一個(gè)基本的構(gòu)建模塊，因此，目前已研究出了多種產(chǎn)生高斯隨機(jī)變量的方法。高斯隨機(jī)變量的CDF為：四、產(chǎn)生不相關(guān)的高斯隨機(jī)變量其中Q(x)是高斯Q函數(shù)：68（一）均勻變量求和法

中心極限定理（CLT）為產(chǎn)生具有高斯pdf的隨機(jī)變量提供了一條很好的途徑。CLT表明，在很廣義的條件下，當(dāng)N->∞時(shí)，N個(gè)獨(dú)立隨機(jī)變量之和的pdf收斂為高斯隨機(jī)變量的pdf。高斯Q函數(shù)不能寫成閉合形式，無法使用逆變換方法，而用舍棄法的效率不高，因此需要尋找其他方法來產(chǎn)生高斯隨機(jī)變量。假設(shè)有N個(gè)獨(dú)立的均勻隨機(jī)變量Ui，i＝1，2，…，N。由這N個(gè)均勻隨機(jī)變量可以構(gòu)成如下隨機(jī)變量69其中B是常量，決定Y的方差。由CLT可知，N->∞時(shí)Y收斂為高斯隨機(jī)變量。因?yàn)镋｛Ui｝=1/2,所以Y均值等于注意到Ui-1/2的方差為可以由此求得Y的方差。因?yàn)榧僭O(shè)了各隨機(jī)變量Ui是相互獨(dú)立的于是有70因此，給定N值，通過選擇合適的B值，就可以將Y的方差設(shè)定為所需的任意值。選擇N時(shí)要折中考慮速度與所得pdf拖尾的精度，通常N取12，因?yàn)檫@會給出簡單結(jié)果B＝y(tǒng)。雖然根據(jù)中心極限定理產(chǎn)生高斯隨機(jī)變量的過程簡單直接。然而，在試圖將其應(yīng)用到數(shù)字通信的實(shí)際問題中時(shí)，會遇到一些較大的困難。首先，因?yàn)閁i-1/2的變化范圍是從-1/2到1/2，由式（7-63）可知Y值在-BN/2到BN/2之間變化，因此，盡管式（7-63）可以在均值附近很精確地近似高斯隨機(jī)變量的pdf,但pdf的拖尾被截短到±BN/2。如果仿真數(shù)字通信系統(tǒng)的目的是確定誤碼率，此時(shí)pdf的拖尾是十分重要的，因?yàn)閜df的拖尾表示導(dǎo)致傳輸差錯(cuò)的大噪聲。若N取得足夠大，對pdf進(jìn)行截尾所帶來的影響可以減到很小。71因此“近似高斯”pdf被截短后，只在以下區(qū)間內(nèi)取非零值具體而言，從式（7-67）可以推出B的值為當(dāng)N=100時(shí)，高斯隨機(jī)變量在±17.32y處被截短。對有些應(yīng)用來說，在17倍標(biāo)準(zhǔn)偏差處截短，仍然可能產(chǎn)生嚴(yán)重的誤差，因此要根據(jù)具體應(yīng)用來確定合適N值。72圖7-13說明了數(shù)字通信系統(tǒng)中對噪聲的概率密度函數(shù)進(jìn)行截短拖尾所帶來的困難，圖7-13畫出了在輸入端低SNR和高SNR時(shí)匹配濾波接收器輸出的條件概率密度函數(shù)。（pdf的拖尾實(shí)際上是連續(xù)的，圖7-13之所以這樣畫是為了強(qiáng)調(diào)截短的影響。）在傳輸二進(jìn)制0的條件下，pdf表示為fV(v|0);在傳輸二進(jìn)制1的條件下，pdf表示為fV(v|1)。在圖7-13中，假設(shè)噪聲方差為常數(shù)，通過改變信號功率來調(diào)節(jié)SNR。對接收器輸入端SNR足夠低這種情況，如圖7-13(a)所示，條件概率密度函數(shù)之間會有較多的重疊，確定的差錯(cuò)概率具有合理的精度。當(dāng)信號功率提高，條件概率密度函數(shù)進(jìn)一步分開，仿真精度會下降。由于對條件概率密度函數(shù)進(jìn)行了截尾，提高信號功率最終使得條件概率密度函數(shù)不再重疊，如圖7-13(b)所示。如果條件概率密度函數(shù)不重疊，無論SNR為何值，差錯(cuò)概率都為零，顯然不符合實(shí)際情況。73圖7-13接收器輸入端高SNR和低SNR時(shí)對應(yīng)的匹配濾波器輸出端的pdf若選擇較大的N值，用式（7-63）來近似高斯隨機(jī)變量又會造成另一個(gè)困難。因?yàn)楫a(chǎn)生一個(gè)X值要調(diào)用均勻隨機(jī)變量發(fā)生器N次，使用式（7-63）給出的算法可能需要過多的CPU時(shí)間。式（7-63）有兩個(gè)優(yōu)點(diǎn)：一個(gè)是它在Y的均值附近能很好地近似高斯隨機(jī)變量；另一個(gè)優(yōu)點(diǎn)是即使組成Y的隨機(jī)變量Ui的pdf不是均勻分布的，Y照樣可以近似為高斯的。高信噪比需要取較大N值，但會增大計(jì)算量，因此需要折中考慮！74雖然我們關(guān)心的主要是用單塊CPU進(jìn)行串行處理所作的仿真，應(yīng)該指出那些不適合于傳統(tǒng)串行處理應(yīng)用的算法可能會相當(dāng)適合用并行處理機(jī)來處理。例如，如果某個(gè)并行處理機(jī)上有100個(gè)CPU，那么在串行機(jī)上產(chǎn)生單個(gè)均勻隨機(jī)變量值所需的時(shí)間內(nèi)，并行機(jī)能夠產(chǎn)生100個(gè)均勻隨機(jī)變量值。大多數(shù)情況下，對100個(gè)均勻隨機(jī)變量求和可得到高斯隨機(jī)變量的一個(gè)很好的近似。如果采用并行處理，可以很快地完成這些工作。這僅僅是一個(gè)例子，說明算法的選擇取決于計(jì)算環(huán)境。75(二)瑞利隨機(jī)變量到高斯隨機(jī)變量的映射由例7-8可知，通過變換式R=SQRT(-22lnU)可由均勻隨機(jī)變量U產(chǎn)生瑞利隨機(jī)變量R?，F(xiàn)在來考慮將瑞利隨機(jī)變量映射為高斯隨機(jī)變量的問題。設(shè)有兩個(gè)獨(dú)立的高斯隨機(jī)變量X和Y,具有相同的方差2。由于X和Y相互獨(dú)立，聯(lián)合概率密度函數(shù)等于邊緣概率密度函數(shù)的乘積76令x＝rcos,y＝rsin，有通過以下變換可由fXY(x,y)求得聯(lián)合概率密度函數(shù)fR(r,):其中dAR

和dAXY分別表示在R,平面和X,Y平面的微分面積。由式（7-73）有77微分面積之比就是該變換的雅可比行列式，即化簡得78現(xiàn)在來考察R和的邊緣概率密度函數(shù)。R的pdf為的pdf為因此，可以看出R是一個(gè)瑞利隨機(jī)變量，是一個(gè)均勻隨機(jī)變量。由于瑞利隨機(jī)變量可由兩個(gè)正交的高斯隨機(jī)變量產(chǎn)生，可以推出，瑞利隨機(jī)變量的正交投影可以產(chǎn)生一對高斯隨機(jī)變量。因此，假設(shè)R是瑞利隨機(jī)變量，在（0，2）上均勻分布，則通過下面式子可產(chǎn)生高斯隨機(jī)變量X和Y79X和Y都是均值為零，方差為2的隨機(jī)變量。根據(jù)它們是高斯變量且不相關(guān)，可以推出X和Y為統(tǒng)計(jì)獨(dú)立的。因此，一對獨(dú)立的高斯隨機(jī)變量可由一對均勻分布的隨機(jī)變量U1和U2通過如下算式產(chǎn)生出來80這里我們使用了式（7-46）來求R。實(shí)現(xiàn)Box-muller算法的MATLAB程序如下：%File:c7_boxmul.mfunction[out1,out2]=c7_boxmul(N)u1=rand(1,N); %generatefirstuniformRVu2=rand(1,N); %generateseconduniformRVray=sqrt(-2*log(u1)); %generateRayleighRVout1=ray.*cos(2*pi*u2); %firstGaussianoutputout2=ray.*sin(2*pi*u2); %secondGaussianoutput%Endoffunctionfile.81（三）極坐標(biāo)法產(chǎn)生一對非相關(guān)、零均值的高斯隨機(jī)變量的另一個(gè)算法是極坐標(biāo)法。極坐標(biāo)法分為以下幾步：1.產(chǎn)生兩個(gè)獨(dú)立的隨機(jī)變量U1和U2，在區(qū)間（0,1)上均勻分布。2.令V1=2U1-1，V2＝2U2-1，使得V1和V2相互獨(dú)立且在區(qū)間（-1，1)上為均勻分布。3.令S=SQRT(V12+V22)。如果S＜1，進(jìn)入第四步；如果S>1,舍棄S值并跳回第一步。4.令A(yù)(S)＝SQRT(（-22lnS）/S)。5．令X＝A(S)V1，Y＝A(S)V2。

用極坐標(biāo)法產(chǎn)生一對高斯隨機(jī)向量的MATLAB代碼如下頁：82%File:c7_polar.mfunction[out1,out2]=c7_polar(N)u1=rand(1,N);u2=rand(1,N); %generateuniformRVsv1=2*u1-1;v2=2*u2-1; %makeuniformin-1to+1outa=zeros(1,N);outb=zeros(1,N); %allocatememoryj=1; %initializecounterfori=1:Ns(i)=v1(i)^2+v2(i)^2; %generatesifs(i)<1 %testj=j+1; %incremantcountera(i)=sqrt((-2*log(s(i)))/s(i));outa(j)=a(i)*v1(i); %firstGaussianRVoutb(j)=a(i)*v2(i); %secondGaussianRVendendout1=outa(1,1:j);out2=outb(1,1:j); %truncatearrays%Endoffunctionfile.83注意在上面的示例代碼中，用MATLAB庫函數(shù)rand產(chǎn)生所需的均勻隨機(jī)變量對。很明顯，這里也可使用用戶提供的基于LCG方法的均勻隨機(jī)數(shù)發(fā)生器。極坐標(biāo)法是舍棄法的一個(gè)例子，因?yàn)榭梢钥吹?，在第三步中舍棄了S的一些值，因此，每次調(diào)用該函數(shù)，產(chǎn)生的隨機(jī)數(shù)少于N。從圖7-14中很容易確定舍棄的概率。隨機(jī)變量V1和V2均勻分布于正方形里，正方形面積Abox＝4，該正方形的內(nèi)切圓的半徑R＝1，圓面積Acircle＝。因此，S被拒絕的概率為：8485通常，由極坐標(biāo)算法所得結(jié)果的相關(guān)特性要比Box-Muller算法好一些，然而，極坐標(biāo)法有一個(gè)缺點(diǎn)。N次調(diào)用Box-Muller算法會產(chǎn)生N對高斯隨機(jī)變量，而N次調(diào)用極坐標(biāo)算法產(chǎn)生高斯隨機(jī)變量對的個(gè)數(shù)卻是個(gè)隨機(jī)變量，其均值為(/4)N。也就是說，要產(chǎn)生一定數(shù)目的高斯隨機(jī)變量，需要調(diào)用極坐標(biāo)算法的次數(shù)未知，這個(gè)缺點(diǎn)往往會使仿真程序變復(fù)雜。86(四)MATLAB實(shí)現(xiàn)

MATLAB5版本以前，MATLAB庫中的高斯隨機(jī)數(shù)發(fā)生函數(shù)randn是應(yīng)用式（7-24）給出的最低標(biāo)準(zhǔn)隨機(jī)數(shù)發(fā)生器以及極坐標(biāo)法，將均勻分布的隨機(jī)數(shù)映射成高斯分布的隨機(jī)數(shù)。從MATLAB5開始，一個(gè)完全不同的算法取代了原來使用的算法，這個(gè)算法沒有用到乘法運(yùn)算或除法運(yùn)算。新的隨機(jī)數(shù)發(fā)生器不需要對均勻分布的隨機(jī)數(shù)進(jìn)行變換，而是直接產(chǎn)生具有高斯pdf的隨機(jī)數(shù)。因?yàn)樗惴ㄖ袥]用到乘法運(yùn)算和除法運(yùn)算，也不需要均勻分布到高斯分布的變換，速度非?？?。后來的MATLAB版本中用到的算法是基本Ziggurt算法的改進(jìn)版本。87到目前為止，產(chǎn)生的隨機(jī)變量都是不相關(guān)。我們現(xiàn)在將注意力轉(zhuǎn)移到另一種情況，即產(chǎn)生具有高斯pdf并且相關(guān)的隨機(jī)變量。首先考察一種簡單的方法，以產(chǎn)生具有給定相關(guān)系數(shù)（一階相關(guān)）的兩個(gè)序列。接著來考慮更一般的情況，即產(chǎn)生具有給定功率譜密度（PSD）的序列，當(dāng)然，確定功率譜密度等價(jià)于確定一個(gè)給定的自相關(guān)函數(shù)。（一）確定給定的相關(guān)系數(shù)我們已學(xué)了兩種方法來產(chǎn)生一對（近似）不相關(guān)的高斯隨機(jī)變量。要把一對不相關(guān)的高斯隨機(jī)變量，記作X和Y，映射成一對具有給定相關(guān)性的高斯隨機(jī)變量是一件很容易的事。假設(shè)X和Y不相關(guān)且具有零均值，接下來產(chǎn)生第三個(gè)隨機(jī)變量Z，定義如下五、產(chǎn)生相關(guān)的高斯隨機(jī)變量88這里是參數(shù)，滿足||≤1。在這個(gè)算法中，X與Z都是零均值隨機(jī)變量且具有相同方差。接下來將證明X與Z的相關(guān)系數(shù)是。

證明：由于Z是高斯隨機(jī)變量線性組合，所以Z也是高斯隨機(jī)變量。如果X和Y是零均值的，則Z也是零均值的。Z的方差為因?yàn)镋{XY}＝E{X}E{Y}=0,X2=Y2=2，上式變?yōu)椋?9協(xié)方差E{XZ}為上面最后一步是根據(jù)X和Y獨(dú)立且有零均值而推出。相關(guān)系數(shù)XZ為恰好等于所需數(shù)值。（二）確定任意的功率譜密度或自相關(guān)函數(shù)要確定一個(gè)隨機(jī)數(shù)序列，使之具有給定自相關(guān)函數(shù)或等價(jià)地具有給定的功率譜密度，一個(gè)通用的方法是，對一組不相關(guān)的樣本進(jìn)行適當(dāng)?shù)臑V波，從而使之具有目標(biāo)功率譜密度。根據(jù)定義，不相關(guān)樣本的功率譜密度在仿真帶寬|f|<fs/2上是常數(shù)，而方差就是Sn(f)曲線下的面積，如圖7-15所示，90因此，為了產(chǎn)生功率譜密度為N0的噪聲，隨機(jī)數(shù)（噪聲）發(fā)生器的方差和采樣頻率必須滿足下式可見，給定功率譜密度所要求的噪聲發(fā)生器的方差是采樣頻率的函數(shù)。要得到一個(gè)具有期望功率譜密度的隨機(jī)序列，可以相對直接地使用隨機(jī)過程基礎(chǔ)理論中的基本結(jié)果。我們知道，如果線性系統(tǒng)的輸入是功率譜密度為Sx(f)的隨機(jī)過程，則系統(tǒng)輸出端的功率譜密度為9192其中，H(f)是線性系統(tǒng)的傳遞函數(shù)，如圖7-16所示。如果噪聲的樣本獨(dú)立，則輸入的功率譜密度為常數(shù)K＝No/2Watts/Hz，于是有為得到目標(biāo)功率譜密度，所需的H(f)為因此，產(chǎn)生符合要求的功率譜密度的問題，簡化為尋找一個(gè)傳遞函數(shù)為H(f)的濾波器，使得H2(f)給出所需的頻譜形狀。93是給定傳遞函數(shù)的最小均方誤差擬合。

例7-10要設(shè)計(jì)具有給定傳遞函數(shù)的濾波器，一種簡便的方法是確定傳遞函數(shù)在最小均方誤差意義下的最佳擬合?？赏ㄟ^MATLAB函數(shù)yulewalk求解Yule-Walker方程來完成這項(xiàng)任務(wù)。具體而言，函數(shù)yulewalk求出濾波器的系數(shù)bk和ak，使得傳遞函數(shù)94為說明這項(xiàng)方法，假設(shè)給定的功率譜密度具有S(f)=K/f的形式。這稱作閃爍或1/f噪聲，經(jīng)常用來建模振蕩器的相位噪聲。要產(chǎn)生這種噪聲，必須產(chǎn)生一個(gè)濾波器，使其傳遞函數(shù)為H(f)=K/SQRT(f)。當(dāng)白噪聲通過濾波器，就產(chǎn)生所需的閃爍噪聲過程，因?yàn)楫?dāng)f->0時(shí)，H(f)->∞，所以傳遞函數(shù)定義為下面的MATLAB程序產(chǎn)生所要求的H(f)的近似值，其中頻率用奈奎斯特頻率fN作了歸一化，并且f0=fN/20。95運(yùn)行程序產(chǎn)生的結(jié)果如圖7-17所示，期望的傳遞函數(shù)用虛線表示，上面窗口表示的是它的3階近似，而下面的窗口表示的是20階近似，后者的效果明顯比前者好。注意，在H(f)相對平坦的頻率區(qū)域，用低階的近似值就足夠了；若H(f)在頻域內(nèi)變化很快，則需要高階的近似。這種方法是產(chǎn)生具有任意功率譜密度的一個(gè)強(qiáng)有力的工具?！鰣D7-17閃爍噪聲的產(chǎn)生96

例7-11一個(gè)運(yùn)動中的無線系統(tǒng)的仿真需要產(chǎn)生具有以下功率譜密度的隨機(jī)過程：以表示多普勒效應(yīng)。式（7-97）中fd表示最大多普勒頻率。如式（7-94）所示，所求濾波器的傳遞函數(shù)為對式（7-98）的樣本值作反DFT變化，得到一個(gè)FIR濾波器的沖激響應(yīng)，這個(gè)FIR濾波器即為所求的濾波器。97產(chǎn)生Jakes濾波器沖激響應(yīng)的程序以及用該濾波器濾除白噪聲的程序見附錄A。執(zhí)行代碼得到的結(jié)果如圖7-18和圖今-19所示。圖7-18表示其沖激響應(yīng)和傳遞函數(shù)H(f)（圖的下部），圖7-19表示復(fù)白噪聲通過該濾波器后的效果，其中圖的上部表示濾波器輸出端估計(jì)的功率譜密度，圖的下部描述的是對數(shù)尺度下的包絡(luò)幅度。在無線通信系統(tǒng)中，這對應(yīng)的是衰減包絡(luò)?！?8圖7-18脈沖響應(yīng)和目標(biāo)功率譜密度99圖7-19估計(jì)的功率譜密度和包絡(luò)函數(shù)100通常來說，要產(chǎn)生一個(gè)噪聲波形同時(shí)滿足給定的功率譜密度和pdf是一件比較困難的事，然而，在一種特定的情況下這個(gè)問題不難解決。如果系統(tǒng)輸出的pdf是高斯型的，我們可以簡單地產(chǎn)生濾波器輸入端的一個(gè)樣本序列，其中的樣本為獨(dú)立的高斯型。由于輸入的樣本值是相互獨(dú)立的，其功率譜密度是常數(shù)Kwatts/Hz。如前一節(jié)所述，可以使用一個(gè)線性濾波器來形成所需的功率譜密度，因?yàn)楣β首V密度成形濾波器是線性的，所以輸出的pdf為高斯型的。我們能得到這一簡單結(jié)果，是因?yàn)楦咚惯^程的任意線性變換會產(chǎn)生另一個(gè)高斯過程。幸運(yùn)的是，許多實(shí)際問題都能用這種方式來解決。六、同時(shí)確定pdf和PSD101

如果同時(shí)給定目標(biāo)波形的pdf和功率譜密度，而要求的pdf是高斯型以外的其他形式，這時(shí)問題就變難了許多。Sondhi針對這個(gè)問題提出了一個(gè)解決方法。實(shí)現(xiàn)Sondhi算法的基本框圖如圖7-20所示。跟慣常一樣，首先產(chǎn)生一列在（0,1)區(qū)間均勻分布的樣本{u[n]}。此外，假設(shè)序列｛u[n]｝是相關(guān)的，即其功率譜密度是白色的。流程中第一個(gè)無記憶的非線性器（記作MNL1)將序列{u[n]}變換成白高斯序列{x[n]}，我們學(xué)過的多種變換方法都可完成這種變換。濾波器進(jìn)行頻譜預(yù)矯正，使功率譜密度等于SY(f)。由于濾波器是線性的，所以序列{Y[n]}仍然具有高斯型pdf。流程中第二個(gè)無記憶非線性器(記作MNL2)將序列{y[n]}的高斯型pdf變換成最終所需的pdf。然而，第二個(gè)非線性器還影響SY(f)，因此，濾波器還必須改變預(yù)矯正過的功率譜密度SY(f)，使樣本序列在經(jīng)過第二個(gè)非線性器后，得出的序列{z(n)}同時(shí)具有所要求的功率譜密度和pdf。102圖7-20Sondhi算法103（七）PN序列發(fā)生器在很多應(yīng)用中尤其是在同步領(lǐng)域中會用到偽噪聲（PN）序列發(fā)生器。舉個(gè)例子，PN序列常用來近似具有均勻概率密度函數(shù)的隨機(jī)變量。PN序列發(fā)生器可以有很多種表示形式，最常見也是我們要集中討論的形式，如圖7-21所示。104在仿真中，使用PN序列最重要的一個(gè)原因是為了建立數(shù)據(jù)源模型。通過使用PN序列發(fā)生器，能用盡可能短的仿真時(shí)間，產(chǎn)生具有給定長度的所有可能的比特組合。在第10章討論半解析仿真時(shí)，我們將用這一特性來研究碼間干擾（ISI）的影響。

PN序列發(fā)生器由三個(gè)基本部分組成：N級移位寄存器、模二加法器和連接向量。這個(gè)連接向量具體定義了移位寄存器各級