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

通信系统的计算机模拟第十一讲1第1页，课件共78页，创作于2023年2月第十讲回顾将均匀分布随机变量映射成任意pdf逆变换法直方图法舍弃法产生不相关的高斯随机变量均匀变量求和法瑞利随机变量到高斯随机变量的映射极坐标法MATLAB实现确定任意的功率谱密度或自相关函数2第2页，课件共78页，创作于2023年2月7.7PN序列发生器在很多应用中尤其是在同频领域中会用到伪噪声（PN）序列发生器。PN序列常用来近似具有均匀概率密度函数的随机变量。PN序列发生器可以有很多表示形式，最常见也是我们要集中讨论的形式，如图7-21所示。3第3页，课件共78页，创作于2023年2月作用在仿真中，使用PN序列最重要的一个原因是为了建立数据源模型。通过使用PN序列发生器，能用尽可能短的仿真时间，产生具有给定长度的所有可能的比特组合。4第4页，课件共78页，创作于2023年2月组成PN序列发生器由三个基本部分组成：N级移位寄存器、模二加法器和连接向量。这个连接向量具体定义了移位寄存器各级与模二加法器之间的连接。它确定了发生器的性能特点，并由以下多项式定义（7-99）5第5页，课件共78页，创作于2023年2月性质可以证明，PN序列发生器输出序列的最大周期是（7-100）当且仅当多项式是本原多项式时，可以达到这个周期。PN序列发生器输出的自相关函数如图7-22所示，这里我们假设数据值（符号）为±1。若或是L的整数倍，则自相关函数等于1；若，则自相关函数-1/L，并且在L足够大时，趋近于零。可见，具有长周期的PN序列，其自相关函数近似为冲激函数。所以，跟期望的一样，功率谱密度近似为白色的。PN序列有很多有趣的特性，在仿真中我们关心的特性包括以下几个：注：本原多项式：一个N阶多项式g(D)能整除Dk+1的最小整数k是k=2N-16第6页，课件共78页，创作于2023年2月7第7页，课件共78页，创作于2023年2月仿真中我们关心的特性·序列几乎是平衡的。也就是说，在序列的一个周期中，1的个数比0的个数只多出一个。在适当的位置上额外添加一个0，则序列平衡。·所有可能的比特组合都会出现一个周期中，这里全0序列除外，但包括全1序列。（注意：如果所有的移位寄存器包含二进制0，则序列发生器会陷于全0的死循环。若已有N-1个0，再添加一个0，则序列平衡，得到的是deBruijn序列）。·虽然自相关函数是周期的，但其波形与随机的二进制波形几乎完全相同。8第8页，课件共78页，创作于2023年2月基于N级移位寄存器的PN序列发生器的设计要找到一个N阶的本原多项式。在N较大时，很难找到这样的本原多项式，幸运的是，文献[13]提供了大量本原多项式的表格。9第9页，课件共78页，创作于2023年2月开发步骤10第10页，课件共78页，创作于2023年2月例11第11页，课件共78页，创作于2023年2月pntaps=[0010000001]; %shiftregistertapspninitial=[0000000001]; %initialshiftregisterstatepndata=zeros(1,1023); %initializeoutputvectorsamp_per_sym=5; %samplespersymbolpnregister=pninitial; %initializeshiftregistern=0; %initializecounterkk=0; %setterminatorindicator12第12页，课件共78页，创作于2023年2月whilekk==0n=n+1; %incrementcounterpndata(1,n)=pnregister(1,1); %saveoutputfeedback=rem((pnregister*pntaps'),2);%calculatefeedbackpnregister=[feedback,pnregister(1,1:9)]; %incrementregisterifpnregister==pninitial;kk=1;end %resetterminationEnd13第13页，课件共78页，创作于2023年2月text=['Theperiodis',num2str(n,15),'.'];disp(text) %displayperiodpndata=replicate(pndata,samp_per_sym);%replicatedatakn=n*samp_per_sym; %lengthofoutputvector pndata=2*pndata-1; %makeoutput+/-onea=fft(pndata);b=a.*conj(a); %PSDofdataRm=real(ifft(b))/kn; %autocorrelationx1=(0:length(Rm)-1)/samp_per_sym;x2=0:100;14第14页，课件共78页，创作于2023年2月%Nowwecanplottheresults.subplot(3,1,1)plot(x1,Rm,'.k');ylabel('R[m]')subplot(3,1,2)stem(x2,Rm(1:101),'.k');ylabel('PartialR[m]')subplot(3,1,3)stem(x2,pndata(1:101),'.k');ylabel('First100outputs')axis([0100-1.51.5]);15第15页，课件共78页，创作于2023年2月注意，该程序作了一个测试，以保证PN序列发生器确实是全周期。（对本例，执行程序后返回的序列周期正好等于N=10时的全周期1023。）该程序还产生了自相关函数图、自相关函数的前101个采样以及发生器输出端的前101个采样。对采样率为1采样/符号时的结果如图7-24所示，若将采样率改为5采样/符号，对应的结果如图7-25所示。16第16页，课件共78页，创作于2023年2月17第17页，课件共78页，创作于2023年2月18第18页，课件共78页，创作于2023年2月7.8信号处理在系统输入是随机的从而输出也是随机的这个前提下，线性系统的几个输入输出关系，我们对以下几个基本结论感兴趣：·系统输入均值和输出均值之间的关系·系统输入方差和输出方差之间的关系·输入-输出互相关·系统输入和系统输出的自相关和功率谱密度之间的关系这些对仿真研究很有用的关系，在系统输入是一个随机过程样本函数的情况下，给我们提供了信号处理的基本工具。接下来的整个分析过程中，假设所研究的系统是线性的和固定不变的，并假设系统的输入具有广义平稳性。19第19页，课件共78页，创作于2023年2月7.8.1输入-输出均值20第20页，课件共78页，创作于2023年2月21第21页，课件共78页，创作于2023年2月7.8.2输入-输出互相关22第22页，课件共78页，创作于2023年2月7.8.3输出自相关函数23第23页，课件共78页，创作于2023年2月24第24页，课件共78页，创作于2023年2月25第25页，课件共78页，创作于2023年2月7.8.4输入-输出方差26第26页，课件共78页，创作于2023年2月第八章后处理将仿真产生的数据处理成有用的形式后处理器通常是图形密集型的视觉显示比数字列表更容易理解，而后者是仿真程序最常见的数据输出。不同系统的误比特率曲线比包含同样信息的数据表格能更快地传达信息后处理程序可能引入（也可能不引入）相当的计算复杂度。有些后处理器只是简单地取得仿真产生的数据，对其进行适当的格式编排后，产生合适的图形输出。产生误比特概率或误符号概率作为函数的曲线，作最少处理来产生图形输出的后处理器的其他例子还包括显示信号波形、眼图和散点图的后处理器。有些后处理子程序会涉及相当多的数据处理，它们中多数要涉及到某种估计，一个简单的例子就是直方图的产生，它是对概率密度函数的估计更复杂的例子有延迟、信码的系统中，从信道误符号率到译码后的误比特率的映射。27第27页，课件共78页，创作于2023年2月基本图形方法系统实例——π/4DQPSK传输28第28页，课件共78页，创作于2023年2月发射信号的相位偏移由d(k)和q(k)的值以及前一符号周期的相位偏移θ（t－1）决定例如001011011129第29页，课件共78页，创作于2023年2月例如0010110111－3π/4,-π,3π/4,0,π/4,……30第30页，课件共78页，创作于2023年2月波形、眼图和散点图如果同相信道信号xd(t)画在平在(D,t)上，正交信号xq(t)画在平面上，就产生了为参数的三维信号。Q,D,t31第31页，课件共78页，创作于2023年2月眼图眼图给出系统性能的一种定性量度。明晰张开的眼孔通常表明好的性能而模糊的眼孔往往表明差的性能眼孔的大小跟符号同步器所要求的精度有关尽管眼图不能对系统性能提供定量的量度，但很难想像一个性能优异的系统会有不明确的眼图。32第32页，课件共78页，创作于2023年2月假设用示波器显示波形，并且示波器在垂直虚线处触发，结果就是如图8-4（b）所示的三段眼图33第33页，课件共78页，创作于2023年2月例8-1产生并显示DQPSK系统中的一些重要信号，仿真这个系统并产生图形输出。34第34页，课件共78页，创作于2023年2月m=200; bits=2*m; %numberofsymbolsandbitssps=10; %samplespersymboliphase=0; %initialphaseorder=5;%filterorderbw=0.2;%normalizedfilterbandwidth%initializevectorsdata=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);%setdirectandquadraturebitstreamsdata=round(rand(1,bits));dd=data(1:2:bits-1);qq=data(2:2:bits);35第35页，课件共78页，创作于2023年2月%mainprogramstheta(1)=iphase;%setinitialphasethetaout(1:sps)=theta(1)*ones(1,sps);fork=2:mifdd(k)==1phi_k=(2*qq(k)-1)*pi/4;elsephi_k=(2*qq(k)-1)*3*pi/4;

end36第36页，课件共78页，创作于2023年2月theta(k)=phi_k+theta(k-1);fori=1:spsj=(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第37页，课件共78页，创作于2023年2月%postprocessorforplottingkk=0;%setexitcounterwhilekk==0%testexitcounterk=menu('pi/4QPSKPlotOptions',...'Unfilteredpi/4QPSKSignalConstellation',...'Unfilteredpi/4QPSKEyeDiagram',...'Filteredpi/4QPSKSignalConstellation',...'Filteredpi/4OQPSKEyeDiagram',...'UnfilteredDirectandQuadratureSignals',...'FilteredDirectandQuadratureSignals',...'ExitProgram');

38第38页，课件共78页，创作于2023年2月1．滤波前的π/4DQPSK信号星座图2．滤波前的π/4DQPSK信号眼图3．滤波后的π/4DQPSK信号星府图4．滤波后的π/4DQPSK信号眼图5．滤波前的同相和正交信号6．滤波后的同相和正交信号7．退出程序（返回MATLAB命令提示符）39第39页，课件共78页，创作于2023年2月ifk==1sigcon(d,q) %plotunfilteredsignalcon.pauseelseifk==2dqeye(d,q,4*sps) %plotunfilteredeyediagrampauseelseifk==3sigcon(df,qf) %plotfilteredsignalcon.pauseelseifk==4dqeye(df,qf,4*sps) %plotfilteredeyediagrampauseelseifk==5numbsym=10; %numberofsymbolsplotteddt=d(1:numbsym*sps); %truncatedvectorqt=q(1:numbsym*sps); %truncateqvectordqplot(dt,qt) %plottruncateddandqsignals40第40页，课件共78页，创作于2023年2月pauseelseifk==6numbsym=10; %numberofsymbolstobeplotteddft=df(1:numbsym*sps); %truncatedftodesiredvalueqft=qf(1:numbsym*sps); %truncateqftodesiredvaluedqplot(dft,qft) %plottruncatedsignalspauseelseifk==7kk=1; %setexitcountertoexitvalueendend41第41页，课件共78页，创作于2023年2月支撑子程序function[]=sigcon(x,y)plot(x,y)axis('square')axis('equal')xlabel('DirectChannel')ylabel('QuadratureChannel')%Endoffunctionfile.42第42页，课件共78页，创作于2023年2月function[]=dqeye(xd,xq,m)lx=length(xd);%samplesindatasegmentkcol=floor(lx/m);%numberofcolumnsxda=[0,xd];xqa=[0,xq];%appendzerosforj=1:kcol%columnindexfori=1:(m+1)%rowindexkk=(j-1)*m+i;%sampleindexy1(i,j)=xda(kk);y2(i,j)=xqa(kk);endend43第43页，课件共78页，创作于2023年2月subplot(211)%directchannelplot(y1);title('D/QEYEDIAGRAM');xlabel('SampleIndex');ylabel('Direct');subplot(212)%quadraturechannelplot(y2);xlabel('SampleIndex');ylabel('Quadratute');subplot(111)%restore44第44页，课件共78页，创作于2023年2月unction[]=dqplot(xd,xq)lx=length(xd);t=0:lx-1;nt=t/(lx-1);nxd=xd(1,1:lx);nxq=xq(1,1:lx);subplot(211)plot(nt,nxd);a=axis;axis([a(1)a(2)1.5*a(3)1.5*a(4)]);title('DirectandQuadratureChannelSignals');xlabel('NormalizedTime');ylabel('Direct');subplot(212)plot(nt,nxq);a=axis;axis([a(1)a(2)1.5*a(3)1.5*a(4)]);xlabel('NormalizedTime');ylabel('Quadratute');subplot(111)45第45页，课件共78页，创作于2023年2月46第46页，课件共78页，创作于2023年2月47第47页，课件共78页，创作于2023年2月48第48页，课件共78页，创作于2023年2月估计-直方图如果可以获得随机过程的一组样本（就像在仿真环境下那样），由这组样本得出的直方图经常可用于估计内在的概率密度函数（pdf）直方图将总共N个样本的数据分组成B个直方（bin）或单元（cell），假定每个直方的宽度为W，其中心记为bi。如果下式成立，给定的样本落在第个直方49第49页，课件共78页，创作于2023年2月Ni表示落在第个直方中的样本数目。Count{N：R}来表示全部N个样本中进入由R所定义的直方中的样本数目50第50页，课件共78页，创作于2023年2月然后画出条状图，每条的高度和Ni成正比，中心在bi为了用作pdf的估计器，我们对直方图进行缩放使其总面积为1，将Ni除以NW即可实现这一点，于是，每条的高度为Ni/NW。代表第个直方的条的面积Ai，等于条的高度乘以宽度。因此，注意到Ai代表第个直方的相对频率注意到每个直方是用一个常数来表示在有限宽度上的pdf。对给定直方中的某个点，pdf的估计是无偏的，然而对给定直方所定义的大多数区域，估计是有偏的51第51页，课件共78页，创作于2023年2月subplot(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,100000);hist(x,50)ylabel('N_i');xlabel('(d)')52第52页，课件共78页，创作于2023年2月53第53页，课件共78页，创作于2023年2月功率密度估计我们感兴趣的波形是随机过程的样本函数，这就使得在给定频率值f1处的PSD是一个随机变量需要最小化谱估计的方差基于快速傅里叶变换（FFT），而且经常用于仿真环境54第54页，课件共78页，创作于2023年2月周期图最简单、最快速也最常用的PSD估计算法是周期图。其定义如下N是数据记录里的样本总数，IN(Kf∆)是用于计算在频率f=Kf∆处PSD估计的数据的N点FFT。周期图的问题在于它是有偏差的对于足够大的N，偏差可以忽略。主要的问题在于大方差55第55页，课件共78页，创作于2023年2月带数据窗的周期图默认的矩形窗每个样本值x[n]乘以w[n]=1,0<n<N矩形窗的影响是将数据样本x[n]与w[n]的傅里叶变换进行卷积从频域看，这种数据窗的旁瓣结构，会导致相当严重的谱泄露。谱泄露会减少估计谱的动态范围或使其失真每个样本值

56第56页，课件共78页，创作于2023年2月使用任意的数据窗U--是数据窗里的能量，理想的数据窗在时域上必须有有限的持续时间，这样就可以用有限的数据记录来精确估计数据的傅里叶变换。此外，对于估计数据记录的傅里叶变换，窗函数一定不能有不良影响。因为时域相乘就是频率卷积，而且只有与冲激函数卷积才不会改变变换前后的函数，所以理想的窗函数是频域冲激函数。然而由于冲激函数的傅里叶变换不是有限伸展的，这些要求互相冲突，于是我们想寻求这样一种数据窗：它在频域处附近有一个窄的主瓣衰减显著。57第57页，课件共78页，创作于2023年2月分段周期图减少周期图方差的常用方法是，将长度为N的采样数据记录分成K段，每一段由M个样本组成，计算每一段的FFT，对结果进行平均，平均过程减少了谱估计的方差，这些段之间可重叠。也可不重叠，如果不重叠K＝MN；否则K>MN。第i段数据记录的周期图由下式给出：58第58页，课件共78页，创作于2023年2月对K个周期图进行平均该估计值是有偏差的，因为数据记录是有限长的。假设K个周期图是互相独立的，则由由于对于固定的我们希望值尽可能大，所以各段往往互相交叠，一般采用50%交叠。当采用了两，值由增加到。如果数据段交叠，那么，个周期图之间不再独立，PSD估计器方差的减小就比式（8-18）预计的要少。数据窗的使用，至少在某种程度上有助于恢复段之间的独立性。59第59页，课件共78页，创作于2023年2月常用窗函数Hanning窗定义如下：60第60页，课件共78页，创作于2023年2月例8-3将独立（白噪声）样本通过有5dB通带纹波的切比雪夫滤波器。试估计滤波器输出端的PSD。61第61页，课件共78页，创作于2023年2月Matlabprogramsettle=100; %ignoretransientfs=1000; %samplingfrequencyN=50000; %sizeofdatarecordf=(0:(N-1))*fs/N;%frequencyscale[b,a]=cheby1(5,5,0.1); %filterNN=N+settle; %allowtransienttodiein=randn(1,NN); %randominputout=filter(b,a,in); %filteroutputout=out((settle+1):NN); %stripoffinitialsampleswindow=hanning(N)'; %setwindowfunction 62第62页，课件共78页，创作于2023年2月winout=out.*window; %windowedfilteroutputfout=abs(fft(winout,N)).^2; %transformandsquaremagU=sum(window.*window); %windowenergyf1out=fout/U; %scalespectrumpsd1=10*log10(abs(f1out)); %logscalesubplot(2,1,1)plot(f(1:5000),psd1(1:5000))grid;axis([0100-7010]);63第63页，课件共78页，创作于2023年2月xlabel('Frequency,Hz')ylabel('PSD')%K=25; %numberofsegmentsM=N/K; %blocksizefK=(0:(M-1))*fs/M; %frequencyscaled=zeros(1,M); %initializevectorpsdk=zeros(1,M); %initializevectorwindow=hanning(M)'; %setwindowfunctionU=sum(window.*window); %windowenergy64第64页，课件共78页，创作于2023年2月fork=1:Kforj=1:Mindex=(k-1)*M+j;d(j)=out(index);enddwin=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([0100-7010]);xlabel('Frequency,Hz')ylabel('PSD')65第65页，课件共78页，创作于2023年2月Result66第66页，课件共78页，创作于2023年2月MATLABPSD估计MATLAB信号处理工具箱包含大量用于PSD估计的子程序其中两个是psd和pwelch67第67页，课件共78页，创作于2023年2月例8-4估计QPSK信号的PSD。假设采用矩形脉冲成型，并以每符号16个点的速率采样QPSK信号的同相和正交分量。68第68页，课件共78页，创作于2023年2月fs=16;x=random_binary(1024,fs)+i*random_binary(1024,fs);fornwin=1:4nwindow=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);ptheory=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('PSDindB')xx=['windowlength=',num2str(nwindow)];xlabel(xx)axis([08-50,10]);grid;end69第69页，课件共78页，创作于2023年2月Result70第70页，课件共78页，创作于2023年2月增益，延迟和信噪比信噪比是一个常用于评估通信系统性能的指标这里介绍的SNR估计方法源于宽带遥测系统中一种测量信道失真误差的方法定义该系统中某点信号上的噪声为该点理想信号和实测信号之间的均方误差（MSE），将理想信号定义为系统输入端的信息承载信号按幅度比例缩放并加入延迟后所得的信号，以此来估计SNR在过去，这种方法的应用包括监视数字脉码调制数据的可靠传输和非线性信道下载波交调比（carrier-to-intermodulationratio）的估计。71第71页，课件共78页，创作于2023年2月实低通信号的理论推导对线性时不变无失真系统，系统中任意点信号是对输入参考信号按幅度比例缩放并加入延迟后所得的信号。因此，我们可以将无失真信号写为：观测信号N（t)代表外部加性噪声，d(t)是系统引入的跟信号相关的内部失真，这种内部失真可能由码间干扰或非线性造成。

72第72页，课件共78页，创作于2023年2月73第73页，课件共78页，创作于2023年2月噪声功率定义A和的理想估计值最小化对于平稳信号，各阶矩独立于时间起点。此外，和的期望等于期望的和

74第74页，课件共78页，创作