matlab 大实验

t=-10:0.01:10;y=stepfun(t,0);plot(t,y) 10.90.80.70.60.50.40.30.20.10-10-8-6-4-20246810 sint(t)

t=0:.00001:10;y=stepfun(t,0);z=sin(t).*y;plot(t,z)

10.80.60.40.20-0.2-0.4-0.6-0.8-1012345678910

et(t)

t=linspace(0,5,51);y=stepfun(t,0);z=exp(-t).*y;plot(t,z) 10.90.80.70.60.50.40.30.20.1000.511.522.533.544.55 etcost(t)

t=linspace(0,5,51);y=stepfun(t,0);z=exp(-t).*y.*cos(t);plot(t,z)

1.210.80.60.40.20-0.200.511.522.533.544.55 g2(t)

t=linspace(-2,2,1001);y=(t>-1)&(t<1);plot(t,y) 10.90.80.70.60.50.40.30.20.10-2-1.5-1-0.500.511.52

Sa(3t)

t=-3*pi:0.001:3*pi;y=sinc(t);plot(t,y);grid on

10.80.60.40.20-0.2-0.4-10-8-6-4-20246810 (2) cos 2k

k=-5:0.2:5;a=pi/2;y=cos(a*k);stem(k,y) 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345

cos

8k

k=-5:0.2:5;a=pi/8;y=cos(a*k);stem(k,y) 10.80.60.40.20-0.2-0.4-5-4-3-2-1012345 cos4k

k=-5:0.2:5;a=pi/4;y=cos(a*k);stem(k,y) 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345

cosk

k=-5:0.2:5;a=pi;y=cos(a*k);stem(k,y) 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345

cos3k 2k=-5:0.2:5;a=1.5*pi;y=cos(a*k);stem(k,y) 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345

cos7k 4k=-5:0.2:5;a=7/4*pi;y=cos(a*k);stem(k,y) 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345 cos15k 8

k=-5:0.2:5;a=15/8*pi;y=cos(a*k);stem(k,y)

10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345 cos2k

k=-5:0.2:5;a=2*pi;y=cos(a*k);stem(k,y) 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345

cos(

7k) 43k=-5:0.2:5;a=7/4*pi;y=cos(a*k+pi/3);stem(k,y) 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345 (3) cos 3kjsin3k

k=-5:0.2:5;a=pi/3;y=cos(a*k);x=sin(a*k);z=y+x*i;stem(k,z) 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345

()(cos12k3kjsin3k)

k=-5:0.2:5;a=pi/3;b=(1/2).^k;y=cos(a*k);x=sin(a*k);z=b.*(y+a*i);stem(k,z) 20151050-5-10-5-4-3-2-1012345 e

k=-5:0.2:5;a=pi/3;z=exp(-i*a*k);stem(k,z) 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345jk3 

2、 熟悉ezplot、polar、mesh等指令

(1) 用ezplot绘ecost(t)的图。

ezplot3('sin(3*t)*cos(t)','sin(3*t)*sin(t)','t','animate'); 10.80.60.40.20-0.2-0.4-0.6-0.8-1-5-4-3-2-1012345t

ezplot('exp(-1*t).*cos(t)',[0,10])

exp(-1 t) cos(t)0.20.150.10.050-0.05012345t678910 ezplot('sin(-3*t)',[0,10]) sin(-3 t)

654z321010.50-0.5y-120t46108

t=0:0.001:5;ezplot('sin(3*t)','cos(2*t)')

x = sin(3 t), y = cos(2 t)10.80.60.40.20-0.2-0.4-0.6-0.8y-1-0.8-0.6-0.4-0.2

ezplot3('sin(3*t)','cos(3*t)*sin(t)','t','animate') x = sin(3 t), y = cos(3 t) sin(t), z = t0x0.20.40.60.817654z321010.50-0.5y-1-1-0.5x0.501

Theta=0:2*pi:361;Beita=pi*sin(Theta);A=abs(sin(3*Beita)./(6*sin(0.5*Beita)));polar(Theta,A)

90120 160 0.8 0.6150 0.4 0.2301800210330240270300

zsin(x2y2)x2y2x,y的取值范围是[8,8]

(2) 用mesh绘的三维曲面。

x=-8:0.5:8; y=x'; X=ones(size(y))*x; Y=y*ones(size(x)); R=sqrt(X.^2+Y.^2)+eps; Z=sin(R)./R; mesh(X,Y,Z); colormap(hot)

xlabel('x'),ylabel('y'),zlabel('z')

3、 计算数值积分：

10t

(1)

e0sintdt

10Sa2tdt

t=linspace(0,pi,1001); x=exp(-t).*sin(t); y1=sum(x)*t(2), y2=trapz(x)*t(2),

f=inline('exp(-t).*sin(t)'), y3=quad(f,0,pi),

y1 =

0.5216 y2 =

0.5216 f =

Inline function:

f(t) = exp(-t).*sin(t) y3 =

0.5216

t=0:0.001:pi;x=exp(-t).*sin(t);y=trapz(x)*t(2)

y =

0.5216

1010Sa2tdt

t=linspace(-10*pi,10*pi,1001); f=inline('sinc(2*t)'), y3=quad(f,-10*pi,10*pi)

f =

Inline function: f(t) = sinc(2*t) y3 =

0.5014

(2) f1(t)cos3t(t)，f2(t)(t)(t4)，求

dt=0.01; k1=0:dt:4; f1=cos(3*k1); k2=k1;

f2=stepfun(k2,0)-stepfun(k2,4); f=dt*conv(f1,f2); k0=k1(1)+k2(1);

k3=length(f1)+length(f2)-2; k=k0:dt:k0+k3*dt;

subplot(2,2,1);plot(k1,f1); subplot(2,2,2);plot(k2,f2); subplot(2,2,3);plot(k,f); 10.500.4-0.5-10.201234001210.80.6f1(t)f2(t)

340.40.20-0.2-0.4-0.602468

一、 系统的描述、时域与变换域分析、求解

1,解方程 (1) yk2yk12yk2fkfk1,y(-1)=1,y(-2)=2, f(k)=ε(k),求系统的h(k)、

g(k)、y(k)

b=[1 1 0];a=[1 2 2]; k=0:19;

y1=dimpulse(b,a,20); y2=dstep(b,a,20); subplot(2,1,1); stem(k,y1);

xlabel('k');ylabel('h(k)');title('dimpulse'); subplot(2,1,2); stem(k,y2);

xlabel('k');ylabel('g(k)');title('dstep') dimpulse1000500h(k)0-5000246810kdstep1214161820600400g(k)2000-2000246810k1214161820 n=1:20+3;

F=stepfun(n-3,0); F(1)=0;F(2)=0; Y(1)=2;Y(2)=1; for k=3:23

Y(k)=F(k)-F(k-1)-2*Y(k-1)-2*Y(k-2); end

n1=n-3;

stem(n1,Y);xlabel('k');ylabel('f(k)');grid on;

6000500040003000f(k)200010000-1000-2000-505k101520

(2) (a)求解微分方程y''4y'3yf'3f，初始条件为：y(0+)=y'(0+)=1，f(t)(t)。理论计算全响应，用ezplot绘结果图。范围为区间[0,4]。

(b)用数值法求解需先化为差分方程。化出对应的差分方程。步长取T=0.1。 (c)确定此差分方程的初始条件。

(d)迭代法计算机求解。请编写程序，运行调试。比较数值解的准确程度。

(a)把微分方程转化为差分方程为

143y(k)240y(k1)100y(k2)13f(k)10f(k1) (b)初始条件需要求解方程组：

y(1)1 y(1)y(2)10.1得到：y(-1)=1；y(-2)=0.9。

Ts=0.1; %时间间隔

n=1:4/Ts+2; %坐标1,2代表初始条件坐标-2,-1 F(n)=1;F(1)=0;F(2)=0; %激励信号初始条件

Y(1)=0.9;Y(2)=1; %输入初始条件y(-2)=0.9;y(-1)=1 for k=3:length(n) %迭代求解

Y(k)=(13*F(k)-10*F(k-1)-100*Y(k-2)+240*Y(k-1))/143; end

n1=(n-3)*Ts;subplot(2,1,1);plot(n1,Y);

xlabel('t');ylabel('y(t)');title('全响应(数值计算)');grid on; subplot(2,1,2);ezplot('-exp(-3*t)+exp(-t)+1',[0,4]); xlabel('t');ylabel('y(t)');title('全响应');grid on;

全响应(数值计算)1.61.4y(t)1.210.8-0.500.511.5t全响应1.41.322.533.54y(t)1.21.1100.511.52t2.533.54

d2vcdvLCRCcvcetdtdt(3) 电压源激励的RLC串联回路微分方程为：。求当R=1Ω，C=1F，

L=1H时二阶电路的冲激响应和阶跃响应并绘波形图。

a=[1 1 1]; b=[1];

y=impulse(b,a,0:0.1:5); z=step(b,a,0:0.1:5);

subplot(2,1,1);impulse(b,a);title('冲激响应') subplot(2,1,2);step(b,a);title('阶跃响应')

冲激响应0.60.4Amplitude0.20-0.20246Time (sec)81012阶跃响应1.5Amplitude10.500246Time (sec)81012

(4) 求解微分方程y''+3y'=3f(t)，初始条件y'(0)=1.5，y(0)=0，f(t)=ε(t)

y=dsolve('(D2y)+Dy*3=3','y(0)=0','Dy(0)=1.5')

y =

t - 1/(6*exp(3*t)) + 1/6

2(p4p3)r(t)(p2)e(t)的冲激响应。 (5) 求微分方程

a=[1 4 3 ];b=[1 2];

impulse(b,a)

title('冲激响应')

冲激响应10.90.80.70.6Amplitude0.50.40.30.20.1000.511.522.5Time (sec)33.544.55.

(6)对连续系统y''+3y'+2y=3(ε(t)-ε(t-2))建立SIMULINK模型，并观察系统的零状态响应。

3Gain2e(t)Add1Addy''1sInt1Gain1e(t-2)3Gain2y_tTo WorkspaceClocky'1sIntyScope

(7)求下面流图的传输函数。

131f2z1x2-11z1x12y

syms z;

Q=vpa(zeros(7,7)); Q(2,1)=2;Q(2,3)=-1;

Q(2,5)=-2;Q(3,2)=1/z;Q(4,1)=3;Q(4,3)=1;Q(5,4)=1/z;Q(6,1)=1;Q(6,5)=2; Q(7,6)=1;

P=[1;0;0;0;0;0;0]; I=eye(size(Q)); H=(I-Q)\\P;

HS=simple(H(7))

HS =

(6*z + 10)/(z^2 + z + 2) + 1

-21、对如下系统函数绘制零极点图、判断稳定性、对稳定的绘制频率特性曲线。

z22z1(1)Hz

2z31b=[1 -2 -1];a=[2 0 0 -1];zplane(b,a)

10.5Imaginary Part0-0.5-1-1-0.500.51Real Part1.522.5

w=0:0.01:10;freqz(b,a,w); 10Magnitude (dB)

50-500.511.522.5Normalized Frequency ( rad/sample)3200Phase (degrees)0-200-40000.511.522.5Normalized Frequency ( rad/sample)3 、Hzz1、。 z31

b=[1 1];a=[1 0 0 -1];zplane(b,a)

10.80.60.4Imaginary Part0.20-0.2-0.4-0.6-0.8-1-1-0.50Real Part0.512

b=[1 1];a=[1 0 0 -1];w=0:0.001:10;freqz(b,a,w); 100Magnitude (dB)500-50-10000.511.522.5Normalized Frequency ( rad/sample)3500Phase (degrees)0-500-100000.511.522.5Normalized Frequency ( rad/sample)3

z22Hz3、 2z2z4z1b=[1 0 2];a=[1 2 -4 1]; zplane(b,a) 1.510.50-0.5-1-1.5-3.5-3-2.5-2-1.5-1Real Part-0.500.51Imaginary Part z3 Hz3 2z0.2z0.3z0.4

b=[1 0 0 0];a=[1 0.2 0.3 0.4];zplane(b,a) 10.80.60.4Imaginary Part0.20-0.2-0.4-0.6-0.8-1-1-0.50Real Part0.513

w=0:0.01:10;freqz(b,a,w); 10Magnitude (dB)50-5-1000.511.522.5Normalized Frequency ( rad/sample)350Phase (degrees)0-5000.511.522.5Normalized Frequency ( rad/sample)3

(2)H(s) s3、。 32s3s6s4b=[1 3];a=[1 3 6 4];pzmap(b,a)

Pole-Zero Map21.510.5Imaginary Axis0-0.5-1-1.5-2-3-2.5-2-1.5Real Axis-1-0.50

b=[1 3];a=[1 3 6 4];w=0:.01:10;freqs(b,a,w) 10Magnitude0 10-110-210-210-110Frequency (rad/s)01010Phase (degrees)-50-100-150-200-210-1011010Frequency (rad/s)10

H(s)s2、 3(s1)(s3)b=[1 2];a=[1 9 27 27];pzmap(b,a)

2x 10-5Pole-Zero Map1.510.5Imaginary Axis0-0.5-1-1.5-2-3.5-3-2.5-2-1.5Real Axis-1-0.50 w=0:0.01:10;freqs(b,a,w); 10Magnitude0 10-110-210-210-110Frequency (rad/s)01010Phase (degrees)-20-40-60-80-210-1011010Frequency (rad/s)10

H(s)

s1 32s2ss1b=[1 1];a=[1 2 1 1];pzmap(b,a) Pole-Zero Map0.80.60.40.2Imaginary Axis0-0.2-0.4-0.6-0.8-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-0.20Real Axis

w=0:0.01:10;freqs(b,a,w); 10Magnitude1 1010100-1-210-210-110Frequency (rad/s)010150Phase (degrees)0-50-100-21010-110Frequency (rad/s)0101

2、 用系统转换的指令得到上述系统函数的零极点增益形式、状态方程形式和多项式形式

z22z1(1)Hz、 32z1

b=[1 -2 -1];a=[2 0 0 -1]; [Z,P,K]=tf2zp(b,a) [A,B,C,D]=tf2ss(b,a) [b,a]=ss2tf(A,B,C,D)

Z =

-0.4142 2.4142 P =

-0.3969 + 0.6874i -0.3969 - 0.6874i 0.7937 K =

0.5000 A =

0 0 0.5000 1.0000 0 0 0 1.0000 0 B = 1 0 0 C =

0.5000 -1.0000 -0.5000 D = 0 b =

0 0.5000 -1.0000 -0.5000 a =

1.0000 -0.0000 -0.0000 -0.5000

Hzz1z31、

b=[1 1];

a=[1 0 0 -1];

[Z,P,K]=tf2zp(b,a) [A,B,C,D]=tf2ss(b,a) [b,a]=ss2tf(A,B,C,D)

Z = -1 P =

-0.5000 + 0.8660i -0.5000 - 0.8660i 1.0000 K = 1 A =

0 0 1 1 0 0 0 1 0 B = 1 0 0 C =

0 1 1 D = 0 b =

0 0.0000 1.0000 1.0000 a =

1.0000 -0.0000 0.0000 -1.0000

Hzz22z32z24z1、

b=[1 0 2];

a=[1 2 -4 1];

[Z,P,K]=tf2zp(b,a) [A,B,C,D]=tf2ss(b,a) [b,a]=ss2tf(A,B,C,D)

Z =

0 + 1.4142i 0 - 1.4142i P =

-3.3028 1.0000 0.3028 K = 1 A =

-2 4 -1 1 0 0 0 1 0 B = 1 0 0 C =

1 0 2 D = 0 b =

0 1.0000 -0.0000 2.0000 a =

1.0000 2.0000 -4.0000 1.0000

zz3Hz30.2z20.3z0.4 b=[1 0 0 0];

a=[1 0.2 0.2 0.4]; [Z,P,K]=tf2zp(b,a) [A,B,C,D]=tf2ss(b,a) [b,a]=ss2tf(A,B,C,D)

Z = 0 0 0 P =

0.2553 + 0.7055i 0.2553 - 0.7055i -0.7106 K = 1 A =

-0.2000 -0.2000 -0.4000 1.0000 0 0 0 1.0000 0 B = 1 0 0 C =

-0.2000 -0.2000 -0.4000 D = 1 b =

1 0 0 0 a =

1.0000 0.2000 0.2000 0.4000

(2)H(s)s3s33s26s4、

b=[1 3];

a=[1 3 6 4];

[Z,P,K]=tf2zp(b,a)

[A,B,C,D]=tf2ss(b,a)

[b,a]=ss2tf(A,B,C,D)

Z = -3 P =

-1.0000 + 1.7321i -1.0000 - 1.7321i -1.0000 K = 1 A =

-3 -6 -4 1 0 0 0 1 0 B = 1 0 0 C =

0 1 3 D = 0 b =

0 0.0000 1.0000 3.0000 a =

1.0000 3.0000 6.0000 4.0000

H(s)s2(s1)(s3)3、

b=[1,2];

a=[1 10 36 54 27]; [Z,P,K]=tf2zp(b,a) [A,B,C,D]=tf2ss(b,a) [b,a]=ss2tf(A,B,C,D)

Z = -2 P =

-3.0000 + 0.0000i -3.0000 - 0.0000i -3.0000 -1.0000 K = 1 A =

-10 -36 -54 -27 1 0 0 0 0 1 0 0

0 0 1 0 B = 1 0 0 0 C =

0 0 1 2 D = 0 b =

0 -0.0000 -0.0000 1.0000 2.0000 a =

1.0000 10.0000 36.0000 54.0000 27.0000

H(s)s1s32s2s1。

b=[1 1];

a=[1 2 1 1];

[Z,P,K]=tf2zp(b,a) [A,B,C,D]=tf2ss(b,a)

[b,a]=ss2tf(A,B,C,D)

Z = -1 P =

-1.7549 -0.1226 + 0.7449i -0.1226 - 0.7449i K = 1 A =

-2 -1 -1 1 0 0 0 1 0 B = 1 0 0 C =

0 1 1 D = 0 b =

0 -0.0000 1.0000 1.0000 a =

1.0000 2.0000 1.0000 1.0000

3、 用Simulink搭建上述系统

z22z1)Hz、 32z1Hz

z1、 z31

Hz

z2 32z2z4z12

、Hz

z。

z30.2z20.3z0.43

(2)H(s)

s3、

s33s26s4

H(s)

s2、

(s1)(s3)3

H(s)

s1。 32s2ss1

二、 Fourier变换

1、 计算Fourier级数

(1)求下图周期三角波y(t)的傅立叶级数系数，绘制振幅谱和相位谱图。

y(t)1π

t=linspace(-pi,pi,201);

n=0:9;n1=linspace(1,10,10);

gt= (n1'*((-1/pi)*abs(t)+1)).*exp(-i*n'*t); Fn=1/(2*pi)*trapz(gt,2)*0.02; sig=abs(Fn)>0.00001; Fn=Fn.*sig; figure;

subplot(2,1,1); stem(n,abs(Fn));

xlabel('n'),ylabel('Fn');title('振幅谱');grid; subplot(2,1,2); stem(n,angle(Fn));

xlabel('n'),ylabel('Theta');title('相位谱');grid;

t

振幅谱0.40.3Fn0.20.100x 10-141234n相位谱5678932Theta10-101234n56789 (2) 计算

，并绘制幅频、相f(t)cos(20t)g6(t)按照T=10周期延拓的周期信号的傅立叶系数Fn频曲线。

t=linspace(-3,3,101);

n=0:9;n1=[1 1 1 1 1 1 1 1 1 1];

gt=(n1'*(cos(20*t).*rectpuls(t,6))).*exp(-i*n'*t); Fn=1/10*trapz(gt,2)*0.02; sig=abs(Fn)>0.00001; Fn=Fn.*sig; figure;

subplot(2,1,1); stem(n,abs(Fn));

xlabel('n'),ylabel('Fn');title('振幅谱');grid; subplot(2,1,2); stem(n,angle(Fn));

xlabel('n'),ylabel('Theta');title('相位谱');grid;

32Fnx 10-3振幅谱1001234n相位谱56789420-2Theta01234n56789

2、 求Fourier变换与反变换（用数值方法作，计算DFT的函数fft()）： f1(t)g4(t)，f2(t)cos(20t)g6(t)，f3(t)Sa(2t)g20(t)，

绘制幅频和相频特性曲线。然后计算傅立叶逆变换计算时域波形并绘图，和原始函数进行对比，观察频谱混叠与频谱泄露现象

f1(t)g4(t) T=10; N=512; Ts=T/N; PT=N/16;

t=linspace(-T/2,T/2-T/N,N); f=0*t;

f(t>-2&t<2)=1;

subplot(3,1,1);plot(t,f);grid; w=[0:N-1]*2*pi/(Ts*N);

Fw=fft(f)*Ts./exp(-i*0.5*T*w); Fw=fftshift(Fw);

w=w-2*pi*(1/(Ts*2)); subplot(3,1,2);

plot(w(N/2-PT:N/2+PT),abs(Fw(N/2-PT:N/2+PT)));grid; subplot(3,1,3);

plot(w(N/2-PT:N/2+PT),angle(Fw(N/2-PT:N/2+PT)));grid;

10.50-56420-2550-5-25-20-15-10-505101520-4-3-2-1012345-20-15-10-505101520 f2(t)cos(20t)g6(t)

T=6; N=200;

t=linspace(-4,4,201);

f=cos(20*t).*stepfun(t,-3)-cos(20*t).*stepfun(t,3); W=16*pi; K=100;

w=linspace(-W/2,W/2-W/K,K); F=0*w;

for k=1:K for n=1:N

F(k)=F(k)+T/N*f(n)*exp(-j*w(k)*t(n)); end end

f1=0*t;

for n=1:N for k=1:K

f1(n)=f1(n)+W/(2*pi*K)*F(k)*exp(j*w(k)*t(n)); end end

subplot(2,1,1);

plot(t,f,'-k',t,f1,':k');

xlabel('t');ylabel('f(t)');legend('f(t)'); subplot(2,1,2); plot(w,F,'-k');

xlabel('w');ylabel('F(w)');grid;

1f(t)0.5f(t) 0-0.5-1 -4-3-2-10t123432F(w)10-1-30-20-100w102030

f3(t)Sa(2t)g20(t)

T=8; N=200;

t=linspace(-12*pi,12*pi,201);

f=sinc(2*t).*stepfun(t,-10*pi)-sinc(20*t).*stepfun(t,10*pi); W=16*pi; K=100;

w=linspace(-W/2,W/2-W/K,K); F=0*w;

for k=1:K for n=1:N

F(k)=F(k)+T/N*f(n)*exp(-j*w(k)*t(n)); end end

f1=0*t;

for n=1:N for k=1:K

f1(n)=f1(n)+W/(2*pi*K)*F(k)*exp(j*w(k)*t(n)); end end

subplot(2,1,1);

plot(t,f,'-k',t,f1,':k');

xlabel('t');ylabel('f(t)');legend('f(t)'); subplot(2,1,2); plot(w,F,'-k');

xlabel('w');ylabel('F(w)');grid;

1f(t)0.5f(t) 0-0.5 -40-30-20-100t102030400.060.04F(w)0.020-0.02-30-20-100w102030

实验总结：不知不觉中，这学期已经接近尾声了 ，MATLAB这么课也已经学习了九周，第一次试着以这种方式学习一门课程，基本上都是自己在看书，然后在不断的实践，不断地从实践中学习到新的知识 ， 再继续前进，每天利用闲暇的时间看书，然后会看来后上机实验，从中学习再进步，这种教学方式不再和以前一样，但我们有自觉 ，了解到它的重要性，然后再自 发的从心里对他感兴趣，兴趣是最好的老师，最后才到了今天，

MATLAB作为一种工具 ，被广泛的应用于各种行业之中，我们学习的都是一些最基础的东西 ，简单的画图，福利叶变换，这两个多月下来，我不仅仅学到了这么多 更多的是一种学习方法，一种对学习态度。我相信，在以后的工作生活中，我可以更好。