Matlab通信系统建模与仿真例题源代码-第三章
Matlab平台下的通信系统建模与仿真
Matlab平台下的通信系统建模与仿真2.1 Matlab平台简介及MATLAB的特点2.1.1 Matlab平台简介MATLAB软件的最初版本是由Cleve Moler博士等组成的名为MathWorks 公司在1992年推向市场的，软件名称MATLAB为矩阵（matrix）和实验室（laboratory）两个英文单词的前三个字母的组合（即矩阵实验室）。
所以，MATLAB非常适合矩阵运算,这也是MATLAB软件的最大特点之一。
2.1.2 MATLAB的特点（1）语言简洁，使用方便，运算符丰富，库函数极其丰富。
利用其丰富的库函数避开繁杂的子程序编程任务。
（2）语法限制不严格，程序设计自由。
例如，在MATLAB里，用户无需对矩阵预定义就可以使用。
（3）MATLAB既具有结构化的控制语句（如for循环、while循环、break 语句和if语句），又有面向对象编程的特性。
（4）程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。
（5）MATLAB的图形功能强大。
数据的可视化非常简单,同时MATLAB 还具有较强的编辑图形界面的能力。
（6）MATLAB的缺点是软件成本高；其次MATLAB的程序不用编译等预处理，也不生成可执行文件，程序为解释执行，和其他高级程序相比，程序的执行速度较慢。
（7）MATLAB由核心部分和各种可选的工具箱两部分组成。
核心部分中包含数百个核心内部函数。
MATLAB的另一重大特色是功能强劲的工具箱。
其工具箱又可分为两类：功能性工具箱和学科性工具箱。
功能性工具箱能用于多种学科，主要用来扩充其符号计算功能、图示建模仿真功能、文字处理功能以及与硬件实时交互功能。
而学科性工具箱是专业性比较强的，如control、toolbox、signal processing toolbox、communication toolbox等。
（8）源程序的开放性。
2.2 Matlab平台下的系统建模与仿真方法在上一章我们介绍了OFDM系统的基本原理，阐明了如果在通信系统中一如OFDM可以提高系统的抗多径干扰能力，从而提高系统的稳定性和可靠性。
第三章_MATLAB与基本PID控制系仿真讲述.讲述
试建立系统的TF模型。
s G( s ) 2 s 2s 10
MATLAB与系统仿真
• 零极点模型
【调用格式】
sys =zpk(z,p,k)
【说明】
z、p、k分别为系统的零点、极点和增益。zpk函数的返回值是一个对象，
【说明】
num和den是离散系统脉冲传递函数的分子和分母多项式系数。 z,p,k是离散系统脉冲传递函数的零点、极点和增益。 Ts是离散系统的采样周期。
MATLAB与系统仿真
数学模型之间的转换
• LTI对象之间的转换 【调用格式】
sys = tf(sys) sys = zpk(sys) % 将sys对象转换为TF模型 % 转换为ZPK模型
• LTI对象属性之间的转换 【调用格式】
[z,p,k] = tf2zp(num, den) [num,den ] = zp2tf(z,p,k) % 将TF对象属性转换为ZPK对象属性 % 将ZPK对象属性转换为TF对象属性
• 连续系统和离散系统之间的转换 【调用格式】
sysd = c2d(sysc,Ts) %将连续系统转换为采样周期为Ts的离散系统
采样周期Ts＝0.1秒,试将其进行离散化处理。
程序： num=10; den=[1,7,10]; ts=0.1; sysc=tf(num,den); sysd=c2d(sysc,ts)
MATLAB与系统仿真
3.2 PID控制概述
比例
rin(k)
+ 微分 _ 积分 + 被控对象
yout(k)
PID控制器是一种线性控制器，它根据给定值rin(t)与实 际输出值yout(t)构成控制偏差：
MATLAB语言与控制系统仿真-参考答案-第3章
3.5 MATLAB 绘图实训3.5.1 实训目的1.学会MA TLAB 绘图的基本知识；2.掌握MA TLAB 子图绘制、图形注释、图形编辑等基本方法；3.学会通过MA TLAB 绘图解决一些实际问题；4.练习二维、三维绘图的多种绘图方式，了解图形的修饰方法；5.学会制作简单的MATLAB 动画。
图3-46 炮弹发射示意图3.5.2 实训内容1. 炮弹发射问题（1）炮弹发射的基础知识炮弹以角度α射出的行程是时间的函数，可以分解为水平距离)(t x 和垂直距离)(t y 。
)cos()(0αtv t x = %水平方向的行程；205.0)sin()(gt tv t y -=α %垂直方向的行程； 其中，0v 是初速度；g 是重力加速度，为9.82m/s ；t 是时间。
（2）炮弹发射程序举例：分析以下程序以及图3-47各个图形的实际意义。
a=pi/4; v0=300; g=9.8;t=0:0.01:50; x=t*v0*cos(a);y=t*v0*sin(a)-0.5*g*t.^2;subplot(221);plot(t,x);grid;title(‘时间-水平位移曲线'); subplot(222);plot(t,y);grid;title(‘时间-垂直位移曲线');subplot(223);plot(x,y);grid;title(‘水平位移-垂直位移曲线'); subplot(224);plot(y,x);grid;title(‘垂直位移-水平位移曲线');图3-4745角发射曲线 （3）编程解决炮弹发射问题①假设在水平地面上以垂直于水平面的角度向上发射炮弹，即发射角90=α，假设初速度分别为[310，290，270]m/s,试绘制时间-垂直位移曲线，编程求取最高射程；绘图要求：◆ 标题设为“炮弹垂直发射问题”； ◆ 在图上通过添加文本的方式表明初速度； ◆ 在x 轴标注“时间”； ◆ 在y 轴上标注“垂直距离”； ◆ 添加网格线； ◆ 将310m/s 的曲线改为线粗为2的红色实线； ◆ 将290m/s 的曲线改为线粗为3的绿色点划线； ◆ 将270m/s 的曲线改为线粗为2的蓝色长点划线；a=pi/2; v1=310; g=9.8;t=0:0.01:50; x1=t*v1*cos(a);y1=t*v1*sin(a)-0.5*g*t.^2;plot(t,y1);grid; title('炮弹垂直发射问题'); xlabel('时间'); ylabel('垂直距离'); hold on; v2=290;x2=t*v2*cos(a);y2=t*v2*sin(a)-0.5*g*t.^2; plot(t,y2); v3=270;x3=t*v3*cos(a);y3=t*v3*sin(a)-0.5*g*t.^2; plot(t,y3);zgsc=[max(y1); max(y2); max(y3)] %三次发射的最高射程 运行结果如下：zgsc =1.0e+003 * 4.9031 4.29083.7194最高射程分别为:4903.1米,4290.8米,3719.4米。
Matlab通信系统建模与仿真例题源代码-第三章
% ch3example1A.mclear;f_p=2400; f_s=5000; R_p=3; R_s=25; % 设计要求指标[n, fn]=buttord(f_p,f_s,R_p,R_s, 's'); % 计算阶数和截止频率Wn=2*pi*fn; % 转换为角频率[b,a]=butter(n, Wn, 's'); % 计算H(s)f=0:100:10000; % 计算频率点和频率范围s=j*2*pi*f; % s=jw=j*2*pi*fH_s=polyval(b,s)./polyval(a,s); % 计算相应频率点处H(s)的值figure(1);subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性axis([0 10000 -40 1]);xlabel('频率Hz');ylabel('幅度dB');subplot(2,1,2); plot(f, angle(H_s)); % 相频特性xlabel('频率Hz');ylabel('相角rad');figure(2); freqs(b,a); % 也可用指令freqs直接画出H(s)的频率响应曲线。
% ch3example1B.mclear;f_p=2400; f_s=5000; R_p=3; R_s=25; % 设计要求指标[n, fn]=ellipord(f_p,f_s,R_p,R_s,'s'); % 计算阶数和截止频率Wn=2*pi*fn; % 转换为角频率[b,a]=ellip(n,R_p,R_s,Wn,'s'); % 计算H(s)f=0:100:10000; % 计算频率点和频率范围s=j*2*pi*f; % s=jw=j*2*pi*fH_s=polyval(b,s)./polyval(a,s); % 计算相应频率点处H(s)的值figure(1);subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性axis([0 10000 -40 1]);xlabel('频率Hz');ylabel('幅度dB');subplot(2,1,2); plot(f, angle(H_s)); % 相频特性xlabel('频率Hz');ylabel('相角rad');figure(2); freqs(b,a); % 也可用指令freqs直接画出H(s)的频率响应曲线。
MATLAB实现通信系统仿真实例
补充内容：模拟调制系统的MATLAB 仿真1。
抽样定理为了用实验的手段对连续信号分析,需要先对信号进行抽样(时间上的离散化），把连续数据转变为离散数据分析。
抽样（时间离散化）是模拟信号数字化的第一步.Nyquist 抽样定律：要无失真地恢复出抽样前的信号，要求抽样频率要大于等于两倍基带信号带宽.抽样定理建立了模拟信号和离散信号之间的关系,在Matlab 中对模拟信号的实验仿真都是通过先抽样，转变成离散信号，然后用该离散信号近似替代原来的模拟信号进行分析的。
【例1】用图形表示DSB 调制波形)4cos()2cos(t t y ππ= 及其包络线。
clf％％计算抽样时间间隔fh=1;%%调制信号带宽（Hz ）fs=100*fh ；％%一般选取的抽样频率要远大于基带信号频率，即抽样时间间隔要尽可能短. ts=1/fs ；％%根据抽样时间间隔进行抽样，并计算出信号和包络t=(0：ts:pi/2)'；%抽样时间间隔要足够小,要满足抽样定理.envelop=cos （2＊pi ＊t)；%%DSB 信号包络y=cos(2＊pi ＊t).＊cos(4*pi ＊t ）;％已调信号%画出已调信号包络线plot （t ，envelop ，'r ：','LineWidth',3）；hold onplot(t ，-envelop,'r:’,’LineWidth'，3）;%画出已调信号波形plot(t ，y,’b'，'LineWidth ’,3）;axis(［0,pi/2,-1,1］)%hold off%xlabel （'t ’);％写出图例【例2】用图形表示DSB 调制波形)6cos()2cos(t t y ππ= 及其包络线。
clf％％计算抽样时间间隔fh=1;％%调制信号带宽(Hz ）fs=100＊fh;%抽样时间间隔要足够小，要满足抽样定理。
第三章 matlab的simulink建模与仿真
3、其它子系统
可配置子系统，代表用 户定义库中的任意模块， 只能在用户定义库中使用。 函数调用子系统。
for循环
3）在一个仿真时间步长内，simulink可以多次进出一 个子系统。 原子子系统：
1）子系统作为一个“实际”的模块，需顺序连续执行。
2）子系统作为整体进行仿真。
3）子系统中的模块在子系统中被排序执行。
建立原子子系统：
1）先建立一空的原子子系统。
2）先建立子系统，再强制转换成原子子系统。
Edit/block parameters
在enabled subsystem
triggered subsystem
enabled and triggered subsystem中。
1）早期simulink版本中，enable和triggered信号需要从 signal&system中调用。
2）simulink后期版本中，在上述模块中含这两个信号。 3）一个系统中不能含多个enable和triggered信号。 4）其它子系统可看成某种形式的条件执行子系统。
3.4创建simulink模型（简单入门）
一、启用simulink并建立系统模型 启动simulink： （1）用命令方式：simulink （2）
二、simulink模块库简介 1、simulink公共模块库 Continuous（连续系统）
连续信号数值积分 输入信号连续时间积分
单步积分延迟，输出为前一输入
动态模型：描述系统动态变化过程
静态模型：平衡状态下系统特性值之间的关系
二、计算机仿真
1、仿真的概念
以相似性原理、控制理论、信息技术及相关领域 的有关知识为基础，以计算机和各种专用物理设备为工 具，借助系统模型对真实系统进行实验研究的一门综合 性技术。 2、仿真分类 实物仿真：建造实体模型 数学模型：将数学语言编制成计算机程序 半实体模型：数学物理仿真
matlab通信仿真实例
matlab通信仿真实例在Matlab中进行通信系统的仿真，可以涉及到多种不同的通信技术和协议，包括调制解调、信道编码、多址接入等。
以下以OFDM系统为例，介绍Matlab 中通信仿真的实例。
OFDM（正交频分复用）是一种常用于现代通信系统中的技术，它将高速数据流分割成多个较低速的子流，并将每个子流分配到不同的子载波上。
优点是能够抵抗多径效应和频率选择性衰落，并提供高数据速率。
首先，我们需要创建一个包含OFDM系统参数的结构体。
例如：ofdmParam.M = 16; % 子载波数量ofdmParam.K = 4; % 用于混合多路复用的用户数量ofdmParam.N = ofdmParam.M * ofdmParam.K; % 总子载波数量ofdmParam.CP = 16; % 循环前缀长度接下来，我们可以生成用于OFDM仿真的数据流。
例如，我们可以使用随机整数生成器生成一系列整数，并将其转换为复数形式的调制符号：data = randi([0, ofdmParam.M-1], 1, ofdmParam.N);dataMod = qammod(data, ofdmParam.M);然后，我们可以创建一个包含OFDM信号的函数。
在OFDM系统中，生成的数据符号将分配到不同的子载波上，然后在时域中通过插入循环前缀进行叠加：function[ofdmSignal] = createOFDMSignal(dataMod, ofdmParam) ofdmSignal = [];for k = 0:ofdmParam.K-1% 提取相应的数据符号，并进行IFFTofdmData =ifft(dataMod(k*ofdmParam.M+1:(k+1)*ofdmParam.M));% 添加循环前缀ofdmDataWithCP = [ofdmData(end-ofdmParam.CP+1:end), ofdmData];% 将OFDM符号添加到OFDM信号中ofdmSignal = [ofdmSignal, ofdmDataWithCP];endend将OFDM信号传输到信道中，我们可以使用加性高斯白噪声（AWGN）信道模型来模拟实际通信环境：EbNo = 10; % 信噪比snr =10*log10(ofdmParam.N*ofdmParam.M/(ofdmParam.N*ofdmParam.M+1 )*(10^(EbNo/10)));ofdmSignalNoisy = awgn(ofdmSignal, snr, 'measured');最后，我们可以对接收到的OFDM信号进行解调和信号恢复。
matlab与通信系统仿真程序及图
实验一2：for i=1:50,j=1:50;A(i,j)=i+j-1;endC=flipud(A),B=fliplr(A),A(1:10,1:10)=zeros;D=A3：clear;a=255*rand(50);disp(a),a(find(a<128))=0;a(find(a>128))=255;w=a;disp(w)实验二：1：（1）clear;n=0:23;N=12;y1=sin(2*4*pi*n/N);y2=sin(2*5*pi*n/N);y3=sin(2*7*pi*n/N);y4=sin(2*10*pi*n/N);plot(n,y1,'+r-',n,y2,'-.',n,y3,'k-',n,y4,'*g-');axis([0 25 -1 2 ]);xlabel('横轴n');ylabel('纵轴y');text(0.8,0.9,'y1');text(19.6,0.9,'y2');text(14,0.9,'y3');text(5,0.9,'y4');title('plot tu xing');legend('y1=sin(2*4*pi*n/N)','y2=sin(2*5*pi*n/N)','y3=sin(2*7*pi*n/N)','y4=sin(2*10*pi*n/N)')（2）clear;n=0:23;N=12;y1=sin(2*4*pi*n/N);y2=sin(2*5*pi*n/N);y3=sin(2*7*pi*n/N);y4=sin(2*10*pi*n/N);stem(n,y1,'+r-');hold on;stem(n,y2,'b-');hold on;stem(n,y3,'k-');hold on;stem(n,y4,'*g-');axis([0,25,-2,4]);xlabel('横轴n');ylabel('纵轴y');text(0.8,1.2,'y1');text(3,1.4,'y2');text(14,1.2,'y3');text(10.9,1.3,'y4');title('n from 0 to 23');legend('y1=sin(2*4*pi*n/N)','y2=sin(2*5*pi*n/N)','y3=sin(2*7*pi*n/N)','y4=sin(2*10*pi*n/N)') 2：clear all;N=input('input N=');T=1;fs=1/T;N_sample=128;dt=1/N_sample;t=0:dt:10*T-dt;n=-N:N;Fn=(-1/2)*sinc(n/2).*exp(-j*n*pi/2);Fn(N+1)=0;ft=zeros(1,length(t));for m=-N:N;ft=ft+ Fn(m+N+1)*exp(j*m*pi*2*fs*t);endplot(t,ft);legend('时域曲线ft=ft+ Fn(m+N+1)*exp(j*m*pi*2*fs*t)',4)3：N_sample=128;dt=1/N_sample;t=0:dt:1-dt;st=[ones(1,N_sample/2),-ones(1,N_sample/2)];subplot(311);plot(t,st);axis([0 1 -2 2]);xlabel('t');ylabel('s(t）'); subplot(312);[f,sf]=T2F(t,st);plot(f,abs(sf));hold on;axis([-10 10 0 1]);xlabel('f');ylabel('|S(f)|');sff=(-1/2)*exp(-j*f*pi/2).*sinc(f/2);plot(f,abs(sff));[t,st]=F2T(f,sf);subplot(313);axis([0 1 -2 2]);xlabel('t');ylabel('恢复的s(t)');plot(t,st)子函数：（1）：function[t,st]=F2T(f,sf)df=f(2)-f(1);Fmx=f(end)-f(1)+df;dt=1/Fmx;N=length(sf);T=dt*N;t=0:dt:T-dt;sff=fftshift(sf);st=Fmx*ifft(sff);（2）：function[f,sf]=T2F(t,st)dt=t(2)-t(1);T=t(end);df=1/T;N=length(st);f=-N/2*df:df:N/2*df-df;sf=fft(st);sf=T/N*fftshift(sf);2图：M=4 M=5M=7 M=10M=4 M=5M=7 M=10N=10 N=10002468100246810实验三：1.（2）clear;p=2;N=5000;u=randn(5,4);c=sqrt(p);u=u*c;power_u=var(u);subplot(211)plot(u(1:20));grid on;ylabel('u(n)');xlabel('n');subplot(212)hist(u,50);grid on;ylabel('histogram of u(n)');2：clear all;N=500;disounce=(sign(rand(1,N)-0.5+eps)+1)/2; n=8;temp1=ones(1,n);temp0=zeros(1,n);new_disounce=[];for i=1:length(disounce)if disounce(i)==0new_disounce=[new_disounce temp0];elsenew_disounce=[new_disounce temp1];endendT=0.01;t=0:T/n:T/n*(length(new_disounce)-1);subplot(211);plot(t,new_disounce)axis([min(t)-0.1,max(t)+0.1,min(new_disounce)-0.1,max(new_disounce)+0.1]) y=new_disounce;Nf=length(t);Y=fft(y,Nf);Pyy=abs(Y).^2/Nf;dt=t(2)-t(1);df=1/dt;f=df*(0:(Nf-1)/2)/Nf;subplot(212);plot(f,10*log10(Pyy(1:((Nf-1)/2+1))));title('功率谱密度');xlabel('频率(HZ)');grid3：clear all;SNRindB1=0:1:12;SNRindB2=0:0.01:12;for i=1:length(SNRindB1)smld_err_prb(i)= exa0308_Fun_singlepe(SNRindB1(i));endsemilogy(SNRindB1,smld_err_prb,'r*');for i=1:length(SNRindB2)SNR=10^(SNRindB2(i)/10);theo_err_prb(i)=(1/2)*erfc(sqrt(SNR/2));endhold on;semilogy(SNRindB2,theo_err_prb);grid on子函数：function [p]=exa0308_Fun_singlepe(snr_in_dB)E=1;SNR=10^(snr_in_dB/10);sgma=E/sqrt(SNR*2);N=10000;dsource=(sign(rand(1,N)-0.5+eps)+1)/2;numoferr=0;for i=1:Nif (dsource(i)==0)r=-E+sgma*randn;elser=E+sgma*randn;endif(r<0.5*E)decis=0;elsedecis=1;endif (decis~=dsource(i)) numoferr=numoferr+1; end end p=numoferr/N;:实验3的图 1：>> x=2*rand(5,4) x =0.7375 1.9921 0.9260 0.5703 1.3593 1.8842 0.2846 1.5663 0.0472 1.6797 0.9457 1.9128 0.8051 0.0293 0.9237 0.82740.8319 1.4452 1.7069 0.359802468101214161820-4-2024u (n )nu =-1.2765 1.1282 0.4296 -1.4450 -1.3768 0.7019 1.3848 -3.8243 -1.9091 -1.3563 1.7813 -2.9299 -0.8431 0.1011 0.4232 0.9075 0.4459 0.3491 -1.6545 0.89682：0.20.40.60.811.21.41.61.820.51050100150200250300350400-400-200200功率谱密度频率(HZ)3：0246810121010101010100。
matlab无线通信仿真代码
matlab无线通信仿真代码
本文分享了一份基于MATLAB平台的无线通信仿真代码。
该代码包含了常用的调制解调算法、信道模型、编码解码等模块，通过调用这些模块可以实现多种无线通信系统的仿真。
具体来说，该代码可以用于仿真以下几种无线通信系统：
1. 传统的模拟调制系统，包括AM、FM、PM等模式的调制解调。
2. 数字调制系统，包括PSK、QAM、FSK等模式的调制解调。
3. 多径衰落信道模型，包括AWGN信道、瑞利衰落信道、莱斯衰落信道等。
4. 信道编码解码，包括卷积码、Turbo码、LDPC码等。
除此之外，该代码还包含了多个演示程序，可以直观地展示无线通信系统的仿真结果。
这些演示程序涉及到信道容量、误码率、信噪比等指标，可以帮助用户深入理解无线通信系统的性能。
总之，该MATLAB无线通信仿真代码简单易用，功能齐全，可用于学术研究、教学实验、工程应用等多个领域。
- 1 -。
基于Matlab(Simulink)《通信原理》实验仿真(模拟部分)
基于Matlab（Simulink）《通信原理》实验仿真（模拟部分）基于Matlab(Simulink)《通信原理》实验仿真(模拟部分)摘要模拟通信在通信系统中的使用非常广泛，而MATLAB(Simulink)是用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，它可以解决包括信号和图像处理、通讯等众多应用领域中的问题。
利用MATLAB集成环境下的M文件和Simulink工具箱可以完成通讯系统设计与仿真，本文主要是利用MATLAB集成环境下的M文件，编写程序来实现AM、FM、VSB调制与解调过程，并分别绘制出其信号波形。
再通过Simulink工具箱对模拟通信系统进行建模仿真。
Simulation of communication in a communication system is very extensive, and the use of MATLAB (Simulink) is used to algorithm development, data visualization, data analysis and numerical calculation of the senior technical calculation language and interactive environment, it can solve the including signal and image processing, communicationetc many applications in question. MATLAB integration environment Mfiles and Simulink tool box can complete communication system design and simulation, this paper is mainly use of MATLAB integration environment, programming of the M files to achieve AM, FM, VSB modulation and demodulation process, and separately plot its signal waveform. Again through Simulink communication system toolbox of simulation modeling simulation .关键词模拟信号;AM;FM;VSB;调制解调;MATLAB(Simulink)目录:第一章绪论第二章理论与方法2.1 matlab简介2.2 Simulink简介2.3 通信原理概述第三章设计方案3.1用 MATLAB的M文件进行模拟调制3.1.1 基于MATLAB的AM信号调制与解调3.1.2 基于MATLAB的FM信号调制与解调3.1.3 基于MATLAB的VSB信号调制与解调3.2用Simulink对模拟通信系统进行建模仿真3.2.1频分复用和超外差接收机的仿真模型3.2.2调频立体声接收机模型第四章小结参考文献致谢附录1:程序清单第一章绪论调制在通信系统中的作用至关重要。
Matlab 通信系统建模与仿真例题源代码-第二章
% ch2example1prg1.mdt=1e-4; % 仿真采样间隔T=3*1e-3; % 仿真终止时间t=0:dt:T;input=2*cos(2*pi*1000*t); % 输入被调信号carrier=5*cos(2*pi*1e4*t); % 载波output=(2+0.5*input).*carrier; % 调制输出% 作图: 观察输入信号, 载波, 以及调制输出subplot(3,1,1); plot(t,input,'LineWidth',3);xlabel('时间 t');ylabel('被调信号');subplot(3,1,2); plot(t,carrier,'LineWidth',3);xlabel('时间t');ylabel('载波');subplot(3,1,3); plot(t,output,'LineWidth',3);xlabel('时间t');ylabel('调幅输出');% ch2example1prg2.mclear; % 清空内存变量，以避免以往运行的结果影响本程序dt=1e-5; % 仿真采样间隔T=3*1e-3; % 仿真终止时间t=0:dt:T;for k=1:length(t) % 基于时间流的仿真计算input(k)=2*cos(2*pi*1000*t(k)); % 第k个仿真步进时的输入被调信号carrier(k)=5*cos(2*pi*1e4*t(k)); % 第k个仿真步进时的载波output(k)=(2+0.5*input(k)).*carrier(k);% 第k个仿真步进时的调制输出end% 作图: 观察输入信号, 载波, 以及调制输出subplot(3,1,1); plot(t,input,'LineWidth',3);xlabel('时间 t');ylabel('被调信号');subplot(3,1,2); plot(t,carrier,'LineWidth',3);xlabel('时间t');ylabel('载波');subplot(3,1,3); plot(t,output,'LineWidth',3);xlabel('时间t');ylabel('调幅输出');% ch2example1prg3.mdt=1e-6; % 仿真采样间隔T=2*1e-3; % 仿真的帧周期for N=0:500 % 总共仿真的帧数t=N*T+(0:dt:T); % 帧中的取样时刻input=2*cos(2*pi*1005*t); % 输入被调信号carrier=5*cos(2*pi*(1e4)*t+0.1*randn); % 载波output=(2+0.5*input).*carrier; % 调制输出noise=randn(size(t)); % 噪声r=output+noise; % 调制信号通过加性噪声信道% 作图: 观察输入信号, 载波, 以及调制输出subplot(3,1,1); plot([0:dt:T],input,'LineWidth',3);xlabel('时间t');ylabel('被调信号');text(T*2/3,1.5,['当前帧数： N=',num2str(N)]); subplot(3,1,2); plot([0:dt:T],carrier,'LineWidth',3);xlabel('时间 t');ylabel('载波');subplot(3,1,3); plot([0:dt:T],r,'LineWidth',3);xlabel('时间 t');ylabel('调幅输出');set(gcf,'DoubleBuffer','on'); % 双缓冲避免作图闪烁drawnow;end% ch2example1prg25.m% 采用较低速率仿真，并用样条插值使得结果精细化simurate=3000; % 较低仿真速率mysimopts = simset('Solver','ode5'); % 设置仿真求解器为ode5的mysimopts = simset(mysimopts,'FixedStep',1./simurate); % 设置仿真步进timespan=[0, 0.01]; % 仿真时间段参数sim('ch2example25',timespan,mysimopts); % 设置仿真参数，仿真时间段并执行仿真plot(simout.time,simout.signals.values,'o-k');hold on;axis([0 0.005 -2 2]);t=timespan(1):1/30000:timespan(2); % 插值时间序列intp_y=interp1(simout.time,simout.signals.values,t,'spline');% 样条插值plot(t,intp_y,'xr'); % 画出插值结果% 采用较高速率仿真simurate=30000; % 较高仿真速率mysimopts = simset(mysimopts,'FixedStep',1./simurate); % 设置仿真步进sim('ch2example25',timespan,mysimopts); % 设置仿真参数，仿真时间段并执行仿真plot(simout.time,simout.signals.values,'b');% 画出仿真结果legend('较低速率仿真输出波形','样条插值精细化的结果','较高速率仿真输出波形');% ch2example2prg1.mdt=1e-5; % 仿真采样间隔R=1e3; % 电阻值C=1e-6; % 电容量T=5*1e-3; % 仿真区间从 -T 到 +Tt=-T:dt:T; % 计算的离散时刻序列y(1)=0; % 电容电压初始值, 在时间小于零区间将保持不变% 如果要仿真零输入响应, 可设置 y(1)=1 等非零值.% ----输入信号设定:可选择: 零输入,阶跃输入,正弦输入,方波输入等----x=zeros(size(t)); % 初始化输入信号存储矩阵x=1*(t>=0); % 在0时刻的输入信号跃变为1, 即输入为阶跃信号.% 如果要仿真零输入响应, 这里可设 x=0 即可% x=sin(2*pi*1000*t).*(t>=0); % 这是从0时刻开始的1000Hz的正弦信号% x=square(2*pi*500*t).*(t>=0); % 这是从0时刻开始的500Hz的方波信号% 仿真开始, 注意: 设零时刻之前电路不工作, 系统状态保持不变for k=1:length(t)time=-T+k*dt;if time>=0y(k+1)=y(k)+1./(R*C)*(x(k)-y(k))*dt; %递推求解下一个仿真时刻的状态值elsey(k+1)=y(k); % 在时间小于零时设电路断开，系统不工作endendsubplot(2,1,1);plot(t,x(1:length(t)),'LineWidth',3);axis([-T T -1.1 1.1]);xlabel('t');ylabel('input');subplot(2,1,2);plot(t,y(1:length(t)),'LineWidth',3);axis([-T T -1.1 1.1]);xlabel('t');ylabel('output');% ch2example3prg1.mdt=0.0001; % 仿真步进T=15; % 仿真时间长度t=0:dt:T; % 仿真计算时间序列g=9.8; % 重力加速度L=1; % 摆线长度m=10; % 摆锤质量k=5; % 空气阻力比例系数theta0=3.1; % 初始摆角设置v0=0; % 初始摆速设置v=zeros(size(t)); % 程序存储变量预先初始化,可提高执行速度theta=zeros(size(t));v(1)=v0; % 初始值赋值theta(1)=theta0;for n=1:length(t) % 仿真求解开始v(n+1)=v(n)+(g*sin(theta(n))-k./m.*v(n)).*dt;theta(n+1)=theta(n)-1./L.*v(n).*dt;end% 使用双坐标系统来作图, 注意作图和图形标注的技巧[AX,H1,H2] = plotyy(t,v(1:length(t)),t,theta(1:length(t)),'plot'); set(H1,'LineStyle','-'); % 设置作图线型set(H2,'LineStyle','-.');set(get(AX(1),'Ylabel'),'String','线速度 v(t) m/s');% 坐标标注set(get(AX(2),'Ylabel'),'String','角位移 \theta(t) rad ');xlabel('时间 t ');legend(H1,'线速度 v(t)',2);legend(H2,'角位移 \theta(t)',1);% ch2example5prg1.mdt=0.01; % 仿真步进T=15; % 仿真时间长度t=0:dt:T; % 仿真计算时间序列g=9.8; % 重力加速度L=1; % 摆线长度m=10; % 摆锤质量k=5; % 空气阻力比例系数theta0=3.1; % 初始摆角设置v0=0; % 初始摆速设置x0=[v0;theta0]; % 初始状态赋值par=[g;k;m;L]; % 系统参数赋值% 以欧拉算法计算并作图[t_out, x_out]=eulerode('pendulumstateeq',t', x0, [], par);plot(t_out, x_out(:,1),'-.b');hold on;% 以ode45算法计算并作图对比[t_out, x_out]=ode45('pendulumstateeq',t', x0, [], par);plot(t_out, x_out(:,1),'-k');xlabel('时间 t ');ylabel('线速度 m/s');legend('欧拉算法','ode45算法');% ch2example6main.mclear;v0=0; y0=1; % 球的初始状态x_state=[v0,y0];% 将初始状态赋值到状态变量中dt=0.01; % 仿真步进t=0:dt:5; % 仿真时间序列K=0.85; % 碰撞衰减系数for k=1:length(t) % 仿真开始，每次循环向前推进一个仿真步进dtx(k,:)=x_state; % 记录并保存当前状态的计算结果[t_out,x_out]=ode45('ch2example6statefun',[t(k),t(k)+dt],x_state);% 计算下一个时刻的新状态% 可换用ode23、ode113、ode23t、ode15s、ode23s以及ode23tb求解器 x_state=x_out(length(x_out),:); % 更新状态if (x_state(2)<=0) & (x_state(1)<0) % 当速度为负(球向下运动)且已经接触碰撞面x_state(1)=-K*x_state(1); % 处理碰撞瞬间情况：速度反向并衰减Kend% 动画作图：显示小球弹跳过程y=x_state(2); % 小球当前位置subplot(2,1,1);plot(0,y,'o');axis([-2 2 -0.1 1]); % 坐标范围固定set(gcf,'DoubleBuffer','on'); % 双缓冲避免作图闪烁drawnow; % 立即显示作图end% 仿真结束。
Matlab通信系统建模与仿真例题源代码-第五章
% ch5example2prg1.mSNR_in_dB=-10:2:30;SNR_in=10.^(SNR_in_dB./10); % 信道信噪比m_a=0.3; % 调制度P=0.5+(m_a^2)/4; % 信号功率for k=1:length(SNR_in)sigma2=P/SNR_in(k); % 计算信道噪声方差并送入仿真模型sim('ch5example2.mdl');% 执行仿真SNRdemod(k,:)=SNR_out; % 记录仿真结果endplot(SNR_in_dB, SNRdemod);xlabel('输入信噪比dB');ylabel('解调输出信噪比dB');legend('包络检波','相干解调');% ch5example7prog1.mclear;Fs=10000; % 仿真的采样率t=1/Fs:1/Fs:1; % 仿真时间点m_t(Fs*1)=0; % 基带信号变量初始化for F=150:400 % 基带信号发生：频率150Hz～400Hzm_t=m_t+0.003*sin(2*pi*F*t)*(400-F); % 幅度随线性递减endm_t90shift=imag(hilbert(m_t)); % 基带信号的希尔伯特变换carriercos=cos(2*pi*1000*t); % 1000Hz载波coscarriersin=sin(2*pi*1000*t); % 1000Hz正交载波sins_SSB1=m_t.*carriercos-m_t90shift.*carriersin; % 上边带SSBs_SSB2=m_t.*carriercos+m_t90shift.*carriersin; % 下边带SSB% 下面作出各波形以及频谱figure(1);subplot(4,2,1); plot(t(1:100),carriercos(1:100),...t(1:100),carriersin(1:100),'--r'); % 载波subplot(4,2,2); plot([0:9999],abs(fft(carriercos))); % 载波频谱axis([0 2000 -500 6000]);subplot(4,2,3); plot(t(1:100),m_t(1:100)); % 基带信号subplot(4,2,4); plot([0:9999],abs(fft(m_t))); % 信号频谱axis([0 2000 -500 6000]);subplot(4,2,5); plot(t(1:100),s_SSB1(1:100)); % SSB波形上边带subplot(4,2,6); plot([0:9999],abs(fft(s_SSB1))); % SSB频谱上边带axis([0 2000 -500 6000]);subplot(4,2,7); plot(t(1:100),s_SSB2(1:100)); % SSB波形下边带subplot(4,2,8); plot([0:9999],abs(fft(s_SSB2))); % SSB频谱下边带axis([0 2000 -500 6000]);% ch5example8prog1.mclear;Fs=10000; % 仿真的采样率t=1/Fs:1/Fs:1; % 仿真时间点m_t(Fs*1)=0; % 基带信号变量初始化for F=150:400 % 基带信号发生：频率150Hz～400Hzm_t=m_t+0.003*sin(2*pi*F*t)*(400-F); % 幅度随线性递减endm_t90shift=imag(hilbert(m_t)); % 基带信号的希尔伯特变换carriercos=cos(2*pi*1000*t); % 1000Hz载波coscarriersin=sin(2*pi*1000*t); % 1000Hz正交载波sins_SSB1=m_t.*carriercos-m_t90shift.*carriersin; % 上边带SSBout=s_SSB1.*carriercos; % 相干解调[a,b]=butter(4, 500/(Fs/2)); % 低通滤波器设计4阶,截止频率为500Hz demodsig=filter(a,b,out); % 解调输出% 下面作出各波形以及频谱figure(1);subplot(3,2,1); plot(t(1:100),s_SSB1(1:100)); % SSB波形subplot(3,2,2); plot([0:9999],abs(fft(s_SSB1))); % SSB频谱axis([0 5000 -500 6000]);subplot(3,2,3); plot(t(1:100),out(1:100)); % 相干解调波形subplot(3,2,4); plot([0:9999],abs(fft(out))); % 相干解调频谱axis([0 3000 -500 6000]);subplot(3,2,5); plot(t(1:100),demodsig(1:100)); % 低通输出信号subplot(3,2,6); plot([0:9999],abs(fft(demodsig)));% 低通输出信号的频谱axis([0 3000 -500 6000]);% ch5example8prog2.mclear;Fs=10000; % 仿真的采样率t=1/Fs:1/Fs:1; % 仿真时间点m_t(Fs*1)=0; % 基带信号变量初始化for F=150:400 % 基带信号发生：频率150Hz～400Hzm_t=m_t+0.003*sin(2*pi*F*t)*(400-F); % 幅度随线性递减endm_t90shift=imag(hilbert(m_t)); % 基带信号的希尔伯特变换carriercos=cos(2*pi*1000*t); % 1000Hz载波coscarriersin=sin(2*pi*1000*t); % 1000Hz正交载波sins_SSB1=m_t.*carriercos-m_t90shift.*carriersin; % 上边带SSBout=s_SSB1.*cos(2*pi*1018*t+1);% 存在频率误差的相位误差时的相干解调[a,b]=butter(4, 500/(Fs/2)); % 低通滤波器设计4阶,截止频率为500Hz demodsig=filter(a,b,out); % 解调输出% 下面作出各波形以及频谱figure(1);subplot(3,2,1); plot(t(1:100),s_SSB1(1:100)); % SSB波形subplot(3,2,2); plot([0:9999],abs(fft(s_SSB1))); % SSB频谱axis([0 5000 -500 6000]);subplot(3,2,3); plot(t(1:100),out(1:100)); % 相干解调波形subplot(3,2,4); plot([0:9999],abs(fft(out))); % 相干解调频谱axis([0 3000 -500 6000]);subplot(3,2,5); plot(t(1:100),demodsig(1:100)); % 低通输出信号subplot(3,2,6); plot([0:9999],abs(fft(demodsig)));% 低通输出信号的频谱axis([0 3000 -500 6000]);% ch5example9prog1.m[wav,fs]=wavread('GDGvoice8000.wav');t_end=1/fs *length(wav); % 计算声音的时间长度Fs=50000; % 仿真系统采样率t=1/Fs:1/Fs:t_end; % 仿真系统采样时间点% 设计300Hz～3400Hz的带通预滤波器H(z)[fenzi,fenmu]=butter(3,[300 3400]/(fs/2));% 对音频信号进行预滤波wav=filter(fenzi,fenmu,wav);% 利用插值函数将音频信号的采样率提升为Fs=50KHzwav=interp1([1/fs:1/fs:t_end],wav,t,'spline');wav_hilbert=imag(hilbert(wav)); % 音频信号的希尔伯特变换fc=10000; % 载波频率HzSSB_OUT=wav.*cos(2*pi*fc*t)-wav_hilbert.*sin(2*pi*fc*t); % 单边带调制figure(1); % 观察调制前后频谱subplot(2,2,1); plot(wav(53550:53750)); axis([0 200 -0.3 0.3]);subplot(2,2,2); psd(wav, 10000, Fs); axis([0 25000 -20 10]);subplot(2,2,3); plot(SSB_OUT(53550:53750)); axis([0 200 -0.3 0.3]);subplot(2,2,4); psd(SSB_OUT, 10000, Fs); axis([0 25000 -20 10]); wavwrite(0.5*SSB_OUT,Fs,'SSB_OUT.wav'); % 将SSB调制输出存盘备用% ch5example9prog2.mfunction out=ch5example9prog2(in, SNRdB)% SNR_dB设定信噪比% in 输入信号序列% out 信道输出序列Fs=50000; %系统采样率Power_of_in=var(in);Power_of_noise=Power_of_in/(10.^(SNRdB/10));bandwidth=13500-10000; % 信道带宽N0=Power_of_noise/bandwidth; % 噪声功率谱密度值W/HzGauss_noise=sqrt(N0*Fs/2).*randn(size(in));[num, den]=butter(4, [10000 13500]/(Fs/2)); % 带通信道10KHz到13.5KHz signal_of_filter_out=filter(num, den, in);noise_of_filter_out=filter(num, den, Gauss_noise);SNR_dB=10*log10(var(signal_of_filter_out)/var(noise_of_filter_out))% 测量得出的信噪比out=signal_of_filter_out+noise_of_filter_out; % 信道输出% ch5example9prog3.mclear;[recvsignal, Fs]=wavread('Channel_OUT.wav');% 读入信道输出信号数据t=(1/Fs:1/Fs: length(recvsignal)/Fs)';fc_local=10000-200; % 本地载波频率9.8KHzlocal_carrier=cos(2*pi*fc_local.*t); % 本地载波xianggan_out=recvsignal.*local_carrier; % 相干解调[fenzi,fenmu]=butter(3,[300 3400]/(Fs/2)); % 设计300Hz～3400Hz的带通滤波器demod_out=filter(fenzi,fenmu,xianggan_out); % 对相干输出信号进行滤波sound(demod_out/max(demod_out), Fs); % 播放解调音频wavwrite(demod_out,Fs,'SSBDemod_OUT.wav'); % 保存输出信号% ch5example10prog1.mclear;Y=[1:-1/7:0].*255; % 8级灰度N=8; % 8*N 列数S=ones(N,1);Y=S*Y;Y=reshape(Y,1,N*8);% 将矩阵扩大指定列数的一行M=40; % M 行数S=ones(M,1);Y=S*Y; % 重复M次矩阵行，得出图像矩阵I=uint8(Y); % 转换为uint8格式imshow(I); % 显示imwrite(I,'graybar.tif','tif'); % 并保存为图像文件% ch5example11prog1.mclear;I=imread('graybar.tif');% 读入图像（由上例生成）img=double(I')/255; % 转换为双精度数以便计算[M,N]=size(img);for K=1:2for column=K:2:N % K=1扫描奇数行for row=1:Mdot=img(row,column); % 取当前扫描点的灰度值greydot=[dot,dot,dot]; % R,G,B均设置为相同值，即显示灰度plot(row,N-column,'s','MarkerFaceColor',greydot,...'MarkerEdgeColor',greydot,'MarkerSize',6);hold on; % 保持扫描点axis([1,M, 1, N]); % 作图范围set(gcf,'DoubleBuffer','on');% 双缓冲避免作图闪烁drawnow; % 立即作图endendend% ch5example12prg1.mclear;R=zeros(100,100);G=zeros(100,100);B=zeros(100,100);R(1:50,30:100)=255;G(30:100,50:100)=255;B(1:70,1:70)=255; % 三基色产生I(:,:,1)=R;I(:,:,2)=G;I(:,:,3)=B; % 合成为100*100图像矩阵imshow(uint8(I));% 显示% ch5example13prg1.m% 生成400*480彩条图像clear;R=[1 1 0 0 1 1 0 0];G=[1 1 1 1 0 0 0 0];B=[1 0 1 0 1 0 1 0];N=60; % N*8列S=ones(N,1);R=S*R; R=reshape(R,1,N*8);G=S*G; G=reshape(G,1,N*8);B=S*B; B=reshape(B,1,N*8);M=400; % 行数S=ones(M,1);R=S*R; G=S*G; B=S*B; % 得出三基色矩阵I(:,:,1)=R; I(:,:,2)=G; I(:,:,3)=B; % 组合I=uint8(I*255); % 类型转换imshow(I); % 显示并保存imwrite(I,'colorbar.tif','tif');% ch5example14prg1.mclear;I=imread('colorbar.tif');% 或用autumn.tif，sydney.JPG等图像文件figure(1);imshow(I);I=double(I);[m,n,p]=size(I);I=-(0.75-0.125)./(255).*I+0.75; % 换算为0.125到0.75电平R=I(:,:,1);G=I(:,:,2);B=I(:,:,3);% 三基色分离R=reshape(R',1,m*n); % 转换为1维（红）G=reshape(G',1,m*n); % 转换为1维（绿）B=reshape(B',1,m*n); % 转换为1维（蓝）Y=0.30*R+0.59*G+0.11*B;% 亮度R_Y=R-Y; % 色差R0B_Y=B-Y; % 色差B0V=0.877*R_Y; % 色差电平压缩U=0.493*B_Y; % 色差电平压缩f_c=283.5*15625; % 副载波频率tvY=[interp1(Y,(1:n/980:m*n),'nearest')];% 插值:行采样980点tvV=[interp1(V,(1:n/980:m*n),'nearest')];% 插值tvU=[interp1(U,(1:n/980:m*n),'nearest')];% 插值% 全电视信号产生for h=1:(m-2)s(1280*h+(1:102))=0.75; % 消隐脉冲s(1280*h+(103:197))=1; % 行同步头s(1280*h+(198:300))=0.75;% 消隐脉冲s(1280*h+(301:1280))=tvY(h*980+(-1:978));% 图像一行像素－亮度% 色差VsV(1280*h+(1:300))=0; % 消隐部分sV(1280*h+(220:265))=0.12;% 色同步选通sV(1280*h+(301:1280))=tvV(h*980+(-1:978));% 图像信号部分% 色差UsU(1280*h+(1:300))=0;sU(1280*h+(301:1280))=tvU(h*980+(-1:978));endt=0:0.5e-7:0.5e-7*(length(s)-1);% 计算时间点序列，采样率20MHzF_v=sV.*cos(2*pi*f_c.*t); % 色差信号正交调制F_u=sU.*sin(2*pi*f_c.*t);c=F_v+F_u; % 调制输出的色差信号TVsignal=s+c; % 合成彩色电视信号figure(2); plot(t,TVsignal); xlabel('时间');% 彩色电视信号figure(3); plot(t,s); xlabel('时间'); % 亮度信号% 功率谱计算[Pxxcc,F]=psd(c,1e5,2e7);figure(4); plot(F,10*log10(Pxxcc),'g');hold on;% 调制后的色差信号频谱[Pxx,F]=psd(s,1e5,2e7); % FFT长度为1e5点，故频率分辨率100Hz figure(4); plot(F,10*log10(Pxx),'k'); % 亮度信号频谱xlabel('频率Hz');ylabel('功率谱密度dB');% ch5example15prg1.mclear;ch5example14prg1; % 用实例5.14程序产生彩色全电视信号%---带限高斯信道---TVsignal=TVsignal+0.1*randn(size(TVsignal)); % 加入噪声[b,a]=butter(2,5e6/(20e6/2)); % 传输信道带宽5MHzTVsignal=filter(b,a,TVsignal); % 通过带限信道%---接收部分---[b,a]=butter(2,0.5e6/(20e6/2)); % 行同步滤波h_pulse=filter(b,a,TVsignal);h_pulse=h_pulse>0.87; % 行同步脉冲选出。
MATLAB仿真实例(通信原理)
值。仿真时我们用这个样值集合来表示信号 s( t ) 。显然 t 反映了仿真系统对信号波形的分 辨率， t 越小则仿真的精确度越高。据通信原理所学，信号被取样以后的频谱是频率的周 1 1 期函数，其重复周期是 。如果信号的最高频率为 f H ，那么必须有 f H 才能保证不 t 2 t 1 发生混叠失真。我们称 Bs 为仿真系统的系统带宽。如果我们的仿真程序中设定的采 2 t 样间隔是 t ，那么我们不能用此仿真程序来研究带宽大于 Bs 的信号或系统。 此外，信号 s( t ) 的频谱 S f 通常来说也是定义在频率区间 , 上的连续函数，所以仿 真频域特性时，我们也必须把 S f 截短并取样。考虑到系统带宽是 Bs ，我们把频域的截短
三、 实验内容
3.1 傅里叶变换与傅里叶反变换 对于确定信号 f (t ) ，其傅里叶变换为：
F( f )
F ( f ) 傅里叶反变换为：
f (t )e 2 j ft dt
f (t ) F ( f )e2 j ft df
在通信原理仿真中，傅里叶变换与傅里叶反变换会经常用到，我们可以利用 MATLAB 的快速傅里叶变换函数 fft 与快速傅里叶反变换函数 ifft 编写傅里叶变换子程序与傅里叶反 变换子程序。其 程 序 代 码 如 下 ：
傅里叶变换子程序：
%傅里叶变换子程序 function X=t2f(x) global dt df N t f T %X=t2f(x) %x 为时域的取样值矢量 %X 为 x 的傅氏变换 %X 与 x 长度相同,并为 2 的整幂。
基于MATLAB的移动通信信道建模与仿真
基于MATLAB的移动通信信道建模与仿真基于MATLAB的移动通信信道建模与仿真1.引言1.1 项目背景1.2 研究目的1.3 研究内容1.4 研究方法2.移动通信信道模型2.1 信道模型概述2.2 多径信道模型2.2.1 多径传播原理2.2.2 多径信道建模方法2.3 表面波信道模型2.3.1 表面波传播原理2.3.2 表面波信道建模方法2.4 阴影衰落信道模型2.4.1 阴影衰落传播原理2.4.2 阴影衰落信道建模方法 2.5 小尺度信道模型2.5.1 小尺度信道传播原理 2.5.2 小尺度信道建模方法 2.6 大尺度信道模型2.6.1 大尺度信道传播原理2.6.2 大尺度信道建模方法3.移动通信信道仿真3.1 仿真平台及工具介绍3.2 仿真参数设置3.3 仿真方法和步骤3.4 仿真结果及分析4.结果与讨论4.1 仿真结果分析4.2 结果对比与验证4.3 结果的实际应用价值5.总结与展望5.1 研究总结5.2 存在问题与改进方向5.3 研究展望6.附件附件1：MATLAB程序代码附件2：仿真数据结果法律名词及注释：1.信道模型：信道模型是对移动通信信道的数学描述，用于模拟信道传输过程。
2.多径信道：多径信道指信号在传播过程中由于反射、折射等现象导致的多条信号路径同时存在的信道。
3.表面波信道：表面波信道是指信号在大地表面的传播过程中，通过地表产生的表面波信道。
4.阴影衰落信道：阴影衰落信道是指信号在传播过程中由于物体阻挡等原因导致信号强度发生波动的信道。
5.小尺度信道：小尺度信道指信道中存在的快速变化的信道衰落现象，如多径信道引起的快速衰落。
6.大尺度信道：大尺度信道指信道的整体特征，如路径损耗、平均衰减等。
3《MATLAB Simulink与控制系统仿真(第3版)》的课件 第3章 Simulink仿真
3.4.3 Simulink模块间的连线处理
（1）改变粗细：线所以有粗细是因为线引出的信号可以是标量信号或向量信号， 当选中Format菜单下的Wide Vector Lines时，线的粗细会根据线所引出的信号是 标量还是向量而改变，如果信号为标量则为细线，若为向量则为粗线。选中 Vector Line Widths则可以显示出向量引出线的宽度，即向量信号由多少个单一信 号合成。 （2）设定标签：只要在线上双击鼠标，即可输入该线的说明标签。也可以通过 选中线，然后打开Edit菜单下的Signal Properties进行设定，其中Signal name属性 的作用是标明信号的名称，设置这个名称反映在模型上的直接效果就是与该信 号有关的端口相连的所有直线附近都会出现写有信号名称的标签。 （3）线的折弯：按住Shift键，再用鼠标在要折弯的线处单击一下，就会出现圆 圈，表示折点，利用折点就可以改变线的形状。 （4）线的分支：按住鼠标右键，在需要分支的地方拉出即可，或者按住Ctrl键 并在要建立分支的地方用鼠标拉出即可。
3.3.1 Simulink模块库分类
Simulink模块库按功能分为16类子模块库
3.4 Simulink功能模块的处理
图3.8 “功能模块参数设置”对话框
图3.10 “示波器属性”对话窗框
3.4.2 Simulink模块的基本操作
（1）移动：选中模块，按住鼠标左键将其拖曳到所需的位置即可。若要脱离线而移动，可按住Shift键再进行拖 曳。 （2）复制：选中模块，按住鼠标右键进行拖曳即可复制同样的一个功能模块。 （3）删除：选中模块，按Delete键即可。若要删除多个模块，可以同时按住Shift键，再用鼠标选中多个模块，按 Delete键即可；也可以用鼠标选取某区域，再按Delete键就可以把该区域中的所有模块和线等全部删除。 （4）转向：为了能够顺序连接功能模块的输入和输出端，功能模块有时需要转向。在菜单Format中选择Flip Block旋转180°，选择Rotate Block顺时针旋转90°；或者直接按Ctrl+F组合键执行Flip Block，按Ctrl+R组合键执 行Rotate Block。 （5）改变大小：选中模块，对模块出现的4个黑色标记进行拖曳即可。 （6）模块命名：先用鼠标在需要更改的名称上单击一下，然后直接更改即可。名称在功能模块上的位置也可以 变换180°，可以用Format菜单中的Flip Name来实现，也可以直接通过鼠标进行拖曳。Hide Name可以隐藏模块名 称。 （7）颜色设定：Format菜单中的Foreground Color可以改变模块的前景颜色，Background Color可以改变模块的背 景颜色，而模型窗口的颜色可以通过Screen Color来改变。 （8）参数设定：用鼠标双击模块就可以进入模块的参数设定窗口，从而对模块进行参数设定。参数设定窗口包 含了该模块的基本功能帮助，为获得更详尽的帮助，可以单击其上的“Help”按钮。通过对模块的参数设定，就 可以获得需要的功能模块。 （9）属性设定：选中模块，打开Edit菜单的Block Properties可以对模块进行属性设定，包括对Description、 Priority、Tag、Open function、Attributes format string等属性的设定。其中Open function属性是一个很有用的属性， 通过它指定一个函数名，当模块被双击之后，Simulink就会调用该函数并执行，这种函数在MATLAB中称为回调 函数。 （10）模块的输入/输出信号：模块处理的信号包括标量信号和向量信号。标量信号是一种单一信号，而向量信号 为一种复合信号，是多个信号的集合，它对应着系统中几条连线的合成。默认情况下，大多数模块的输出都为标 量信号，对于输入信号，模块都具有一种“智能”的识别功能，能自动进行匹配。某些模块通过对参数的设定， 可以使模块输出向量信号。
MatlabSimulink通信系统建模与仿真课程设计
MatlabSimulink通信系统建模与仿真课程设计MatlabSimulink通信系统建模与仿真课程设计电子信息课程设计题目：Matlab/Simulink通信系统建模与仿真班级：2008级电子（X）班学号：姓名：电子信息课程设计Matlab/Simulink通信系统建模与仿真一、设计目的：学习Matlab/Simulink的功能及基本用法，对给定系统进行建模与仿真。
二、基本知识：Simulink是用来对动态系统进行建模、仿真和分析的软件包，依托于MATLAB丰富的仿真资源，可应用于任何使用数学方式进行描述的动态系统，其最大优点是易学、易用，只需用鼠标拖动模块框图就能迅速建立起系统的框图模型。
三、设计内容：1、基本练习：（1）启动*****K：先启动MATLAB，在命令窗口中键入：simulink，回车；或点击窗口上的*****K图标按钮。
图（1）建立simulink （2）点击File＼new＼Model或白纸图标，打开一个创建新模型的窗口。
（3）移动模块到新建的窗口，并按需要排布。
（4）连接模块：将光标指向起始模块的输出口，光标变为“+”，然后拖动鼠标到目标模块的输入口；或者，先单击起始模块，按下Ctrl键再单击目标模块。
（5）在连线中插入模块：只需将模块拖动到连线上。
（6）连线的分支与改变：用鼠标单击要分支的连线，光标变为“+”，然后拖动到目标模块；单击并拖动连线可改变连线的路径。
（7）信号的组合：用Mux模块可将多个标量信号组合成一个失量信号，送到另一模块（如示波器Scope）。
（8）生成标签信号：双击需要加入标签的信号线，会出现标签编辑框，键入标签文本即可。
或点击Edit＼Signal Properties。
传递：选择信号线并双击，在标签编辑框中键入，并在该尖括号内键入信号标签即可。
四、建立模型1. 建立仿真模型（1）在simulink library browser 中查找元器件，并放置在创建的新模型的窗口中，连接元器件，得到如下的仿真模型。
通信原理基于matlab的计算机仿真源代码
例错误!文档中没有指定样式的文字。
-1%周期信号（方波）的展开,fb_jinshi.mclose all;clear all;N=100; %取展开式的项数为2N＋1项T=1;fs=1/T;N_sample=128; %为了画出波形，设置每个周期的采样点数dt = T/N_sample;t=0:dt:10*T-dt;n=-N:N;Fn = sinc(n/2).*exp(-j*n*pi/2);Fn(N+1)=0;ft = zeros(1,length(t));for m=-N:Nft = ft + Fn(m+N+1)*exp(j*2*pi*m*fs*t);endplot(t,ft)例错误!文档中没有指定样式的文字。
-4利用FFT计算信号的频谱并与信号的真实频谱的抽样比较。
脚本文件T2F.m定义了函数T2F，计算信号的傅立叶变换。
function [f,sf]= T2F(t,st)%This is a function using the FFT function to calculate a signal's Fourier %Translation%Input is the time and the signal vectors,the length of time must greater %than 2%Output is the frequency and the signal spectrumdt = t(2)-t(1);T=t(end);df = 1/T;N = length(st);f=-N/2*df:df:N/2*df-df;sf = fft(st);sf = T/N*fftshift(sf);脚本文件F2T.m定义了函数F2T，计算信号的反傅立叶变换。
function [t st]=F2T(f,sf)%This function calculate the time signal using ifft function for the input %signal's spectrumdf = f(2)-f(1);Fmx = ( f(end)-f(1) +df);dt = 1/Fmx;N = length(sf);T = dt*N;%t=-T/2:dt:T/2-dt;t = 0:dt:T-dt;sff = fftshift(sf);st = Fmx*ifft(sff);另写脚本文件fb_spec.m如下：%方波的傅氏变换, fb_spec.mclear all;close all;T=1;N_sample = 128;dt=T/N_sample;t=0:dt:T-dt;st=[ones(1,N_sample/2), -ones(1,N_sample/2)]; %方波一个周期subplot(211);plot(t,st);axis([0 1 -2 2]);xlabel('t'); ylabel('s(t)');subplot(212);[f sf]=T2F(t,st); %方波频谱plot(f,abs(sf)); hold on;axis([-10 10 0 1]);xlabel('f');ylabel('|S(f)|');%根据傅氏变换计算得到的信号频谱相应位置的抽样值sff= T^2*j*pi*f*0.5.*exp(-j*2*pi*f*T).*sinc(f*T*0.5).*sinc(f*T*0.5);plot(f,abs(sff),'r-')例错误!文档中没有指定样式的文字。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
% ch3example1A.mclear;f_p=2400; f_s=5000; R_p=3; R_s=25; % 设计要求指标[n, fn]=buttord(f_p,f_s,R_p,R_s, 's'); % 计算阶数和截止频率Wn=2*pi*fn; % 转换为角频率[b,a]=butter(n, Wn, 's'); % 计算H(s)f=0:100:10000; % 计算频率点和频率范围s=j*2*pi*f; % s=jw=j*2*pi*fH_s=polyval(b,s)./polyval(a,s); % 计算相应频率点处H(s)的值figure(1);subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性axis([0 10000 -40 1]);xlabel('频率Hz');ylabel('幅度dB');subplot(2,1,2); plot(f, angle(H_s)); % 相频特性xlabel('频率Hz');ylabel('相角rad');figure(2); freqs(b,a); % 也可用指令freqs直接画出H(s)的频率响应曲线。
% ch3example1B.mclear;f_p=2400; f_s=5000; R_p=3; R_s=25; % 设计要求指标[n, fn]=ellipord(f_p,f_s,R_p,R_s,'s'); % 计算阶数和截止频率Wn=2*pi*fn; % 转换为角频率[b,a]=ellip(n,R_p,R_s,Wn,'s'); % 计算H(s)f=0:100:10000; % 计算频率点和频率范围s=j*2*pi*f; % s=jw=j*2*pi*fH_s=polyval(b,s)./polyval(a,s); % 计算相应频率点处H(s)的值figure(1);subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性axis([0 10000 -40 1]);xlabel('频率Hz');ylabel('幅度dB');subplot(2,1,2); plot(f, angle(H_s)); % 相频特性xlabel('频率Hz');ylabel('相角rad');figure(2); freqs(b,a); % 也可用指令freqs直接画出H(s)的频率响应曲线。
% ch3example2A.mf_N=8000; % 采样率f_p=2100; f_s=2500; R_p=3; R_s=25; % 设计要求指标Ws=f_s/(f_N/2); Wp=f_p/(f_N/2); % 计算归一化频率[n, Wn]=buttord(Wp,Ws,R_p,R_s); % 计算阶数和截止频率[b,a]=butter(n, Wn); % 计算H(z)figure(1);freqz(b,a, 1000, 8000) % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率)subplot(2,1,1); axis([0 4000 -30 3])figure(2); % 第二种作图方法f=0:40:4000; % 计算频率点和频率范围z=exp(j*2*pi*f./(f_N)); %H_z=polyval(b,z)./polyval(a,z); % 计算相应频率点处H(s)的值subplot(2,1,1); plot(f, 20*log10(abs(H_z))); % 幅频特性axis([0 4000 -40 1]);xlabel('频率Hz');ylabel('幅度dB');subplot(2,1,2); plot(f, angle(H_z)); % 相频特性xlabel('频率Hz');ylabel('相角rad');% ch3example3A.mf_N=8000; % 采样率f_p=1000; f_s=700; R_p=3; R_s=20; % 设计要求指标Ws=f_s/(f_N/2); Wp=f_p/(f_N/2); % 计算归一化频率[n, Wn]=cheb1ord(Wp,Ws,R_p,R_s); % 计算阶数和截止频率[b,a]=cheby1(n, R_p, Wn, 'high'); % 计算H(z)freqz(b,a, 1000, 8000) % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率) subplot(2,1,1); axis([0 4000 -30 3])% ch3example4A.mf_N=10000; % 采样率f_p=[1000, 1500]; f_s=[600, 1900]; R_p=3; R_s=20; % 设计要求指标Ws=f_s/(f_N/2); Wp=f_p/(f_N/2); % 计算归一化频率[n, Wn]=ellipord(Wp,Ws,R_p,R_s); % 计算阶数和截止频率[b,a]=ellip(n, R_p, R_s, Wn); % 计算H(z)freqz(b,a, 1000, 10000) % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率) subplot(2,1,1); axis([0 5000 -30 3])% ch3example5A.mf_N=10000; %采样率f_p=[1000, 1500]; f_s=[1200, 1300]; R_p=3; R_s=30; %设计要求指标Ws=f_s/(f_N/2); Wp=f_p/(f_N/2); %计算归一化频率[n, Wn]=cheb2ord(Wp,Ws,R_p,R_s); %计算阶数和截止频率[b,a]=cheby2(n,R_s, Wn, 'stop'); %计算带阻H(z)系数freqz(b,a, 1000, 10000) %作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率) subplot(2,1,1); axis([0 5000 -35 3])% ch3example6A.mFs=8000; Ts=1/8000; % 采样率f0=500; % 梳状滤波器开槽基频率bw=60/(Fs/2); % 归一化开槽带宽ab=-3 ; % 计算开槽带宽位置处的衰减分贝值n=Fs/f0; % 计算滤波器阶数[num,den] = iircomb(n,bw,ab, 'notch'); % 计算H(z)freqz(num,den, 4000, 8000); % 作出H(z)的幅频相频图axis([0 4000 -30 5]);% ch3example10.msim('ch3example10.mdl');bar([-10:1/5:10-1/5],histdata);axis([-10 10 0 0.05]);% 作出归一化统计直方图% ch3example10prg.msim('ch3example10.mdl');bar([-10:1/5:10-1/5],histdata*5);axis([-10 10 0 0.25]);% 作出归一化统计直方图% ch3example11prg1.mt=-1: 0.01 :10 ; % 计算时间范围f_t=exp(-t).*(t>0);subplot(3,1,1);plot(t,f_t); % 时域波形axis([-1 10 -0.1 1.1]); xlabel('time (sec)');w=-40: 0.1 : 40; % 计算角频率范围F_w=1./(1+j*w); % 频谱理论结果subplot(3,1,2);plot(w, abs(F_w)); % 频域幅度谱axis([-40 40 0 1.1]);xlabel('freq (rad/s)');subplot(3,1,3);plot(w, angle(F_w)); % 频域相位谱axis([-40 40 -pi/2 pi/2]);xlabel('freq (rad/s)');% ch3example11prg2.mw_m=40; % 截断频率T=pi/w_m; % 采样间隔L=5;t=0: T :L ; % 时域截断x_t=exp(-t).*(t>0); %信号序列N=length(x_t); %序列长度(点数)%-----------------------------------X_k=fft(x_t); % FFT计算%-----------------------------------w0=2*pi/(N*T); % 离散频率间隔kw=2*pi/(N*T) .*[0: N-1]; % 离散频率样点X_kw=T.*X_k; % 乘以T得到连续傅利叶变换频谱的样值%-----------------------------------subplot(2,1,1);plot(kw, abs(X_kw),'.','MarkerSize', 14); % 作出数值计算的幅度谱点axis([-40 90 -0.1 1.1]);xlabel('freq (rad/s)');hold on; % 保持前面作图曲线不被擦除w=-40: 0.1 : 40;X_w=1./(1+j*w); % 理论计算频谱表达式plot(w, abs(X_w)); % 作图对比subplot(2,1,2);plot(kw, angle(X_kw),'.','MarkerSize', 14); % 作出数值计算的相位谱点hold on;plot(w, angle(X_w)); % 频域相位谱axis([-40 90 -pi pi]);xlabel('freq (rad/s)');% ch3example12prg1.mDf=5; % 频率间隔f_s=2*500; % 采样率N=f_s/Df; % 序列点数t=0:1./f_s:(N-1)./f_s; % 计算时间段freq=0:Df:(N-1)*Df; % 计算频率段f_t=2*sin(2*pi*100*t)+cos(2*pi*180*t); % 信号F_f=1/f_s*fft(f_t,N); % 用FFT计算频谱plot(freq-f_s/2,abs(fftshift(F_f))); % 将零频率移动到FFT中心xlabel('频率Hz'); ylabel('幅度谱'); % 并作出幅度频谱% ch3example14prg1.mfs=200; % 采样率Delta_f = 1; % 频率分辨率T = 1/fs; % 时间分辨率L=1/ Delta_f; % 时域截取长度N= floor(fs/Delta_f)+1; % 计算截断信号的采样点数t=0:T:L; % 截取时间段和采样时间点freq=0: Delta_f :fs; % 分析的频率范围和频率分辨率f_t=(sin(2*pi*50*t)+0.7*sin(2*pi*75*t))';% 在截取范围内的分析的信号时域波形f_t_rectwin= rectwin(N).*f_t./sqrt(sum(abs(rectwin(N).^2))./N);% 矩形窗（功率归一化）f_t_hamming= hamming(N) .*f_t./sqrt(sum(abs(hamming(N).^2))./N);% 海明窗（功率归一化）f_t_hann = hann(N) .*f_t./sqrt(sum(abs(hann(N).^2))./N);% 汉宁窗（功率归一化）F_w_rectwin =T.* fft(f_t_rectwin, N);% 进行N点FFT，并乘以采样时间间隔T得到频谱F_w_hamming =T.* fft(f_t_hamming, N); % 加海明窗的频谱F_w_hann =T.* fft(f_t_hann, N); % 加汉宁窗的频谱figure(1); subplot(2,2,1);plot(t,f_t);title('Original Signal');subplot(2,2,2);plot(t, f_t_rectwin);title('Rectwin Windowing');subplot(2,2,3);plot(t, f_t_hamming);title('hamming Windowing');subplot(2,2,4);plot(t, f_t_hann);title('hanning Windowing');figure(2);subplot(3,1,1);plot(freq, 10*log10(abs(F_w_rectwin)));title('Rectwin Windowing Spectrum');ylabel('幅度dB');axis([0,200,-50,0]);grid on;subplot(3,1,2);plot(freq, 10*log10(abs(F_w_hamming)));title('hamming Windowing Spectrum');ylabel('幅度dB');axis([0,200,-50,0]);grid on;subplot(3,1,3);plot(freq, 10*log10(abs(F_w_hann)));title('hanning Windowing Spectrum');ylabel('幅度dB');axis([0,200,-50,0]);grid on;p_original_signal=var(f_t) % 计算原始信号的平均功率p_hannwindowed=sum((abs(F_w_hann/T)).^2)/(N^2) % 计算加窗后的信号功率% ch3example16prg1.mfs=500; % 采样率Df=1; % 频率分辨率N=floor(fs/Df)+1; % 计算的序列点数t=0:1/fs:(N-1)/fs; % 截取信号的时间段F=0:Df:fs; % 功率谱估计的频率分辨率和范围xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));% 截取时间段上的离散信号样点序列Pxx=abs(fft(xk)).^2/(N^2); % 功率谱估计Pav_timedomaim=sum(xk.^2)/N % 在时域计算信号功率Pav_freqdomain=sum(Pxx) % 通过功率谱计算信号功率plot(F,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB') % 作出功率谱密度图% ch3example16prg2.mfs=500; % 采样率Df=1; % 频率分辨率N=floor(fs/Df)+1; % 计算的序列点数t=0:1/fs:(N-1)/fs; % 截取信号的时间段F=0:Df:fs; % 功率谱估计的频率分辨率和范围xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));% 截取时间段上的离散信号样点序列Pxx=(abs(fft(xk(1:167))).^2+...abs(fft(xk(168:334))).^2+...abs(fft(xk(335:501))).^2)/3/((N/3)^2);Pav_timedomaim=sum(xk.^2)/N % 在时域计算信号功率Pav_freqdomain=sum(Pxx) % 通过功率谱计算信号功率plot(0:3:fs,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB');%作出功率谱密度图% ch3example16prg3.mfs=500; % 采样率Df=1; % 频率分辨率N=floor(fs/Df)+1; % 计算的序列点数t=0:1/fs:(N-1)/fs; % 截取信号的时间段F=0:Df:fs; % 功率谱估计的频率分辨率和范围xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));% 截取时间段上的离散信号样点序列Pxx=(abs(fft(xk(1:167))).^2+...abs(fft(xk(83:249))).^2+...abs(fft(xk(168:334))).^2+...abs(fft(xk(250:416))).^2+...abs(fft(xk(335:501))).^2)/5/((N/3)^2);Pav_timedomaim=sum(xk.^2)/N % 在时域计算信号功率Pav_freqdomain=sum(Pxx) % 通过功率谱计算信号功率plot(0:3:fs,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB');%作出功率谱密度图% ch3example16prg4.mfs=500; % 采样率Df=1; % 频率分辨率N=floor(fs/Df)+1; % 计算的序列点数t=0:1/fs:(N-1)/fs; % 截取信号的时间段F=0:Df:fs; % 功率谱估计的频率分辨率和范围xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));% 截取时间段上的离散信号样点序列w=hamming(167)';%海明窗w=w*sqrt(167/sum(w.*w));% 使得海明窗与矩形窗等能量，即加窗后不对信号功率产生影响Pxx=(abs(fft(w.*xk(1:167))).^2+...abs(fft(w.*xk(83:249))).^2+...abs(fft(w.*xk(168:334))).^2+...abs(fft(w.*xk(250:416))).^2+...abs(fft(w.*xk(335:501))).^2)/5/((N/3)^2);Pav_timedomaim=sum(xk.^2)/N % 在时域计算信号功率Pav_freqdomain=sum(Pxx) % 通过功率谱计算信号功率plot(0:3:fs,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB');%作出功率谱密度图% ch3example16prg5.mfs=500; % 采样率Df=1; % 频率分辨率N=floor(fs/Df)+1; % 计算的序列点数t=0:1/fs:(N-1)/fs; % 截取信号的时间段F=0:Df:fs; % 功率谱估计的频率分辨率和范围xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));% 截取时间段上的离散信号样点序列[Pxx,F] = psd(xk,512,500,hamming(256),128);plot(F,10*log10(Pxx/(512/2))); xlabel('freq Hz');ylabel('PSD dB')Pav=sum(Pxx/(512/2)) % 通过功率谱计算信号功率% ch3example23prg1.mkc=150e3; % Hz/V VCO控制灵敏度omega_n=2*pi*15e3/2.5; % PLL自然角频率K=2*pi*(0.5*kc); % 估算环路增益zeta=1; % 临界阻尼tau_1=K/((omega_n).^2);tau_2=2*zeta/omega_n-1/K;freq=0:10:100e3; % 计算频率范围0到100KHzs=j*2*pi*freq;G_s=(1+tau_2*s)./(1+tau_1*s); % 环路滤波器传递函数figure(1);semilogx(freq,(abs(G_s)));% 作出环路滤波器的频率响应title('环路滤波器幅频响应');xlabel('Hz');ylabel('|G(s)|');grid on;b=[tau_2,1]; % 环路滤波器分子系数向量a=[tau_1,1]; % 环路滤波器分母系数向量H_s=(G_s*K./s)./(1+G_s*K./s);figure(2);semilogx(freq,20*log10(abs(H_s)));% 作出闭环频率响应title('PLL线性相位模型闭环频率响应');xlabel('Hz');ylabel('20log|H(s)|(dB)'); grid on;。