数字信号处理合成地震记录fft.docx

1、数字信号处理合成地震记录fft数字信号处理实验报告实验一、地震子波波形显示及一维地震记录合成一、实验目的1、 认识地震子波（以雷克子波为例），对子波有直观的认识。2、 利用线性褶积公式合成一维地震记录。二、实验内容1、 雷克子波：（零相位子波）、（最小相位子波），其中代表子波的中心频率，代表子波宽度,随着的增大，子波能量后移，当=7时，最小相位子波可视为混合相位子波，这里取 = 25 Hz，= 3;2、 根据公式编程实现零相位子波、最小相位子波的波形显示；3、 设计反射系数(n=500)，其中，其它为0;4、 应用褶积公式合成一维地震记录，并图形显示；5、 根据所学知识对实验结果进行分析。三、

2、实验结果： 1、零相位子波：（1）程序源代码：%编写零相位子波t=0.002;r=3;fm=25;for n=1:51 w(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(2*pi*fm*t*n);endplot(w)xlabel(n)ylabel(w)title(零相位子波)（2）图像：2、最小相位子波：（1）程序源代码：%最小相位子波t=0.002;r=3;fm=25;for n=1:51 w(n)=exp(-(2*pi*fm/r)2*(t*n)2)*sin(2*pi*fm*t*n);endplot(w)xlabel(n);ylabel(w);title(最小相位子波)（

3、2）图像：3、对比不同时的波形图 （1）程序：t=0.002;r=3;fm=25;for n=1:51 w1(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(2*pi*fm*t*n);endr=4;for n=1:51 w2(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(2*pi*fm*t*n);endr=5;for n=1:51 w3(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(2*pi*fm*t*n);endsubplot(1,3,1),plot(w1);axis(0,55,-0.7,1);xlabel(n);title(r=3时

4、);subplot(1,3,2),plot(w2);axis(0,55,-0.7,1);xlabel(n);title(r=4时);subplot(1,3,3),plot(w3);axis(0,55,-0.7,1);xlabel(n);title(r=5时); （2）图像：（3）分析： 代表子波宽度，随着的增大，子波能量后移。4、一维地震记录：（1）零相位子波程序：t=0.002;r=3;fm=25;for n=1:51 w(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(2*pi*fm*t*n);end%设置反射系数r=zeros(500);r(100)=1.0;r(200

5、)=-0.7;r(300)=0.5;r(400)=0.4;r(500)=0.6;%编写褶积公式f=zeros(1,550);for n=1:550 for m=1:500 if(1=(n-m)&(n-m)=51) f(n)=f(n)+r(m)*w(n-m); end endendplot(f)（2）零相位子波图像： （3）最小相位子波程序：t=0.002;r=3;fm=25;for n=1:51 w(n)=exp(-(2*pi*fm/r)2*(t*n)2)*sin(2*pi*fm*t*n);end%设置反射系数r=zeros(1,500);r(100)=1.0;r(200)=-0.7;r(30

6、0)=0.5;r(400)=0.4;r(500)=0.6;%编写褶积公式f=zeros(1,550);for n=1:550 for m=1:500 if(1=(n-m)&(n-m)=51) f(n)=f(n)+r(m)*w(n-m); end endendplot(f) (4)最小相位子波图像：（5）对比零相位子波和最小相位子波的一维地震记录： 放大如下图：局部放大可发现，最小相位子波比零相位子波的地震记录要滞后。最小相位子波的能量要稍小于零相位子波的能量。程序：t=0.002;r=3;fm=25;for n=1:51 w1(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(

7、2*pi*fm*t*n); w2(n)=exp(-(2*pi*fm/r)2*(t*n)2)*sin(2*pi*fm*t*n);endr=zeros(1,500);r(100)=1.0;r(200)=-0.7;r(300)=0.5;r(400)=0.4;r(500)=0.6;f1=zeros(1,550);f2=zeros(1,550);for n=1:550 for m=1:500 if(0(n-m)&(n-m)=51) f1(n)=f1(n)+r(m)*w1(n-m); f2(n)=f2(n)+r(m)*w2(n-m); end endendplot(f1)holdplot(f2,r)tit

8、le(对比最小相位和零相位子波的一维地震记录)四、结果分析： 1、离散雷克子波，取采样间隔一般是2毫秒或者1毫秒，在这里我取了2毫秒。采样点数取51个。 2、代表子波宽度，随着的增大，子波能量后移。 3、褶积公式： 中，r的长度为M=500，w的长度为N=51，则f的长度为N+M-1=550。所以计算过程中应该用二重循环，外层是n从1到550，内层是m从1到500.同时循环过程要满足1=n-m=51。实验二、带通滤波及频谱分析一、实验内容对复合频率信号进行频谱分析，并根据其振幅谱设计带陷滤波器，滤掉某些频率成分。二、实验步骤1、设计某一信号，包含多种频率成分（可用雷克子波）；2、将离散，并应用

9、fft进行频谱分析，绘出振幅谱；3、分析振幅谱有什么特点，在频率域设计带陷滤波器（可加斜坡），以消除某频段（大于62.5Hz）的频率成分，并显示滤波后的振幅谱。（要求绘出滤波器图形）4、将滤波后的信号反变换回时间域，并绘出信号曲线，观察其与原信号的差别。5、根据所学知识对实验结果进行分析。三、实验过程 以主频为25Hz的雷克子波为例。设置的子波长度为200，滤波器长度为51。（1） 不加斜坡的滤波器分析：关于N/2=128对称。但是吉卜斯现象波动明显。程序如下：clft=0.002;r=3;fm=25;for n=1:200 x(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos

10、(2*pi*fm*t*n);endfigure(1)plot(x)title(主频为25Hz的雷克子波)for n=1:256 for n=1:200 x(n)=exp(-(2*pi*fm/r)2*(t*n)2)*cos(2*pi*fm*t*n); end for n=201:256 x(n)=0; endends=fft(x);real_s=real(s);imag_s=imag(s);p=sqrt(real_s.2+imag_s.2);figure(2)plot(p)title(快速傅立叶变换后的振幅谱)ylabel(振幅)fx=62.5; %消除大于62.5Hz的成分N=256;ff=1

11、/(N*t);k1=fx/ff;k2=N-k1;for n=1:256 for n=1:k1 r(n)=1; end for n=k1-1:k2 r(n)=0; end for n=k2+1:N r(n)=1; end xx(n)=r(n)*p(n);endfigure(3)plot(r);title(滤波器图形)xxx=abs(xx);figure(4)plot(xxx);title(滤波后的振幅谱)ylabel(振幅)s_s=ifft(xxx);figure(5)plot(real(s_s)title(滤波后的信号曲线)（2） 加斜坡的滤波器程序：除滤波器处的程序不同外，其余相同，不做赘述

12、。下面附上加斜坡的滤波器的程序如下：for n=1:256 for n=1:k1-5 r(n)=1; end for n=k1-4:k1-1 r(n)=8-k1/4; end for n=k1:k2 r(n)=0; end for n=k2:k2+4 r(n)=k2/4-56; end for n=k2+5:N r(n)=1; end xx(n)=r(n)*p(n);end分析，在0和1过渡时没有直接跳跃，而是设计了一个线性函数，即加了一个斜坡过渡。吉卜斯效应：当频率特性曲线是不连续函数而对滤波因子去有限项时，导致对应的频率特征曲线是一波动的曲线，频率域滤波算子反生畸变。（3） 对比两种滤波器（4）对比两种滤波后的信号曲线（5）对比不同斜率的斜坡的滤波作用：