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;

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.m

dt=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.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的方波信号% 仿真开始, 注意: 设零时刻之前电路不工作, 系统状态保持不变

for k=1:length(t)

time=-T+k*dt;

if time>=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([-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.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;

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.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');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.m

clear;

v0=0; y0=1; % 球的初始状态

x_state=[v0,y0];% 将初始状态赋值到状态变量中

dt=0.01; % 仿真步进

t=0:dt:5; % 仿真时间序列

K=0.85; % 碰撞衰减系数

for k=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([-2 2 -0.1 1]); % 坐标范围固定

set(gcf,'DoubleBuffer','on'); % 双缓冲避免作图闪烁

drawnow; % 立即显示作图

end

% 仿真结束。最后输出计算时间序列上的结果：随时间变化的速度和位移曲线subplot(2,1,2);plotyy(t,x(:,1),t,x(:,2),'plot');

function xdot=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]; % 初始状态赋值

for n=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([-1 N 0 1.5]); % 输入信号波形subplot(4,1,2);stem(t_n,x(1,1:N));axis([-1 N -0.6 0.6]);% 状态1的波形subplot(4,1,3);stem(t_n,x(2,1:N));axis([-1 N 0 1.5]); % 状态2的波形

subplot(4,1,4);stem(t_n,y);axis([-1 N -0.5 1.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);

grid on; 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)]);

hold on;

%--理论结果

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('冲激响应仿真数值结果','冲激响应理论计算结果');

grid on; xlabel('时间 t ');ylabel('电容电压 ');

function xdot=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); % 状态方程

function input=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-');

grid on; 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-'); hold on;

% 计算幅度频率相应的理论曲线并作图与扫频输出结果相互对比

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('系统扫频响应仿真结果','系统幅度频率响应理论值');

function xdot=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=[0 1 0;0 0 1;-4 -6 -3];

B=[0;0;1];

C=[0 4 0];

D=0;

switch flag,

case 0 % flag=0 初始化

[sys,x0,str,ts]=mdlInitializeSizes(A,B,C,D);

% 可将A,B,C,D矩阵送入初始化函数

case 1 % flag=1 计算连续系统状态方程(导数)

sys=mdlDerivatives(t,x,u,A,B,C,D);

case 3 % flag=3 计算输出

sys=mdlOutputs(t,x,u,A,B,C,D);

case { 2, 4, 9 } % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',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 = [0 0]; % 表示连续取样时间的仿真

% 初始化函数结束

% 子函数实现（2）系统状态方程函数----------------------------- function sys=mdlDerivatives(t,x,u,A,B,C,D) % 系统状态方程函数sys = A*x+B*u ; % 这里写入系统的状态方程矩阵形式即可

% 系统状态函数结束

% 子函数实现（3）系统输出方程函数----------------------------- function sys = mdlOutputs(t,x,u,A,B,C,D)

sys = C*x; % 这里写入系统的输出方程矩阵形式即可

% 注意，如果使用语句sys = C*x+D*u ;代替上句，即使D=0，

% 也要将初始化函数中的sizes.DirFeedthrough 设为1

% 即系统存在输入输出之间的直通项, 否则执行将出现错误

% 系统输出方程函数结束

function [sys,x0,str,ts] =ch2example15SfunB(t,x,u,flag)

% 连续系统状态方程;

% x'=Ax+Bu

% y =Cx+Du

% 定义A,B,C,D矩阵

A=[0 1 0;0 0 1;-4 -6 -3];

B=[0;0;1];

C=[0 4 0];

D=0;

%------显示调用该S函数时，调用标志,仿真时间, 系统状态, 输入信号flag

t

x

u

% -----以下程序代码和ch2example15Sfun.m中的对应部分相同（略）switch flag,

case 0 % flag=0 初始化

[sys,x0,str,ts]=mdlInitializeSizes(A,B,C,D);

% 可将A,B,C,D矩阵送入初始化函数

case 1 % flag=1 计算连续系统状态方程(导数)

sys=mdlDerivatives(t,x,u,A,B,C,D);

case 3 % flag=3 计算输出

sys=mdlOutputs(t,x,u,A,B,C,D);

case { 2, 4, 9 } % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',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 = [0 0]; % 表示连续取样时间的仿真

% 初始化函数结束

% 子函数实现（2）系统状态方程函数----------------------------- function sys=mdlDerivatives(t,x,u,A,B,C,D) % 系统状态方程函数

sys = A*x+B*u ; % 这里写入系统的状态方程矩阵形式即可

% 系统状态函数结束

% 子函数实现（3）系统输出方程函数----------------------------- function sys = mdlOutputs(t,x,u,A,B,C,D)

sys = C*x; % 这里写入系统的输出方程矩阵形式即可

% 注意，如果使用语句sys = C*x+D*u ;代替上句，即使D=0，

% 也要将初始化函数中的sizes.DirFeedthrough 设为1

% 即系统存在输入输出之间的直通项, 否则执行将出现错误

% 系统输出方程函数结束

function [sys,x0,str,ts] =ch2example16Sfun(t,x,u,flag)

% 连续系统状态方程;

% x'=Ax+Bu

% y =Cx+Du

% 定义A,B,C,D矩阵

A=[0 1 0;0 0 1;-4 -6 -3];

B=[0;0;1];

C=[0 4 0];

D=0;

switch flag,

case 0 % flag=0 初始化

[sys,x0,str,ts]=mdlInitializeSizes(A,B,C,D);

% 可将A,B,C,D矩阵送入初始化函数

case 1 % flag=1 计算连续系统状态方程(导数)

sys=mdlDerivatives(t,x,u,A,B,C,D);

case 3 % flag=3 计算输出

sys=mdlOutputs(t,x,u,A,B,C,D);

case { 2, 4, 9 } % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',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 = [0 0]; % 表示连续取样时间的仿真

% 初始化函数结束

% 子函数实现（2）系统状态方程函数-----------------------------

function sys=mdlDerivatives(t,x,u,A,B,C,D) % 系统状态方程函数

sys = A*x+B*u ; % 这里写入系统的状态方程矩阵形式即可

% 系统状态函数结束

% 子函数实现（3）系统输出方程函数-----------------------------

function sys = mdlOutputs(t,x,u,A,B,C,D)

sys = C*x; % 这里写入系统的输出方程矩阵形式即可

% 注意，如果使用语句sys = C*x+D*u ;代替上句，即使D=0，

% 也要将初始化函数中的sizes.DirFeedthrough 设为1

% 即系统存在输入输出之间的直通项, 否则执行将出现错误

% 系统输出方程函数结束

function [sys,x0,str,ts] =ch2example17Sfun(t,x,u,flag,Amp,Freq,Phase)

% 正弦波信号源

switch flag,

case 0 % flag=0 初始化

[sys,x0,str,ts]=mdlInitializeSizes;

case 3 % flag=3 计算输出

sys=mdlOutputs(t,Amp,Freq,Phase);

case {1, 2, 4, 9 } % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',num2str(flag)]);

end

% 主函数结束

% 子函数实现（1）初始化函数----------------------------------

function [sys,x0,str,ts] = mdlInitializeSizes %

sizes = simsizes; % 获取SIMULINK仿真变量结构

sizes.NumContStates = 0; % 连续系统的状态数为0

sizes.NumDiscStates = 0; % 离散系统的状态数为0

sizes.NumOutputs = 1; % 输出信号数目是1

sizes.NumInputs = 0; % 输入信号数目是0

sizes.DirFeedthrough = 0; % 该系统不是直通的

sizes.NumSampleTimes = 1; % 这里必须为1

sys = simsizes(sizes);

str = []; % 通常为空矩阵

x0 = []; % 初始状态矩阵x0 （零状态情况）

ts = [0 0]; % 表示连续取样时间的仿真

% 初始化函数结束

% 子函数实现（2）系统输出方程函数-----------------------------

function sys = mdlOutputs(t,Amp,Freq,Phase)

sys = Amp*sin(2*pi*Freq*t+Phase);% 这里写入系统的输出方程矩阵形式即可% 修改这个函数可以得到任意的波形输出

% 系统输出方程函数结束

function [sys,x0,str,ts] =ch2example18Sfun(t,x,u,flag,numofshowsteps) % 简单的示波器

global T Y N; % 声明使用的全局变量

switch flag,

case 0 % flag=0 初始化

sizes = simsizes; % 获取SIMULINK仿真变量结构

sizes.NumContStates = 0; % 连续系统的状态数为0

sizes.NumDiscStates = 0; % 离散系统的状态数为0

sizes.NumOutputs = 0; % 输出信号数目是0

sizes.NumInputs = 1; % 输入信号数目是1

sizes.DirFeedthrough = 1; % 该系统是直通的

sizes.NumSampleTimes = 1; % 这里必须为1

sys = simsizes(sizes);

str = []; % 通常为空矩阵

x0 = []; % 初始状态矩阵x0 （零状态情况）

ts = [0 0]; % 表示连续取样时间的仿真

N=0; % 缓冲区全局变量初始化

T=zeros(1,numofshowsteps-1);

Y=zeros(1,numofshowsteps-1);

case 3 % flag=3 计算输出

if N < numofshowsteps-1 % 将输入暂存到缓冲区中

N=N+1; % 缓冲区数组跑标

Y(N)=u; % 记录当前信号

T(N)=t; % 记录当前时刻

else % 缓冲区满（一帧完成）则作图

figure(1); plot(T-min(T), Y); % 作出一帧信号波形

axis([0 max(T)-min(T) 1.1*min(Y) 1.1*max(Y)]);% 坐标范围

set(gcf,'DoubleBuffer','on'); % 双缓冲避免作图闪烁

drawnow;% 作图

N=0; % 缓冲区跑标复位

end

case {1, 2, 4, 9 } % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',num2str(flag)]);

end

function [sys,x0,str,ts] =ch2example18SfunB(t,x,u,flag,sampletime) % 简单的示波器2

% 输出信号u是一个取自缓存区的矢量

% 直接对u作图即可

switch flag,

case 0 % flag=0 初始化

sizes = simsizes; % 获取SIMULINK仿真变量结构

sizes.NumContStates = 0; % 连续系统的状态数为0

sizes.NumDiscStates = 0; % 离散系统的状态数为0

sizes.NumOutputs = 0; % 输出信号数目是0

sizes.NumInputs = -1; % 输入信号数目是自适应的

sizes.DirFeedthrough = 1; % 该系统是直通的

sizes.NumSampleTimes = 1; % 这里必须为1

sys = simsizes(sizes);

str = []; % 通常为空矩阵

x0 = []; % 初始状态矩阵x0 （零状态情况）

ts = [-1 0]; % 采样时间由外部模块给出

case 3 % flag=3 计算输出

% 计算缓存区信号的对应时间段长度，然后作图

T=0:sampletime:sampletime*(length(u)-1);

plot(T,u); % 作出一帧信号波形

axis([0 max(T) -5 5]); % 坐标范围

set(gcf,'DoubleBuffer','on'); % 双缓冲避免作图闪烁

drawnow;% 作图

case {1, 2, 4, 9 } % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',num2str(flag)]);

end

function [sys,x0,str,ts] =ch2example19Sfun(t,x,u,flag)

% 输入信号的代数运算实例：调幅

% 输出u是3行1列矩阵，u(1),u(2),u(3)分别表示基带信号，直流偏置和载波% size(u) % 取消本句注释可观察输入信号u的矩阵维数

switch flag,

case 0 % flag=0 初始化

sizes = simsizes; % 获取SIMULINK仿真变量结构

sizes.NumContStates = 0; % 连续系统的状态数是0

sizes.NumDiscStates = 0; % 离散系统的状态数是0

sizes.NumOutputs = 1; % 输出信号数目是1

sizes.NumInputs = -1; % 输入信号数目是自适应的

sizes.DirFeedthrough = 1; % 该系统是直通的

sizes.NumSampleTimes = 1; % 这里必须为1

sys = simsizes(sizes);

str = []; % 通常为空矩阵

x0 = []; % 初始状态矩阵x0

ts = [-1 0]; % 采样时间由外部模块给出

case 3 % flag=3 计算输出

sys=(u(1)+u(2))*u(3);% 调幅输出计算

case {1, 2, 4, 9 } % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',num2str(flag)]);

end

% ch2example20main.m

%执行以下命令，然后启动仿真

b=[2 1]; % H(z)的分子

a=[1 0.5 0.8]; % 分母

[A,B,C,D]=tf2ss(b,a);% 转换为状态方程

function [sys,x0,str,ts] =ch2example20Sfun(t,x,u,flag,b,a)

% 离散系统传递函数的S-函数实现

% 参数b,a 分别为H(z)的分母、分子多项式的系数向量

b=[2 1]; % H(z)的分子

a=[1 0.5 0.8]; % 分母

[A,B,C,D]=tf2ss(b,a);% 将H(z)转换为状态空间方程系数矩阵

switch flag,

case 0 % flag=0 初始化

sizes = simsizes;% 获取SIMULINK仿真变量结构

sizes.NumContStates = 0;% 连续系统的状态数是0

sizes.NumDiscStates = size(A,1);% 离散系统的状态数

sizes.NumOutputs = size(D,1);% 输出信号数目

sizes.NumInputs = size(D,2);% 输入信号数目

sizes.DirFeedthrough = 1; % 该系统是直通的

sizes.NumSampleTimes = 1; % 这里必须为1

sys = simsizes(sizes);

str = []; % 通常为空矩阵

x0 = zeros(sizes.NumDiscStates,1);% 零状态

ts = [-1 0]; % 采样时间由外部模块给出

case 2 % flag=2 离散状态方程计算

sys = A*x+B*u;

case 3 % flag=3 输出方程计算

sys = C*x+D*u;

case {1,4, 9 }% 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',num2str(flag)]);

end

function [sys,x0,str,ts] =ch2example21Sfun(t,x,u,flag,threshold) % 两路切换开关的实现

% u(2)为切换控制信号

switch flag,

case 0 % flag=0 初始化

sizes = simsizes; % 获取SIMULINK仿真变量结构

sizes.NumContStates = 0; % 连续系统的状态数是0

sizes.NumDiscStates = 0; % 离散系统的状态数是0

sizes.NumOutputs = 1; % 输出信号数目是1

sizes.NumInputs = -1; % 输入信号数目是自适应的

sizes.DirFeedthrough = 1; % 该系统是直通的

sizes.NumSampleTimes = 1; % 这里必须为1

sys = simsizes(sizes);

str = []; % 通常为空矩阵

x0 = []; % 初始状态矩阵x0

ts = [-1 0]; % 采样时间由外部模块给出

case 3 % flag=3 计算输出

sys=u(1)*(u(2)>threshold)+u(3)*(u(2)<=threshold);% 选择输出

case {1, 2, 4, 9 } % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',num2str(flag)]);

end

function [sys,x0,str,ts] =ch2example22Sfun(t,x,u,flag,Amp,f0, k, phi0) % VCO的实现

switch flag,

case 0 % flag=0 初始化

sizes = simsizes; % 获取SIMULINK仿真变量结构

sizes.NumContStates = 1; % 连续系统的状态数是1

sizes.NumDiscStates = 0; % 离散系统的状态数是0

sizes.NumOutputs = 1; % 输出信号数目是1

sizes.NumInputs = -1; % 输入信号数目是自适应的

sizes.DirFeedthrough = 0; % 该系统不是直通的

sizes.NumSampleTimes = 1; % 这里必须为1

sys = simsizes(sizes);

str = []; % 通常为空矩阵

x0 = [0]; % 初始状态矩阵x0

ts = [-1 0]; % 采样时间由外部模块给出

case 1 % flag=1 计算连续系统状态方程

sys=2*pi*f0+2*pi*k*u; % VCO 状态方程

case 3 % flag=3 计算输出

sys= Amp*sin(x+phi0); % VCO 输出方程

case {2,4, 9 } % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',num2str(flag)]);

end

function [sys,x0,str,ts] =ch2example23Sfun(t,x,u,flag,v0, theta0,L,m,g,k)

% 单摆的实现

switch flag,

case 0 % flag=0 初始化

sizes = simsizes; % 获取SIMULINK仿真变量结构

sizes.NumContStates = 2; % 连续系统的状态数是2

sizes.NumDiscStates = 0; % 离散系统的状态数是0

sizes.NumOutputs = 1; % 输出信号数目是1

sizes.NumInputs = 0; % 输入信号数目是0

sizes.DirFeedthrough = 0; % 该系统不是直通的

sizes.NumSampleTimes = 1; % 这里必须为1

sys = simsizes(sizes);

str = []; % 通常为空矩阵

x0 = [v0, theta0]; % 初始状态矩阵x0

ts = [0 0]; % 连续系统

case 1 % flag=1 计算连续系统状态方程

vdot=g*sin(x(2))-k/m*x(1);

thetadot=-1/L*x(1);

sys=[vdot; thetadot]; %状态方程

case 3 % flag=3 计算输出

sys=x(2) ; % 输出方程

case {2,4,9} % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵

otherwise % 异常处理

error(['Unhandled flag = ',num2str(flag)]);

end

function

[sys,x0,str,ts]=ch2example23Sfunshow(t,x,u,flag,v0,theta0,L,m,g,k) % 单摆的动画显示

switch flag,

case 0 % flag=0 初始化

sizes = simsizes; % 获取SIMULINK仿真变量结构

sizes.NumContStates = 0; % 连续系统的状态数是0

sizes.NumDiscStates = 0; % 离散系统的状态数是0

sizes.NumOutputs = 0; % 输出信号数目是0

sizes.NumInputs = 1; % 输入信号数目是0

sizes.DirFeedthrough = 1; % 该系统是直通的

sizes.NumSampleTimes = 1; % 这里必须为1

sys = simsizes(sizes);

str = []; % 通常为空矩阵

x0 = []; % 初始状态矩阵x0

ts = [-1 0];

case 3 % flag=3 计算输出

sita=u;

x=L*sin(sita); % 计算摆锤的位置

y=-L*cos(sita);

plot([0,x],[0,y],'-o');% 动画作图

text(.5,-.5,['t=',num2str(t)]);% 显示参数

text(.5,-.6,['L=',num2str(L)]);

text(.5,-.7,['m=',num2str(m)]);

text(.5,-.8,['k=',num2str(k)]);

axis equal; % 坐标纵横比例相同

axis([-L L -L L]);% 坐标范围固定

set(gcf,'DoubleBuffer','on');% 双缓冲避免作图闪烁 drawnow; % 立即执行作图命令 case {1,2,4,9} % 其他作不处理的flag

sys=[]; % 无用的flag时返回sys为空矩阵otherwise % 异常处理

error(['Unhandled flag = ',num2str(flag)]);

end

clear

t=[0;1;2];

f1=[1;2;3];

f2=[4;2;5];

simin=[t,f1,f2];

mysimopts=simset('Solver','ode5');

mysimopts=simset(mysimopts,'FixedStep',0.5);

timespan=[0,4];

sim('ch2example24',timespan,mysimopts);

load chirp

wavplay(y,Fs,'sync')

load gong

wavplay(y,Fs)

pause

load chirp

wavplay(y,Fs,'async')

load gong

wavplay(y,Fs)

% ch2example25prg1.m

% 采用较低速率仿真，并用样条插值使得结果精细化

simurate=3000; % 较低仿真速率

mysimopts = simset('Solver','ode5'); % 设置仿真求解器为ode5的mysimopts = simset(mysimopts,'FixedStep',1./simurate); % 设置仿真步进timespan=[0, 0.02]; % 仿真时间段参数

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('较低速率仿真输出波形','样条插值精细化的结果','较高速率仿真输出波形');

function [t_out, x_out]=eulerode(odefuncname,tspan, x0, options, par) % 欧拉法求状态方程odefuncname的数值解

% odefuncname为描述状态方程函数名称（也即文件名）

% tspan为求解的时间序列（要求是列向量）

% x0为系统初始状态列向量

% options 此处不使用

% par 是传入odefuncname的参数，与odefuncname 函数的参数par设计相匹配

% 返回值：t_out 计算的时间点序列, 等于输入的tspan

% x_out是对应于t_out的状态变量解矩阵

t_out=tspan;

x_out=zeros(length(t_out),length(x0));% 初始化输出解矩阵, 以加快程序计算速度

% x_out的第k行状态值对应着t_out中的第k个时刻

for k=1:length(t_out)

x_out(k,:)=x0'; % 将第k步计算结果存入输出解矩阵(注意要转置)

if k

h=t_out(k+1)-t_out(k); % 到下一时刻的计算步长

else

h=0; % 最后一个步长

end

x0=x0+h*eval([odefuncname, '(t_out(k),x0,[], par)']);% 欧拉算法 % 用eval调用状态方程函数, 其参数flag不使用, 置为空矩阵或任意数值End

function xdot=pendulumstateeq(t, x, flag, par)

% 考虑空气阻力的单摆系统的状态方程函数

% 输入: t 当前计算时刻, flag 此处不用

% x为2#1矩阵, x(1)为当前摆锤的切向速度; x(2)为当前角位移

% par为4#1系统参数矩阵

% par(1) 为重力加速度参数 g

% par(2) 为空气阻力比例系数 k

% par(3) 为摆锤质量 m

% par(4) 为摆线长度 L

xdot=zeros(2,1); % 状态变量矩阵初始化

g=par(1);k=par(2);m=par(3);L=par(4);% 系统参数设置

xdot(1)=g*sin(x(2))-k/m*x(1); % 状态方程组

xdot(2)=-1/L*x(1);

function xdot=pendulumstateeqlinear(t, x, flag, par)

% 考虑空气阻力的单摆系统的状态方程函数: 近似线性模型

% 输入: t 当前计算时刻, flag 此处不用

% x为2#1矩阵, x(1)为当前摆锤的切向速度; x(2)为当前角位移

% par为4#1系统参数矩阵

% par(1) 为重力加速度参数 g

% par(2) 为空气阻力比例系数 k

% par(3) 为摆锤质量 m

% par(4) 为摆线长度 L

xdot=zeros(2,1); % 状态变量矩阵初始化

g=par(1);k=par(2);m=par(3);L=par(4);% 系统参数设置

A=[ -k/m, g; -1/L, 0]; % 状态矩阵赋值

B=0;

f=0; % 无输入信号

xdot=A*x+B*f; % 通用的状态方程

function xdot=rcstateequation(t, x, flag, par)

现代信号处理Matlab仿真——例611

例6.11 利用卡尔曼滤波估计一个未知常数 题目： 设已知一个未知常数x 的噪声观测集合，已知噪声v(n)的均值为零， 方差为 ，v(n)与x 不相关，试用卡尔曼滤波估计该常数 题目分析： 回忆Kalman 递推估计公式 由于已知x 为一常数，即不随时间n 变化，因此可以得到： 状态方程： x(n)=x(n-1) 观测方程： y(n)=x(n)+v(n) 得到A(n)=1,C(n)=1, , 将A(n)=1,代入迭代公式 得到：P(n|n-1)=P(n-1|n-1) 用P(n-1)来表示P(n|n-1)和P(n-1|n-1),这是卡尔曼增益表达式变为 从而 2v σ1??(|1)(1)(1|1)(|1)(1)(1|1)(1)()()(|1)()[()(|1)()()]???(|)(|1)()[()()(|1)](|)[()()](|1)H w H H v x n n A n x n n P n n A n P n n A n Q n K n P n n C n C n P n n C n Q n x n n x n n K n y n C n x n n P n n I K n C n P n n --=----=----+=--+=-+--=--2()v v Q n σ=()0w Q n =(|1)(1)(1|1)(1)()H w P n n A n P n n A n Q n -=----+21 ()(|1)[(|1)]v K n P n n P n n σ-=--+22(1)()[1()](1)(1)v v P n P n K n P n P n σσ-=--=-+

matlab课后习题解答第二章

第2章符号运算 习题2及解答 1 说出以下四条指令产生的结果各属于哪种数据类型，是“双精度” 对象，还是“符号”符号对象 3/7+; sym(3/7+; sym('3/7+'); vpa(sym(3/7+) 〖目的〗 不能从显示形式判断数据类型，而必须依靠class指令。 〖解答〗 c1=3/7+ c2=sym(3/7+ c3=sym('3/7+') c4=vpa(sym(3/7+) Cs1=class(c1) Cs2=class(c2) Cs3=class(c3) Cs4=class(c4) c1 = c2 = 37/70 c3 = c4 = Cs1 = double Cs2 = sym Cs3 = sym Cs4 = sym 2 在不加专门指定的情况下，以下符号表达式中的哪一个变量被认 为是自由符号变量. sym('sin(w*t)'),sym('a*exp(-X)'),sym('z*exp(j*th)') 〖目的〗 理解自由符号变量的确认规则。 〖解答〗 symvar(sym('sin(w*t)'),1) ans = w symvar(sym('a*exp(-X)'),1)

ans = a symvar(sym('z*exp(j*th)'),1) ans = z 5求符号矩阵???? ??????=3332 31 232221 131211 a a a a a a a a a A 的行列式值和逆，所得结果应采用“子表达式置换”简洁化。 〖目的〗 理解subexpr 指令。 〖解答〗 A=sym('[a11 a12 a13;a21 a22 a23;a31 a32 a33]') DA=det(A) IA=inv(A); [IAs,d]=subexpr(IA,d) A = [ a11, a12, a13] [ a21, a22, a23] [ a31, a32, a33] DA = a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31 IAs = [ d*(a22*a33 - a23*a32), -d*(a12*a33 - a13*a32), d*(a12*a23 - a13*a22)] [ -d*(a21*a33 - a23*a31), d*(a11*a33 - a13*a31), -d*(a11*a23 - a13*a21)] [ d*(a21*a32 - a22*a31), -d*(a11*a32 - a12*a31), d*(a11*a22 - a12*a21)] d = 1/(a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31) 8（1）通过符号计算求t t y sin )(=的导数 dt dy 。（2）然后根据此结果，求- =0t dt dy 和2 π = t dt dy 。 〖目的〗 diff, limit 指令的应用。 如何理解运行结果。 〖解答〗 syms t

Matlab仿真实例-卫星轨迹

卫星轨迹 一．问题提出 设卫星在空中运行的运动方程为： 其中是k 重力系数(k=401408km3/s)。卫星轨道采用极坐标表示，通过仿真，研究发射速度对卫星轨道的影响。实验将作出卫星在地球表面（r=6400KM ，θ=0）分别以v=8KM/s,v=10KM/s,v=12KM/s 发射时，卫星绕地球运行的轨迹。 二．问题分析 1．卫星运动方程一个二阶微分方程组，应用Matlab 的常微分方程求解命令ode45求解时，首先需要将二阶微分方程组转换成一阶微分方程组。若设，则有： 2．建立极坐标如上图所示，初值分别为：卫星径向初始位置，即地球半径：y(1,1)=6400；卫星初始角度位置：y(2,1)=0；卫星初始径向线速度：y(3,1)=0；卫星初始周向角速度：y(4,1)=v/6400。 3．将上述一阶微分方程及其初值带入常微分方程求解命令ode45求解，可得到一定时间间隔的卫星的径向坐标值y(1)向量；周向角度坐标值y(2)向量；径向线速度y(3)向量；周向角速度y(4)向量。 4．通过以上步骤所求得的是极坐标下的解，若需要在直角坐标系下绘制卫星的运动轨迹，还需要进行坐标变换，将径向坐标值y(1)向量；周向角度坐标值y(2)向量通过以下方程转换为直角坐标下的横纵坐标值X,Y 。 5．卫星发射速度速度的不同将导致卫星的运动轨迹不同，实验将绘制卫星分别以v=8KM/s ，v=10KM/s ，v=12KM/s 的初速度发射的运动轨迹。 三．Matlab 程序及注释 1．主程序 v=input('请输入卫星发射速度单位Km/s ：\nv=');%卫星发射速度输入。 axis([-264007000-1000042400]);%定制图形输出坐标范围。 %为了直观表达卫星轨迹，以下语句将绘制三维地球。 [x1,y1,z1]=sphere(15);%绘制单位球。 x1=x1*6400;y1=y1*6400;???????-=+-=dt d dt dr r dt d dt d r r k dt r d θ θθ2)(2 22222θ==)2(,)1(y r y ?????????????**-=**+*-===)1(/)4()3(2)4()4()4()1()1()1()3()4()2() 3()1(y y y dt dy y y y y y k dt dy y dt dy y dt dy ???*=*=)] 2(sin[)1(Y )]2(cos[)1(X y y y y

matlab课后答案完整版

ones表示1矩阵 zeros表示0矩阵 ones(4)表示4x4的1矩阵 zeros(4)表示4x4的0矩阵 zeros(4,5)表示4x5的矩阵 eye(10,10)表示10x10的单位矩阵rand(4,5)表示4x5的伴随矩阵 det(a)表示计算a的行列式 inv(a)表示计算a的逆矩阵 Jordan(a)表示求a矩阵的约当标准块rank(a)表示求矩阵a的秩 [v,d]=eig(a)对角矩阵 b=a’表示求a矩阵的转置矩阵 sqrt表示求平方根 exp表示自然指数函数 log自然对数函数 abs绝对值 第一章 一、5（1） b=[97 67 34 10;-78 75 65 5;32 5 -23 -59]; >> c=[97 67;-78 75;32 5;0 -12]; >> d=[65 5;-23 -59;54 7]; >> e=b*c e = 5271 11574 -11336 664 1978 3112 （2）a=50:1:100 二、1 、x=-74; y=-27; z=(sin(x.^2+y.^2))/(sqrt(tan(abs(x+y)))+pi) z = 2、a=::; >> b=exp*a).*sin(a+ 3、x=[2 4; 5]; y=log(x+sqrt(1+x.^2))/2 y =4、a*b表示a矩阵和b矩阵相乘 a.*b表示a矩阵和b矩阵单个元素相乘A(m,n)表示取a矩阵第m行，第n列 A(m,:)表示取a矩阵第m行的全部元素 A(:,n)表示取a矩阵的第n列全部元素 A./B表示a矩阵除以b矩阵的对应元素， B.\A等价于A./B A.^B表示两个矩阵对应元素进行乘方运算A.^2表示a中的每个元素的平方 A^2表示A*A 例：x=[1,2,3]; y=[4,5,6]; z=x.^y z= 1 3 2 729 指数可以是标量（如y=2）.底数也可以是标量（如x=2） 5、a=1+2i; >> b=3+4i; >> c=exp((pi*i)/6) c = + d=c+a*b/(a+b) d = + 第二章 二、4、（1） y=0;k=0; >> while y<3 k=k+1; y=y+1/(2*k-1); end >> display([k-1,y-1/(2*k-1)]) ans = 第三章 二1（1）

matlab第二章习题答案

第一大题： （1） a = 7/3 b = sym(7/3) c = sym(7/3,'d') d = sym('7/3') v1=vpa(abs(a-d)) v2=vpa(abs(b-d)) v3=vpa(abs(c-d)) a = 2.3333 b = 7/3 c = 2.3333333333333334813630699500209 d = 7/3 v1 = 0.0 v2 = 0.0 v3 = 0.00000000000000014802973661668756666666667788716（2） a = pi/3 b = sym(pi/3) c = sym(pi/3,'d') d = sym('pi/3') v1=vpa(abs(a-d)) v2=vpa(abs(b-d)) v3=vpa(abs(c-d)) a = 1.0472 b = pi/3 c = 1.047197551196597631317786181171 d = pi/3 v1 = 0.0 v2 = 0.0 v3 = 0.00000000000000011483642827992216762806615818554 (3) a = pi*3^(1/3) b = sym(pi*3^(1/3)) c = sym(pi*3^(1/3),'d') d = sym('pi*3^(1/3)') v1=vpa(abs(a-d)) v2=vpa(abs(b-d)) v3=vpa(abs(c-d))

a = 4.5310 b = 1275352044764433/281474976710656 c = 4.5309606547207899041040946030989 d = pi*3^(1/3) v1 = 0.00000000000000026601114166290944374842393221638 v2 = 0.00000000000000026601114166290944374842393221638 v3 = 0.0000000000000002660111416629094726767991785515 第二大题： （1） c1=3/7+0.1 c1 = 0.5286 双精度 （2） c2=sym(3/7+0.1) c2 = 37/70 符号 （3） c3=vpa(sym(3/7+0.1)) c3 = 0.52857142857142857142857142857143 完整显示精度 第三大题： （1） findsym(sym('sin(w*t)'),1) ans = w (2) findsym(sym('a*exp(-X)' ) ,1) ans = a (3) findsym(sym('z*exp(j*theta)'),1) ans = z

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 on plot(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;%抽样时间间隔要足够小，要满足抽样定理。 ts=1/fs; %%根据抽样时间间隔进行抽样

MATLAB基础教程 薛山第二版 课后习题答案

《MATLAB及应用》实验指导书《MATLAB及应用》实验指导书 班级：T1243-7 姓名：柏元强 学号：20120430724 总评成绩： 汽车工程学院 电测与汽车数字应用中心

目录 实验04051001 MATLAB语言基础 (1) 实验04051002 MATLAB科学计算及绘图 (18) 实验04051003 MATLAB综合实例编程 (31)

实验04051001 MATLAB语言基础 1实验目的 1)熟悉MATLAB的运行环境 2)掌握MATLAB的矩阵和数组的运算 3)掌握MATLAB符号表达式的创建 4)熟悉符号方程的求解 2实验内容 第二章 1.创建double的变量，并进行计算。 (1)a=87，b=190，计算 a+b、a-b、a*b。 clear,clc a=double(87); b=double(190); a+b,a-b,a*b (2)创建 uint8 类型的变量，数值与(1)中相同，进行相同的计算。 clear,clc a=uint8(87); b=uint8(190); a+b,a-b,a*b 2．计算：

(1) () sin 60 (2) e3 (3) 3cos 4??π ??? clear,clc a=sind(60) b=exp(3) c=cos(3*pi/4) 3．设2u =，3v =，计算： (1) 4 log uv v (2) () 2 2 e u v v u +- (3) clear,clc u=2;v=3; a=(4*u*v)/log(v) b=((exp(u)+v)^2)/(v^2-u) c=(sqrt(u-3*v))/(u*v) 4．计算如下表达式： (1) ()() 3542i i -+ (2) () sin 28i - clear,clc (3-5*i)*(4+2*i) sin(2-8*i)

MATLAB运算基础(第2章)答案

实验01讲评、参考答案 讲评 未交实验报告的同学名单 批改情况： 问题1： 不仔细，式子中出错。 问题2： 提交的过程不完整。 问题3： 使用语句尾分号(;)不当，提交的过程中不该显示的结果显示。 问题4： 截屏窗口没有调整大小。

附参考答案： 实验01 MATLAB 运算基础 （第2章 MATLAB 数据及其运算） 一、实验目的 1. 熟悉启动和退出MATLAB 的方法。 2. 熟悉MATLAB 命令窗口的组成。 3. 掌握建立矩阵的方法。 4. 掌握MATLAB 各种表达式的书写规则以及常用函数的使用。 二、实验内容 1. 数学表达式计算 先求下列表达式的值，然后显示MATLAB 工作空间的使用情况并保存全部变量。 1.1 计算三角函数 12 2sin 851z e =+（注意：度要转换成弧度，e 2如何给出） 示例：点击Command Window 窗口右上角的，将命令窗口提出来成悬浮窗口，适 当调整窗口大小。 命令窗口中的执行过程： 《MATLAB 软件》课内实验 王平

1.2 计算自然对数 221 ln(1)2z x x =++，其中2120.45 5i x +??=??-??（提示：clc 命令擦除命令窗口，clear 则清除工作空间中的所有变量，使用时注意区别，慎用clear 命令。 应用点乘方） 命令窗口中的执行过程： 1.3 求数学表达式的一组值 0.30.330.3sin(0.3)ln , 3.0, 2.9,,2.9,3.022 a a e e a z a a --+=++=--L 提示：利用冒号表达式生成a 向量，求各点的函数值时用点乘运算。 命令窗口中的执行过程：

matlab答案

第一章 2.在命令窗口输入‘’w=3+2‘’，然后依次使用clear和clc命令，分别观察命令窗口、工作空间窗口和历史命令窗口的变化。 使用clear命令时，命令窗口无变化，工作空间窗口中的内容被删除，历史命令窗口多出一条命令记录。 使用clc命令时，命令窗口中的内容被删除，工作空间窗口无变化，历史命令窗口中多出一条命令记录 3.将硬盘上一已有目录，加入到搜索路径，并将其设置为当前工作目录。 File-set path-add folder-save 第二章 1.计算复数3+4i与5-6i的乘积。 a=3+4i b=5-6i c=a*b 2.构建结构体Students，属性包含Name、age和Email，数据包括｛’Zhang’，18，*‘Zhang@16 https://www.sodocs.net/doc/a64506194.html,’，’Zhang@https://www.sodocs.net/doc/a64506194.html,’+｝、｛’Wang’，21，[]｝和｛’Li’，[],[]｝，构建后读取所有Name属性值，并且修改’Zhang’的Age属性值为19。 Students(1).Name='Zhang' Students(1).Age=18 Students(1).Email='Zhang@https://www.sodocs.net/doc/a64506194.html,','Zhang@https://www.sodocs.net/doc/a64506194.html,' Students(2).Name='Wang' Students(2).Age=21 Students(2).Email=[] Students(3).Name='Li' Students(3).Age=[] Students(3).Email=[] https://www.sodocs.net/doc/a64506194.html, Student(1).Age=19

Student(1).Age 3.用满矩阵和稀疏矩阵存储方式分别构造下属矩阵： A=[0 1 0 0 0;1 0 0 0 0;0 0 0 0 0;0 0 0 1 0] ; F=full(A) S=sparse(A) S=sparse([2,1,4],[1,2,4],[1,1,1],4,5) 4.采用向量构造符得到向量[1,5,9....,41]. A=1:4:41 5.按水平和竖直方向分别合并下述两个矩阵： A=[1 0 0;1 1 0;0 0 1],B=[2 3 4;5 6 7;8 9 10] A=[1 0 0;1 1 0;0 0 1] ; B=[2 3 4;5 6 7;8 9 10] ; C=[A B] D=[A;B] 6.分别删除第五题两个结果的第2行。 A=[1 0 0;1 1 0;0 0 1] B=[2 3 4;5 6 7;8 9 10] C=[A B] D=[A;B] C(2,:)=[] D(2,:)=[] 7.分别将第5题两个结果的第2行最后3列的数值改为[11 12 13]。 A=[1 0 0;1 1 0;0 0 1] B=[2 3 4;5 6 7;8 9 10] C=[A B] D=[A;B] C(2,4:6)=[11 12 13] D(2,:)=[11 12 13] 8.分别查看第5题两个结果的各方向长度

matlab课后答案

第二章习题 1．观察一个数（在此用@记述）在以下四条不同指令作用下的异同： a=@ b=sym(@) c=sym(@,’d’) d=sym(‘@’) 在此，@分别代表具体数值7/3，pi/3，pi*3^(1/3)；而异同通过vpa(abs(a-d)),vpa(abs(b-d)),vpa(abs(c-d))等来观察。 a=7/3 b=sym(7/3) c=sym(7/3,'d') d=sym('7/3') vpa(abs(a-d)) vpa(abs(a-b)) vpa(abs(c-d)) a = 2.3333 b = 7/3 c = 2.3333333333333334813630699500209 d = 7/3 ans = 0.0 ans = 0.0 ans = 0.00000000000000014802973661668756666666667788716 a=pi/3 b=sym(pi/3) c=sym(pi/3,'d') d=sym('pi/3') vpa(abs(a-d)) vpa(abs(a-b)) vpa(abs(c-d)) a = 1.0472 b = pi/3 c = 1.047197551196597631317786181171

d = pi/3 ans = 0.0 ans = 0.0 ans = 0.00000000000000011483642827992216762806615818554 2．说出以下三条指令产生的结果各属于那种数据类型，是“双精度”对象，还是“符号”对象？ 3/7+0.1，sym(3/7+0.1),vpa(sym(3/7+0.1)) a=3/7+0.1 b=sym(3/7+0.1) c=vpa(sym(3/7+0.1)) class(a) class(b) class(c) a = 0.5286 b = 37/70 c = 0.52857142857142857142857142857143 ans = double ans = sym ans = sym 3．在不加专门指定的情况下，以下符号表达中的那一个变量是独立自由变量。 sym(‘sin(w*t)’),sym(‘a*exp(-x)’),sym(‘z*exp(j*theta)’) syms t w x a j z theta a=sin(w*t) b=a*exp(-x) c=z*exp(j*theta) findsym(a) findsym(b) findsym(c) a = sin(t*w) b =

第2章MATLAB数据及其运算习题答案.doc

第2章 MATLAB数据及其运算 习题2 一、选择题 1．下列可作为MA TLAB合法变量名的是（）。D A．合计B．123 C．@h D．xyz_2a 2．下列数值数据表示中错误的是（）。C A．+10 B．1.2e-5 C．2e D．2i 3．使用语句t=0:7生成的是（）个元素的向量。A A．8 B．7 C．6 D．5 4．执行语句A=[1,2,3;4,5,6]后，A(3)的值是（）。B A．1 B．2 C．3 D．4 5．已知a为3×3矩阵，则a(:,end)是指（）。D A．所有元素B．第一行元素 C．第三行元素D．第三列元素 6．已知a为3×3矩阵，则运行a (1)=[]后（）。A A．a变成行向量B．a变为2行2列 C．a变为3行2列D．a变为2行3列 7．在命令行窗口输入下列命令后，x的值是（）。B >> clear >> x=i*j A．不确定B．-1 C．1D．i*j 8．fix(354/100)+mod(354,10)*10的值是（）。D A．34 B．354 C．453D．43 9．下列语句中错误的是（）。B A．x==y==3 B．x=y=3 C．x=y==3 D．y=3,x=y 10．find(1:2:20>15)的结果是（）。C A．19 20 B．17 19 C．9 10 D．8 9 11．输入字符串时，要用（）将字符括起来。C A．[ ] B．{ } C．' ' D．" " 12．已知s='显示"hello"'，则s的元素个数是（）。A A．9 B．11 C．7 D．18

13．eval('sqrt(4)+2')的值是（）。B A．sqrt(4)+2 B．4 C．2 D．2，2 14．有3×4的结构矩阵student，每个结构有name（姓名）、scores（分数）两个成员，其中scores是以1×5矩阵表示的5门课的成绩，那么要删除第4个学生的第2门课成绩，应采用的正确命令是（）。D A．rmfield(student(1,2).scores(2)) B．rmfield(student(4).scores) C．student(4).scores(2)=0 D．student(1,2).scores(2)=[] 15．有一个2行2列的单元矩阵c，则c(2)是指（）。B A．第一行第二列的元素内容B．第二行第一列的元素内容 C．第一行第二列的元素D．第二行第一列的元素 二、填空题 1．从键盘直接输入矩阵元素来建立矩阵时，将矩阵的元素用括起来，按矩阵行的顺序输入各元素，同一行的各元素之间用分隔，不同行的元素之间用分隔。中括号，逗号或空格，分号 2．设A=[1,2;3,4]，B=[5,6;7,8]，则A*B= ，A.*B= 。 A*B=[19,22;43,50]，A.*B=[5,12;21,32] 3．有3×3矩阵，求其第5个元素的下标的命令是，求其第三行、第三列元素的序号的命令是。[i,j]=ind2sub([3 3],5)，ind=sub2ind([3 3],3,3) 4．下列命令执行后的输出结果是。20 >> ans=5; >> 10; >> ans+10 5．下列命令执行后，new_claim的值是。This is a great example. claim= 'This is a good example.'; new_claim=strrep(claim,'good','great') 三、应用题 1．命令X=[]与clear X有何不同？请上机验证结论。 Clear X是将X从工作空间中删除，而X=[]是给X赋空矩阵。空矩阵存在于工作空间中，只是没有任何元素。 2．在一个MA TLAB命令中，6 + 7i和6 + 7*i有何区别？i和I有何区别？ 3．设A和B是两个同大小的矩阵，试分析A*B和A.*B、A./B和B.\A、A/B和B\A 的区别？如果A和B是两个标量数据，结论又如何？ 4．写出完成下列操作的命令。 （1）将矩阵A第2～5行中第1，3，5列元素赋给矩阵B。 （2）删除矩阵A的第5号元素。 （3）求矩阵A的大小和维数。 2

MATLAB 实用教程 课后习题答案

第二章 1.计算复数3+4i与5-6i的乘积。 a=3+4i b=5-6i c=a*b 2.构建结构体Students，属性包含Name、age和Email，数据包括｛’Zhang’，18，*‘Zhang@16 https://www.sodocs.net/doc/a64506194.html,’，’Zhang@https://www.sodocs.net/doc/a64506194.html,’+｝、｛’Wang’，21，[]｝和｛’Li’，[],[]｝，构建后读取所有Name属性值，并且修改’Zhang’的Age属性值为19。 Students(1).Age=18 Students(1).Email='zhang@https://www.sodocs.net/doc/a64506194.html,','zhang@https://www.sodocs.net/doc/a64506194.html,' Students(2).Name='Wang' Students(2).Age=21 Students(2).Email=[] Students(3).Name='Li' Students(3).Age=[] Students(3).Email=[] https://www.sodocs.net/doc/a64506194.html, Student(1).Age(1)=19 Student.Age 3.用满矩阵和稀疏矩阵存储方式分别构造下属矩阵： A=[0 1 0 0 0;1 0 0 0 0;0 0 0 0 0;0 0 0 1 0] A=[0 1 0 0 0;1 0 0 0 0;0 0 0 0 0;0 0 0 1 0] S=sparse(A) S=sparse([2,1,4],[1,2,4],[1,1,1],4,5) 4.采用向量构造符得到向量[1,5,9....,41]. A=1:4:41 5.按水平和竖直方向分别合并下述两个矩阵：A=[1 0 0;1 1 0;0 0 1],B=[2 3 4;5 6 7;8 9 10] A=[1 0 0;1 1 0;0 0 1] B=[2 3 4;5 6 7;8 9 10] C=[A B] D=[A;B] 6.分别删除第五题两个结果的第2行。 A=[1 0 0;1 1 0;0 0 1] B=[2 3 4;5 6 7;8 9 10] C=[A B] D=[A;B] C(2,:)=[] D(2,:)=[] 7.分别将第5题两个结果的第2行最后3列的数值改为[11 12 13]。 A=[1 0 0;1 1 0;0 0 1] B=[2 3 4;5 6 7;8 9 10] C=[A B] D=[A;B]

计算方法及其MATLAB实现第二章作业

作者：夏云木子 1、 >> syms re(x) re(y) re(z) >> input('计算相对误差：'),re(x)=10/1991,re(y)=0.0001/1.991,re(y)=0.0000001/0.0001991 所以可知re(y)最小，即y精度最高 2、 >> format short,A=sqrt(2) >> format short e,B=sqrt(2) >> format short g,C=sqrt(2)

>> format long,D=sqrt(2) >> format long e,E=sqrt(2) >> format long g,F=sqrt(2) >> format bank,H=sqrt(2) >> format hex,I=sqrt(2) >> format +,J=sqrt(2) >> format,K=sqrt(2)

3、 >> syms A >> A=[sqrt(3) exp(7);sin(5) log(4)];vpa(pi*A,6) 4、1/6251-1/6252=1/6251*6252 5、（1）1/（1+3x）-(1-x)/(1+x)=x*(3*x-1)/[(1+3*x)*(1+x)] (2) sqrt(x+1/x)-sqrt(x-1/x)=2/x/[sqrt(x-1/x)+sqrt(x+1/x)] (3) log10(x1)-log(x2)=log10(x1/x2) (4) [1-cos(2*x)]/x =x^2/factorial(2)-x^4/factorial(4)+x^6/factorial(6)-…

MATLAB及其在大学物理中的应用——第二章习题答案

MATLAB 及其在大学物理中的应用——第二章习题答案 作者：荆楚理工吴世华 2.1试求下列极限： （1）x x x 1)93(lim +∞→（2）5 232)5()3()2(lim +++∞→+++x x x x x x x 2.2求下列函数的导数： （1） x e x x x y -=1sin )(（2）) 4)(3()2)(1()(----=x x x x x y （3）)ln(tan 22y x x y a +=（4）0,ln 1)(>+-=n x a x na x y n n 2.3已知参数方程22,sin cos cos ln dx y d dx dy t t t y t x 和求? ??-==。2.4设22,00y u xv yu yv xu ?????=+=+求。2.5设已知函数矩阵 ?? ????++=z y x z e x y x f y sin 3),(23，试求其雅可比矩阵。2.6求下列不定积分：

（1）dx a x x a x x f ?++=2222) (3)(（2） dx x x x x x f ?+++=1)1()(（3） dx bx xe x f ax ?=cos )(（4）cxdx bx e x f ax sin sin )(?=2.7求下列函数的泰勒幂级数展开。（1）dt t t x ?0sin （2）x x -+11ln （3） )3/3sin(5π+-x e x 分别关于x=0,x=a 的幂级数展开。2.8分别用roots 函数和多项式伴随矩阵的特征值求根法求解方程 6251234=++-x x x 的所有根。 2.9分别用矩阵除法和linsolve 函数法求解下列方程组。 ???????=++-+=+--=+-+--=+-+-0 2660835502492633092458432142143214321x x x x x x x x x x x x x x x 2.10用符号法求解方程组：? ??=-=+151926628y x y x

matlab第2章作业

深圳大学实验报告课程名称：可视化计算语言 实验项目名称：Matlab矩阵及其运算 学院：医学院 专业：生物医学工程 指导教师： 报告人：学号：班级： 实验时间：2012-9-26 实验报告提交时间： 教务部制

实验目的与要求： 1．掌握MATLAB 数据对象的特点以及数据的运算规则。 2．掌握MATLAB 中建立矩阵的方法以及矩阵处理的方法。 3．掌握MATLAB 分析的方法。 实验过程及内容： 1．求下列表达式的值。 （1）w=2×(1+0.34245×10-6) （2）a c b abc c b a x ++-+++ =)tan(e 22ππ，其中a=3.5，b=5，c=-9.8。 （3），])4 8333.0()41[(22απβπαy ---=π其中α=3.32，β=-7.9。 （4）)1ln(2122t t e z t ++=，其中t =??????--65.05i 312。 2．已知 A = ??????????--7613870451， B =???? ??????--023352138 求下列表达式的值： （1）A+6B 和A 2-B+I （其中I 为单位矩阵）。 （2）A*B 、A.*B 和B*A 。 （3）A/B 及B\A 。 （4）[A,B]和[A([1,3],:);B^2]。 3．已知 ????????????---=14.35454 .9632053256545410778.01023Α 完成下列操作： （1）输出A 在[10，25]范围内的全部元素。 （2）取出A 前3行构成矩阵B ，前两列构成矩阵C ，右下角3×2子矩阵构成矩阵D ，B 与C 的乘积构成矩阵E 。 （3）分别求表达式E

第2章--MATLAB数据及其运算-习题答案教学内容

第2章--M A T L A B数据及其运算-习题答案

第2章 MATLAB数据及其运算 习题2 一、选择题 1．下列可作为MATLAB合法变量名的是（）。D A．合计 B．123 C．@h D．xyz_2a 2．下列数值数据表示中错误的是（）。C A．+10 B．1.2e-5 C．2e D．2i 3．使用语句t=0:7生成的是（）个元素的向量。A A．8 B．7 C．6 D．5 4．执行语句A=[1,2,3;4,5,6]后，A(3)的值是（）。B A．1 B．2 C．3 D．4 5．已知a为3×3矩阵，则a(:,end)是指（）。D A．所有元素 B．第一行元素 C．第三行元素 D．第三列元素 6．已知a为3×3矩阵，则运行a (1)=[]后（）。A A．a变成行向量 B．a变为2行2列 C．a变为3行2列 D．a变为2行3列 7．在命令行窗口输入下列命令后，x的值是（）。B >> clear >> x=i*j A．不确定 B．-1 C．1 D．i*j 8．fix(354/100)+mod(354,10)*10的值是（）。D A．34 B．354 C．453 D．43 9．下列语句中错误的是（）。B A．x==y==3 B．x=y=3 C．x=y==3 D．y=3,x=y 10．find(1:2:20>15)的结果是（）。C A．19 20 B．17 19 C．9 10 D．8 9 11．输入字符串时，要用（）将字符括起来。C A．[ ] B．{ } C．' ' D．" " 12．已知s='显示"hello"'，则s的元素个数是（）。A A．9 B．11 C．7 D．18

MATLAB2014b第二章部分答案

4. (1) syms x positive assumeAlso(x,'integer') f1=x^4-5.1*x^3+58.04*x^2-264.384*x+321.408; r1=solve(f1) (2) syms x clear syms a positive solve(x^2-a*x+a^2) 5. A=sym('[a11,a12,a13;a21,a22,a23;a31,a32,a33]') 或 A=sym('a%d%d',[3,3]) Dt=det(A) Iv=inv(A) [Ivs,d]=subexpr(Iv,'d') 6. syms x k f=x^k; s=symsum(f,k,0,inf) subs(s,x,{sym('-1/3'),sym('1/pi'),sym('3')}) 7. syms k syms x positive f=2/(2*k+1)*((x-1)/(x+1))^(2*k+1); s=simple(symsum(f,k,0,inf)) 8. syms t y=abs(sin(t)) d=diff(y) d0=limit(d,t,0,'left') d1=limit(d,t,pi/2) (or subs(d,t,sym('pi/2'))) 9. syms x clear syms x f=exp(-abs(x))*abs(sin(x)) vpa(int(f,x,-5*pi,1.7*pi),64) 10. syms x y

int(int(f,y,1,x^2),x,1,2) 11. syms t x f=sin(t)/t y=int(f,t,0,x) y5=subs(y,x,sym('4.5')) ezplot(y,[0,2*pi]) 12. syms x syms n positive f=sin(x)^n y=int(f,x,0,pi/2) ys=vpa(subs(y,n,sym('1/3'))) yd=vpa(subs(y,n,1/3)) 13. syms a b k n x=a^k; h=b^k; f=subs(h,k,n)*subs(x,k,k-n) y=symsum(f,n,0,k) 15. syms A t w syms a positive ft=A*exp(-a*abs(t)); Fw=fourier(ft,t,w) 16. syms t A w syms tao positive ft=A*((1+t/tao)*(heaviside(t+tao)-heaviside(t))+(1-t/tao)*(heaviside(t)-heaviside(t-tao))); Fw=simple(fourier(ft,t,w)) Fw1=simple(subs(Fw,[A,tao],[2,2])) ezplot(Fw1) 17. syms s t Fs=(s+3)/(s^3+3*s^2+6*s+4) ft=ilaplace(Fs, s,t) 18. syms t s y=sym('f(t)')

MATLAB电路仿真实例

题14.14 图(a)所示电路，已知 V )2cos(15S t u =二端口网络阻抗参数矩阵 Ω?? ????=46j 6j 10Z 求ab 端戴维南等效电路并计算电压o u 。 u -+o u 图题14.14 (一)手动求解： 将网络N 用T 型电路等效，如图(b)所示 S U +-o U 等效阻抗 Ω=-+-?+ -=4.6j615j6j6)15(6j 6j 4i Z 开路电压 V 2j302 15j6j6105j6OC =?∠?+-+=U V 1482 18.3j46.42j3j4j4Z j4OC o ?∠=+?=?+=U U i

所以 )1482cos(18.3o ?+=t u V (二)Matlab 仿真： ⒈分析：本次仿真需输入各阻抗Zl 、Z1、Z2、Z3、Z4以及激励源Us 的参数值，仿真结果需输出开路电压Uoc 、等效阻抗Zi 以及电感两端电压U0的幅值和相位信息，并绘制Uoc ，U0的值随时间变化的波形曲线。其中各元件与原图的对应关系如下图所示： ⒉编辑M 文件的源程序如下： clear %清空自定义变量 z1=4-6j;z2=6j;z3=10-6j;z4=5;us=15*exp(j*0);zl=4j;%输入各元件参数 zi=z1+(z2*(z3+z4)/(z2+z3+z4));%等效阻抗zi 的计算表达式uoc=us*z2/(z2+z3+z4);%开路电压uoc 的计算表达式u0=zl/(zi+zl)*uoc;%电感两端电压uo 的计算表达式disp('The magnitude of zi is'); %在屏幕上显示“The magnitude of zi is ”disp(abs(zi)) %显示等效阻抗zi 的模disp('The phase of zi is'); %在屏幕上显示“The phase of zi is ”disp(angle(zi)*180/pi)%显示等效阻抗zi 的辐角 disp('The magnitude of uoc is'); %在屏幕上显示“The magnitude of uoc is ” disp(abs(uoc))%显示开路电压uoc 的模