地震波观测系统MATLAB仿真报告课程设计.doc

地震波观测系统MATLAB仿真报告课程设计地震波观测系统的MATLAB仿真 课程名称 数字信号处理 实验工程 题目6 地震波观测系统的MATLAB仿真 指导老师 学 院 光电信息与通信工程 _专 业 电子信息工程 班级/学号 学生姓名 课设时间 2022-12- 28至20_-1-5 “数字信号处理课程设计”任务书 题目6 地震波观测系统的MATLAB仿真 主要 内容 掌握地震波观测系统的数字信号处理方法。实现宽频带系统的输出仿真到窄频带输出及地面运动恢复。设计 要求 要求 以某地震台站记录的地震观测文件为例，选择适宜滤波器提醒地面运动恢复和仿真的概念 步骤 1读取地震波观测文件数据，做出时域、频域图形。设计一个包含所有频率成分的宽频带滤波器，假定为宽频带地震仪，恢复地面运动。绘出滤波器频率特性、地面运动时域图。2短周期窄带仪器的阻带边界频率为0.01 4.5Hz，通带边界频率为0.1 3.8Hz，通带波纹为1dB,阻带衰减20dB; 将宽频带仪器的输出仿真到短周期窄带仪器上；并与窄带仪器的输出进展比拟画图。绘出窄带仪器的频谱图。3长周期地震仪的窄带仪器用低通滤波器表示，其阻带边界频率为0.1Hz，通带边界频率为0.02Hz，通带波纹为1dB，阻带衰减为30dB，将宽频带仪器的输出仿真到长周期窄带仪器上；并与窄带仪器的输出比拟。同步骤2作图。主要仪 器设备 1、计算机1台，安装MATLAB软件 主要参 考文献 美数字信号处理使用MATLABM.西安：西安交通大学出版社，20_2.课程设计进度方案起止时间、工作内容本课程设计共安排6个题目，这是其中题目之一。整个课程设计共24学时，分1.5周安排，详细进度如下：4学时 复习题目相关知识，掌握实现的原理；12学时 用MATLAB语言实现题目要求；4学时 进一步完善功能，现场检查、辩论；4学时 完成课程设计报告。课程设计开场日期 2022.12.26 课程设计完成日期 20_.1.6 课程设计实验室名称 信号与信息处理实验室 地 点 实验楼3-603、605 资料下载地址 目录 【摘要】：p 】： - 4 - 正文 - 4 - 一、目的 - 4 - 二、原理 - 4 - 三、要求 - 5 - 四、步骤 - 5 - 五、程序实现 - 6 - 实验结果 - 12 - 六、体会 - 15 - 【参考文献】：p 】： - 15 - 【摘要】：p 】： 本文的目的是实现地震波观测系统的MATLAB仿真。一个线性系统y(t)=h(t)_(t)，_(t)为地面运动，h(t)为系统的冲击响应，y(t)为系统输出。根据卷积定理，有Y=H()_()。由地震波观测文件数据yt，再设计一个宽频带滤波器ht，就可以恢复地面运动_t。对于短周期地震仪，其系统函数为H1w，对于输入地面运动_t，有Y1()=H1()_，我们可以推导出Y1w=H1(w)Y(w)/H(w)，再对Y1w作ifft就可以实现宽频带仪器到短周期窄带仪器的仿真。同样，对长周期地震仪，其系统函数为H2w，我们也可以得到Y2w=H1(w)Y(w)/H(w)，然后对Y2w作ifft实现仿真。椭圆滤波器、巴特沃斯滤波器和切比雪夫滤波器的设计都很简单，只要滤波器的指标没问题，调用相应的函数就能实现。仿真的结果请参考本文的正文局部。正文 一、目的 运用所学数字信号处理的根本知识，掌握地震波观测系统的数字信号处理方法。实现宽频带系统的输出仿真到窄频带输出及地面运动恢复。二、原理 对于一个线性系统，可以用它的系统函数或脉冲响应来表示 y(t)=h(t)_(t) 式中，_(t)为输入信号，相当于地震观测系统的地面运动；y(t)为系统的输出，相当于地震观测系统的地震记录；h(t)为系统的冲击响应。在频率域内，根据卷积定理，该式可以表示为 Y=H()_() 式中，H()为系统的传递函数， _()、 Y为_(t)、y(t)的傅里叶变换。设想一个频带范围很宽的线性系统，如宽频带地震仪，其系统函数为H();另一个频带较窄的系统，如短周期地震仪，其系统函数为H1()，对于同样的输入_()有 Y()=H()_(), Y1()=H1()_ 式中，Y1()为频带较窄的系统记录的频谱；H1(为频带较窄系统的传递函数。由式可得 H1() Y() Y1()= H() 将上式变换到时间域就得到频带较窄系统的输出y1(t)。也就是说，假如知道宽频带和窄频带系统的传递函数H()和 H1()，原那么上可以从宽频带系统的输出推测出窄频带系统的输出。但假如我们知道窄频带系统输出及其两种系统的传递函数，却无法得到宽频带系统的输出。这样就使得我们在记录某种信号时采用宽频带记录，然后仿真到各种窄频带的记录仪器上对信号进展分析p 。假如地震仪的输出和地震仪的传递函数，我们可以求出地面运动为 _= Y()/ H() 三、要求 以某地震台站记录的地震观测文件为例，选择适宜滤波器提醒地面运动恢复和仿真的概念 四、步骤 1、 读取地震波观测文件数据，做出时域、频域图形。设计一个包含所有频率成分的宽频带滤波器，假定为宽频带地震仪，恢复地面运动。绘出滤波器频率特性、地面运动时域图。2、 短周期窄带仪器的阻带边界频率为0.01 4.5Hz，通带边界频率为0.1 3.8Hz，通带波纹为1dB,阻带衰减20dB；将宽频带仪器的输出仿真到短周期窄带仪器上；并与窄带仪器的输出进展比拟画图。绘出窄带仪器的频谱图。长周期地震仪的窄带仪器用低通滤波器表示，其阻带边界频率为0.1Hz，通带边界频率为0.02Hz，通带波纹为1dB，阻带衰减为30dB，将宽频带仪器的输出仿真到长周期窄带仪器上；并与窄带仪器的输出比拟。同步骤2作图。五、程序实现 close all,clear all,clc load hns.dat ; 读取数据序列 _t=hns; 把数据赋值给变量 Fs=50; 设定采样率 单位(Hz) dt=1/Fs; 求采样间隔 单位(s) N=length(_t); 得到序列的长度 t=0:N-1_dt; 时间序列 Yf=fft(_t); 对信号进展快速Fourier变换(FFT) figure(1); subplot(2,1,1),plot(0:N-1/Fs,_t); 绘制原始值序列 _label(时间/s),title(时间域); grid on; subplot(2,1,2),plot(0:N-1/N_Fs,abs(Yf);绘制信号的振幅谱 _label(频率/Hz),title(幅频图); ylabel(振幅); _lim(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,_t); subplot(2,1,1),plot(f,20_log10(abs(H); 画出宽带滤波器的幅频特性 _label(lambda);ylabel(A(lambda)/db); title(宽频带滤波器幅频特性);grid on; subplot(2,1,2),plot(f,angle(H) 画出宽带滤波器的相频特性 _label(频率/Hz);ylabel(相位/o);title(宽频带滤波器相频特性);grid on; 宽频带地震仪的频率特性,恢复地面运动 H,f=freqz(b,a,N,Fs,whole); 得到地震仪的特性 _f=zeros(1,N); for i=1:N if (H(i)>1.0e-4) _f(i)=Yf(i)./H(i); 得到地面运动的频率域表示 end end figure(3); _t=real(ifft(_f); 得到地面运动 subplot(2,1,1); plot(t,_t,r); _label(时间/s); ylabel(振幅); title(地面运动时域图); grid on; subplot(2,1,2); plot(t,_t,g); _label(时间/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) _label(频率/Hz);ylabel(振幅/dB);grid on; subplot(2,1,2),plot(f,180/pi_unwrap(angle(H) _label(频率/Hz);ylabel(相位/o);grid on; y=filtfilt(b,a,_t); 在宽带滤波器上的输出 figure(5) subplot(2,1,1),plot(t,_t) _label(时间/s),title(输入信号); ylabel(振幅); grid on; subplot(2,1,2),plot(t,y) _label(时间/s),title(椭圆宽带滤波器输出信号); ylabel(振幅); grid on; figure(6) subplot(2,1,1),plot(t,y1,g); _label(时间/s),title(切比雪夫1型宽频带滤波器输出信号); ylabel(振幅); grid on; subplot(2,1,2),plot(t,y,r) _label(时间/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); _label(lambda);ylabel(A(lambda)/db);title(长周期窄带滤波器幅频特性);grid on; subplot(2,1,2),plot(f,angle(H2); _label(频率/Hz);ylabel(相位/o);title(长周期窄带滤波器相频特性);grid on; figure(8); y2=filtfilt(b,a,_t); 在窄带滤波器上的输出 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 _2=ifft(Yf2); subplot(2,1,1); plot(t,y2,g); 绘制实际输出信号 _label(时间/s); ylabel(振幅); title(长周期地震仪实际输出); grid on; subplot(2,1,2); plot(t,real(_2),r); 绘制仿真输出信号 title(长周期地震仪仿真输出); _label(时间/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,_t); 在窄带滤波器上的输出 H1,f=freqz(b,a,N,Fs,whole); 得到地震仪的特性 _1=zeros(1,N); for ii=1:N if (abs(H1(ii)>1.0e-4) 为了防止H值太小将该频率的信号放大 _1(ii)=Yf(ii)._H1(ii)./H(ii); 得到仿真结果 end end _1=ifft(_1); subplot(1,2,1); plot(t,y1); title(实际输出); _label(时间/s); ylabel(振幅); grid on; subplot(1,2,2); plot(t,real(_1); title(仿真输出); _label(时间/s); ylabel(振幅); grid on; figure(10); subplot(2,1,1),plot(t,y2,g); _label(时间/s),title(巴特沃思滤波器滤波器输出信号); ylabel(振幅); grid on; subplot(2,1,2),plot(t,y1,r); _label(时间/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) _label(频率/Hz);ylabel(振幅/dB);grid on; subplot(2,1,2),plot(f,180/pi_unwrap(angle(H1) _label(频率/Hz);ylabel(相位/o);grid on; figure(12) y1=filtfilt(b,a,_t); 在窄带滤波器上的输出 H1,f=freqz(b,a,N,Fs,whole); 得到地震仪的特性 _1=zeros(1,N); for ii=1:N 得到仿真结果 if (abs(H1(ii)>1.0e-4) _1(ii)=Yf(ii)._H1(ii)/H(ii); end end _1=ifft(_1); plot(t,y1,t,real(_1),r) 绘制输入信号 legend(实际输出,仿真输出,1) _label(时间/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,_t); subplot(2,1,1),plot(f,20_log10(abs(H);画出宽带滤波器的幅频特性 _label(lambda);ylabel(A(lambda)/db); title(宽频带滤波器幅频特性); grid on; subplot(2,1,2),plot(f,angle(H) 画出宽带滤波器的相频特性 _label(频率/Hz);ylabel(相位/o); title(宽频带滤波器相频特性); grid on; figure(14); subplot(2,1,1),plot(t,y3,g); _label(时间/s),title(切比雪夫2型滤波器滤波器输出信号); ylabel(振幅); grid on; subplot(2,1,2),plot(t,y1,r); _label(时间/s),title(椭圆宽带滤波器输出信号); ylabel(振幅); grid on; close all,clear all,clc load hns1.dat ; 读取数据序列 _t=hns1; 把数据赋值给变量 Fs=50; 设定采样率 单位(Hz) dt=1/Fs; 求采样间隔 单位(s) N=length(_t); 得到序列的长度 t=0:N-1_dt; 时间序列 Yf=fft(_t); 对信号进展快速Fourier变换(FFT) figure(1); subplot(2,1,1),plot(0:N-1/Fs,_t); 绘制原始值序列 title(P波); _label(时间/s),title(时间域); title(P波); grid on; subplot(2,1,2),plot(0:N-1/N_Fs,abs(Yf); 绘制信号的振幅谱 _label(频率/Hz),title(幅频图); ylabel(振幅); _lim(0 2); 频率轴只画出2Hz频率之前的局部 grid on; load hns2.dat ; 读取数据序列 _t=hns2; 把数据赋值给变量 Fs=50; 设定采样率 单位(Hz) dt=1/Fs; 求采样间隔 单位(s) N=length(_t); 得到序列的长度 t=0:N-1_dt; 时间序列 Yf=fft(_t); 对信号进展快速Fourier变换(FFT) figure(2); subplot(2,1,1),plot(0:N-1/Fs,_t); 绘制原始值序列 title(S波); _label(时间/s),title(时间域); title(S波); grid on; subplot(2,1,2),plot(0:N-1/N_Fs,abs(Yf); 绘制信号的振幅谱 _label(频率/Hz),title(幅频图); ylabel(振幅); _lim(0 2); 频率轴只画出2Hz频率之前的局部 grid on; load hns3.dat ; 读取数据序列 _t=hns3; 把数据赋值给变量 Fs=50; 设定采样率 单位(Hz) dt=1/Fs; 求采样间隔 单位(s) N=length(_t); 得到序列的长度 t=0:N-1_dt; 时间序列 Yf=fft(_t); 对信号进展快速Fourier变换(FFT) figure(3); subplot(2,1,1),plot(0:N-1/Fs,_t); 绘制原始值序列 title(面波); _label(时间/s),title(时间域); title(面波); grid on; subplot(2,1,2),plot(0:N-1/N_Fs,abs(Yf);绘制信号的振幅谱 _label(频率/Hz),title(幅频图); ylabel(振幅); _lim(0 2); 频率轴只画出2Hz频率之前的局部 grid on; 实验结果 图1 切比雪夫1型宽频带滤波器与椭圆宽带滤波器输出信号比照 地面运动时域与原始信号比照 图2 输入信号与输出信号 图3 宽频带振幅与相位 图4 短周期窄带振幅与相位 图5 实际输出与仿真输出比照 图6 巴特沃夫长周期实际输出与仿真输出比照 图7 地震波面波、P波、S波幅频图 图8 长周期窄带滤波器幅频特性 长周期窄带滤波器相频特性 六、实验体会 通过这次实验，我进一步复习了数字信号处理关于滤波器的根底，也理解了理论和实际的不同。在我们身边处处都能看到数字信号处理的相关知识的应用，从语音的识别采集处理到地震波观测，这直观的证实了数字信号处理这门课程的重要性。在这次实验中，我们在实际操作中加强理论才能，稳固了数字信号处理理论知识，培养了我们解决实际问题的才能，在设计过程中，进步我们的考虑才能、动手才能。让我们在学习理论知识的同时，明白如何把这些应用于实际。这次的课程设计让我认识到了自己的缺乏，也认识到了我们学习的根底知识终究能运用于什么领域，如何运用。在老师和同学的耐心指导下我发现了自己在选择巴特沃斯、切比雪夫滤波器上的问题，经过修改和调试，终于得到了应有的效果，这让我看到了理论与理论相结合的优势与用途，让我受益匪浅。【参考文献】：p 】： 1焦瑞莉 罗倩 汪毓铎 顾奕.数字信号处理M.机械工艺出版社.pp:184-195 2/auto/old/protel99/yyong1 袁宇波.用Matlab和Protel设计微机保护中Butterworth模拟低通滤波器 3曾庆禹.电力系统数字光电量测系统的原理及技术J.电网技术.20_1.25(4) pp:1-5 第 10 页 共 10 页