<信号与系统MATLAB实践>
练习一
实验一
二. 熟悉简单的矩阵输入
1.实验代码
>>A=[1,2,3;4,5,6;7,8,9]
实验结果
A =
1 2 3
4 5 6
7 8 9
3．实验代码
>>B=[9,8,7;6,5,4;3,2,1]
C=[4,5,6;7,8,9;1,2,3]
实验结果：
B =
9 8 7
6 5 4
3 2 1
C =
4 5 6
7 8 9
1 2 3
4．>> A
A =
1 2 3
4 5 6
7 8 9
>> B
B =
9 8 7
6 5 4
3 2 1
>> C
C =
4 5 6
7 8 9
1 2 3
三. 基本序列运算
1.>>A=[1,2,3],B=[4,5,6]
A =
1 2 3
B =
4 5 6
>> C=A+B
C =
5 7 9
>> D=A-B
D =
-3 -3 -3
>> E=A.*B
E =
4 10 18
>> F=A./B
F =
0.2500 0.4000 0.5000
>> G=A.^B
G =
1 32 729
>> stem(A)
>> stem(B)
>> stem(C)
>> stem(D)
>> stem(E)
>> stem(F)
>> stem(G)
再举例:
>> a=[-1,-2,-3]
a =
-1 -2 -3
>> b=[-4,-5,-6]
b =
-4 -5 -6
>> c=a+b
c =
-5 -7 -9
>> d=a-b
d =
3 3 3
>> e=a.*b
e =
4 10 18
>> f=a./b
f =
0.2500 0.4000 0.5000
>> g=a.^b
g =
1.0000 -0.0313 0.0014
>> stem(a)
>> stem(b)
>> stem(c)
>> stem(d)
>> stem(e)
>> stem(f)
>> stem(g)
2. >>t=0:0.001:10
f=5*exp(-t)+3*exp(-2*t);
plot(t,f)
ylabel('f(t)');
xlabel('t');
title('(1)');
>> t=0:0.001:3;
f=(sin(3*t))./(3*t);
plot(t,f)
ylabel('f(t)');
xlabel('t');
title('(2)');
>> k=0:1:4;
f=exp(k);
stem(f)
四. 利用MATLAB求解线性方程组
2.
>>A=[1,1,1;1,-2,1;1,2,3]
b=[2;-1;-1]
x=inv(A)*b
A =
1 1 1
1 -2 1
1 2 3
b =
2
-1
-1
x =
3.0000
1.0000
-2.0000
4.
>> A=[2,3,-1;3,-2,1;1,2,1]
b=[18;8;24]
x=inv(A)*b
A =
2 3 -1
3 -2 1
1 2 1
b =
18
8
24
x =
4
6
8
实验二
二.
1.
>> k=0:50
x=sin(k);
stem(x)
xlabel('k');
ylabel('sinX');
title('sin(k)ε(k)');
2.
>> k=-25:1:25
x=sin(k)+sin(pi*k);
stem(k,x)
xlabel('k');
ylabel('f(k)');
title('sink+sinπk');
3.
>> k=3:50
x=k.*sin(k);
stem(k,x)
xlabel('k');
ylabel('f(k)');
title('ksinkε(k-3)');
4.
%函数
function y=f1(k)
if k<0
y=(-1)^k;
else y=(-1)^k+(0.5)^k;
end
%运行代码
for k=-10:1:10;
y4(k+11)=f1(k);
end
k=-10:1:10;
stem(k,y4);
xlabel('k');
ylabel('f(k)');
title('4');
七．
2．>> f1=[1 1 1 1];
f2=[3 2 1];
conv(f1,f2)
ans =
3 5 6 6 3 1
3.
函数定义：
function [r]= pulse( k )
if k<0
r=0;
else
r=1;
end
end
运行代码
for k=1:10
f1(k)=pulse(k);
f2(k)=(0.5^k)*pulse(k);
end
conv(f1,f2)
结果
ans =
Columns 1 through 10
0.5000 0.7500 0.8750 0.9375 0.9688 0.9844 0.9922 0.9961 0.9980 0.9990
Columns 11 through 20
0.9995 0.9998 0.9999 0.9999 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
Columns 21 through 30
0.5000 0.2500 0.1250 0.0625 0.0312 0.0156 0.0078 0.0039 0.0020 0.0010
Columns 31 through 39
0.0005 0.0002 0.0001 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000
4
for i=1:10
f1(i)=pulse(i);
f2(i)=((-0.5)^i)*pulse(i);
end
conv(f1,f2)
结果
ans =
Columns 1 through 10
-0.5000 -0.2500 -0.3750 -0.3125 -0.3438 -0.3281 -0.3359 -0.3320 -0.3340 -0.3330
Columns 11 through 20
-0.3325 -0.3323 -0.3322 -0.3321 -0.3321 -0.3320 -0.3320 -0.3320 -0.3320 -0.3320
Columns 21 through 30
0.1680 -0.0820 0.0430 -0.0195 0.0117 -0.0039 0.0039 -0.0000 0.0020 0.0010
Columns 31 through 39
0.0005 0.0002 0.0001 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000
实验三
2．
clear;
x=[1,2,3,4,5,6,6,5,4,3,2,1];
N=0:11;
w=-pi:0.01:pi;
m=length(x);
n=length(w);
for i=1:n
F(i)=0;
for k=1:m
F(i)=F(i)+x(k)*exp(-1j*w(i)*k);
end
end
F=F/10;
subplot(2,1,1);
plot(w,abs(F),'b-');xlabel('w');ylabel('F');title('幅度频谱');grid
subplot(2,1,2);
plot(w,angle(F),'b-');xlabel('w');
X=fftshift(fft(x))/10;
subplot(2,1,1);
hold on;
plot(N*2*pi/12-pi,abs(X),'r.');
legend('DIFT算法','DFT算法');
subplot(2,1,2);hold on;
plot(N*2*pi/12-pi,angle(X),'r.');
xlabel('w');ylabel('相位');title('相位频谱');grid
三．
1.
%fun1.m
function y=fun1(x)
if((-pi
y=pi+x;
elseif ((0
y=pi-x;
else
y=0
end
%new.m
clear all
clc
for i=1:1000
g(i)=fun1(2/1000*i-1);
w(i)=(i-1)*0.2*pi;
end
for i=1001:10000
g(i)=0;
w(i)=(i-1)*0.2*pi;
end
G=fft(g)/1000;
subplot(1,2,1);
plot(w(1:50),abs(G(1:50)));
xlabel('w');ylabel('G');title('DFT幅度频谱');
subplot(1,2,2);
plot(w(1:50),angle(G(1:50)))
xlabel('w');ylabel('Fi');title('DFT相位频谱');
2.
%fun2.m
function y=fun2(x)
if x<1 && x>-1
y=cos(pi*x/2);
else
y=0;
end
%new2.m
for i=1:1000
g(i)=fun2(2/1000*i-1);
w(i)=(i-1)*0.2*pi;
end
for i=1001:10000
g(i)=0;
w(i)=(i-1)*0.2*pi;
end
G=fft(g)/1000;
subplot(1,2,1);
plot(w(1:50),abs(G(1:50)));
xlabel('w');ylabel('G');title('幅度频谱');
subplot(1,2,2);
plot(w(1:50),angle(G(1:50)))
xlabel('w');ylabel('Fi');title('相位频谱');
3.
%fun3.m
function y=fun3(x)
if x<0 && x>-1
y=1;
elseif x>0 && x<1
y=-1;
else
y=0
end
%new.m
for i=1:1000
g(i)=fun3(2/1000*i-1);
w(i)=(i-1)*0.2*pi;
end
for i=1001:10000
g(i)=0;
w(i)=(i-1)*0.2*pi;
end
G=fft(g)/1000;
subplot(1,2,1);
plot(w(1:50),abs(G(1:50)));
xlabel('w');ylabel('G');title('DFT幅度频谱');
subplot(1,2,2);
plot(w(1:50),angle(G(1:50)))
xlabel('w');ylabel('Fi');title('DFT相位频谱');
练习二
实验六
一．用MATLAB语言描述下列系统，并求出极零点、
1.
>> Ns=[1];
Ds=[1,1];
sys1=tf(Ns,Ds)
实验结果：
sys1 =
1
-----
s + 1
>> [z,p,k]=tf2zp([1],[1,1])
z =
Empty matrix: 0-by-1
p =
-1
k =
1
2.
>>Ns=[10]
Ds=[1,-5,0]
sys2=tf(Ns,Ds)
实验结果：
Ns =
10
Ds =
1 -5 0
sys2 =
10
---------
s^2 - 5 s
>>[z,p,k]=tf2zp([10],[1,-5,0])
z =
Empty matrix: 0-by-1
p =
0
5
k =
10
二．已知系统的系统函数如下，用MATLAB描述下列系统。
1．
>> z=[0];
p=[-1,-4];
k=1;
sys1=zpk(z,p,k)
实验结果：
sys1 =
s
-----------
(s+1) (s+4)
Continuous-time zero/pole/gain model.
2.
>> Ns=[1,1]
Ds=[1,0,-1]
sys2=tf(Ns,Ds)
实验结果：
Ns =
1 1
Ds =
1 0 -1
sys2 =
s + 1
-------
s^2 - 1
Continuous-time transfer function.
3．
>> Ns=[1,6,6,0];
Ds=[1,6,8];
sys3=tf(Ns,Ds)
实验结果：
Ns =
1 6 6 0
Ds =
1 6 8
sys3 =
s^3 + 6 s^2 + 6 s
-----------------
s^2 + 6 s + 8
Continuous-time transfer function.
六．已知下列H(s)或H(z)，请分别画出其直角坐标系下的频率特性曲线。
1.
>> clear;
for n = 1:400
w(n) = (n-1)*0.05;
H(n) = (1j*w(n))/(1j*w(n)+1);
end
mag = abs(H);
phase = angle(H);
subplot(2,1,1)
plot(w,mag);title('幅频特性')
subplot(2,1,2)
plot(w,phase);title('相频特性')
实验结果：
2.
>> clear;
for n = 1:400
w(n) = (n-1)*0.05;
H(n) = (2*j*w(n))/((1j*w(n))^2+sqrt(2)*j*w(n)+1);
end
mag = abs(H);
phase = angle(H);
subplot(2,1,1)
plot(w,mag);title('幅频特性')
subplot(2,1,2)
plot(w,phase);title('相频特性')
实验结果：
3.
>>clear;
for n = 1:400
w(n) = (n-1)*0.05;
H(n) = (1j*w(n)+1)^2/((1j*w(n))^2+0.61);
end
mag = abs(H);
phase = angle(H);
subplot(2,1,1)
plot(w,mag);title('幅频特性')
subplot(2,1,2)
plot(w,phase);title('相频特性')
实验结果：
4.
>>clear;
for n = 1:400
w(n) = (n-1)*0.05;
H(n) =3*(1j*w(n)-1)*(1j*w(n)-2)/(1j*w(n)+1)*(1j*w(n)+2);
end
mag = abs(H);
phase = angle(H);
subplot(2,1,1)
plot(w,mag);title('幅频特性')
subplot(2,1,2)
plot(w,phase);title('相频特性')
实验结果：
实验七
三．已知下列传递函数H(s)或H(z)，求其极零点，并画出极零图。
1.
>> z=[1,2]';
p=[-1,-2]';
zplane(z,p)
实验结果：
2.
>> z=[1,2];
p=[-1,-2];
zplane(z,p)
>> num=[1];
den=[1,0];
[z,p,k]=tf2zp(num,den);
zplane(z,p)
>> num=[1];
den=[1,0];
[z,p,k]=tf2zp(num,den)
zplane(z,p)
实验结果：
z =
Empty matrix: 0-by-1
p =
0
k =
1
3.
>> num=[1,0,1];
den=[1,2,5];
[z,p,k]=tf2zp(num,den)
zplane(z,p)
实验结果：
z =
0 + 1.0000i
0 - 1.0000i
p =
-1.0000 + 2.0000i
-1.0000 - 2.0000i
k =
1
4.
>> num=[1.8,1.2,1.2,3];
den=[1,3,2,1];
[z,p,k]=tf2zp(num,den)
zplane(z,p)
实验结果：
z =
-1.2284
0.2809 + 1.1304i
0.2809 - 1.1304i
p =
-2.3247
-0.3376 + 0.5623i
-0.3376 - 0.5623i
k =
1.8000
5．
>> clear;
A=[0,1,0; 0,0,1; -6,-11,-6];
B=[0;0;1];
C=[4,5,1];
D=0;
sys5=ss(A,B,C,D);
pzmap(sys5)
实验结果：
五．求出下列系统的极零点，判断系统的稳定性。
1.
>> clear;
A=[5,2,1,0; 0,4,6,0; 0,-3,-6,-1;1,-2,-1,3];
B=[1;2;3;4];
C=[1,2,5,2];
D=0;
sys=ss(A,B,C,D);
[z,p,k]=ss2zp(A,B,C,D,1)
pzmap(sys)
实验结果：
z =
4.0280 + 1.2231i
4.0280 - 1.2231i
0.2298
p =
-3.4949
4.4438 + 0.1975i
4.4438 - 0.1975i
0.6074
k =
28
由求得的极点，该系统不稳定。
4.
>>z=[-3]
P=[-1,-5,-15]
所以该系统为稳定的。
5.
>>num=100*conv([1,0],conv([1,2],conv([1,2],conv([1,3,2],[1,3,2]))));
den=conv([1,1],conv([1,-1],conv([1,3,5,2],conv([1,0,2,0,4],[1,0,2,0,4]))));
[z,p,k]=tf2zp(num,den)
实验结果：
z =
0
-2.0005 + 0.0005i
-2.0005 - 0.0005i
-1.9995 + 0.0005i
-1.9995 - 0.0005i
-1.0000 + 0.0000i
-1.0000 - 0.0000i
p =
1.0000
0.7071 + 1.2247i
0.7071 - 1.2247i
0.7071 + 1.2247i
0.7071 - 1.2247i
-1.2267 + 1.4677i
-1.2267 - 1.4677i
-0.7071 + 1.2247i
-0.7071 - 1.2247i
-0.7071 + 1.2247i
-0.7071 - 1.2247i
-1.0000
-0.5466
>>zplane(z,p)
所以该系统不稳定。
七．已知反馈系统开环转移函数如下，试作其奈奎斯特图，并判断系统是否稳定。
1.
>> b=[1];
a=[1,3,2];
sys=tf(b,a);
nyquist(sys);
实验结果：
由于奈奎斯特图并未围绕上-1点运动，同时其开环转移函数也是稳定的，由此，该线性负反馈系统也是稳定的。
2．
>> b=[1];
a=[1,4,4,0];
sys=tf(b,a);
nyquist(sys);
实验结果：
由于奈奎斯特图并未围绕上-1点运动，同时其开环转移函数也是稳定的，由此，该线性负反馈系统也是稳定的。
3.
>> b=[1];
a=[1,2,2];
sys=tf(b,a);
nyquist(sys);
实验结果：
由于奈奎斯特图并未围绕上-1点运动，同时其开环转移函数也是稳定的，由此，该线性负反馈系统也是稳定的。
练习三
实验三
五．
1．
>>help window
WINDOW Window function gateway.
WINDOW(@WNAME,N) returns an N-point window of type specified
by the function handle @WNAME in a column vector. @WNAME can
be any valid window function name, for example:
@bartlett - Bartlett window.
@barthannwin - Modified Bartlett-Hanning window.
@blackman - Blackman window.
@blackmanharris - Minimum 4-term Blackman-Harris window.
@bohmanwin - Bohman window.
@chebwin - Chebyshev window.
@flattopwin - Flat Top window.
@gausswin - Gaussian window.
@hamming - Hamming window.
@hann - Hann window.
@kaiser - Kaiser window.
@nuttallwin - Nuttall defined minimum 4-term Blackman-Harris window.
@parzenwin - Parzen (de la Valle-Poussin) window.
@rectwin - Rectangular window.
@tukeywin - Tukey window.
@triang - Triangular window.
WINDOW(@WNAME,N,OPT) designs the window with the optional input argument
specified in OPT. To see what the optional input arguments are, see the help
for the individual windows, for example, KAISER or CHEBWIN.
WINDOW launches the Window Design & Analysis Tool (WinTool).
EXAMPLE:
N = 65;
w = window(@blackmanharris,N);
w1 = window(@hamming,N);
w2 = window(@gausswin,N,2.5);
plot(1:N,[w,w1,w2]); axis([1 N 0 1]);
legend('Blackman-Harris','Hamming','Gaussian');
See also bartlett, barthannwin, blackman, blackmanharris, bohmanwin,
chebwin, gausswin, hamming, hann, kaiser, nuttallwin, parzenwin,
rectwin, triang, tukeywin, wintool.
Overloaded functions or methods (ones with the same name in other directories)
help fdesign/window.m
Reference page in Help browser
doc window
2.
>>N = 128;
w = window(@rectwin,N);
w1 = window(@bartlett,N);
w2 = window(@hamming,N);
plot(1:N,[w,w1,w2]); axis([1 N 0 1]);
legend('矩形窗','Bartlett','Hamming');
3.
>>wvtool(w,w1,w2)
六．
ts=0.01;
N=20;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(1):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
ts=0.01;
N=30;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(2):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
ts=0.01;
N=50;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(3):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
ts=0.01;
N=100;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(4):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
ts=0.01;
N=150;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(5):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
实验八
1．
%冲激响应
>> clear;
b=[1,3];
a=[1,3,2];
sys=tf(b,a);
impulse(sys);
结果：
%求零输入响应
>> A=[1,3;0,-2];
B=[1;2];
Q=A\B
Q =
4
-1
>> clear
B=[1,3];
A=[1,3,2];
[a,b,c,d]=tf2ss(B,A)
sys=ss(a,b,c,d);
x0=[4;-1];
initial(sys,x0);
grid;
a =
-3 -2
1 0
b =
1
0
c =
1 3
d =
0
2.
%冲激响应
>> clear;
b=[1,3];
a=[1,2,2];
sys=tf(b,a);
impulse(sys)
%求零输入响应
>> A=[1,3;1,-2];
B=[1;2];
Q=A\B
Q =
1.6000
-0.2000
>> clear
B=[1,3];
A=[1,2,2];
[a,b,c,d]=tf2ss(B,A)
sys=ss(a,b,c,d);
x0=[1.6;-0.2];
initial(sys,x0);
grid;
a =
-2 -2
1 0
b =
1
0
c =
1 3
d =
0
3.
%冲激响应
>> clear;
b=[1,3];
a=[1,2,1];
sys=tf(b,a);
impulse(sys)
%求零输入响应
>> A=[1,3;1,-1];
B=[1;2];
Q=A\B
Q =
1.7500
-0.2500
>> clear
B=[1,3];
A=[1,2,1];
[a,b,c,d]=tf2ss(B,A)
sys=ss(a,b,c,d);
x0=[1.75;-0.25];
initial(sys,x0);
grid;
a =
-2 -1
1 0
b =
1
0
c =
1 3
d =
0
二．
>> clear;
b=1;
a=[1,1,1,0];
sys=tf(b,a);
subplot(2,1,1);
impulse(sys);title('冲击响应');
subplot(2,1,2);
step(sys);title('阶跃响应');
t=0:0.01:20;
e=sin(t);
r=lsim(sys,e,t);
figure;
subplot(2,1,1);
plot(t,e);xlabel('Time');ylabel('A');title('激励信号');
subplot(2,1,2);
plot(t,r);xlabel('Time');ylabel('A');title('响应信号');
三．
1.
>> clear;
b=[1,3];
a=[1,3,2];
t=0:0.08:8;
e=[exp(-3*t)];
sys=tf(b,a);
lsim(sys,e,t);
2.
>> clear;
b=[1,3];
a=[1,2,2];
t=0:0.08:8;
sys=tf(b,a);
step(sys)
3．
>> clear;
b=[1,3];
a=[1,2,1];
t=0:0.08:8;
e=[exp(-2*t)];
sys=tf(b,a);
lsim(sys,e,t);
Doc:
1.
>> clear;
B=[1];
A=[1,1,1];
sys=tf(B,A,-1);
n=0:200;
e=5+cos(0.2*pi*n)+2*sin(0.7*pi*n);
r=lsim(sys,e);
stem(n,r);
2.
>> clear;
B=[1,1,1];
A=[1,-0.5,-0.5];
sys=tf(B,A,-1);
e=[1,zeros(1,100)];
n=0:100;
r=lsim(sys,e);
stem(n,r);
练习三
实验三
五．
1．
>>help window
WINDOW Window function gateway.
WINDOW(@WNAME,N) returns an N-point window of type specified
by the function handle @WNAME in a column vector. @WNAME can
be any valid window function name, for example:
@bartlett - Bartlett window.
@barthannwin - Modified Bartlett-Hanning window.
@blackman - Blackman window.
@blackmanharris - Minimum 4-term Blackman-Harris window.
@bohmanwin - Bohman window.
@chebwin - Chebyshev window.
@flattopwin - Flat Top window.
@gausswin - Gaussian window.
@hamming - Hamming window.
@hann - Hann window.
@kaiser - Kaiser window.
@nuttallwin - Nuttall defined minimum 4-term Blackman-Harris window.
@parzenwin - Parzen (de la Valle-Poussin) window.
@rectwin - Rectangular window.
@tukeywin - Tukey window.
@triang - Triangular window.
WINDOW(@WNAME,N,OPT) designs the window with the optional input argument
specified in OPT. To see what the optional input arguments are, see the help
for the individual windows, for example, KAISER or CHEBWIN.
WINDOW launches the Window Design & Analysis Tool (WinTool).
EXAMPLE:
N = 65;
w = window(@blackmanharris,N);
w1 = window(@hamming,N);
w2 = window(@gausswin,N,2.5);
plot(1:N,[w,w1,w2]); axis([1 N 0 1]);
legend('Blackman-Harris','Hamming','Gaussian');
See also bartlett, barthannwin, blackman, blackmanharris, bohmanwin,
chebwin, gausswin, hamming, hann, kaiser, nuttallwin, parzenwin,
rectwin, triang, tukeywin, wintool.
Overloaded functions or methods (ones with the same name in other directories)
help fdesign/window.m
Reference page in Help browser
doc window
2.
>>N = 128;
w = window(@rectwin,N);
w1 = window(@bartlett,N);
w2 = window(@hamming,N);
plot(1:N,[w,w1,w2]); axis([1 N 0 1]);
legend('矩形窗','Bartlett','Hamming');
3.
>>wvtool(w,w1,w2)
六．
ts=0.01;
N=20;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(1):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
ts=0.01;
N=30;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(2):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
ts=0.01;
N=50;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(3):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
ts=0.01;
N=100;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(4):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
ts=0.01;
N=150;
t=0:ts:(N-1)*ts;
x=2*sin(4*pi*t)+5*cos(6*pi*t);
g=fft(x,N);
y=abs(g)/100;
figure(5):plot(0:2*pi/N:2*pi*(N-1)/N,y);
grid;
实验八
1．
%冲激响应
>> clear;
b=[1,3];
a=[1,3,2];
sys=tf(b,a);
impulse(sys);
结果：
%求零输入响应
>> A=[1,3;0,-2];
B=[1;2];
Q=A\B
Q =
4
-1
>> clear
B=[1,3];
A=[1,3,2];
[a,b,c,d]=tf2ss(B,A)
sys=ss(a,b,c,d);
x0=[4;-1];
initial(sys,x0);
grid;
a =
-3 -2
1 0
b =
1
0
c =
1 3
d =
0
2.
%冲激响应
>> clear;
b=[1,3];
a=[1,2,2];
sys=tf(b,a);
impulse(sys)
%求零输入响应
>> A=[1,3;1,-2];
B=[1;2];
Q=A\B
Q =
1.6000
-0.2000
>> clear
B=[1,3];
A=[1,2,2];
[a,b,c,d]=tf2ss(B,A)
sys=ss(a,b,c,d);
x0=[1.6;-0.2];
initial(sys,x0);
grid;
a =
-2 -2
1 0
b =
1
0
c =
1 3
d =
0
3.
%冲激响应
>> clear;
b=[1,3];
a=[1,2,1];
sys=tf(b,a);
impulse(sys)
%求零输入响应
>> A=[1,3;1,-1];
B=[1;2];
Q=A\B
Q =
1.7500
-0.2500
>> clear
B=[1,3];
A=[1,2,1];
[a,b,c,d]=tf2ss(B,A)
sys=ss(a,b,c,d);
x0=[1.75;-0.25];
initial(sys,x0);
grid;
a =
-2 -1
1 0
b =
1
0
c =
1 3
d =
0
二．
>> clear;
b=1;
a=[1,1,1,0];
sys=tf(b,a);
subplot(2,1,1);
impulse(sys);title('冲击响应');
subplot(2,1,2);
step(sys);title('阶跃响应');
t=0:0.01:20;
e=sin(t);
r=lsim(sys,e,t);
figure;
subplot(2,1,1);
plot(t,e);xlabel('Time');ylabel('A');title('激励信号');
subplot(2,1,2);
plot(t,r);xlabel('Time');ylabel('A');title('响应信号');
三．
1.
>> clear;
b=[1,3];
a=[1,3,2];
t=0:0.08:8;
e=[exp(-3*t)];
sys=tf(b,a);
lsim(sys,e,t);
2.
>> clear;
b=[1,3];
a=[1,2,2];
t=0:0.08:8;
sys=tf(b,a);
step(sys)
3．
>> clear;
b=[1,3];
a=[1,2,1];
t=0:0.08:8;
e=[exp(-2*t)];
sys=tf(b,a);
lsim(sys,e,t);
Doc:
1.
>> clear;
B=[1];
A=[1,1,1];
sys=tf(B,A,-1);
n=0:200;
e=5+cos(0.2*pi*n)+2*sin(0.7*pi*n);
r=lsim(sys,e);
stem(n,r);
2.
>> clear;
B=[1,1,1];
A=[1,-0.5,-0.5];
sys=tf(B,A,-1);
e=[1,zeros(1,100)];
n=0:100;
r=lsim(sys,e);
stem(n,r);
本文来源：https://www.2haoxitong.net/k/doc/2d72724f02768e9951e738be.html
文档为doc格式