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

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

本章首先介紹隨機過程樣本函數(shù)的產(chǎn)生，在仿真環(huán)境下考察平穩(wěn)的概念，然后簡要介紹一下數(shù)字調(diào)制器的仿真模型。有了這些預備知識后，我們再轉向本章的重點，來討論如下問題：4※在（0，1)上產(chǎn)生均勻分布且不相關的隨機數(shù)※將不相關且均勻分布的隨機數(shù)映射成不相關的、具有任意（需要的）概率密度函數(shù)（pdf）的隨機數(shù)※產(chǎn)生不相關的、具有高斯型pdf的隨機數(shù)※產(chǎn)生相關的、具有高斯型pdf的隨機數(shù)※產(chǎn)生相關的、具有任意（需要的）pdf的隨機數(shù)※偽噪聲（PN）序列的產(chǎn)生以及幾種應用于隨機數(shù)序列的計算方法。5一、平穩(wěn)與遍歷性過程

在仿真通信系統(tǒng)時，用于表示信號、噪聲和干擾而產(chǎn)生的樣本函數(shù)通常假設為各態(tài)歷經(jīng)的。作這一假設，是因為我們通常依次處理系統(tǒng)中波形的時域樣本，而且系統(tǒng)中的每個點上只有一個波形（樣本函數(shù)）。我們還假設仿真所處理的波形是由其內(nèi)在的統(tǒng)計模型定義的總體(Ensemble）中的一個典型成員，這時各種統(tǒng)計量如各階矩、信噪比和誤比特率，就可以當作時間平均來計算。將仿真結果與對應的理論結果進行比較時，常有一個潛在的假設：仿真計算得到的時間平均等于總體均值。因此，隱含的假設條件是對應的隨機過程為遍歷性的。6遍歷性過程總是平穩(wěn)的，所以，總是可以假設仿真中產(chǎn)生的樣本函數(shù)屬于平穩(wěn)隨機過程。由基本的隨機過程理論可知，平穩(wěn)性是這樣定義的：所有統(tǒng)計量與時間原點無關。為了說明其中的幾點思想，先來看一個簡單的例子。例7-1假設隨機過程的樣本函數(shù)表達式為式中，i是對應隨機試驗的樣本空間中的一個輸出，每一個i映射為一個相位i。再假定對應的隨機試驗就是從均勻隨機數(shù)發(fā)生器輸出端抽取一個數(shù)，抽取得到的結果i=ui，這里ui在(0,1)區(qū)間上均勻分布。然后ui映射成相位i＝kui。當幅度A和頻率f固定，i的值便決定了波形。7這個例子中我們感興趣的是k的兩個值：k=2，對應的相位均勻分布在（0，2）區(qū)間，和k＝/2，對應的相位均勻分布在（0，/2）區(qū)間。作為第二個例子，假設隨機過程由以下表達式描述：在這種情況下，幅度在（A，2A）上均勻分布。下面MATLAB程序產(chǎn)生一個隨機過程的三組樣本函數(shù)。1）波形表示為x(t)，對應式（7-1)中k＝2的情況。2）波形表示為y(t)，對應式（7-1)中k＝/2的情況。3）波形表示為z(t)，由式（7-2）所定義。對于所有的波形有A＝1，f＝1。每個仿真過程產(chǎn)生兩秒鐘長度的數(shù)據(jù)以及20個樣本函數(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ù)的時間平均都等于零。很容易便可證明：對于大量滿足0≤ti≤2的ti,計算出來的x(t)的總體均值近似為零；若樣本函數(shù)的個數(shù)趨近∞，時間平均收斂于零。然而，對于y(t)，當t在0.875和1.875附近時總體均值近似為1，當t在0.375和1.375附近時總體均值近似為-1，而當t在0.125、0.625、1.125和1.625附近時總體均值近似為零。這是周期平穩(wěn)過程的一個例子，這種過程的矩是周期的。

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

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

例7-2在后面的工作中，我們經(jīng)常要用到數(shù)字調(diào)制器的模型。這些調(diào)制器基本構建模塊是函數(shù)randombinary，它產(chǎn)生電平值為+1或-1的二進制波形，產(chǎn)生的比特數(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可以仿真多個數(shù)字調(diào)制器，例如，可用如下MATLAB語句仿真一個QPSK調(diào)制器：x=randombinary(nbits,nsamples)+i*randombinary(nbits,nsamples);下列MATLAB程序產(chǎn)生一個長10比特的QPSK信號，采樣頻率為每比特8個采樣點：%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運行該程序，得到QPSK信號的同相分量和正交分量，如圖7-2所示。圖7-2QPSK信號的同相分量和正交分量（10比特）15圖7-2-QPSK信號的同相分量和正交分量（100比特）16具有均勻概率密度函數(shù)的隨機變量很容易轉換成具有其他所需pdf的隨杭變量，因此，要產(chǎn)生一個具有特定pdf的隨機變量，首先需要產(chǎn)生一個在（0，1)區(qū)間均勻分布的隨機變量。通常，先產(chǎn)生一列介于0和M之間的數(shù)（整數(shù)），然后將序列中每個元素除以M。實現(xiàn)隨機數(shù)發(fā)生器最常用的方法是線性同余（linearcongruence).（一）線性同余

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

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

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

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

例7-4前一例子中設計的LCG確實具有周期m＝50000在以下MATLAB程序中，輸入種子數(shù)，并運行程序直到種子數(shù)再次出現(xiàn)。設產(chǎn)生n個整數(shù)，如果n＞m而種子數(shù)沒有再次出現(xiàn)，則認為發(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.可以看到，跟期望的一樣，周期確實是5000。■20方法B：具有素模數(shù)的乘性算法·乘性發(fā)生器的定義式為

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

．m是素數(shù)（通常要求m取較大值）

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

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

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

．乘子amod（8）結果為3或5

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

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

xi＋1＝(1229xi＋1)mod(2048）（7-18）用例7-4中的程序c7_LCGperiod.m驗證，這兩個發(fā)生器都是全周期的。每個發(fā)生器產(chǎ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)生的散點圖見圖7-3。人們希望得到的散點圖中縱坐標xi+1和橫坐標xi的各種組合都26會出現(xiàn)，這種情形下的散點圖缺乏結構性，從圖7-3看來，a＝65的發(fā)生器比a＝1229的發(fā)生器具有更小的序列相關性。在例7-6中我們將看到情況確實如此。■圖7-3a1=65（左圖）和a2=1229（右圖）時的散點圖27Durbin-Watson測試Durbin-Watson獨立性測試可以通過計算Durbin參數(shù)來完成：其中X[n]是一個零均值(zero-mean)隨機變量。下面會講到當D值在2附近時，X[n]和X[n-1］相關性很小。為了說明Durbin-Watson測試的性質，假設X[n]和X[n-1]相關并且X[n]是一個遍歷性過程。為了簡化記號，假設N足夠大，使得N-1≈N,式（7-19)可寫成：28式中X表示X[n]，Y表示X[n-1]，E{·｝表示數(shù)學期望。由于假定X[n]和X[n-1]相關，可令其中X和Z不相關，是X和Y的相關系數(shù)。注意，X、Y和Z有相同的方差，用2表示。將式（7-21）代入式（7-20）得：X和Z不相關且均值為零，所以中間項等于零。由于X和Z方差相同，所以有29因為-1≤≤1，Durbin參數(shù)D在O與4之間，且有＝O時D＝2。D＜2表示正相關，而D＞2表示值為負。以下的MATLAB函數(shù)計算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本例將計算例7-5中討論到的兩個噪聲發(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對應的相關值約為0;a2＝1229對應的相關值約為0.2。因此，從Durbin測試結果得出a1=65給出的結果優(yōu)于a2＝1229，這個結論與圖7-3所示的散點圖一致?！觯ㄈ┳畹蜆藴时砻鹘o定LCG能通過各種隨機性的統(tǒng)計測試，以便徹底測試它的質量，這是一項首要任務。對于產(chǎn)生長的序列來說更是如此。為了部分地解決這個問題，多種算法被確定為最低標準（minimumstandard）算法。最低標準算法具有以下性質：·全周期/·能通過所有適用的隨機性統(tǒng)計測試、·易于從一臺計算機移植到另一臺計算機32一旦算法經(jīng)確定并適當?shù)赜涗浵聛恚统闪艘粋€最低標準算法，于是這個算法可以在別處放心地使用而無須額外的測試。如果使用一個最低標準算法，不用擔心算法本身的正確性，但必須保證在給定計算環(huán)境下算法得到正確的實現(xiàn)編程時要特別注意的是，由算法產(chǎn)生的所有的數(shù)都要能唯一地表示出來。33Lewis，Goodman和Miller最低標準

Lewis，Goodman和Miller最低標準定義為：

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

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

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

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

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

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

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

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

在最簡單的形式下，gx(x)均勻分布在(0,a)區(qū)間。如果目標pdffX(x),在(0,a)以外為零，則有由于Mgx(x)界定fX(x)，在上式中有58這一點如圖7-10中所示。用于產(chǎn)生pdf為fX(x)的隨機變量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，從第一步開始重復以上過程。下面證明，以上算法產(chǎn)生的隨機變量X具有目標pdffx(x)。5960因為V1和V2都是均勻分布，所以樣本對（V1，V2）產(chǎn)生的點等概率地分布于ab區(qū)域內(nèi)的任何位置。V1的接受概率等于fX(x)所覆蓋的部分在區(qū)域ab所占比重，即圖7-10陰影部分面積跟總面積ab的比值，因此有由概率密度函數(shù)的定義，上式中分子等于1，于是有：因此有61這就定義了目標pdf。例7-9在本例中，我們將剛討論的方法應用于以下pdf:令a＝R，b＝M/R，就可獲得產(chǎn)生pdf為fx(x)的隨機變量X的算法。對這種特殊情況，圖7-10變?yōu)閳D7-11。MATLAB代碼如下：62取N＝3000點的情況執(zhí)行該程序，所得結果如圖7-12所示。3000點中，接受2301點，舍棄699點，效率達到76.70%，接近理論值（N->∞時）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運用舍棄法得到的結果N=300065圖7-12-2運用舍棄法得到的結果N=3000066盡管我們前面介紹的舍棄法是針對界定pdf為均勻分布這種情況，這里所講的方法只要稍加修改，便可適用于界定pdf非均勻的情況。理想情況下，應該緊密界定目標pdf,從而最小化舍棄概率。如果目標pdf具有有限支撐（只在有限范圍內(nèi)有非零值），舍棄法的應用效果良好。然而，有限支撐不是必需的，這里闡述的舍棄法可以很容易地擴展到無限支撐的情況。67從學習中可知，在通信系統(tǒng)中經(jīng)常會碰到高斯隨機變量，它恰當?shù)亟o出了熱噪聲和許多其他現(xiàn)象的模型。在很多仿真中，高斯噪聲發(fā)生器都是一個基本的構建模塊，因此，目前已研究出了多種產(chǎn)生高斯隨機變量的方法。高斯隨機變量的CDF為：四、產(chǎn)生不相關的高斯隨機變量其中Q(x)是高斯Q函數(shù)：68（一）均勻變量求和法

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

用極坐標法產(chǎn)生一對高斯隨機向量的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)生所需的均勻隨機變量對。很明顯，這里也可使用用戶提供的基于LCG方法的均勻隨機數(shù)發(fā)生器。極坐標法是舍棄法的一個例子，因為可以看到，在第三步中舍棄了S的一些值，因此，每次調(diào)用該函數(shù)，產(chǎn)生的隨機數(shù)少于N。從圖7-14中很容易確定舍棄的概率。隨機變量V1和V2均勻分布于正方形里，正方形面積Abox＝4，該正方形的內(nèi)切圓的半徑R＝1，圓面積Acircle＝。因此，S被拒絕的概率為：8485通常，由極坐標算法所得結果的相關特性要比Box-Muller算法好一些，然而，極坐標法有一個缺點。N次調(diào)用Box-Muller算法會產(chǎn)生N對高斯隨機變量，而N次調(diào)用極坐標算法產(chǎn)生高斯隨機變量對的個數(shù)卻是個隨機變量，其均值為(/4)N。也就是說，要產(chǎn)生一定數(shù)目的高斯隨機變量，需要調(diào)用極坐標算法的次數(shù)未知，這個缺點往往會使仿真程序變復雜。86(四)MATLAB實現(xiàn)

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

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

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

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

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

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