Matlab 通信系统建模与仿真例题源代码第二章.docx
- 文档编号:5750032
- 上传时间:2022-12-31
- 格式:DOCX
- 页数:29
- 大小:24.25KB
Matlab 通信系统建模与仿真例题源代码第二章.docx
《Matlab 通信系统建模与仿真例题源代码第二章.docx》由会员分享,可在线阅读,更多相关《Matlab 通信系统建模与仿真例题源代码第二章.docx(29页珍藏版)》请在冰豆网上搜索。
Matlab通信系统建模与仿真例题源代码第二章
%ch2example1prg1.m
dt=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.m
clear;%清空内存变量,以避免以往运行的结果影响本程序
dt=1e-5;%仿真采样间隔
T=3*1e-3;%仿真终止时间
t=0:
dt:
T;
fork=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.m
dt=1e-6;%仿真采样间隔
T=2*1e-3;%仿真的帧周期
forN=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');
holdon;axis([00.005-22]);
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.m
dt=1e-5;%仿真采样间隔
R=1e3;%电阻值
C=1e-6;%电容量
T=5*1e-3;%仿真区间从-T到+T
t=-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的方波信号
%仿真开始,注意:
设零时刻之前电路不工作,系统状态保持不变
fork=1:
length(t)
time=-T+k*dt;
iftime>=0
y(k+1)=y(k)+1./(R*C)*(x(k)-y(k))*dt;%递推求解下一个仿真时刻的状态值
else
y(k+1)=y(k);%在时间小于零时设电路断开,系统不工作
end
end
subplot(2,1,1);plot(t,x(1:
length(t)),'LineWidth',3);axis([-TT-1.11.1]);
xlabel('t');ylabel('input');
subplot(2,1,2);plot(t,y(1:
length(t)),'LineWidth',3);axis([-TT-1.11.1]);
xlabel('t');ylabel('output');
%ch2example3prg1.m
dt=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;
forn=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.m
dt=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');holdon;
%以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.m
clear;
v0=0;y0=1;%球的初始状态
x_state=[v0,y0];%将初始状态赋值到状态变量中
dt=0.01;%仿真步进
t=0:
dt:
5;%仿真时间序列
K=0.85;%碰撞衰减系数
fork=1:
length(t)%仿真开始,每次循环向前推进一个仿真步进dt
x(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);%处理碰撞瞬间情况:
速度反向并衰减K
end
%动画作图:
显示小球弹跳过程
y=x_state
(2);%小球当前位置
subplot(2,1,1);plot(0,y,'o');
axis([-22-0.11]);%坐标范围固定
set(gcf,'DoubleBuffer','on');%双缓冲避免作图闪烁
drawnow;%立即显示作图
end
%仿真结束。
最后输出计算时间序列上的结果:
随时间变化的速度和位移曲线
subplot(2,1,2);plotyy(t,x(:
1),t,x(:
2),'plot');
functionxdot=ch2example6statefun(t,x,flag)
%乒乓球弹跳模型的标准状态方程
%x
(1)为小球速度,x
(2)为小球位移
xdot=zeros(2,1);%状态变量矩阵初始化
xdot
(1)=-9.8;%速度加速度方程
xdot
(2)=x
(1);%位移速度方程
%ch2example7prg1.m
clear;
N=5;%仿真计算的时间序列点数
f=[1,zeros(1,N-1)];%输入:
单位数字冲激信号
x=zeros(2,N+1);%状态变量存储矩阵初始化
x(:
1)=[0;0];%初始状态赋值
forn=1:
N%开始递推计算
x(1,n+1)=-0.5.*(x(1,n)+f(n));%状态方程1
x(2,n+1)=f(n);%状态方程2
y(n)=x(1,n)+x(2,n)+f(n);%输出方程
end
t_n=0:
N-1;%得到序列对应的离散时间点并作出波形
subplot(4,1,1);stem(t_n,f);axis([-1N01.5]);%输入信号波形
subplot(4,1,2);stem(t_n,x(1,1:
N));axis([-1N-0.60.6]);%状态1的波形
subplot(4,1,3);stem(t_n,x(2,1:
N));axis([-1N01.5]);%状态2的波形
subplot(4,1,4);stem(t_n,y);axis([-1N-0.51.3]);%输出信号波形
%ch2example10prg1.m
clear;
R=100;L=2e-3;C=1e-7;%设置电路元件的参数
ts=2e-6;
t_start=-1e-4;
t_end=4e-4;
t=t_start:
ts:
t_end;%设置求解的离散时间点序列
i_L0=0;u_C0=0;%系统初始状态为零
x0=[i_L0,u_C0];%系统状态变量初始赋值
tic
[t_out,x_out]=ode45('ch2example10statefun',t,x0,[],R,L,C);%仿真计算
toc
s_t_simu=x_out(:
2);%阶跃响应仿真结果
h_t_simu=x_out(:
1)./C;%等价的冲激响应仿真结果
figure
(1);plot(t_out,s_t_simu,'k-','LineWidth',3);
gridon;xlabel('时间t');ylabel('电容电压');
axis([t_start,t_end,1.1*min(s_t_simu),1.1*max(s_t_simu)]);
legend('单位阶跃响应仿真结果');
figure
(2);plot(t_out,h_t_simu,'r.','MarkerSize',20);
axis([t_start,t_end,1.1*min(h_t_simu),1.1*max(h_t_simu)]);
holdon;
%--理论结果
alfa=R/(2*L);
beta=sqrt(1/(L*C)-(R^2)/(4*L^2));
h_t=(L*C*beta)^(-1)*exp(-alfa*t).*sin(beta*t).*(t>=0);%冲激响应
plot(t,h_t,'k','LineWidth',2);legend('冲激响应仿真数值结果','冲激响应理论计算结果');
gridon;xlabel('时间t');ylabel('电容电压');
functionxdot=ch2example10statefun(t,x,flag,R,L,C)
%考虑RLC串联环路的状态方程函数
%输入:
t当前计算时刻,flag此处不用
%x为2#1矩阵,x
(1)为电感上的电流;x
(2)为电容电压
%R:
电阻值
%L:
电感值
%C:
电容值
xdot=zeros(2,1);%状态变量矩阵初始化
A=[-R/L,-1/L;1/C,0];%状态方程系数矩阵
B=[1/L;0];
xdot=A*x+B*f(t);%状态方程
functioninput=f(t)
input=(t>=0);%输入信号为单位阶跃
%ch2example11prg1.m
clear;
R=100;L=2e-3;C=1e-7;%设置电路元件的参数
ts=1e-6;%仿真步进1微秒
t_start=0;%起始时间
t_end=10e-3;%终止时间
t=t_start:
ts:
t_end;%设置求解的离散时间点序列
i_L0=0;u_C0=0;%系统初始状态为零
x0=[i_L0;u_C0];%系统状态变量初始赋值
inputtimespan=t;%输入信号的计算时间点
input=sin(2*pi*(2e3*t+3e6/2*(t.^2)));%扫频信号10ms内从2KH线性变化到32KH
%仿真计算开始
[t_out,x_out]=ode45('ch2example11statefun',t,x0,[],R,L,C,...
input,inputtimespan);
%仿真完毕,进入仿真后处理,显示仿真结果
s_t_simu=x_out(:
2);%阶跃响应仿真结果
figure
(1);plot(t_out,s_t_simu,'k-');
gridon;xlabel('时间t');ylabel('电容电压');
axis([0,t_end,1.1*min(s_t_simu),1.1*max(s_t_simu)]);
%计算扫频信号相应时间点对应的瞬时频率值,并以频率为横坐标作图
freq=2e3+3e6*t_out;
figure
(2);plot(freq,s_t_simu,'k-');holdon;
%计算幅度频率相应的理论曲线并作图与扫频输出结果相互对比
f1=0;df=1e3;f2=32e3;%计算的频率范围和步进
f=f1:
df:
f2;
w=2*pi*f;%相应的角频率序列
s=j*w;
H_s=1./(L*C*s.^2+R*C.*s+1);%代入传递函数计算
plot(f,abs(H_s),'ko-');%求幅度频率响应并作图
xlabel('频率Hz');ylabel('电容电压');
axis([0,max(freq),1.1*min(s_t_simu),1.1*max(s_t_simu)]);
legend('系统扫频响应仿真结果','系统幅度频率响应理论值');
functionxdot=ch2example11statefun(t,x,flag,R,L,C,input,inputtimespan)
%考虑RLC串联环路的状态方程函数
%输入:
t当前计算时刻,flag此处不用
%x为2#1矩阵,x
(1)为电感上的电流;x
(2)为电容电压
%R:
电阻值
%L:
电感值
%C:
电容值
%input为在时间点序列inputtimespan上给定输入信号
xdot=zeros(2,1);%状态变量矩阵初始化
A=[-R/L,-1/L;1/C,0];%状态方程系数矩阵
B=[1/L;0];
f_t=interp1(inputtimespan,input,t);
%利用插值来计算任意时刻的信号值
%因为求解器调用本函数时时间变量t是任意连续值
xdot=A*x+B*f_t;%状态方程
function[sys,x0,str,ts]=ch2example15Sfun(t,x,u,flag)
%连续系统状态方程;
%x'=Ax+Bu
%y=Cx+Du
%定义A,B,C,D矩阵
A=[010;001;-4-6-3];
B=[0;0;1];
C=[040];
D=0;
switchflag,
case0%flag=0初始化
[sys,x0,str,ts]=mdlInitializeSizes(A,B,C,D);
%可将A,B,C,D矩阵送入初始化函数
case1%flag=1计算连续系统状态方程(导数)
sys=mdlDerivatives(t,x,u,A,B,C,D);
case3%flag=3计算输出
sys=mdlOutputs(t,x,u,A,B,C,D);
case{2,4,9}%其他作不处理的flag
sys=[];%无用的flag时返回sys为空矩阵
otherwise%异常处理
error(['Unhandledflag=',num2str(flag)]);
end
%主函数结束
%子函数实现
(1)初始化函数----------------------------------
function[sys,x0,str,ts]=mdlInitializeSizes(A,B,C,D)%
sizes=simsizes;%获取
sizes.NumContStates=3;%连续系统的状态数为3
sizes.NumDiscStates=0;%离散系统的状态数,对于本系统此句可不用
sizes.NumOutputs=1;%输出信号数目是1
sizes.NumInputs=1;%输入信号数目是1
sizes.DirFeedthrough=0;%因为该系统不是直通的
sizes.NumSampleTimes=1;%这里必须为1
sys=simsizes(sizes);
str=[];%通常为空矩阵
x0=[0;0;0];%初始状态矩阵x0(零状态情况)
ts=[00];%表示连续取样时间的仿真
%初始化函数结束
%子函数实现
(2)系统状态方程函数-----------------------------
functionsys=mdlDerivatives(t,x,u,A,B,C,D)%系统状态方程函数
sys=A*x+B*u;%这里写入系统的状态方程矩阵形式即可
%系统状态函数结束
%子函数实现(3)系统输出方程函数---------------
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 通信系统建模与仿真例题源代码第二章 通信 系统 建模 仿真 例题 源代码 第二
![提示](https://static.bdocx.com/images/bang_tan.gif)