数字信号处理合成地震记录簿fft.docx
《数字信号处理合成地震记录簿fft.docx》由会员分享，可在线阅读，更多相关《数字信号处理合成地震记录簿fft.docx（23页珍藏版）》请在冰豆网上搜索。
数字信号处理合成地震记录簿fft
数字信号处理实验报告
实验一、地震子波波形显示及一维地震记录合成
一、实验目的
1、认识地震子波（以雷克子波为例），对子波有直观的认识。
2、利用线性褶积公式合成一维地震记录。
二、实验内容
1、雷克子波：
（零相位子波）、
（最小相位子波），
其中
代表子波的中心频率，
代表子波宽度,随着
的增大，子波能量后移，当
=7时，最小相位子波可视为混合相位子波，这里取
=25Hz，
=3;
2、根据公式编程实现零相位子波、最小相位子波的波形显示；
3、设计反射系数
（n=500），其中
，
，
，
，
，其它为0;
4、应用褶积公式
合成一维地震记录，并图形显示；
5、根据所学知识对实验结果进行分析。
三、实验结果：
1、零相位子波：
（1）程序源代码：
%编写零相位子波
t=0.002;
r=3;
fm=25;
forn=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;
forn=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;
forn=1:
51
w1（n）=exp（-（2*pi*fm/r）^2*（t*n）^2）*cos（2*pi*fm*t*n）;
end
r=4;
forn=1:
51
w2（n）=exp（-（2*pi*fm/r）^2*（t*n）^2）*cos（2*pi*fm*t*n）;
end
r=5;
forn=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;
forn=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）;
forn=1:
550
form=1:
500
if（1<=（n-m）&&（n-m）<=51）
f（n）=f（n）+r（m）*w（n-m）;
end
end
end
plot（f）
（2）零相位子波图像：
（3）最小相位子波程序：
t=0.002;
r=3;
fm=25;
forn=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）;
forn=1:
550
form=1:
500
if（1<=（n-m）&&（n-m）<=51）
f（n）=f（n）+r（m）*w（n-m）;
end
end
end
plot（f）
（4）最小相位子波图像：
（5）对比零相位子波和最小相位子波的一维地震记录：
放大如下图：
局部放大可发现，最小相位子波比零相位子波的地震记录要滞后。
最小相位子波的能量要稍小于零相位子波的能量。
程序：
t=0.002;
r=3;
fm=25;
forn=1:
51
w1（n）=exp（-（2*pi*fm/r）^2*（t*n）^2）*cos（2*pi*fm*t*n）;
w2（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;
f1=zeros（1,550）;
f2=zeros（1,550）;
forn=1:
550
form=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
end
end
plot（f1）
hold
plot（f2,'r'）
title（'对比最小相位和零相位子波的一维地震记录'）
四、结果分析：
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、将
离散，并应用fft进行频谱分析，绘出振幅谱；
3、分析振幅谱有什么特点，在频率域设计带陷滤波器（可加斜坡），以消除某频段（大于62.5Hz）的频率成分，并显示滤波后的振幅谱。
（要求绘出滤波器图形）
4、将滤波后的信号反变换回时间域，并绘出信号曲线，观察其与原信号的差别。
5、根据所学知识对实验结果进行分析。
三、实验过程
以主频为25Hz的雷克子波为例。
设置的子波长度为200，滤波器长度为51。
（1）不加斜坡的滤波器
分析：
关于N/2=128对称。
但是吉卜斯现象波动明显。
程序如下：
clf
t=0.002;
r=3;
fm=25;
forn=1:
200
x（n）=exp（-（2*pi*fm/r）^2*（t*n）^2）*cos（2*pi*fm*t*n）;
end
figure
（1）
plot（x）
title（'主频为25Hz的雷克子波'）
forn=1:
256
forn=1:
200
x（n）=exp（-（2*pi*fm/r）^2*（t*n）^2）*cos（2*pi*fm*t*n）;
end
forn=201:
256
x（n）=0;
end
end
s=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/（N*t）;
k1=fx/ff;
k2=N-k1;
forn=1:
256
forn=1:
k1
r（n）=1;
end
forn=k1-1:
k2
r（n）=0;
end
forn=k2+1:
N
r（n）=1;
end
xx（n）=r（n）*p（n）;
end
figure（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）加斜坡的滤波器
程序：
除滤波器处的程序不同外，其余相同，不做赘述。
下面附上加斜坡的滤波器的程序如下：
forn=1:
256
forn=1:
k1-5
r（n）=1;
end
forn=k1-4:
k1-1
r（n）=8-k1/4;
end
forn=k1:
k2
r（n）=0;
end
forn=k2:
k2+4
r（n）=k2/4-56;
end
forn=k2+5:
N
r（n）=1;
end
xx（n）=r（n）*p（n）;
end
分析，在0和1过渡时没有直接跳跃，而是设计了一个线性函数，即加了一个斜坡过渡。
吉卜斯效应：
当频率特性曲线是不连续函数而对滤波因子去有限项时，导致对应的频率特征曲线是一波动的曲线，频率域滤波算子反生畸变。
（3）对比两种滤波器
（4）对比两种滤波后的信号曲线
（5）对比不同斜率的斜坡的滤波作用：
升级会员 