地震波观测系统的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资料下载地址11/实践环节/数字信号处理课程设计目录摘要- 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) figur