首页 > 条据书信 > 其他信函 / 正文
地震波观测系统的MATLAB仿真报告课程设计
2020-12-14 10:04:21 ℃
地震波观测系统的MATLAB仿真
课程名称
数字信号处理
实验项目
题目6 地震波观测系统的MATLAB仿真
指导教师
学
院
光电信息与通信工程
_
专
业
电子信息工程
班级/学号
学生姓名
课设时间
2011-12- 28至2012-1-5
“数字信号处理课程设计”任务书 题目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 资料下载地址
目录
摘要 - 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
- 上一篇:预备党员个人入党思想报告
- 下一篇:X区居委会换届选举工作汇报
猜你喜欢
- 2021-10-05 赴延安革命基地学党史感悟
- 2021-10-05 美丽乡村先进个人事迹
- 2021-10-02 实施人才强国战略神党课讲稿【五篇】
- 2021-05-06 学好高一数学需要注意方面
- 2021-05-05 “防风险、守底线”专题民主生活会组织生活会个对照检查剖析材料
- 2021-04-28 推进“人才强县”战略工作经验总结
- 2021-04-28 县委机要局局长民主生活会个人对照检查材料
- 2021-04-26 建党百年学习-“长征精神”-中国精神范本五篇,(图文)
- 2021-04-26 大型企业集团党委2021党建工作经验材料7篇
- 2021-04-26 2021年公安局领导班子教育整顿专题民主生活会对照检查材料
- 搜索
-
- 2020第三次中央新疆工作座谈会学习心得 09-30
- 民主生活会谈心谈话记录(3篇) 06-15
- 机关、事业单位工作人员工作季度记实个 06-15
- 五年级道德与法治上学期期末质量分析 08-11
- 人才工作存在的问题 06-12
- 关于协助做好相关工作的函 07-21
- 党支部委员会会议记录(15篇) 07-12
- 党政领导干部选拔任用工作中存在的主要 03-24
- 组织生活会征求意见表 06-12
- 地理教师简短精彩自我介绍 地理教师自 01-10
- 11-25国庆70周年庆典晚会 庆典晚会串词
- 11-25办公室礼仪的十大原则 浅谈办公室的电话礼仪
- 01-17用心灵轻轻地歌唱_心灵的歌唱
- 01-17也许你不是我一生的唯一|也许不是我
- 01-17爱了,请珍惜;不爱,趁早放手|爱就珍惜不爱就放手
- 01-17岁月带走的是记忆,但回忆会越来越清晰|有趣又有深意的句子
- 01-17曾经的美好只是曾经,我只想珍惜身边的人|我只想珍惜你
- 01-18从容不惊 [学会笑眼去看世界,不惊不乍,淡定从容]
- 02-03当代大学生学习态度调查报告
- 02-03常用护患英语会话
- 标签列表