系统仿真实验报告

发布时间:2021-10-23 14:38:19

系统仿真实验报告 本文关键词：仿真，实验，报告，系统

系统仿真实验报告 本文简介：系统仿真实验报告班级：电气工程及其自动化1301班学号：姓名：指导老师：完成时间：2015年4月19日目录实验一MATLAB中矩阵与多项式的基本运算-1-实验二MATLAB绘图命令8实验三MATLAB程序设计10实验四MATLAB的符号计算与SIMULINK的使用11实验五MATLAB在控制系统分析

系统仿真实验报告 本文内容：

系统仿真实验报告

班

级：电气工程及其自动化1301班

学

号：

姓

名：

指导老师：

完成时间：2015年4月19日

目

录

实验一

MATLAB中矩阵与多项式的基本运算-

实验二

MATLAB绘图命令8

实验三

MATLAB程序设计10

实验四

MATLAB的符号计算与SIMULINK的使用11

实验五

MATLAB在控制系统分析中的应用13

实验六

连续系统数字仿真的基本算法22

实验一

MATLAB中矩阵与多项式的基本运算

实验任务

1．了解MATLAB命令窗口和程序文件的调用。

2．熟悉如下MATLAB的基本运算：

①

矩阵的产生、数据的输入、相关元素的显示；

②

矩阵的加法、乘法、左除、右除；

③

特殊矩阵：单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算；

④

多项式的运算：多项式求根、多项式之间的乘除。

【命令】与【运行结果】

eye(3)

ans

ones(3)

ans

zeros(3,4)

ans

rand(3,4)

ans

0.8147

0.9134

0.2785

0.9649

0.9058

0.6324

0.5469

0.1576

0.1270

0.0975

0.9575

0.9706

diag(3)

ans

A=rand(3),diag(A)

0.9572

0.1419

0.7922

0.4854

0.4218

0.9595

0.8003

0.9157

0.6557

ans

0.9572

0.4218

0.6557

A=[1

3;4

6;7

9];.

B=[1

5;7

2;4

8];.

A＼B,B/A,inv(A)*B,B*inv(A),inv(B)*A

Warning:

Matrix

singular

badly

scaled.

Results

may

inaccurate.

RCOND

1.541976e-018.

ans

1.0e+016

4.0532

-4.0532

-8.1065

8.1065

4.0532

-4.0532

Warning:

Matrix

singular

working

precision.

ans

NaN

-Inf

Inf

NaN

-Inf

Inf

NaN

-Inf

Inf

Warning:

Matrix

singular

badly

scaled.

Results

may

inaccurate.

RCOND

1.541976e-018.

ans

1.0e+016

4.0532

-4.0532

-8.1065

8.1065

4.0532

-4.0532

Warning:

Matrix

singular

badly

scaled.

Results

may

inaccurate.

RCOND

1.541976e-018.

ans

1.0e+016

0.0000

4.0532

-8.1065

4.0532

ans

3.5000

3.0000

2.5000

-2.5000

-2.0000

-1.5000

1.0000

【结果分析】

题目中的矩阵A是没有逆的，但由于计算机精度的问题，使得上面的结果比较奇怪：

首先，B/A应该和B*inv（A）都出错显示结果为无定义，但是事实上却是前者结果无定义，而后者却又结果。另外A/B和inv(A)*B

应该也是没有定义的，但还是有结果，因此，我们有必要检查运算结果，鉴于软件在精度方面不可避免的误差。

>>syms

f=3*x^5+4*x^3-5*x^2-7*x+5;

p=[3,0,4,-5,-7,5];

X=roots(p)

poly(X)

-0.3074

1.6176i

-0.3074

1.6176i

-1.0000

1.0000

0.6147

ans

1.0000

0.0000

1.3333

-1.6667

-2.3333

1.6667

A=fix(10*rand(1,3)),B=fix(20*rand(1,4))

conv(A,B)，[Q,r]=deconv(B,A)

ans

123

1.5000

2.9167

1.5833

8.0833

A=fix(10*rand(3)),B=fix(20*rand(3))

A*B，A.*B

ans

110

155

273

342

149

151

ans

135

【结果分析】：A*B是正常的规矩乘法，而A.*B是一种扩展功能，其结果Cij=Aij*Bij。

who

Your

variables

are:

ans

whos

Name

Size

Bytes

Class

Attributes

3x3

double

3x3

double

1x2

double

5x1

double

complex

ans

3x3

double

1x1

170

sym

1x6

double

1x4

double

1x1

126

sym

A=fix(10*rand(1,3)),B=fix(20*rand(3)),C=fix(10*rand(3,2)),D=[3]

diag(A),diag(B),diag(C),diag(D)

ans

disp(

lalalalaalalllallla

);

A=fix(10*rand(3,2)),B=zeros(3),C=ones(3,6)

size(A),length(A),size(B),length(B),size(C),length(C)

lalalalaalalllallla

ans

【结果分析】：disp不受”;分号”影响，size输出矩阵的行数和列数，length输出行列数较大值

实验二

MATLAB绘图命令

实验任务

熟悉MATLAB基本绘图命令，掌握如下绘图方法：

1．坐标系的选择、图形的绘制；

2．图形注解（题目、标号、说明、分格线）的加入；

3．图形线型、符号、颜色的选取。

【程序】

x=linspace(0,100,pi);

y1=x.*x;

y2=x.*x.*x;

subplot(2,2,1);plot(x,y1,-.b,x,y2,r

),legend(

y=x^{2},y=x^{3}

);

grid

off;

subplot(2,2,2);loglog(x,y1,rp-,x,y2,:b

),title(

说啥好呢

),grid

on;

subplot(2,2,3);semilogx(x,y1,b*-,x,y2,r+-

),xlabel(

),ylabel(

);

subplot(2,2,4);semilogy(x,y1,g-+,x,y2,m

),text(50,8000,y=x^{2}

),text(40,100000,y=x^{3}

)

【结果截图】

【程序】

t=0:0.01:2*pi;

r=sin(2*t).*cos(2*t);

subplot(2,2,1);polar(t,r)

subplot(2,2,2);bar(t,r,r

),axis([0,1,-0.5,0.5]),title(

bar

)

subplot(2,1,2);stairs(t,r,c

),axis([0,1,-0.5,0.5]),title(

stairs

)

【结果截图】

【程序】

[x,y,z]=peaks(100);

subplot(1,2,1),mesh(x,y,z)

subplot(1,2,2),contour(x,y,z)

【结果截图】

实验三

MATLAB程序设计

实验任务

用一个MATLAB语言编写一个程序：输入一个自然数，判断它是否是素数，如果是，输出“It

one

prime”，如果不是，输出“It

not

one

prime.”。要求通过调用子函数实现。最好能具有如下功能：①设计较好的人机对话界面，程序中含有提示性的输入输出语句。②能实现循环操作，由操作者输入相关命令来控制是否继续进行素数的判断。如果操作者希望停止这种判断，则可以退出程序。③如果所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。例：输入

6，

因为6不是素数。则程序中除了有“It

not

one

prime”的结论外，还应有：“6=2*3”的说明。

【程序】：

a=input(

请输入一个整数(end

0):

);

while(a~=0)

a==1

disp(

not

prime

composite

number.

)

elseif

isprime(a)==1

fprintf(

prime./n,a)

elseif

isprime(a)==0

f=factor(a);m=length(f);

fprintf(

not

prime,%d=%d,a,a,f(1))

for

i=2:m

fprintf(%d,f(i));

end

fprintf(

);

end

a=input(

请输入一个整数(end

0):

);

end

【运行结果】：

请输入一个整数(end

0):1

not

prime

composite

number.

请输入一个整数(end

0):2

prime.

请输入一个整数(end

0):35

not

prime,35=5*7

请输入一个整数(end

0):0

实验四

MATLAB的符号计算与SIMULINK的使用

实验任务

1．掌握MATLAB符号计算的特点和常用基本命令；

2．掌握SIMULINK的使用。

【程序1】

syms

A=[a

b;c

d];B=[1,4;2,9];

det(A),det(B),eig(A),eig(B)

【运行结果1】

ans

a*d-b*c

ans

1/2*a+1/2*d+1/2*(a^2-2*a*d+d^2+4*b*c)^(1/2)

1/2*a+1/2*d-1/2*(a^2-2*a*d+d^2+4*b*c)^(1/2)

ans

0.1010

9.8990

【程序2】

r=solve(

x^3+5*x^2-9*x+17

)

r4=vpa(r,4)

r10=vpa(r,10)

【运行结果2】

-1/3*(557+3*18849^(1/2))^(1/3)-52/3/(557+3*18849^(1/2))^(1/3)-5/3

1/6*(557+3*18849^(1/2))^(1/3)+26/3/(557+3*18849^(1/2))^(1/3)-5/3+1/2*i*3^(1/2)*(-1/3*(557+3*18849^(1/2))^(1/3)+52/3/(557+3*18849^(1/2))^(1/3))

1/6*(557+3*18849^(1/2))^(1/3)+26/3/(557+3*18849^(1/2))^(1/3)-5/3-1/2*i*3^(1/2)*(-1/3*(557+3*18849^(1/2))^(1/3)+52/3/(557+3*18849^(1/2))^(1/3))

-6.717

.858-1.339*i

.858+1.339*i

r10

-6.716750613

.858375306-1.339469138*i

.858375306+1.339469138*i

【simulink仿真】

【仿真结果】

实验五

MATLAB在控制系统分析中的应用

一、实验任务

1．掌握MATLAB在控制系统时间响应分析中的应用；

2．掌握MATLAB在系统根轨迹分析中的应用；

掌握MATLAB控制系统频率分析中的应用；

掌握MATLAB在控制系统稳定性分析中的应用

求下面系统的单位阶跃响应

【程序】

num=[4]

;

den=[1,1,4]

;step(num,den)

[y,x,t]=step(num,den)

;

tp=spline(y,t,max(y))

%计算峰值时间

ymaxmax(y)

%计算峰值

【结果】

=1.5708

ymax

=1.4419

求下面系统的单位脉冲响应：

【程序】

num=[4]

;

den=[1,1,4]

;

impulse(num,den)

【结果如右图】

3.已知二阶系统的状态方程为：

求系统的零输入响应和脉冲响应。

【程序如下】：

【结果如下图】

a=[0,1

;

-10,-2]

;

b=[0

;

c=[1,0]

;

d=[0]

;

x0=[1,0];

subplot(1,2,1)

;

initial(a,b,c,d,x0)

subplot(1,2,2)

;

impulse(a,b,c,d)

4.系统传递函数为：

输入正弦信号时，观察输出信号的相位差。

【程序如下】：

【结果如下图】：

num=[1]

;

den=[1,1]

;

t=0

0.01

;

u=sin(2*t)

;

hold

plot(t,u,‘r’)

lsim(num,den,u,t)

5.有一二阶系统，求出周期为4秒的方波的输出响应

【程序如下】：

【结果如下图】：

num=[2

1];

den=[1

3];

t=(0:.1:10);

period=4;

u=(rem(t,period)>=period./2);

lsim(num,den,u,t);

6.已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性

【程序如下】：

num=[1

2];den1=[1

3];den=conv(den1,den1);figure(1)

rlocus(num,den)

[k,p]=

rlocfind(num,den)

figure(2)

k=55;num1=k*[1

2];den=[1

3];den1=conv(den,den);[num,den]=cloop(num1,den1,-1);

impulse(num,den)

title(

impulse

response

(k=55)

)

figure(3)

k=56;num1=k*[1

2];den=[1

3];den1=conv(den,den);[num,den]=cloop(num1,den1,-1);

impulse(num,den)

title(

impulse

response(k=56)

)

【结果】：

K=55时系统稳定，k=56时系统发散。

作如下系统的bode图

【程序如下】：

【结果】：

n=[1,1]

;

d=[1,4,11,7]

;

bode(n,d)

8.系统传函如下

求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频率响应

【程序如下】：

num=[1];den=conv([1

2],conv([1

2],[1

2]));

w=logspace(-1,2);

t=0.5;

[m1,p1]=bode(num,den,2);

p1=p1-t*w180/pi;

[n2,d2]=pade(t,4);

numt=conv(n2,num);dent=(conv(den,d2));

[m2,p2]=bode(numt,dent,w);

subplot(2,1,1);

semilogx(w,20*log10(m1),w,20*log10(m2),k--

);

title(

bode

plot

);xlabel(

frequency

);ylabel(

gain

);

subplot(2,1,2);

semilogx(w,p1,w,p2,k--

);

xlabel(

frequency

);ylabel(

phase

);

【结果如下图】：

9.已知系统模型为

求它的幅值裕度和相角裕度

【程序如下】：

n=[3.5];

d=[1

2];

[Gm,Pm,Wcg,Wcp]=margin(n,d)

【结果】：

=1.1433

=7.1688

Wcg

=1.7323

Wcp

=1.6541

10.二阶系统为：

令wn=1,分别作出ξ=2,1,0.707,0.5时的nyquist曲线。

【程序如下】：

n=[1]

;

d1=[1,4,1]

;

d2=[1,2,1]

;

d3=[1,1.414,1];

d4=[1,1,1];

nyquist(n,d1)

;

hold

nyquist(n,d2)

;

nyquist(n,d3)

;

nyquist(n,d4)

;

【结果如下图】：

11.一多环系统，其结构图如下，使用Nyquist频率曲线判断系统的稳定性。

【程序如下】：

k1=16.7/0.0125;z1=[0];

p1=[-1.25

-4

-16];

[num1,den1]=zp2tf(z1,p1,k1);

[num,den]=cloop(num1,den1);

[z,p,k]=tf2zp(num,den);

figure(1)

nyquist(num,den)

figure(2)

[num2,den2]=cloop(num,den);

impulse(num2,den2);

【结果如下图】：

12.

已知系统为：

作该系统的nichols曲线。

【程序如下】：

n=[1]

;

d=[1,1,0]

;

ngrid(‘new’)

;

nichols(n,d)

;

【结果如下图】：

实验六

连续系统数字仿真的基本算法

实验任务

1．理解欧拉法和龙格-库塔法的基本思想；

2．理解数值积分算法的计算精度、速度、稳定性与步长的关系；

取h=0.2，试分别用欧拉法、RK2法和RK4法求解微分方程的数值解，并比较计算精度。

注：解析解：

【程序】

clear

t(1)=0;

置自变量初值

y(1)=1;

y_euler(1)=1;

y_rk2(1)=1;

y_rk4(1)=1;

置解析解和数值解的初值

h=0.2;

步长

求解析解

for

k=1:5

t(k+1)=t(k)+h;

y(k+1)=sqrt(1+2*t(k+1));

end

利用欧拉法求解

for

k=1:5

y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k));

end

利用RK2法求解

for

k=1:5

k1=y_rk2(k)-2*t(k)/y_rk2(k);

k2=(y_rk2(k)+h*k1)-2*(t(k)+h)/(y_rk2(k)+h*k1);

y_rk2(k+1)=y_rk2(k)+h*(k1+k2)/2;

end

利用RK4法求解

for

k=1:5

k1=y_rk4(k)-2*t(k)/y_rk4(k);

k2=(y_rk4(k)+h*k1/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k1/2);

k3=(y_rk4(k)+h*k2/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k2/2);

k4=(y_rk4(k)+h*k3)-2*(t(k)+h)/(y_rk4(k)+h*k3);

y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6;

end

disp(

时间

解析解

欧拉法

RK2法

RK4法

)

yt=[t,y,y_euler,y_rk2,y_rk4

];

disp(yt)

【运行结果】

时间

解析解

欧拉法

RK2法

RK4法

1.0000

0.2000

1.1832

1.2000

1.1867

1.1832

0.4000

1.3416

1.3733

1.3483

1.3417

0.6000

1.4832

1.5315

1.4937

1.4833

0.8000

1.6125

1.6811

1.6279

1.6125

1.0000

1.7321

1.8269

1.7542

1.7321

%H=0.05

时间

解析解

欧拉法

RK2法

RK4法

1.0000

0.0500

1.0488

1.0500

1.0489

1.0488

0.1000

1.0954

1.0977

1.0956

1.0954

0.1500

1.1402

1.1435

1.1403

1.1402

0.2000

1.1832

1.1876

1.1834

1.1832

0.2500

1.2247

1.2301

1.2250

1.2247

H=0.001

时间

解析解

欧拉法

RK2法

RK4法

1.0000

0.0010

1.0010

0.0020

1.0020

0.0030

1.0030

0.0040

1.0040

0.0050

1.0050

%欧拉法需要步长足够小时才逼近解析解，RK4逼近得最快，其次是RK2

考虑如下二阶系统：

在上的数字仿真解（已知：，），

并将不同步长下的仿真结果与解析解进行精度比较。

说明：

已知该微分方程的解析解分别为：

采用RK4法进行计算，选择状态变量：

则有如下状态空间模型及初值条件

采用RK4法进行计算。程序如下：

clear

h=input(

请输入步长h=

);

输入步长

M=round(10/h);

置总计算步数

t(1)=0;

置自变量初值

y_0(1)=100;

y_05(1)=100;

置解析解的初始值(y_0和y_05分别对应于为R=0和R=0.5)

x1(1)=100;

x2(1)=0;

置状态向量初值

y_rk4_0(1)=x1(1);

y_rk4_05(1)=x1(1);

置数值解的初值

求解析解

for

k=1:M

t(k+1)=t(k)+h;

y_0(k+1)=100*cos(t(k+1));

y_05(k+1)=100*exp(-t(k+1)/2).*cos(sqrt(3)/2*t(k+1))+100*sqrt(3)/3*exp(-t(k+1)/2).*sin(sqrt(3)/2*t(k+1));

end

利用RK4法求解

R=0

for

k=1:M

k11=x2(k);

k12=-x1(k);

k21=x2(k)+h*k12/2;

k22=-(x1(k)+h*k11/2);

k31=x2(k)+h*k22/2;

k32=-(x1(k)+h*k21/2);

k41=x2(k)+h*k32;

k42=-(x1(k)+h*k31);

x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6;

x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6;

y_rk4_0(k+1)=x1(k+1);

end

R=0.5

for

k=1:M

k11=x2(k);

k12=-x1(k)-x2(k);

k21=x2(k)+h*k12/2;

k22=-(x1(k)+h*k11/2)-(x2(k)+h*k12/2);

k31=x2(k)+h*k22/2;

k32=-(x1(k)+h*k21/2)-(x2(k)+h*k22/2);

k41=x2(k)+h*k32;

k42=-(x1(k)+h*k31)-(x2(k)+h*k32);

x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6;

x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6;

y_rk4_05(k+1)=x1(k+1);

end

求出误差最大值

err_0=max(abs(y_0-y_rk4_0));

err_05=max(abs(y_05-y_rk4_05));

输出结果

disp(

最大误差（R=0)

最大误差（R=0.5)

)

err_max=[err_0,err_05];

disp(err_max)

步长h取0.0001

到0.5之间变化，并且将数值解直接与解析解比较，求出最大误差数据列入下表中。

步长h

0.0001

0.0005

0.001

0.005

0.01

0.05

0.1

0.5

R=0

最大误差

5.4330×10-10

1.6969×10-10

1.0574×10-10

4.1107×10-9

6.6029×10-8

4.1439×10-5

6.6602×10-4

4.2988×10-1

R=0.5

最大误差

2.7649×10-11

6.8123×10-12

5.3753×10-12

4.0902×10-10

6.5425×10-9

4.1365×10-6

6.7152×10-5

4.5976×10-2

从上表中可以看出，当步长h=0.001时，总误差最小；当步长h小于0.001时，由于舍入误差变大而使总误差增加；当步长h大于0.001时，则由于截断误差的增加也使得总误差加大。另外，当系统的解变化激烈时（如R=0），误差对步长的变化较为敏感；当系统的解变化平稳时，步长的变化对误差的影响就要缓和得多。数值积分算法确定以后，在选择步长时，需要综合考虑。

篇2：操稳仿真实验报告

操稳仿真实验报告 本文关键词：仿真，实验，报告

操稳仿真实验报告 本文简介：西安交通大学汽车试验学实验报告实验名称汽车操稳性仿真学院机械工程学院班级车辆01班学号10015020姓名王放日期2013/6/20汽车操稳性仿真实验报告一、实验目的1.了解和掌握汽车操纵稳定性试验条件、试验规程、数据处理方法及试验仪器设备2.熟悉掌握Adams/Car软件的应用并能实际操作完成汽车

操稳仿真实验报告 本文内容：

西安交通大学

汽车试验学实验报告

实验名称

汽车操稳性仿真

学

院

机械工程学院

班

级

车辆01班

学

号

10015020

姓

名

王

放

日

期

2013/6/20

汽车操稳性仿真实验报告

一、实验目的

了解和掌握汽车操纵稳定性试验条件、试验规程、数据处理方法及试验仪器设备

2.熟悉掌握Adams/Car软件的应用并能实际操作完成汽车操稳性仿真的全过程。

二、实验器材

Adams软件，计算机一台

三、实验步骤

（一）回正性仿真

回正性试验本质是一种方向盘阶跃仿真实验，在回正性仿真试验时，首先执行转弯行驶达到设定的稳态条件，然后将方向盘闭环控制解除，模拟试车员突然松开方向盘以研究汽车在回正力作用下的响应：分析的重点是横摆加速度等变量恢复直线行驶状态时的过渡特性。

1.设置仿真

单击File→Open→Assembly。

在Assembly

Name文本框内单击右键，指向Search（搜索）并选择/assemblies.tbl。

双击MIDI_Demo_Vehicle.asy。

单击OK按钮。

单击Simulate→Full-Vehicle

Analysis→Cornering

Events→Cornering

w/Steer

Release。

按照图示填写对话框。

单击OK，通过Z+左键以及T+左键的调整得到下图。

2.动画处理

按F8进入ADAMS/Postprocessor后处理器界面，选择Animation。

在屏幕上右键，选择Load

Animation。

选择test1_car，单击OK。

点击中的可以播放动画。

录制动画之前要先设定工作目录。从后处理菜单栏选择File>Select

Directory，选定自己的ADAMS

工作目录。

在动画控制栏里选择Record

选项，设置按默认值。

点击开始录制，再点击开始播放，动画开始后选择自己想要的截止的时机，再次按后录制结束。则该动画被保存到自己设定的工作目录里。

3.仿真曲线

在ADAMS/Postprocessor后处理器界面，选择Plotting。

选择test1_car→condition_sensors→steering_wheel_angle。

得到如下曲线。

选择test1_car→chassis_velocities→longitudinal。

得到如下曲线。

选择test1_car→chassis_accelerations→lateral。

得到如下曲线。

（二）收油门转弯仿真

该仿真试验车辆在稳态回转时突然抬起油门踏板，使油门开度为0；然后以设定的角速度继续向回转圆内转向直至汽车的纵向速度小于2.5m/s。仿真分为两个阶段执行，第一阶段在圆环形车道上将车辆加速，使车辆产生设定的侧向加速度，作为收油门转弯试验前的稳态条件；第二阶段关闭油门，发动机怠速或反牵引转速，然后以设定的角速度转动方向盘使汽车的回转半径逐步减小直到仿真结束。

1.设置仿真

单击Simulate→Full-Vehicle

Analysis→Cornering

Events→Lift-Off

Turn-In。

按照图示填写对话框。

单击OK，通过Z+左键以及T+左键的调整得到下图。

2.动画处理

按F8进入ADAMS/Postprocessor后处理器界面，选择Animation。

在屏幕上右键，选择Load

Animation。

选择example_1ti，单击OK。

点击中的可以播放动画。

录制动画之前要先设定工作目录。从后处理菜单栏选择File>Select

Directory，选定自己的ADAMS

工作目录。

在动画控制栏里选择Record

选项，设置按默认值。

点击开始录制，再点击开始播放，动画开始后选择自己想要的截止的时机，再次按后录制结束。则该动画被保存到自己设定的工作目录里。

3.仿真曲线

在ADAMS/Postprocessor后处理器界面，选择Plotting。

选择example_1ti→vas_streering_demand_data→value。

得到如下曲线。

选择example_1it→streering_displacements→angle_front。

得到如下曲线。

选择，单击，对曲线求导，得到如下曲线。

选择example_1it→vas_thorttle_demand_data→value。

得到如下曲线。

（三）弯道收油门仿真

弯道收油门仿真仅试验在稳态回转时油门突然关闭对车辆运动的影响，在关闭油门时，可以设置方向盘锁定或由驱动器调节转向值维持原转弯半径。从这种仿真中可以得到的车辆性能参数有航向偏移、纵向减速度、侧滑角、横摆角、倾斜角等。弯道收油门仿真对回转半径的变化比较敏感；仿真在收油门动作执行到位5s结束。

1.设置仿真

单击Simulate→Full-Vehicle

Analysis→Cornering

Events→Power-Off

Cornering。

按照图示填写对话框。

单击OK，通过Z+左键以及T+左键的调整得到下图。

2.动画处理

按F8进入ADAMS/Postprocessor后处理器界面，选择Animation。

在屏幕上右键，选择Load

Animation。

选择test_poc，单击OK。

点击中的可以播放动画。

录制动画之前要先设定工作目录。从后处理菜单栏选择File>Select

Directory，选定自己的ADAMS

工作目录。

在动画控制栏里选择Record

选项，设置按默认值。

点击开始录制，再点击开始播放，动画开始后选择自己想要的截止的时机，再次按后录制结束。则该动画被保存到自己设定的工作目录里。

3.仿真曲线

在ADAMS/Postprocessor后处理器界面，选择Plotting。

选择test_poc→vas_throttle_demand_data→value。

得到如下曲线。

（四）蛇形穿越

蛇形穿越原指在路面上直线竖立一系列的标竿，让汽车在两侧多次变速穿越行驶。为了在ADAMS中实现，首先必须指定路径和车速。为简化起见，对MDI_Demo_Vehilcle_It只作单速测试：计划使用事件建模器（Event

Builder）建立事件文件后再进行文件驱动仿真。

设车速65km/h，穿越中线间距30m，偏移距3m。

1.在标准界面打开模型MDI_Demo_Vehicle_It，启动Event

Builder。

单击File→Open→Assembly。

右键→Search→/assemblies.tbl。

选择MDI_Demo_Vehicle.asy，单击打开。

单击OK按钮。

单击Simulate→Full-Vehicle

Analysis→Event

Builder。

单击File→New。

新建事件文件wangfang，单击OK按钮激活Event

Builder主界面。

核对单位设置为meter、kilogram、newton、second、degree（默认值）。在全局参数设置部分，输入初始速度为18m/s，档位等于3，static

setup为等速直行（straight，纵向加速度为零），其他默认。

2.定义微操纵

由于试验实际上只需要操作纵转向，用一个微操纵即可定义。按顺序定义转向、油门、制动、离合器、结束条件，转向的控制使用mechine+path

map。

单击Path

Map

Table

Editor进入行驶轨迹编辑器，输入如下XY坐标。

输入完成后单击redraw→fit→OK。

油门、制动控制车速维持初速度，档位和离合器采用开环控制，档位在整个微操纵期间不变，Contorl

Value为3档；离合器保持接合，Contorl

Value为0。

结束条件使用同时组成条件：满足侧向加速度等于0和纵向行驶距离大于250m。

所有选项卡设置完毕后单击保存按钮保存到当前工作目录，退出Event

Builder。

3.执行仿真

使用新建的wangfang.xml运行文件驱动仿真。单击Simulate→Full-Vencle

Analysis→File

Driven

Events。

右键→Browse。

单击OK按钮。得到下图。

按F8进入后处理。

单击wangfang→condition_sensors→yaw_angle，得到下面的曲线。

单击，再单击，单击曲线，得到下面的曲线。

单击wangfang→condition_sensors→roll_angle，得到下面的曲线。

单击，单击曲线，得到下面的曲线。

右键→Load

Animation

选择wangfang。

播放动画。

4.观察仿真结果

按F8键启动后处理模块，观察结果。

四、实验结果

（一）回正性仿真

（二）收油门转弯仿真

（三）弯道收油门仿真

（四）蛇形穿越

篇3：地震波观测系统的MATLAB仿真报告

地震波观测系统的MATLAB仿真报告 本文关键词：地震波，观测，仿真，报告，系统

地震波观测系统的MATLAB仿真报告 本文简介：地震波观测系统的MATLAB仿真课程名称数字信号处理实验项目题目6地震波观测系统的MATLAB仿真指导教师赵双琦学院光电信息与通信工程_专业电子信息工程班级/学号学生姓名课设时间2011-12-28至2012-1-509级“数字信号处理课程设计”任务书题目6地震波观测系统的MATLAB仿真主要内容掌

地震波观测系统的MATLAB仿真报告 本文内容：

地震波观测系统的MATLAB仿真

课程名称

数字信号处理

实验项目

题目6

地震波观测系统的MATLAB仿真

指导教师

赵双琦

学

院

光电信息与通信工程

专

业

电子信息工程

班级/学号

学生姓名

课设时间

2011-12-

28至2012-1-5

09级“数字信号处理课程设计”任务书

题目6

地震波观测系统的MATLAB仿真

主要

内容

掌握地震波观测系统的数字信号处理方法。实现宽频带系统的输出仿真到窄频带输出及地面运动恢复。

设计

要求

以某地震台站记录的地震观测文件为例，选择合适滤波器揭示地面运动恢复和仿真的概念

步骤

1读取地震波观测文件数据，做出时域、频域图形。设计一个包含所有频率成分的宽频带滤波器，假定为宽频带地震仪，恢复地面运动。绘出滤波器频率特性、地面运动时域图。

2已知短周期窄带仪器的阻带边界频率为[0.01

4.5]Hz，通带边界频率为[0.1

3.8]Hz，通带波纹为1dB,阻带衰减20dB;

将宽频带仪器的输出仿真到短周期窄带仪器上；并与窄带仪器的输出进行比较（画图）。绘出窄带仪器的频谱图。

3长周期地震仪的窄带仪器用低通滤波器表示，其阻带边界频率为0.1Hz，通带边界频率为0.02Hz，通带波纹为1dB，阻带衰减为30dB，将宽频带仪器的输出仿真到长周期窄带仪器上；并与窄带仪器的输出比较。同步骤2作图。

主要仪

器设备

1、计算机1台，安装MATLAB软件

主要参

考文献

[美]数字信号处理——使用MATLAB[M].西安：西安交通大学出版社，2002.

课程设计进度计划（起止时间、工作内容）

本课程设计共安排6个题目，这是其中题目之一。整个课程设计共24学时，分1.5周安排，具体进度如下：

4学时

复习题目相关知识，掌握实现的原理；

12学时

用MATLAB语言实现题目要求；

4学时

进一步完善功能，现场检查、答辩；

4学时

完成课程设计报告。

课程设计开始日期

2011.12.26

课程设计完成日期

2012.1.6

课程设计实验室名称

信号与信息处理实验室

地

点

实验楼3-603、605

资料下载地址

http://59.64.74.111/实践环节/数字信号处理课程设计

摘要-

正文-

一、目的-

二、原理-

三、要求-

四、步骤-

五、程序实现-

实验结果-

六、体会-

参考文献-

摘要

本文的目的是实现地震波观测系统的MATLAB仿真。一个线性系统y(t)=h(t)*x(t)，x(t)为地面运动，h(t)为系统的冲击响应，y(t)为系统输出。根据卷积定理，有Y（ω）=H(ω)X(ω)。由地震波观测文件数据y（t），再设计一个宽频带滤波器h（t），就可以恢复地面运动x（t）。对于短周期地震仪，其系统函数为H1（w），对于输入地面运动x（t），有Y1(ω)=H1(ω)X（ω），我们可以推导出Y1（w）=H1(w)Y(w)/H(w)，再对Y1（w）作ifft就可以实现宽频带仪器到短周期窄带仪器的仿真。同样，对长周期地震仪，其系统函数为H2（w），我们也可以得到Y2（w）=H1(w)Y(w)/H(w)，然后对Y2（w）作ifft实现仿真。椭圆滤波器、巴特沃斯滤波器和切比雪夫滤波器的设计都很简单，只要滤波器的指标没问题，调用相应的函数就能实现。仿真的结果请参考本文的正文部分。

正文

一、目的

运用所学数字信号处理的基本知识，掌握地震波观测系统的数字信号处理方法。实现宽频带系统的输出仿真到窄频带输出及地面运动恢复。

二、原理

对于一个线性系统，可以用它的系统函数或脉冲响应来表示

y(t)=h(t)*x(t)

①

式中，x(t)为输入信号，相当于地震观测系统的地面运动；y(t)为系统的输出，相当于地震观测系统的地震记录；h(t)为系统的冲击响应。在频率域内，根据卷积定理，该式可以表示为

Y（ω）=H(ω)X(ω)

②

式中，H(ω)为系统的传递函数，

X(ω)、

Y（ω）为x(t)、y(t)的傅里叶变换。

设想一个频带范围很宽的线性系统，如宽频带地震仪，其系统函数为H(ω);另一个频带较窄的系统，如短周期地震仪，其系统函数为H1(ω)，对于同样的输入X(ω)有

Y(ω)=H(ω)X(ω),Y1(ω)=H1(ω)X（ω）

③

式中，Y1(ω)为频带较窄的系统记录的频谱；H1(ω为频带较窄系统的传递函数。由式③可得

H1(ω)

Y(ω)

Y1(ω)=

H(ω)

④

将上式变换到时间域就得到频带较窄系统的输出y1(t)。也就是说，如果知道宽频带和窄频带系统的传递函数H(ω)和

H1(ω)，原则上可以从宽频带系统的输出推测出窄频带系统的输出。但如果我们知道窄频带系统输出及其两种系统的传递函数，却无法得到宽频带系统的输出。这样就使得我们在记录某种信号时采用宽频带记录，然后仿真到各种窄频带的记录仪器上对信号进行分析。

如果已知地震仪的输出和地震仪的传递函数，我们可以求出地面运动为

X（ω）=

Y(ω)/

H(ω)

⑤

三、要求

以某地震台站记录的地震观测文件为例，选择合适滤波器揭示地面运动恢复和仿真的概念

四、步骤

1、

读取地震波观测文件数据，做出时域、频域图形。设计一个包含所有频率成分的宽频带滤波器，假定为宽频带地震仪，恢复地面运动。绘出滤波器频率特性、地面运动时域图。

2、

已知短周期窄带仪器的阻带边界频率为[0.01

4.5]Hz，通带边界频率为[0.1

3.8]Hz，通带波纹为1dB,阻带衰减20dB；将宽频带仪器的输出仿真到短周期窄带仪器上；并与窄带仪器的输出进行比较（画图）。绘出窄带仪器的频谱图。长周期地震仪的窄带仪器用低通滤波器表示，其阻带边界频率为0.1Hz，通带边界频率为0.02Hz，通带波纹为1dB，阻带衰减为30dB，将宽频带仪器的输出仿真到长周期窄带仪器上；并与窄带仪器的输出比较。同步骤2作图。

五、程序实现

all,clear

all,clc

load

hns.dat

;

%读取数据序列

Xt=hns;

%把数据赋值给变量

Fs=50;

%设定采样率

单位(Hz)

dt=1/Fs;

%求采样间隔

单位(s)

N=length(Xt);

%得到序列的长度

t=[0:N-1]*dt;

%时间序列

Yf=fft(Xt);

%对信号进行快速Fourier变换(FFT)

figure(1);

subplot(2,1,1),plot([0:N-1]/Fs,Xt);

%绘制原始值序列

xlabel(

时间/s

),title(

时间域

);

grid

on;

subplot(2,1,2),plot([0:N-1]/N*Fs,abs(Yf));%绘制信号的振幅谱

xlabel(

频率/Hz

),title(

幅频图

);

ylabel(

振幅

);

xlim([0

2]);

%频率轴只画出2Hz频率之前的部分

grid

on;

%----------设计一个切比雪夫1型宽频带滤波器，假定为宽频带地震仪---------------

ws=[0.00001

25.0]*2/Fs;

%阻带边界频率（归一化频率）

wp=[0.001

25.0]*2/Fs;

%通带边界频率（归一化频率）

Rp=1;Rs=20;Nn=513;

%通带波纹和阻带衰减以及绘制频率特性的数据点数

[Order,Wn]=cheb1ord(wp,ws,Rp,Rs);%求取数字滤波器的最小阶数和归一化截止频率

[b,a]=cheby1(Order,Rp,Wn);

%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器

figure(2);

[H,f]=freqz(b,a,Nn,Fs);

%按传递函数系数、数据点数和采样频率求得滤波器的频率特性

y1=filtfilt(b,a,Xt);

subplot(2,1,1),plot(f,20*log10(abs(H)));

%画出宽带滤波器的幅频特性

xlabel(

/lambda

);ylabel(

A(/lambda)/db

);

title(

宽频带滤波器幅频特性

);grid

on;

subplot(2,1,2),plot(f,angle(H))

%画出宽带滤波器的相频特性

xlabel(

频率/Hz

);ylabel(

相位/^o

);title(

宽频带滤波器相频特性

);grid

on;

%已知宽频带地震仪的频率特性,恢复地面运动

[H,f]=freqz(b,a,N,Fs,whole

);

%得到地震仪的特性

Xf=zeros(1,N);

for

i=1:N

(H(i)>1.0e-4)

Xf(i)=Yf(i)./H(i);

%得到地面运动的频率域表示

end

figure(3);

xt=real(ifft(Xf));

%得到地面运动

subplot(2,1,1);

plot(t,xt,r

);

xlabel(

时间/s

);

ylabel(

振幅

);

title(

地面运动时域图

);

grid

on;

subplot(2,1,2);

plot(t,Xt,g

);

xlabel(

时间/s

);

ylabel(

振幅

);

title(

原始信号

);

grid

on;

%设计一个椭圆宽带滤波器，假定为宽频带地震仪

ws=[0.00001

25.0]*2/Fs;wp=[0.001

25.0]*2/Fs;

%通带和阻带边界频率（归一化频率）

Rp=1;Rs=50;Nn=512;

%通带波纹和阻带衰减以及绘制频率特性的数据点数

[Order,Wn]=ellipord(wp,ws,Rp,Rs);

%求取数字滤波器的最小阶数和归一化截止频率

[b,a]=ellip(Order,Rp,Rs,Wn);

%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器

figure(4)

[H,f]=freqz(b,a,Nn,Fs);

%按传递函数系数、数据点数和采样频率求得滤波器的频率特性

subplot(2,1,1),plot(f,20*log10(abs(H)))

xlabel(

频率/Hz

);ylabel(

振幅/dB

);grid

on;

subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))

xlabel(

频率/Hz

);ylabel(

相位/^o

);grid

on;

y=filtfilt(b,a,Xt);

%在宽带滤波器上的输出

figure(5)

subplot(2,1,1),plot(t,Xt)

xlabel(

时间/s

),title(

输入信号

);

ylabel(

振幅

);

grid

on;

subplot(2,1,2),plot(t,y)

xlabel(

时间/s

),title(

椭圆宽带滤波器输出信号

);

ylabel(

振幅

);

grid

on;

figure(6)

subplot(2,1,1),plot(t,y1,g

);

xlabel(

时间/s

),title(

切比雪夫1型宽频带滤波器输出信号

);

ylabel(

振幅

);

grid

on;

subplot(2,1,2),plot(t,y,r

)

xlabel(

时间/s

),title(

椭圆宽带滤波器输出信号

);

ylabel(

振幅

);

grid

on;

%--------仿真到长周期地震仪上，长周期地震仪用一个巴特沃思滤波器来表示----------

ws=0.1*2/Fs;wp=0.02*2/Fs;

%通带和阻带边界频率（归一化频率）

Rp=1;Rs=30;Nn=512;

%通带波纹和阻带衰减以及绘制频率特性的数据点数

[Order,Wn]=buttord(wp,ws,Rp,Rs);

%求取数字滤波器的最小阶数和归一化截止频率

[b,a]=butter(Order,Wn);

%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器

figure(7);

[H2,f]=freqz(b,a,Nn,Fs);

%按传递函数系数、数据点数和采样频率求得滤波器的频率特性

subplot(2,1,1),plot(f,20*log10(abs(H2)));

xlabel(

/lambda

);ylabel(

A(/lambda)/db

);title(

长周期窄带滤波器幅频特性

);grid

on;

subplot(2,1,2),plot(f,angle(H2));

xlabel(

频率/Hz

);ylabel(

相位/^o

);title(

长周期窄带滤波器相频特性

);grid

on;

figure(8);

y2=filtfilt(b,a,Xt);

%在窄带滤波器上的输出

[H2,f]=freqz(b,a,N,Fs,whole

);

%得到地震仪的特性

Yf2=zeros(1,N);

for

i=1:N

(abs(H2(i))>1.0e-4)

%为了防止H值太小将该频率的信号放大

Yf2(i)=Yf(i).*H2(i)./H(i);

%得到仿真结果

end

x2=ifft(Yf2);

subplot(2,1,1);

plot(t,y2,g

);

%绘制实际输出信号

xlabel(

时间/s

);

ylabel(

振幅

);

title(

长周期地震仪实际输出

);

grid

on;

subplot(2,1,2);

plot(t,real(x2),r

);

%绘制仿真输出信号

title(

长周期地震仪仿真输出

);

xlabel(

时间/s

);

ylabel(

振幅

);

grid

on;

%仿真到长周期地震仪上，长周期地震仪用一个窄带椭圆滤波器来表示

ws=0.1*2/Fs;wp=0.02*2/Fs;

%通带和阻带边界频率（归一化频率）

Rp=1;Rs=30;Nn=512;

%通带波纹和阻带衰减以及绘制频率特性的数据点数

[Order,Wn]=ellipord(wp,ws,Rp,Rs);

%求取数字滤波器的最小阶数和归一化截止频率

[b,a]=ellip(Order,Rp,Rs,Wn);

%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器

figure(9)

y1=filtfilt(b,a,Xt);

%在窄带滤波器上的输出

[H1,f]=freqz(b,a,N,Fs,whole

);

%得到地震仪的特性

XX1=zeros(1,N);

for

ii=1:N

(abs(H1(ii))>1.0e-4)

%为了防止H值太小将该频率的信号放大

XX1(ii)=Yf(ii).*H1(ii)./H(ii);

%得到仿真结果

end

x1=ifft(XX1);

subplot(1,2,1);

plot(t,y1);

title(

实际输出

);

xlabel(

时间/s

);

ylabel(

振幅

);

grid

on;

subplot(1,2,2);

plot(t,real(x1));

title(

仿真输出

);

xlabel(

时间/s

);

ylabel(

振幅

);

grid

on;

figure(10);

subplot(2,1,1),plot(t,y2,g

);

xlabel(

时间/s

),title(

巴特沃思滤波器滤波器输出信号

);

ylabel(

振幅

);

grid

on;

subplot(2,1,2),plot(t,y1,r

);

xlabel(

时间/s

),title(

椭圆宽带滤波器输出信号

);

ylabel(

振幅

);

grid

on;

%仿真到短周期地震仪上，短周期地震仪用一个窄带椭圆滤波器来表示

ws=[0.01

4.5]*2/Fs;wp=[0.1

3.8]*2/Fs;

%通带和阻带边界频率（归一化频率）

Rp=1;Rs=20;Nn=512;

%通带波纹和阻带衰减以及绘制频率特性的数据点数

[order,Wn]=ellipord(wp,ws,Rp,Rs);

%求取数字滤波器的最小阶数和归一化截止频率

[b,a]=ellip(order,Rp,Rs,Wn);

%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器

figure(11)

[H1,f]=freqz(b,a,Nn,Fs);

%按传递函数系数、数据点数和采样频率求得滤波器的频率特性

subplot(2,1,1),plot(f,20*log10(abs(H1)))

xlabel(

频率/Hz

);ylabel(

振幅/dB

);grid

on;

subplot(2,1,2),plot(f,180/pi*unwrap(angle(H1)))

xlabel(

频率/Hz

);ylabel(

相位/^o

);grid

on;

figure(12)

y1=filtfilt(b,a,Xt);

%在窄带滤波器上的输出

[H1,f]=freqz(b,a,N,Fs,whole

);

%得到地震仪的特性

XX1=zeros(1,N);

for

ii=1:N

%得到仿真结果

(abs(H1(ii))>1.0e-4)

XX1(ii)=Yf(ii).*H1(ii)/H(ii);

end

x1=ifft(XX1);

plot(t,y1,t,real(x1),r

)

%绘制输入信号

legend(

实际输出,仿真输出,1)

xlabel(

时间/s

);

ylabel(

振幅

);

grid

on;

%-------仿真到短周期地震仪上，短周期地震仪用一个切比雪夫2型滤波器来表示------

ws=[0.01

4.5]*2/Fs;wp=[0.1

3.8]*2/Fs;

Rp=1;Rs=20;Nn=512;

[Order,Wn]=cheb2ord(wp,ws,Rp,Rs);%求取数字滤波器的最小阶数和归一化截止频率

[b,a]=cheby2(Order,Rp,Wn);%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器

figure(13);

[H,f]=freqz(b,a,Nn,Fs);%按传递函数系数、数据点数和采样频率求得滤波器的频率特性

y3=filtfilt(b,a,Xt);

subplot(2,1,1),plot(f,20*log10(abs(H)));%画出宽带滤波器的幅频特性

xlabel(

/lambda

);ylabel(

A(/lambda)/db

);

title(

宽频带滤波器幅频特性

);

grid

on;

subplot(2,1,2),plot(f,angle(H))

%画出宽带滤波器的相频特性

xlabel(

频率/Hz

);ylabel(

相位/^o

);

title(

宽频带滤波器相频特性

);

grid

on;

figure(14);

subplot(2,1,1),plot(t,y3,g

);

xlabel(

时间/s

),title(

切比雪夫2型滤波器滤波器输出信号

);

ylabel(

振幅

);

grid

on;

subplot(2,1,2),plot(t,y1,r

);

xlabel(

时间/s

),title(

椭圆宽带滤波器输出信号

);

ylabel(

振幅

);

grid

on;

all,clear

all,clc

load

hns1.dat

;

%读取数据序列

Xt=hns1;

%把数据赋值给变量

Fs=50;

%设定采样率

单位(Hz)

dt=1/Fs;

%求采样间隔

单位(s)

N=length(Xt);

%得到序列的长度

t=[0:N-1]*dt;

%时间序列

Yf=fft(Xt);

%对信号进行快速Fourier变换(FFT)

figure(1);

subplot(2,1,1),plot([0:N-1]/Fs,Xt);

%绘制原始值序列

title(

P波

);

xlabel(

时间/s

),title(

时间域

);

title(

P波

);

grid

on;

subplot(2,1,2),plot([0:N-1]/N*Fs,abs(Yf));

%绘制信号的振幅谱

xlabel(

频率/Hz

),title(

幅频图

);

ylabel(

振幅

);

xlim([0

2]);

%频率轴只画出2Hz频率之前的部分

grid

on;

load

hns2.dat

;

%读取数据序列

Xt=hns2;

%把数据赋值给变量

Fs=50;

%设定采样率

单位(Hz)

dt=1/Fs;

%求采样间隔

单位(s)

N=length(Xt);

%得到序列的长度

t=[0:N-1]*dt;

%时间序列

Yf=fft(Xt);

%对信号进行快速Fourier变换(FFT)

figure(2);

subplot(2,1,1),plot([0:N-1]/Fs,Xt);

%绘制原始值序列

title(

S波

);

xlabel(

时间/s

),title(

时间域

);

title(

S波

);

grid

on;

subplot(2,1,2),plot([0:N-1]/N*Fs,abs(Yf));

%绘制信号的振幅谱

xlabel(

频率/Hz

),title(

幅频图

);

ylabel(

振幅

);

xlim([0

2]);

%频率轴只画出2Hz频率之前的部分

grid

on;

load

hns3.dat

;

%读取数据序列

Xt=hns3;

%把数据赋值给变量

Fs=50;

%设定采样率

单位(Hz)

dt=1/Fs;

%求采样间隔

单位(s)

N=length(Xt);

%得到序列的长度

t=[0:N-1]*dt;

%时间序列

Yf=fft(Xt);

%对信号进行快速Fourier变换(FFT)

figure(3);

subplot(2,1,1),plot([0:N-1]/Fs,Xt);

%绘制原始值序列

title(

面波

);

xlabel(

时间/s

),title(

时间域

);

title(

面波

);

grid

on;

subplot(2,1,2),plot([0:N-1]/N*Fs,abs(Yf));%绘制信号的振幅谱

xlabel(

频率/Hz

),title(

幅频图

);

ylabel(

振幅

);

xlim([0

2]);

%频率轴只画出2Hz频率之前的部分

grid

on;

实验结果

图1

切比雪夫1型宽频带滤波器与椭圆宽带滤波器输出信号对比

地面运动时域与原始信号对比

图2

输入信号与输出信号

图3

宽频带振幅与相位

图4

短周期窄带振幅与相位

图5

实际输出与仿真输出对比

图6

巴特沃夫长周期实际输出与仿真输出对比

图7

地震波面波、P波、S波幅频图

图8

长周期窄带滤波器幅频特性

长周期窄带滤波器相频特性

六、实验体会

通过这次实验，我进一步复习了数字信号处理关于滤波器的基础，也了解了理论和实际的不同。在我们身边处处都能看到数字信号处理的相关知识的应用，从语音的识别采集处理到地震波观测，这直观的证实了数字信号处理这门课程的重要性。

在这次实验中，我们在实际操作中加强实践能力，巩固了数字信号处理理论知识，培养了我们解决实际问题的能力，在设计过程中，提高我们的思考能力、动手能力。让我们在学习理论知识的同时，明白如何把这些应用于实际。

这次的课程设计让我认识到了自己的不足，也认识到了我们学习的基础知识究竟能运用于什么领域，如何运用。在老师和同学的耐心指导下我发现了自己在选择巴特沃斯、切比雪夫滤波器上的问题，经过修改和调试，终于得到了应有的效果，这让我看到了理论与实践相结合的优势与用处，让我受益匪浅。

参考文献

[1]焦瑞莉

罗倩

汪毓铎

顾奕.数字信号处理[M].机械工艺出版社.pp:184-195

[2]http://lgb.ougz.com.cn/html/auto/old/protel99/yyong1.htm

袁宇波.用Matlab和Protel设计微机保护中Butterworth模拟低通滤波器

[3]曾庆禹.电力系统数字光电量测系统的原理及技术[J].电网技术.2001.25(4)

pp:1-5