最新范文 方案 计划 总结 报告 体会 事迹 讲话 倡议书 反思 制度 入党

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

日期:2021-01-18  类别:最新范文  编辑:一流范文网  【下载本文Word版

地震波观测系统的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/实践环节/数字信号处理课程设计

目录

摘要-

4

-

正文-

4

-

一、目的-

4

-

二、原理-

4

-

三、要求-

5

-

四、步骤-

5

-

五、程序实现-

6

-

实验结果-

12

-

六、体会-

15

-

参考文献-

15

-

摘要

本文的目的是实现地震波观测系统的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作图。

五、程序实现

close

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

if

(H(i)>1.0e-4)

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

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

end

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

if

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

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

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

%得到仿真结果

end

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

if

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

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

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

%得到仿真结果

end

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

%得到仿真结果

if

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

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

end

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;

close

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

    以上《地震波观测系统的MATLAB仿真报告》范文由一流范文网精心整理,如果您觉得有用,请收藏及关注我们,或向其它人分享我们。转载请注明出处 »一流范文网»最新范文»地震波观测系统的MATLAB仿真报告
‖大家正在看...
设为首页 - 加入收藏 - 关于范文吧 - 返回顶部 - 手机版
Copyright © 一流范文网 如对《地震波观测系统的MATLAB仿真报告》有疑问请及时反馈。All Rights Reserved