地震波观测系统的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