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

通信系统的计算机模拟第十二讲1第一页，共五十九页，编辑于2023年，星期三关于期末报告下下次课为Project报告时间今天请提交纸质报告一份（首页注明姓名、学号信息）电子版和Matlab程序email到:btian01@请将email标题设为Simulation2008＋姓名＋学号请准备好报告PPT,每人报告时间为（10±2）分钟如有特殊时间需求请提前声明2第二页，共五十九页，编辑于2023年，星期三第十一讲回顾将仿真产生的数据处理成有用的形式波形、眼图和散点图，书上的例子要学会，典型的标准模块估计直方图功率密度估计周期图带数据窗的周期图分段周期图增益，延迟和信噪比3第三页，共五十九页，编辑于2023年，星期三蒙特卡罗方法导论一个简单易懂的场院景中考察一些基本方法对置信区间和收敛性等重要问题也作了简单的讨论假定设估计器所采有的观测值是相到独立的4第四页，共五十九页，编辑于2023年，星期三9.1蒙特卡罗方法的基本概念蒙特卡罗仿真建立在几率游戏的基础上“蒙特卡罗”——赌博著称的地中海城市。蒙特卡罗仿真是那些利用蒙特卡罗方法估计系统参数如误比特率（BER）的仿真而蒙特卡罗估计同是指通过内在的随机试验来估计参数值的过程。5第五页，共五十九页，编辑于2023年，星期三9.1.1相对频率蒙特卡罗估计是基于概率的相对频率解释的为定义相对频率，第一步是确定随机试验和一个感兴趣的事件A。随机试验中试验结结果无法准确预测，用统计的方法加以描述。一个最基本的随机试验就是掷硬币，正面、反面随着机试验中的事件可以是一个结果或几个结果的集合以数学字通信系统为例，随机试验可定度为发送一个二进制数1，接收机输出端结果为对所发送二进制符号的估计，它是二进制数0或者二进制数1我们所感兴趣的事件可能在发送1的过程中所产生的差错要确定系统的BER，得估计在发送1的条件下接收到0这一事件的概率。大量的随机试验，试次数为N，以表示事件A发生的次数。将事件A发生的概率近似为相对频率，其定义为。在相对频率的意义下，事件A的概率可以通过重复无限多次随机实验获得N是总的发送比特数或符号数，NA－错误数6第六页，共五十九页，编辑于2023年，星期三讨论

由于试验的随机性，在试验次数N有限的情况下，NA是随机变量，Pr(A)也是随机变量。估计器（仿真）随机变量的统计性质决定了估计器精度，因而也决定了仿真的质量。7第七页，共五十九页，编辑于2023年，星期三9.1.2无偏和一致估计器蒙特卡罗估计器必须满足几个重要性质在实际中才有意义首先，我们希望蒙特卡罗估计器是无偏的。如果A的估计值Ă，我们希望换句话说，在平均意义下可以得到正确的结果。估值有比较小的方差如果估计值是无偏的并具有小的方差，则估计器所产生的估计值会聚集在待估计参数真值的周围，且有较小的散布范围。除非所研究的事件是统计独立的，用解析的方法确定蒙特卡罗估计器的方差通常是困难的，但是几乎可以肯定，随着估计值数量的增加而减少，称为一致的。对于一致估计器，当时，而对于无偏和和一致估计器，误差具有零均值，而其方差在收敛为0。收敛过程通常非常缓慢8第八页，共五十九页，编辑于2023年，星期三9.1.3蒙特卡罗估计考虑确定一个形状不规则的区域的面积作为一个蒙特卡罗估计器的简单例子。侍估计面积的区域完全包含在一个面积已知的方框中，随机试验定义在包围方框中随机抽取采样，事件A定义为采样点落在面积待确定区域内。要得到一个未知面积的无偏估计器，随机采用样点必须均匀分布在面积已知范围区域中，利用具有两个均匀随机数发生器的计算机程序就可以很容易地实现这一点。9第九页，共五十九页，编辑于2023年，星期三Step110第十页，共五十九页，编辑于2023年，星期三下一步是定义感兴趣的事件A。我们要估计如图9-2所示的爆炸形（sunburst）区域的面积，分别定义和落在包围中和爆炸形区域的采样点数因为采样点在包围方框中是均匀分布的，近似为采样点落在钻石形区域内的数目与落在包围方框中的数目之比，即也就是在采样点为均匀分布这一条件下，随着采样点数的增加，近似精度会不断提高。11第十一页，共五十九页，编辑于2023年，星期三

9.1.4π的估计估计数值π的方法之一是用一个具有单位面积的正方形包围一个馅饼状（pie-shaped）的区域，即单位圆的第一象限。如图9-3所示为这种情况，以及总的采样点数Nbox。如果正方形在x轴上所点的区间是（0，1），在y轴上所占的区间也是（0，1），显然馅饼状区域（四分之一圆）的面积为山巅一寺一壶酒，尔乐苦煞吾，把酒吃，酒杀尔，杀不死，乐尔乐。3．1415926535897932384626

3．141592653589793238

12第十二页，共五十九页，编辑于2023年，星期三假设采样点是均匀分布的，则和Nbox

的比构成的无偏估计，所以13第十三页，共五十九页，编辑于2023年，星期三例9-1m=input('EnterM,thenumberofexperiments>');n=input('EnterN,thenumberoftrialsperexperiment>');z=zeros(1,m);data=zeros(n,m);forj=1:mx=rand(1,n);y=rand(1,n);k=0;fori=1:nifx(i)^2+y(i)^2<=1 %Fallinsidepieslice?k=k+1;enddata(i,j)=4*(k/i); %jthestimateofpiendz(j)=data(n,j); %Storedataendplot(data,'k') %Plotcurvesxlabel('NumberofTrials')ylabel('Estimateofpi')14第十四页，共五十九页，编辑于2023年，星期三Result如果对5个结果进行平均，则如果对5个结果进行平均，则

，这样的结果等价于2500次的试验结果15第十五页，共五十九页，编辑于2023年，星期三讨论这个例子具有所有蒙特卡罗仿真的许多重要特点其中都有一个测试条件以及两个计数器随机试验每进行一次，第一个计数器就增加1，如果测试条件每满足一次，第二个计数器就增加1。在数字通信系统的仿真中，仿真目的是估计误比特率。测试条件就是发送给定的比特或数据符号时是否发生了差错。在仿真中，每当发送了一个比特或数据符号，第一计数器就增加1。每当观测到一次差错，第二个计数器就增加1。下一节中的两个例子会演示这个方法。但我们先考察一下AWGN信道的特点。16第十六页，共五十九页，编辑于2023年，星期三9.2在通信系统中的应用——AWGN信道

为了用蒙特卡罗方法估计通信系统的性能，让N个符号通过系统（实际上是系统的计算机仿真模型），并计算发送差错的个数Ne。如果在N次的符号发送中有Ne次差错，则符号差错概率为这一估计是有偏的还是无偏的呢？它是否为一致估计？为了在这一个尽可能简单的场景下研究这些重要的问题，我们假设信道是AWGN（加性高斯白噪声）信道。在AWGN条件下，由信道噪声所产生的差错相互独立，而且在次符号发送中出现的差错次数可以用二项分布来描述17第十七页，共五十九页，编辑于2023年，星期三9.2.1二项式分布我们下面的任务是确定的统计特性。首先是确定的均值和方差。对于相互独立的差错事件，N次符号发送中有Ne次差错的概率服从以下的二项式分布是二项式分布系数，PE是单次发送的差错概率。18第十八页，共五十九页，编辑于2023年，星期三二项式分布的性质

及蒙特卡罗估计器的性质代入则差错概率的蒙特卡罗估计器的均值为差错概率的蒙特卡罗估计器是无偏的。差错概率蒙特卡罗估计器的方差为估计器是一致的

假设了二项式分布，这只有在差错事件相互独立这一条件下才成立

19第十九页，共五十九页，编辑于2023年，星期三讨论希望蒙特卡罗估计器无偏、一致无偏平均意义下正确小的方差实用的一致随着N的增大，方差减小PE未定，给定要求的方差，N不能确定有偏一致的呢？不正确，除非知道如何消除偏差20第二十页，共五十九页，编辑于2023年，星期三尽管知道估计器的性质非常重要，但在许多情况下很难证明估计器是无偏的和一致的本节所得到的所有结果很理想成立的条件是信噪声引起的差错相互独立，从而对应的差错概率是二项分布的如果差错事件是相关的，例如在带限信道的情况下，这里所给出的结果不再正确，从而面临的问题就更加困难不过，即使差错事件不是相互独立的，式（9-11）仍然是一个有效的差错概率估计器。21第二十一页，共五十九页，编辑于2023年，星期三例9-2在差错事件相互独立的情况下，可用掷硬币事件来模拟二进制发送，N个符号的发送可建模为对一枚不均匀的硬币进行N次投掷。可以假设第i次投掷结果为“反面”相当于对第i次发送作出了正确的判决，而第i次投掷结果为“正面”相当于对第i次发送作出了错误的判决。在这个例子中，掷币过程所具有的统计特性由仿真来决定。因为投掷是相互独立的，所以这个试验模拟了AWGN信道中的二进制数据传输。22第二十二页，共五十九页，编辑于2023年，星期三Solution假设投掷结果为“反面”（没有差错）的概率为1-p，结果为“正面”（有差错）的概率为p。我们希望通过投掷次硬币来估计p值。p的蒙特卡罗估计器为其中Nheads表示在连续投掷次硬币中出现“正面”的次数。当然，在N次投掷中，可以是0到N中的任何数，但是出现k次“正面”的概率为因此，为了估计Nheads的统计分布并确定p的估计器，我们必须重复进行这个试验多次，比如是M次。23第二十三页，共五十九页，编辑于2023年，星期三Matlab程序M=2000; %numberofexperimentsN=500; %Numberoftosses/experimentH=zeros(1,M);%InitializearrayH_theor=zeros(1,M); %Initializearrayforj=1:MA=rand(1,N);heads=0;fork=1:NifA(k)<=0.2heads=heads+1;endendH(j)=heads;end24第二十四页，共五十九页，编辑于2023年，星期三H_max=max(H);H_min=min(H);r=H_min:H_max;[Nb]=hist(H,r);％fork=H_min:H_maxH_theor(k)=M*nkchoose(N,k)*((0.2)^k)*((0.8)^(N-k));endsubplot(2,1,1)hist(H,r)xlabel('Numberofheads')ylabel('Numberofoccurences')subplot(2,1,2)plot(r,Nb,'ok',r,H_theor(1,H_min:H_max),'k')xlabel('Numberofheads')ylabel('Numberofoccurences')25第二十五页，共五十九页，编辑于2023年，星期三二项式系数functionout=nkchoose(n,k)%Computesn!/k!/(n-k)!a=sum(log(1:n)); %lnofn!b=sum(log(1:k)); %lnofk!c=sum(log(1:(n-k))); %lnof(n-k)!out=round(exp(a-b-c)); %result%Endoffunctionfile.26第二十六页，共五十九页，编辑于2023年，星期三Result发生次数发生次数27第二十七页，共五十九页，编辑于2023年，星期三简单的蒙特卡罗仿真举例假设*在这射机中没有进行脉冲成形。*假设信道是AWGN信道。*信源输出端的数据符号是相互独立和等概率的。*在系统中没有滤波处理，因此不存在码间干扰。非常简单，存在解析解，不需要仿真仍以此为例28第二十八页，共五十九页，编辑于2023年，星期三Findout:当蒙特卡罗仿真应用于更贴近我们要集中研究的主题（数字通信系统）方面的问题时，其性能究竟如何?这里，时延模块是没必要的提醒：几乎所有的仿真都需要这个重要部分29第二十九页，共五十九页，编辑于2023年，星期三利用这一方法，调制器输出端假设的带通信号可表示为其中表示Ac载波幅值，Km是跟调制有关的常数，d[n]是第n个数据符号，θ是参考相位。X(t,n)的复包络只是符号序数的函数，表示为30第三十页，共五十九页，编辑于2023年，星期三Letθ=0对如图9-6所示的系统，为了确定和画出Eb/N0以为函数的BER，保持Eb恒定，并让噪声功率在感兴趣的范围内递增，这就要求对图9-6中的噪声发生器输出端的噪声功率进行校准。9-249-2531第三十一页，共五十九页，编辑于2023年，星期三由第7章知，噪声方差和噪声功率谱密度（PSD）的关系为信噪比SNR定义为Eb/N0，其中fs为采样频率。因此有如果能量Eb和采样频率都归一化为1，则有用于确定后面仿真中的噪声标准偏差。32第三十二页，共五十九页，编辑于2023年，星期三例9-3（二进制相移键控，BPSK）为了产生二进制PSK信号的同相和正交信号空间分量，在式（9-24）和式（9-25）中令Ac=1和Km=π，由此得到

同时，

对于AWGN条件下，等概率、等能量的信号，接收机的阈值总是零

33第三十三页，共五十九页，编辑于2023年，星期三clearallsnrdB_min=0;snrdB_max=10; %SNR(indB)limitssnrdB=snrdB_min:1:snrdB_max;Nsymbols=input('Enternumberofsymbols>');snr=10.^(snrdB/10); %convertfromdBh=waitbar(0,'SNRIteration');len_snr=length(snrdB);34第三十四页，共五十九页，编辑于2023年，星期三forj=1:len_snr %incrementSNRwaitbar(j/len_snr)sigma=sqrt(1/(2*snr(j))); %noisestandarddeviationerror_count=0;

fork=1:Nsymbols %simulationloopbeginsd=round(rand(1)); %dataifd==0x_d=1; %directtransmitteroutput x_q=0; %quadraturetransmitteroutput elsex_d=0; %directtransmitteroutputx_q=1; %quadraturetransmitteroutputendn_d=sigma*randn(1); %directnoisecomponentn_q=sigma*randn(1); %quadraturenoisecomponenty_d=x_d+n_d; %directreceiverinputy_q=x_q+n_q; %quadraturereceiverinputify_d>y_q%testcondition 35第三十五页，共五十九页，编辑于2023年，星期三

d_est=0; %conditionaldataestimateelsed_est=1; %conditionaldataestimateendif(d_est~=d) error_count=error_count+1; %errorcounterendend %simulationloopendserrors(j)=error_count; %storeerrorcountforplotendclose(h)ber_sim=errors/Nsymbols; %BERestimateber_theor=q(sqrt(snr)); %theoreticalBERsemilogy(snrdB,ber_theor,snrdB,ber_sim,'o')axis([snrdB_minsnrdB_max0.00011])xlabel('SNRindB')ylabel('BER')legend('Theoretical','Simulation')36第三十六页，共五十九页，编辑于2023年，星期三对于每一个SNR值，发送Nsymbols=1000个符号，运行结果如图9-8所示。注意，当SNR增加时，BER估计器的可靠度会变差，这是由于差错发生次数减少的缘故。这一现象表明，仿真所发送的符号数可能与SNR有关，或者是连续运行仿真程序，直到对每一个SNR值都记录到相同数目的差错为止。37第三十七页，共五十九页，编辑于2023年，星期三例9-4（二进制频移键控，BFSK）为了产生二进制FSK信号空间的同相和正交分量，令式（9-24）和式（9-25）中的km=π/2，从而得到（9-33）类似地有，（9-34）38第三十八页，共五十九页，编辑于2023年，星期三信号空间表示如图9-9所示，其中1和2是信号空间的基函数。（根据第4章的内容，对于二维空间，可以认为基函数定义了低通复包络信号的同相和正交分量。）因为二进制FSK信号空间是二维的，仿真必须产生信号和噪声的同相和正交分量。判决区域也在图9-9中表示出来了39第三十九页，共五十九页，编辑于2023年，星期三snrdB_min=-3;snrdB_max=8; %SNR(indB)limitssnrdB=snrdB_min:1:snrdB_max;Nsymbols=input('Enternumberofsymbols>');snr=10.^(snrdB/10); %convertfromdBh=waitbar(0,'SNRIteration');len_snr=length(snrdB);Matlabprogram40第四十页，共五十九页，编辑于2023年，星期三forj=1:len_snr %incrementSNRwaitbar(j/len_snr)sigma=sqrt(1/(2*snr(j))); %noisestandarddeviation

error_count=0;fork=1:Nsymbols %simulationloopbeginsd=round(rand(1)); %datax_d=2*d-1; %transmitteroutputn_d=sigma*randn(1); %noisey_d=x_d+n_d; %receiverinputify_d>0 %testconditiond_est=1; %conditionaldataestimateelsed_est=0; %conditionaldataestimateendif(d_est~=d)error_count=error_count+1; %errorcounterendend 41第四十一页，共五十九页，编辑于2023年，星期三%simulationloopendserrors(j)=error_count; %storeerrorcountforplotendclose(h)ber_sim=errors/Nsymbols; %BERestimateber_theor=q(sqrt(2*snr)); %theoreticalBERsemilogy(snrdB,ber_theor,snrdB,ber_sim,'o')axis([snrdB_minsnrdB_max0.00011])xlabel('SNRindB')ylabel('BER')legend('Theoretical','Simulation')42第四十二页，共五十九页，编辑于2023年，星期三9.3蒙特卡罗积分AWGN信道，通过采样积分清除（integrate-and-dump）检测器的输出所得到的充分统计量V，是一个高斯随机变量，其均值由数据符号决定，方差由信道噪声决定。在d[n]=0和d[n]=1的条件下，条件概率密度函数如图9-11所示，其中KT是接收机阈值。在d[n]=1的条件下的条件差错概率为系统的差错概率估计43第四十三页，共五十九页，编辑于2023年，星期三蒙特卡罗积分基本概念假设要求（9-38）其中g（x）是一个在积分区域内有界的函数。由基本的概率论，可知函数的总体均值为（9-39）其中是fx(x)随机变量的概率密度函数。在区间（0，1）上满足fx(x)=1，而在其他地方为零，则。因此，如果U是在区间上均匀分布的随机变量，那么就有（9-40）利用相对频率的观点可以得到（9-41）因此，我们对被积函数进行仿真，再在区间上对它进行次采样，采样点的平均值可用来估计积分值。系统的蒙特卡罗仿真基本上也是这个思路，因为通常无法获得充分统计量在差错区域内的解析表达式，所以采用系统仿真来产生统计量的采样。44第四十四页，共五十九页，编辑于2023年，星期三如果无法求式（9-41）中的极限（实际应用中往往如此），则得到积分的一个近似。将此近似记作得蒙特卡罗估计器为（9-42）总之，通过对函数在N个均匀分布的采样点上求值再取平均，就可实现积分的估计器。这种方法适用于任何常义积分。通过简单的变量代换，可用蒙特卡罗方法来估算任意区间上的常义积分。例如，积分（9-43）可以通过变量代换变成标准形式并求得（9-44）45第四十五页，共五十九页，编辑于2023年，星期三例9-5为了用蒙特卡罗方法估计π值，只需找到一个定积分值为π的函数。很容易就可想到如下积分因此，我们就可以利用式（9-42）定义的算法，求解积分值，并将所得结果乘以4。显然我们所能做到的顶多只是选用一个很大但有限的，据此获得的不是的准确值，而是其近似值。因此，这个积分值，也就是的估计值，是一个随机变量。如图9-12所示为五次估计值的结果，每一次都基于500个试验。(9-42)46第四十六页，共五十九页，编辑于2023年，星期三Matlab程序M=5; %NumberofexperimentsN=500; %Trialsperexperimentu=rand(N,M); %Generaterandomnumbersuu=1./(1+u.*u);%Definefunctiondata=zeros(N,M); %Initializearray%Thefollowingfourlinesofcodedetermine%Mestimatesasafunctionofj,0<j<=N.data(1,:)=4*uu(1,:);forj=2:Ndata(j,:)=4*sum(uu(1:j,:))/j;endest=data(N,:) %Mestimatesofpiest1=sum(est)/M %Averageestimateplot(data,'k') %Plotresultsxlabel('NumberofTrials')ylabel('Estimateofpi')47第四十七页，共五十九页，编辑于2023年，星期三48第四十八页，共五十九页，编辑于2023年，星期三收敛性假设要估计积分值I，同时取得了N个随机观测或采样值X。I的估计器为（9-48）假设N个观测为独立同分布（IID）。样本的算术平均为（9-49）因此估计值是无偏的。由于已经假设观测的独立性，样本方差为（9-50）该式表明积分估计器是一致的。49第四十九页，共五十九页，编辑于2023年，星期三假设，样本的方差，记作，表达为（9-51）因此，给定被积函数g(u)，可以确定要达到某一给定误差方差所需要的N值。因为估计器是一致的，只要N充分得大，就可以得到精确的积分估计。同时可知如果样本g（ui）具有小方差，也可以得到对I的精确估计。对已给定的N，如果g(u)在积分区域内近似恒定（平滑），积分的蒙特卡罗估计就会非常精确。即使，只要g(u)在积分区域内恒定，积分的估计也是精确的。50第五十页，共五十九页，编辑于2023年，星期三置信区间估计器的性能常常用置信区间来表示。置信区间表明了估计值以概率(1-α)落在数值的给定范围内(±βбi).因此，其中有两个关键参数：由α确定的概率和由β确定的范围。置信区间定义为如下表达式51第五十一页，共五十九页，编辑于2023年，星期三置信区间。

我们把区间称为置信区间。52第五十二页，共五十九页，编辑于2023年，星期三现在来确定参数β。式（9-52）可写成误差的形式（9-53）其中，这里由式（9-51）确定。假设误差是高斯随机变量，则的概率密度函数近似为

因为积分估计的无偏性，是零均值随机变量，所以有（9-54）通过变量代换可得（9-55）其中Q()代表高斯Q-函数。53第五十三页，共五十九页，编辑于2023年，星期三从图9-13可得（9-56）从而（9-57）I的估计落在区间的概率是1-α，其中бx由式（9-51）决定。表达式的数量大小决定了上下置信界。必须确定бx

。54第五十四页，共五十九页，编辑于2023年，星期三例9-6为说明前面的概念，考虑以下积分（9-58）利用数值积分的方法，（如MATLAB函数quad），求得I值为（9-59）用同样的方法得（9-60）将式（9-59）和式（9-60）代入式（9-51）可得（9-61）所以积分估计的标准偏差为（9-62）上下置信界为