1、通信系統(tǒng)的計(jì)算機(jī)模擬第十一講1第十講 回顧將均勻分布隨機(jī)變量映射成任意pdf逆變換法直方圖法舍棄法產(chǎn)生不相關(guān)的高斯隨機(jī)變量均勻變量求和法瑞利隨機(jī)變量到高斯隨機(jī)變量的映射極坐標(biāo)法MATLAB實(shí)現(xiàn)確定任意的功率譜密度或自相關(guān)函數(shù)27.7 PN序列發(fā)生器在很多應(yīng)用中尤其是在同頻領(lǐng)域中會用到偽噪聲（PN）序列發(fā)生器。PN序列常用來近似具有均勻概率密度函數(shù)的隨機(jī)變量。PN序列發(fā)生器可以有很多表示形式，最常見也是我們要集中討論的形式，如圖7-21所示。3作用在仿真中，使用PN序列最重要的一個原因是為了建立數(shù)據(jù)源模型。通過使用PN序列發(fā)生器，能用盡可能短的仿真時間，產(chǎn)生具有給定長度的所有可能的比特組合。4組

2、成PN序列發(fā)生器由三個基本部分組成：N級移位寄存器、模二加法器和連接向量。這個連接向量具體定義了移位寄存器各級與模二加法器之間的連接。它確定了發(fā)生器的性能特點(diǎn)，并由以下多項(xiàng)式定義 （7-99）5性質(zhì)可以證明，PN序列發(fā)生器輸出序列的最大周期是 （7-100）當(dāng)且僅當(dāng)多項(xiàng)式是本原多項(xiàng)式時，可以達(dá)到這個周期。PN序列發(fā)生器輸出的自相關(guān)函數(shù)如圖7-22所示，這里我們假設(shè)數(shù)據(jù)值（符號）為1。若或是L的整數(shù)倍，則自相關(guān)函數(shù)等于1；若，則自相關(guān)函數(shù)-1/L，并且在L足夠大時，趨近于零?？梢姡哂虚L周期的PN序列，其自相關(guān)函數(shù)近似為沖激函數(shù)。所以，跟期望的一樣，功率譜密度近似為白色的。PN序列有很多有趣的特

3、性，在仿真中我們關(guān)心的特性包括以下幾個：注：本原多項(xiàng)式：一個N階多項(xiàng)式g(D)能整除Dk+1的最小整數(shù)k是k=2N-167仿真中我們關(guān)心的特性序列幾乎是平衡的。也就是說，在序列的一個周期中，1的個數(shù)比0的個數(shù)只多出一個。在適當(dāng)?shù)奈恢蒙项~外添加一個0，則序列平衡。所有可能的比特組合都會出現(xiàn)一個周期中，這里全0序列除外，但包括全1序列。（注意：如果所有的移位寄存器包含二進(jìn)制0，則序列發(fā)生器會陷于全0的死循環(huán)。若已有N-1個0，再添加一個0，則序列平衡，得到的是deBruijn序列）。雖然自相關(guān)函數(shù)是周期的，但其波形與隨機(jī)的二進(jìn)制波形幾乎完全相同。8基于N級移位寄存器的PN序列發(fā)生器的設(shè)計(jì) 要找到一

4、個N階的本原多項(xiàng)式。在N較大時，很難找到這樣的本原多項(xiàng)式，幸運(yùn)的是，文獻(xiàn)13提供了大量本原多項(xiàng)式的表格。9開發(fā)步驟10例11pntaps = 0 0 1 0 0 0 0 0 0 1; % shift register tapspninitial = 0 0 0 0 0 0 0 0 0 1; % initial shift register statepndata = zeros(1,1023); % initialize output vectorsamp_per_sym = 5; % samples per symbolpnregister = pninitial;% initialize

5、shift registern = 0;% initialize counterkk = 0;% set terminator indicator12while kk = 0 n = n+1; % increment counter pndata(1,n) = pnregister(1,1);% save output feedback = rem(pnregister*pntaps),2); % calculate feedback pnregister = feedback,pnregister(1,1:9);% increment register if pnregister = pni

6、nitial; kk = 1; end % reset terminationEnd13text = The period is ,num2str(n,15),.;disp(text) % display periodpndata=replicate(pndata,samp_per_sym); % replicate datakn = n*samp_per_sym; % length of output vectorpndata = 2*pndata - 1; % make output +/- onea = fft(pndata);b = a.*conj(a); % PSD of dataR

7、m = real(ifft(b)/kn;% autocorrelationx1 = (0:length(Rm)-1)/samp_per_sym;x2 = 0:100;14% Now we can plot the results.subplot(3,1,1)plot(x1,Rm,.k); ylabel(Rm)subplot(3,1,2)stem(x2,Rm(1:101),.k); ylabel(Partial Rm)subplot(3,1,3)stem(x2,pndata(1:101),.k); ylabel(First 100 outputs)axis(0 100 -1.5 1.5);15注

8、意，該程序作了一個測試，以保證PN序列發(fā)生器確實(shí)是全周期。（對本例，執(zhí)行程序后返回的序列周期正好等于N=10時的全周期1023。）該程序還產(chǎn)生了自相關(guān)函數(shù)圖、自相關(guān)函數(shù)的前101個采樣以及發(fā)生器輸出端的前101個采樣。對采樣率為1采樣/符號時的結(jié)果如圖7-24所示，若將采樣率改為5采樣/符號，對應(yīng)的結(jié)果如圖7-25所示。1617187.8 信號處理在系統(tǒng)輸入是隨機(jī)的從而輸出也是隨機(jī)的這個前提下，線性系統(tǒng)的幾個輸入輸出關(guān)系，我們對以下幾個基本結(jié)論感興趣：系統(tǒng)輸入均值和輸出均值之間的關(guān)系系統(tǒng)輸入方差和輸出方差之間的關(guān)系輸入-輸出互相關(guān)系統(tǒng)輸入和系統(tǒng)輸出的自相關(guān)和功率譜密度之間的關(guān)系這些對仿真研究很

9、有用的關(guān)系，在系統(tǒng)輸入是一個隨機(jī)過程樣本函數(shù)的情況下，給我們提供了信號處理的基本工具。接下來的整個分析過程中，假設(shè)所研究的系統(tǒng)是線性的和固定不變的，并假設(shè)系統(tǒng)的輸入具有廣義平穩(wěn)性。 197.8.1 輸入-輸出均值20217.8.2 輸入-輸出互相關(guān)227.8.3 輸出自相關(guān)函數(shù)2324257.8.4 輸入-輸出方差26第八章 后處理 將仿真產(chǎn)生的數(shù)據(jù)處理成有用的形式后處理器通常是圖形密集型的視覺顯示比數(shù)字列表更容易理解，而后者是仿真程序最常見的數(shù)據(jù)輸出。不同系統(tǒng)的誤比特率曲線比包含同樣信息的數(shù)據(jù)表格能更快地傳達(dá)信息后處理程序可能引入（也可能不引入）相當(dāng)?shù)挠?jì)算復(fù)雜度。有些后處理器只是簡單地取得仿

10、真產(chǎn)生的數(shù)據(jù)，對其進(jìn)行適當(dāng)?shù)母袷骄幣藕?，產(chǎn)生合適的圖形輸出。產(chǎn)生誤比特概率或誤符號概率作為函數(shù)的曲線，作最少處理來產(chǎn)生圖形輸出的后處理器的其他例子還包括顯示信號波形、眼圖和散點(diǎn)圖的后處理器。有些后處理子程序會涉及相當(dāng)多的數(shù)據(jù)處理，它們中多數(shù)要涉及到某種估計(jì)，一個簡單的例子就是直方圖的產(chǎn)生，它是對概率密度函數(shù)的估計(jì)更復(fù)雜的例子有延遲、信碼的系統(tǒng)中，從信道誤符號率到譯碼后的誤比特率的映射。27基本圖形方法 系統(tǒng)實(shí)例/4 DQPSK傳輸28發(fā)射信號的相位偏移由d(k)和q(k)的值以及前一符號周期的相位偏移（t1）決定 例如 001011011129例如 00101101113/4, - ,3/4,

11、0, /4,30波形、眼圖和散點(diǎn)圖 如果同相信道信號xd(t)畫在平在(D,t)上，正交信號xq(t)畫在平面上，就產(chǎn)生了為參數(shù)的三維信號。 Q,D,t31眼圖 眼圖給出系統(tǒng)性能的一種定性量度。明晰張開的眼孔通常表明好的性能而模糊的眼孔往往表明差的性能眼孔的大小跟符號同步器所要求的精度有關(guān)盡管眼圖不能對系統(tǒng)性能提供定量的量度，但很難想像一個性能優(yōu)異的系統(tǒng)會有不明確的眼圖。32假設(shè)用示波器顯示波形，并且示波器在垂直虛線處觸發(fā)，結(jié)果就是如圖8-4（b）所示的三段眼圖 33例8-1產(chǎn)生并顯示DQPSK系統(tǒng)中的一些重要信號，仿真這個系統(tǒng)并產(chǎn)生圖形輸出。34m = 200;bits = 2*m;% nu

12、mber of symbols and bitssps = 10; % samples per symboliphase = 0; % initial phaseorder = 5; % filter orderbw = 0.2; % normalized filter bandwidth % initialize vectorsdata = zeros(1,bits); d = zeros(1,m); q = zeros(1,m);dd = zeros(1,m); qq = zeros(1,m); theta = zeros(1,m);thetaout = zeros(1,sps*m); %

13、 set direct and quadrature bit streamsdata = round(rand(1,bits); dd = data(1:2:bits-1);qq = data(2:2:bits);35% main programstheta(1) = iphase; % set initial phasethetaout(1:sps) = theta(1)*ones(1,sps);for k=2:m if dd(k) = 1 phi_k = (2*qq(k)-1)*pi/4; else phi_k = (2*qq(k)-1)*3*pi/4; end 36theta(k) =

14、phi_k + theta(k-1); for i=1:sps j = (k-1)*sps+i; thetaout(j) = theta(k); endendd = cos(thetaout);q = sin(thetaout);b,a = butter(order,bw);df = filter(b,a,d);qf = filter(b,a,q);37% postprocessor for plottingkk = 0; % set exit counterwhile kk = 0 % test exit counterk = menu(pi/4 QPSK Plot Options,. Un

15、filtered pi/4 QPSK Signal Constellation,. Unfiltered pi/4 QPSK Eye Diagram,. Filtered pi/4 QPSK Signal Constellation,. Filtered pi/4 OQPSK Eye Diagram,. Unfiltered Direct and Quadrature Signals,. Filtered Direct and Quadrature Signals,. Exit Program); 381濾波前的/4DQPSK信號星座圖2濾波前的/4DQPSK信號眼圖3濾波后的/4DQPSK信

16、號星府圖4濾波后的/4DQPSK信號眼圖5濾波前的同相和正交信號6濾波后的同相和正交信號7退出程序（返回MATLAB命令提示符）39if k = 1 sigcon(d,q) % plot unfiltered signal con. pause elseif k =2 dqeye(d,q,4*sps) % plot unfiltered eye diagram pause elseif k = 3 sigcon(df,qf) % plot filtered signal con. pause elseif k = 4 dqeye(df,qf,4*sps) % plot filtered eye

17、 diagram pause elseif k = 5 numbsym = 10; % number of symbols plotted dt = d(1:numbsym*sps); % truncate d vector qt = q(1:numbsym*sps); % truncate q vector dqplot(dt,qt) % plot truncated d and q signals40pause elseif k = 6 numbsym = 10; % number of symbols to be plotted dft=df(1:numbsym*sps); % trun

18、cate df to desired value qft=qf(1:numbsym*sps); % truncate qf to desired value dqplot(dft,qft) % plot truncated signals pause elseif k = 7 kk = 1; % set exit counter to exit value endend41支撐子程序function =sigcon(x,y)plot(x,y)axis(square)axis(equal)xlabel(Direct Channel)ylabel(Quadrature Channel)% End

19、of function file.42function = dqeye(xd,xq,m)lx = length(xd); % samples in data segmentkcol = floor(lx/m); % number of columnsxda = 0,xd; xqa = 0,xq; % append zerosfor j = 1:kcol % column index for i = 1:(m+1) % row index kk = (j-1)*m+i; % sample index y1(i,j) = xda(kk); y2(i,j) = xqa(kk); endend43su

20、bplot(211) % direct channelplot(y1);title(D/Q EYE DIAGRAM);xlabel(Sample Index);ylabel(Direct);subplot(212) % quadrature channelplot(y2);xlabel(Sample Index);ylabel(Quadratute);subplot(111) % restore44unction = dqplot(xd,xq)lx = length(xd);t = 0:lx-1;nt = t/(lx-1);nxd = xd(1,1:lx);nxq = xq(1,1:lx);s

21、ubplot(211)plot(nt,nxd);a = axis;axis(a(1) a(2) 1.5*a(3) 1.5*a(4);title(Direct and Quadrature Channel Signals);xlabel(Normalized Time);ylabel(Direct);subplot(212)plot(nt,nxq);a = axis;axis(a(1) a(2) 1.5*a(3) 1.5*a(4);xlabel(Normalized Time);ylabel(Quadratute);subplot(111)45464748估計(jì)-直方圖 如果可以獲得隨機(jī)過程的一組

22、樣本（就像在仿真環(huán)境下那樣），由這組樣本得出的直方圖經(jīng)常可用于估計(jì)內(nèi)在的概率密度函數(shù)（pdf）直方圖將總共N個樣本的數(shù)據(jù)分組成B個直方（bin）或單元（cell），假定每個直方的寬度為W，其中心記為bi。如果下式成立，給定的樣本落在第個直方49Ni表示落在第個直方中的樣本數(shù)目。 CountN：R來表示全部N個樣本中進(jìn)入由R所定義的直方中的樣本數(shù)目 50然后畫出條狀圖，每條的高度和Ni成正比，中心在bi為了用作pdf的估計(jì)器，我們對直方圖進(jìn)行縮放使其總面積為1，將Ni除以NW即可實(shí)現(xiàn)這一點(diǎn)，于是，每條的高度為Ni/NW。代表第個直方的條的面積Ai，等于條的高度乘以寬度。因此，注意到Ai代表第個直

23、方的相對頻率 注意到每個直方是用一個常數(shù)來表示在有限寬度上的pdf。對給定直方中的某個點(diǎn)，pdf的估計(jì)是無偏的，然而對給定直方所定義的大多數(shù)區(qū)域，估計(jì)是有偏的 51subplot(2,2,1)x = randn(1,100); hist(x,20)ylabel(N_i); xlabel(a)subplot(2,2,2)x = randn(1,100); hist(x,5)ylabel(N_i); xlabel(b)subplot(2,2,3)x = randn(1,1000); hist(x,50)ylabel(N_i); xlabel(c)subplot(2,2,4)x = randn(1,

24、100000); hist(x,50)ylabel(N_i); xlabel(d)5253功率密度估計(jì) 我們感興趣的波形是隨機(jī)過程的樣本函數(shù)，這就使得在給定頻率值f1處的PSD是一個隨機(jī)變量需要最小化譜估計(jì) 的方差基于快速傅里葉變換（FFT），而且經(jīng)常用于仿真環(huán)境 54周期圖最簡單、最快速也最常用的PSD估計(jì)算法是周期圖。其定義如下 N是數(shù)據(jù)記錄里的樣本總數(shù)，IN(Kf)是用于計(jì)算在頻率f=Kf 處PSD估計(jì)的數(shù)據(jù)的N點(diǎn)FFT。周期圖的問題在于它是有偏差的 對于足夠大的N，偏差可以忽略。主要的問題在于大方差 55帶數(shù)據(jù)窗的周期圖默認(rèn)的矩形窗 每個樣本值 xn 乘以 wn=1, 0nMN。第i段

25、數(shù)據(jù)記錄的周期圖由下式給出：58對K個周期圖進(jìn)行平均 該估計(jì)值是有偏差的，因?yàn)閿?shù)據(jù)記錄是有限長的。假設(shè)K個周期圖是互相獨(dú)立的，則由由于對于固定的我們希望值盡可能大，所以各段往往互相交疊，一般采用50%交疊。當(dāng)采用了兩，值由增加到。如果數(shù)據(jù)段交疊，那么，個周期圖之間不再獨(dú)立，PSD估計(jì)器方差的減小就比式（8-18）預(yù)計(jì)的要少。數(shù)據(jù)窗的使用，至少在某種程度上有助于恢復(fù)段之間的獨(dú)立性。 59常用窗函數(shù)Hanning窗定義如下： 60例8-3將獨(dú)立（白噪聲）樣本通過有5dB通帶紋波的切比雪夫?yàn)V波器。試估計(jì)濾波器輸出端的PSD。61Matlab programsettle = 100; % ignore

26、 transientfs = 1000; % sampling frequencyN = 50000; % size of data recordf = (0:(N-1)*fs/N; % frequency scaleb,a = cheby1(5,5,0.1); % filterNN = N+settle; % allow transient to diein = randn(1,NN); % random inputout = filter(b,a,in); % filter outputout = out(settle+1):NN); % strip off initial samples

27、window = hanning(N); % set window function62winout = out.*window; % windowed filter outputfout = abs(fft(winout,N).2; % transform and square magU = sum(window.*window); % window energyf1out = fout/U; % scale spectrumpsd1 = 10*log10(abs(f1out); % log scalesubplot(2,1,1)plot(f(1:5000),psd1(1:5000)grid

28、; axis(0 100 -70 10);63xlabel(Frequency, Hz)ylabel(PSD) %K = 25; % number ofsegments M = N/K; % block sizefK = (0:(M-1)*fs/M; % frequency scaled = zeros(1,M); % initialize vectorpsdk = zeros(1,M); % initialize vectorwindow = hanning(M); % set window functionU = sum(window.*window); % window energy64

29、for k=1:Kfor j=1:M index = (k-1)*M+j; d(j) = out(index); end dwin = d.*window; psdk = (abs(fft(dwin,M).2)/U + psdk;endpsd2 = 10*log10(psdk/K);subplot(2,1,2)plot(fK(1:250),psd2(1:250)grid; axis(0 100 -70 10);xlabel(Frequency, Hz)ylabel(PSD)65Result66MATLAB PSD估計(jì)MATLAB信號處理工具箱包含大量用于PSD估計(jì)的子程序其中兩個是psd和pw

30、elch67例8-4估計(jì)QPSK信號的PSD。假設(shè)采用矩形脈沖成型，并以每符號16個點(diǎn)的速率采樣QPSK信號的同相和正交分量。68fs = 16;x = random_binary(1024,fs)+i*random_binary(1024,fs);for nwin=1:4 nwindow = nwin*1024; pxx,f = pwelch(x,nwindow,fs); pxx = pxx/sum(sum(pxx); n2 = length(f)/2; pxxdB = 10*log10(pxx/pxx(1); ptheory = sin(pi*f+eps)./(pi*f+eps); pth

31、eory = ptheory.*ptheory; ptheorydB = 10*log10(ptheory/ptheory(1); subplot(2,2,nwin) plot(f(1:n2),pxxdB(1:n2),f(1:n2),ptheorydB(1:n2) ylabel(PSD in dB) xx = window length = ,num2str(nwindow); xlabel(xx) axis(0 8 -50, 10); grid;end 69Result70增益，延遲和信噪比 信噪比是一個常用于評估通信系統(tǒng)性能的指標(biāo)這里介紹的SNR估計(jì)方法源于寬帶遙測系統(tǒng)中一種測量信道失真誤

32、差的方法定義該系統(tǒng)中某點(diǎn)信號上的噪聲為該點(diǎn)理想信號和實(shí)測信號之間的均方誤差（MSE），將理想信號定義為系統(tǒng)輸入端的信息承載信號按幅度比例縮放并加入延遲后所得的信號，以此來估計(jì)SNR在過去，這種方法的應(yīng)用包括監(jiān)視數(shù)字脈碼調(diào)制數(shù)據(jù)的可靠傳輸和非線性信道下載波交調(diào)比（carrier-to-intermodulation ratio）的估計(jì)。71實(shí)低通信號的理論推導(dǎo)對線性時不變無失真系統(tǒng)，系統(tǒng)中任意點(diǎn)信號是對輸入?yún)⒖夹盘柊捶缺壤s放并加入延遲后所得的信號。因此，我們可以將無失真信號寫為：觀測信號N（t)代表外部加性噪聲，d(t)是系統(tǒng)引入的跟信號相關(guān)的內(nèi)部失真，這種內(nèi)部失真可能由碼間干擾或非線性造成。 7273噪聲功率定義A和的理想估計(jì)值最小化 對于平穩(wěn)信號，各階矩獨(dú)立于時間起點(diǎn)。此外，和的期望等于期望的和 7475Y(t)中信號分量的功率為因此，信噪比為76實(shí)現(xiàn)系統(tǒng)增益、延遲和SNR的后處理估計(jì)器MATLAB函數(shù) function gain,delay,px,py,rxy,rho,snrdb = snrmse(x,y)ln = le