数字信号处理(合成地震记录&fft)

数字信号处理实验报告

实验一、地震子波波形显示及一维地震记录合成

一、实验目的

1、认识地震子波（以雷克子波为例），对子波有直观的认识。

2、利用线性褶积公式合成一维地震记录。 二、实验内容 1、 雷克子波：

()()t f e t w m t f m πγπ2cos 2

2

/2-=（零相位子波）、 ()()t f e

t w m t f m πγ

π2sin 2

2/2-=（最小相位子波）

， 其中m f 代表子波的中心频率，γ代表子波宽度,随着γ的增大，子波能量后移，当γ=7时，最小相位子波可视为混合相位子波，这里取m f = 25 Hz ，γ= 3; 2、 根据公式编程实现零相位子波、最小相位子波的波形显示；

3、 设计反射系数)(n r (n=500)，其中0.1)100(=r ，7.0)200(-=r ，5.0)300(=r ，

4.0)400(=r ，6.0)500(=r ，其它为0;

4、 应用褶积公式∑=-=

*=N

m m n w m r n w n r n f 1

)()()()()(合成一维地震记录，

并图形显示；

5、 根据所学知识对实验结果进行分析。 三、实验结果： 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); end

plot(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);

end

plot(w)

xlabel('n');ylabel('w');

title('最小相位子波')

（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);

end

r=4;

for n=1:51

w2(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);

end

r=5;

for n=1:51

w3(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);

end

subplot(1,3,1),plot(w1);

axis([0,55,-0.7,1]);xlabel('n');title('r=3时');

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)=-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

end

end

plot(f)

（2）零相位子波图像：

0100

-0.8

-0.6-0.4

-0.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(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 end end

0100200300400500600

-0.8

-0.6-0.4-0.2

0.20.40.60.8零相位子波的一维地震记录