1、Matlab 通信系统仿真实验实验一 熟悉基本的 Matlab 仿真环境一、实验目的1、熟悉 Matlab 仿真环境，编制简单的 matlab 程序，熟悉基本的调试技巧等。认为学生已经掌握 Matlab 的基本语法和基本操作。2、熟悉基本的 Matlab 中通信仿真工具，相关的函数和命令等的基本使用，包括基本的通信模块相关命令函数，plot 相关的命令函数3、计算机通信仿真的基本的技术和方法二、知识要点1、 Matlab 概述Matlab 是由美国的 MathWorks 公司推出的一种科学计算和工程仿真软件。Matlab 将高性能的科学计算、结果可视化和编程集中在一个易于操作的环境中，提供了大
2、量的内置函数，具有强大的矩阵计算和绘图功能，适用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。目前，在世界范围内被科研工作者、工程技术人员和院校师生广泛采用。2、 Matlab 中的通信仿真工具实现基本的 Matlab 通信仿真，有两种基本的途径：第一种，用 matlab 的基本运算和操作实现基本的通信功能模块，当然前提是对这些基本的通信功能模块的概念和原理非常的清晰。另一种途径是，利用 Matlab 中提供的专业通信工具箱中的函数实现。前提是对这些函数功能非常明确，并熟悉其使用的算法和调用的方法，尤其是参数的理解和设置。Matlab工具箱中包括 100 多个 Matalb 函数
3、可用于通信算法的开发、系统分析及设计。通信工具箱能完成如下任务：1) 信源编码及量化2) 高斯白噪声信道模型3) 差错控制编码4) 调制和解调5) 发送和接收滤波器6) 基带和调制信道模型7) 同步，包括模拟和数字锁相环8) 多址接入，包括 CDMA,FDMA,TDMA.9) 分析结果和比较系统误码率的图形用户界面10) 用于通信信号可视化图形分析和绘制，包括眼图，星座表等。11) 新增的信道可视化工具用于进行时变信道的可视化和开发。3、 Matlab 中的绘图功能Matlab 为用户提供了结果可视化功能，只要在命令窗口输入相应的命令，结果就会用图形直接表示出来。Matlab 提供的常用绘图类
4、函数列表如下：plot 绘制二维线性图形subplot 绘制子图figure 创建一个图的窗口legend 图的注释title 图的标题 xlabel 横轴标注ylabel 纵轴标注 grid 图上加网格hold 保持当前图形 clf 清除图形及属性mesh 三维网线图 plot3 绘制三维线性图形surf 三维表面图 stem 绘制脉冲图errorbar 误差条形图4、 基本的计算机通信仿真的技术和方法蒙特卡罗仿真建立在机率游戏的基础上，因此，以赌城蒙特卡罗命名。其含义是利用蒙特卡罗方法估计系统参数（例如，误比特率）的仿真。几个相关的概念和说明： 蒙特卡罗估计：通过随机试验估计参数值。 相对
5、频率：进行大量的随机试验，试验次数为 N，以 NA 表示事件 A 发生的次数。将事件 A 发生的概率近似为相对频率，其定义为 NA/ N 。 在相对频率的意义下，事件 A 的概率可以通过重复无限多次随机试验来求得，即：NA/ N 的值就是 P(A)的估计器。 由于实验是随机的，N 是有限的实验次数，NA 是随机变量，故估计器也是随机变量。该变量的统计形式决定了其精准度三、实验内容及步骤1、 使用下面的例子熟悉描述绘图的基本步骤：【matlab 程序】：x=-5:0.1:5;y1=sin(x);y2=cos(x) ;%准备绘图数据figure(1) ;%打开图形窗口subplot(2,1,1)
6、;%确定第一幅图绘图窗口plot(x,y1) ;%以 x，y1 绘图title(plot(x,y1);%为第一幅图取名为 plot(x,y1)grid on;%为第一幅绘制网格线subplot(2,1,2);%确定第二幅图绘图窗口plot(x,y2);%以 x，y2 绘图xlabel(time);%第二副图横坐标名为timeylabel(y);%第二幅图纵坐标名为yfigure(2) ;%打开图形窗口 2subplot(2,1,1) ;%图形窗口 2 中第一幅图的绘图窗口stem(x,y1, r) ;%绘制红色的脉冲图subplot(2,1,2);%图形窗口 2 中第二幅图的绘图窗口error
7、bar(x,y1, g);%绘制绿色的误差条形图【程序运行结果】 ；2、 运用下面的实例熟悉三维绘图步骤和相关函数【matlab 程序】：%空间曲线绘制x=-2 :0.2 :2 ;y1=sin(x) ;y2=cos(x) ;plot3(y1,y2,x) ;grid on ;%空间曲面绘制x=-2 :0.2 :2 ;y=x ;X,Y=meshgrid(x,y) ;%生成 x-y 坐标“格点”矩阵 Z=2.*exp(-X.2-Y.2) ;subplot(2,2,1);surf(Z);%绘制曲面shading flat;%把曲面上的小格平滑掉subplot(2,2,2);mesh(Z);%绘制网格曲
8、面subplot(2,2,3);meshc(Z);%等高线投影到平面上subplot(2,2,4);surfl(Z);view(20,0);%变换立体图视角【程序运行结果】：实验 2 常用信号源和噪声源的仿真一、实验目的1、 针对教材第二章的教学内容：随机变量的产生仿真、随机过程的功率谱分析仿真、随机过程通过线性系统的输出信号仿真、带通过程的低通等效仿真等，进行相关内容的熟悉，实践和知识的应用。2、 熟悉上述仿真中的重要的 matlab 命令或函数。3、 理解信源和噪声源的随机性及其数学模型，进行不同性质随机过程（例如，包括低通随机过程、带限低通随机过程、带通随机过程）的分析和仿真（主要包括产
9、生的方法，随机过程的数字统计特征，包括：均值，方差矩阵，自相关函数，功率谱等） 。二、知识要点1、 一种产生高斯分布随机变量的方法（步骤）：产生（0,1）内均匀分布的随机变量 A利用 产生瑞利分布随机变量 R产生（0,1）内均匀分布的随机变量 B，则： 是服从均匀分布的随机变量利用产生一对高斯分布的随机变量。注意：参数 是 C 和 D 的方差2、 一个平稳随机过程 的频域特性用功率谱表征，功率谱和随机过程自相关函数互为傅里叶变换和反变换，即：对于离散时间序列 ，计算该序列的的自相关，其定义如下：同样，求自相关函数的离散傅里叶变换（DFT）可得到 的功率谱。DFT 的定义为：可以采用快速傅里叶变
10、换算法高效地计算出结果。3、 随机过程的线性滤波假设一个平稳随机过程 通过某一线性时不变滤波器，该滤波器在时域用它的冲击响应来表征，而在频域则用它的频率响应：来表征，则该线性滤波器的输出随机过程为：的均值是：其中， 是频率响应在 在 的值。的自相关函数：和 的功率谱之间的关系：4、 离散时间情况下，随机序列通过线性滤波器的输出仿真：假设一平稳随机过程 被采样，其样本通过某一离散线性滤波器，该滤波器的脉冲响应为 。该线性滤波器的输出为（输入和冲击响应的卷积和）：其中， 是输入随机过程 的离散时间值， 是离散时间滤波器的输出。输出过程均值：其中， 是频率响应在 在 的值，输出过程的自相关函数：输出
11、过程的功率谱：其中，输入过程的功率谱为：故：5、 低通和带通过程若过程的功率谱在 附近是大的，而在高频域很小（接近于 0） ，这个随机过程就称之为低通随机过程。即：低通随机过程的功率大部分集中在低频域。若过程的功率谱，则称这个低通随机过程是带限的，参数 B 称为该随机过程的带宽。6、 带通过程的低通等效：若随机过程的功率谱在某中心频率 临近的一个频带内是大的，而在该频带以外相对很小，就称这个随机过程是带通过程。若通带 ，则这个过程是一个窄带过程。通信过程中，携带信息的信号通常是一个低通随机过程。用它去调制某一载波，而在一个带通（窄带）通信信道上传输，因此，已调信号时一个带通随机过程。带通随机过
12、程的低通表示：和确定信号一样，带通随机过程 能表示为：三、实验内容、步骤及实验结果1、 理解教材第 46 页阐述的一个产生高斯分布随机变量的方法，结合知识要点 1，调试并运行教材第 47 页的程序 gngauss 程序，输出运行的结果。【matlab 程序】：function gsrv1,gsrv2 = gngauss( m,sgma )%UNTITLED4 Summary of this function goes here% Detailed explanation goes hereif nargin = 0m=0;sgma=1;elseif nargin = 1 sgma=m;m=0;
13、end;u=rand;z=sgma*(sqrt(2*log(1/(1-u);u=rand;gsrv1=m+z*cos(2*pi*u);gsrv2=m+z*sin(2*pi*u);endfunction X = gaus_mar( X0,rho,N )%UNTITLED5 Summary of this function goes here% Detailed explanation goes herefor i=1:2:NWs(i) Ws(i+1)=gngauss;endX(1)=rho*X0+Ws(1);for i=2:NX(i)=rho*X(i-1)+Ws(i);endendrho=0.9
14、5;X0=0;N=1000;X=gaus_mar(X0,rho,N);plot(X);【程序运行结果】：【实验结果分析】：产生了高斯-马尔科夫系列，说明也产生了高斯分布随机变量，进而得到了验证。【说明】：调用了 gngauss 函数和 gaus_mar 函数，gngauss 函数是产生高斯分布随机变量，gaus_mar 产生高斯-马尔科夫系列。2、 调试并运行教材第 55 页解说题 2.4 的程序，输出运行结果（补充程序中有关结果输出的 plot 程序语句） 。观察不同次数运行结果的差异，理解程序中采用 10 次实现求平均自相关的原因。在调试程序的过程中，关注 fftshift，fft ， r
15、and 等 matlab 函数的功能及应用，结合自相关函数估计的方法，理解子程序 Rx_est 的细节内容。function Rx=Rx_est(X,M)% Rx=Rx_est(X,M)% RX_EST estimates the autocorrelation of the sequence of random% variables given in X. Only Rx(0), Rx(1), . , Rx(M) are computed.% Note that Rx(m) actually means Rx(m-1). N=length(X);Rx=zeros(1,M+1); for m=
16、1:M+1, for n=1:N-m+1,Rx(m)=Rx(m)+X(n)*X(n+m-1); end;Rx(m)=Rx(m)/(N-m+1); end;N=1000;M=50;Rx_av=zeros(1,M+1);Sx_av=zeros(1,M+1);for j=1:10;X=rand(1,N)-1/2;Rx=Rx_est(X,M);Sx=fftshift(abs(fft(Rx);Rx_av=Rx_av+Rx;Sx_av=Sx_av+Sx;end;Rx_av=Rx_av/10;Sx_av=Sx_av/10;figure(1) ;plot(Rx_av);%显示自相关函数f=-0.5:0.02:
17、0.5;figure(2) ;plot(f,Sx_av);axis(-0.5,0.5,0,0.14);%给横坐标定范围为 -0.5,0.5、纵坐标为 0，0.14（都是闭区间）【程序运行结果】:自相关函数功率谱【实验结果分析】：获得的结果与理论结果比较接近，本次的实验比较成功。【说明】：调用 Rx_est(X,M)，该函数是根据 X 序列来产生 1 个自相关函数。3、 调试并运行教材第 57 页解说题 2.5 的程序，输出运行结果（补充 plot 相关的语句） 。程序附下：【matlab 程序】% MATLAB script for Illustrative Problem 5, Chapte
18、r 2.echo on% first partSx1=ones(1,32);Rx1=ifft(Sx1,32); % second partSx2=ones(1,16),zeros(1,224),ones(1,16);Rx2=ifft(Sx2,256);% plotting commands followsubplot(2,1,1)t=1:1:32;plot(t,fftshift(Rx1);subplot(2,1,2)t=1:1:256;plot(t,fftshift(Rx2);【程序运行结果】:【实验结果分析】:对比两图，由于 Sx1 由 32 个样本表示，Sx2 由 256 个样本表示。因此
19、得到的 Rx1 是一个粗略的显示，相比之下 Rx2 中的那些中间值在 Rx1 中都为零。4、 调试并运行教材第 60 页解说题 2.6，第 61 页解说题 2.7 的程序，输出运行结果（补充程序中有关结果输出的 plot 程序语句） 。 （知识要点 3）【matlab 程序】：题 2.6 程序：deIta=0.01;F_min=-2;F_max=2;f=F_min:deIta:F_max;Sx=ones(1,length(f);H=1./(1+(2*pi*f).2);Sy=Sx.*H.2;plot(f,Sy);题 2.7 程序：N=256;deItaf=0.1;f=0:deItaf:(N/2)
20、*deItaf,-(N/2-1)*deItaf:deItaf:-deItaf;Sy=1./(1+(2*pi*f).2);Ry=ifft(Sy);plot(fftshift(real(Ry);【程序运行结果】：题 2.6 功率谱题 2.7 自相关函数【实验结果分析】：功率谱 Sx=1(对于全部 f)，而 Sy 只有在 f=0 的时候才为 1，实际上信号是经过了滤波器，输出了滤波信号。从图也可以分析出这一点。5、 调试并运行教材第 63 页解说题 2.8 的程序，输出运行结果（补充程序中有关结果输出的 plot 程序语句） 。 （知识要点 4）【matlab 程序】：delta_w=2*pi/10
21、0;w=-pi:delta_w:pi;Sy=1./(1.9025-1.9*cos(w);plot(w,Sy);【程序运行结果】：题 2.8 功率谱图【实验结果分析】：该实验进行的是白噪音通过一个线性滤波器，获取输出的信号，从实验结果看，白噪声的功率谱为 1（对于所有 f） ，经过滤波器后，获得的功率谱只有在 f=0时才为 1，而相邻比较远的频率区域功率谱为 0，对滤波器进行了验证。6、 调试并运行教材第 65 页解说题 2.9 的程序，输出运行结果（补充程序中有关结果输出的 plot 程序语句）并比较线性滤波前后信号自相关函数与功率谱的差异并作出分析。【matlab 程序】：N=1000;M=
22、50;Rxav=zeros(1,M+1);Ryav=zeros(1,M+1);Sxav=zeros(1,M+1);Syav=zeros(1,M+1);for i=1:10X=rand(1,N)-0.5;Y(1)=0;for n=2:NY(n)=0.9*Y(n-1)+X(n);endRx=Rx_est(X,M);Ry=Rx_est(Y,M);Sx=fftshift(abs(fft(Rx);Sy=fftshift(abs(fft(Ry);Rxav=Rxav+Rx;Ryav=Ryav+Ry;Sxav=Sxav+Sx;Syav=Syav+Sy;endRxav=Rxav/10;Ryav=Ryav/10;
23、Sxav=Sxav/10;Syav=Syav/10;figure(1)plot(Rxav);axis(0,50,-0.01,0.09);figure(2)plot(Ryav);axis(0,50,0,0.45);figure(3)f=-0.5:0.02:0.5;plot(f,Sxav);axis(-0.5,0.5,0,0.14);figure(4)plot(f,Syav);axis(-0.5,0.5,0,6);【程序运行结果】：自相关函数 Rx自相关函数 Ry功率谱 Sx功率谱 Sy【实验结果分析】:这个实验中，从功率谱来看，Sx 的功率谱在 f 范围内为 1，经过了低通滤波器后产出低通过程，
24、获得的 Sy 的功率谱在 f=0 的时候其值最大，离 f=0 处越远的频率对于的功率谱越小，验证了低通过程。从自相关函数来看，Rx 为一个冲击响应，而 Ry 表现出一条单调递减的曲线，说明输出的信号的自相关性很大，也可以验证低通过程。7、 调试并运行教材第 69 页解说题 2.10 的程序，输出运行结果（补充程序中有关结果输出的 plot 程序语句） 。 （知识要点 6）【matlab 程序】：N=1000;for i=1:2:NX1(i),X1(i+1)=gngauss;X2(i),X2(i+1)=gngauss;end;A=1,-0.9;B=1;Xc=filter(B,A,X1);Xs=f
25、ilter(B,A,X2);fc=1000/pi;for i=1:N;band_pass_process(i)=Xc(i)*cos(2*pi*fc*i)-Xs(i)*sin(2*pi*fc*i);end;M=50;bpp_autocorr=Rx_est(band_pass_process,M);bpp_spectrum=fftshift(abs(fft(bpp_autocorr);f=-0.5:0.02:0.5;plot(f,bpp_spectrum);【程序运行结果】：【实验结果分析】：从实验结果看，在中心频率f 0 附近的频带内功率谱达到最大，而在其他频带内越远离f 0 功率谱越来越小。因
26、此验证了带通过程。8、 利用 matlab 函数 rand(1,N)在区间【0,1】上产生 1000 个均匀随机数集合，画出这个序列直方图和概率分布函数。 （教材第 70 页 2.1）【matlab 程序】N=1000;x=rand(1,N);xmin=min(x);xmax=max(x);t=linspace(xmin,xmax,10); y=hist(x,t);figure(1);bar(t,y);y=y/length(x);F=zeros(1,10);F(1)=y(1);for i=2:1:10F(i)=F(i-1)+y(i);endfigure(2);plot(F);【程序运行结果】:
27、直方图概率分布函数9、 利用 matlab 函数 rand(1,N)在区间【-1/2,1/2】上产生 1000 个均匀随机数集合，画出这个序列直方图和概率分布函数。 （教材第 70 页 2.2）【matlab 程序】N=1000;x=rand(1,N)-0.5;xmin=min(x);xmax=max(x);t=linspace(xmin,xmax,10); y=hist(x,t);figure(1);bar(t,y);y=y/length(x);F=zeros(1,10);F(1)=y(1);for i=2:1:10F(i)=F(i-1)+y(i);endfigure(2);plot(F);
28、【程序运行结果】:直方图概率分布函数10、 利用 matlab 函数 rand(1,N)在区间【-2,2】上产生 1000 个均匀随机数集合，画出这个序列直方图和概率分布函数。 （教材第 71 页 2.3）【matlab 程序】N=1000;x=4*(rand(1,N)-0.5);xmin=min(x);xmax=max(x);t=linspace(xmin,xmax,10); y=hist(x,t);figure(1);bar(t,y);y=y/length(x);F=zeros(1,10);F(1)=y(1);for i=2:1:10F(i)=F(i-1)+y(i);endfigure(2
29、);plot(F);【程序运行结果】:直方图概率分布函数11、 利用教材 2.2 节叙述的方法，产生 1000 个具有零均值和单位方差的高斯随机数集合，画出这个序列的直方图和概率分布函数。 （教材第 71 页 2.5）【matlab 程序】：N=1000;for i=1:Nx(i),x(i+1)=gngauss;end;xmin=min(x);xmax=max(x);t=linspace(xmin,xmax,10); y=hist(x,t);figure(1);bar(t,y);y=y/length(x);F=zeros(1,10);F(1)=y(1);for i=2:1:10F(i)=F(i
30、-1)+y(i);endfigure(2);plot(F);【程序运行结果】：直方图概率分布函数12、 利用 matlab 函数 randn(1,N)产生 1000 个具有零均值和单位方差的高斯随机数的集合，画出这个序列的直方图和概率分布函数。输出结果与 11、的结果作比较。 （教材第71 页 2.6）【matlab 程序】：N=1000;x=randn(1,N);xmin=min(x);xmax=max(x);t=linspace(xmin,xmax,10); y=hist(x,t);figure(1);bar(t,y);y=y/length(x);F=zeros(1,10);F(1)=y(
31、1);for i=2:1:10F(i)=F(i-1)+y(i);endfigure(2);plot(F);【程序运行结果】：直方图概率分布图【实验结果分析】：输出结果与 11、的结果作比较，从这两个实验可以看到，图像和直方图的大致轮廓相似，都产生高斯随机数，不过仿真会出现误差，因此两者的图形不会完全相同。13、 当一带限随机过程的功率谱是：由该功率谱计算随机过程的自相关函数。 （教材第 72 页 2.10）a) 在带内 中用 N 个样本代表功率谱 （假设 B=1）b) 比较两种情况下不同的结果。 （只考虑带内的功率谱样本点和考虑带内的功率谱样本点和带外的零点）【solution】：功率谱表达式
32、为只考虑带内的功率谱样本点，N=33时，功率谱表示为：Sx1=0:0.0625:1, 0.9375:-0.0625: 0;考虑带内的功率谱样本点和带外的零点，N=255 时，功率谱表示为Sx2=0:0.0625:1,zeros(1,222),0.9375:-0.0625:0;【matlab 程序】:Sx1=0:0.0625:1, 0.9375:-0.0625: 0;Rx1=ifft(Sx1,33); % second partSx2=0:0.0625:1,zeros(1,222),0.9375:-0.0625:0;Rx2=ifft(Sx2,255);figure(1);t1=1:1:33;pl
33、ot(t1,fftshift(Rx1);figure(2);t2=1:1:255;plot(t2,fftshift(Rx2);【程序运行结果】N=33 的自相关函数N=255 的自相关函数【实验结果分析】:对比两图，由于 Sx1 由 33 个样本表示，Sx2 由 255 个样本表示，增加样本数以包括|f|B 的样本，来得到 Rx2 的中间值。因此得到的 Rx1 是一个粗略的显示，相比之下 Rx2 中的那些中间值在 Rx1 中都为零。14、 假设由一个功率谱为： （对全部 f）的白色随机过程 激励某一线性滤波器，假设线性滤波器的冲激响应为：求该滤波器输出的功率谱。 （教材第 72 页 2.11）
34、【solution】：求出该滤波器的频率响应为所以就可求得功率谱。【matlab 程序】：deIta=0.01;F_min=-2;F_max=2;f=F_min:deIta:F_max;Sx=ones(1,length(f);H=1./(9+(2*pi*f).2);Sy=Sx.*H.2;plot(f,Sy);【程序运行结果】：功率谱图15、 假设一样本为 的白色随机过程通过一线性滤波器，该滤波器的脉冲响应 h(n)是：求输出过程 的功率谱。 （教材第 72 页 2.13）【solution】：求得：因此有：功率谱也为：【matlab 程序】:delta_w=2*pi/100;w=-pi:delta_w:pi;Sy=1./(1.64-1.6*cos(w);plot(w,Sy);【程序运行结果】:【实验结果分析】：该实验进行的是白噪音通过一个线性滤波器，获取输出的信号，从实验结果看，白噪声的功率谱为 1（对于所有 f） ，经过滤波器后，获得的功率谱只有在 f=0时才为 1，而相邻比较远的频率区域功率谱为 0，对滤波器进行了验证。16、 产生一个在 内，N=1 000 均匀分布随机数的独立同分布序列 ，这个序列通过一线性滤波器，其脉冲响应为：
