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

地震波观测系统的MATLAB仿真课程名称 数字信号处理 实验项目 题目6 地震波观测系统的MATLAB仿真 指导教师 赵双琦 学 院 光电信息与通信工程 _ 专 业 电子信息工程 班级/学号 学生姓名 课设时间 2011-12- 28至2012-1-5 09级“数字信号处理课程设计”任务书题目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.西安：西安交通大学出版社，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.5Hz，通带边界频率为0.1 3.8Hz，通带波纹为1dB,阻带衰减20dB；将宽频带仪器的输出仿真到短周期窄带仪器上；并与窄带仪器的输出进行比较（画图）。绘出窄带仪器的频谱图。长周期地震仪的窄带仪器用低通滤波器表示，其阻带边界频率为0.1Hz，通带边界频率为0.02Hz，通带波纹为1dB，阻带衰减为30dB，将宽频带仪器的输出仿真到长周期窄带仪器上；并与窄带仪器的输出比较。同步骤2作图。五、程序实现close all,clear all,clcload 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); %得到地面运动的频率域表示 endendfigure(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:Nif (abs(H2(i)1.0e-4) %为了防止H值太小将该频率的信号放大 Yf2(i)=Yf(i).*H2(i)./H(i); %得到仿真结果endendx2=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); %得到仿真结果 endendx1=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); endendx1=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,clcload 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-1952http:/lgb.ougz.com.cn/html/auto/old/protel99/yyong1.htm 袁宇波.用Matlab和Protel设计微机保护中Butterworth模拟低通滤波器3曾庆禹.电力系统数字光电量测系统的原理及技术J.电网技术.2001.25(4) pp:1-5