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

数字信号处理实验报告实验一、地震子波波形显示及一维地震记录合成一、实验目的1、认识地震子波（以雷克子波为例），对子波有直观的认识。2、利用线性褶积公式合成一维地震记录。二、实验内容1、雷克子波：tfetwmtfm2cos22/2（零相位子波）、tfetwmtfm2sin22/2（最小相位子波），其中mf代表子波的中心频率，代表子波宽度,随着的增大，子波能量后移，当=7时，最小相位子波可视为混合相位子波，这里取mf=25Hz，=3;2、根据公式编程实现零相位子波、最小相位子波的波形显示；3、设计反射系数)(nr(n=500)，其中0.1)100(r，7.0)200(r，5.0)300(r，4.0)400(r，6.0)500(r，其它为0;4、应用褶积公式Nmmnwmrnwnrnf1)()()()()(合成一维地震记录，并图形显示；5、根据所学知识对实验结果进行分析。三、实验结果：1、零相位子波：（1）程序源代码：%编写零相位子波t=0.002;r=3;fm=25;forn=1:51w(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;forn=1:51w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*sin(2*pi*fm*t*n);endplot(w)xlabel('n');ylabel('w');title('最小相位子波')（2）图像：3、对比不同时的波形图（1）程序：t=0.002;r=3;fm=25;forn=1:51w1(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);endr=4;forn=1:51w2(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);endr=5;forn=1:51w3(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时');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:51w(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:550form=1:500if(1=(n-m)&&(n-m)=51)f(n)=f(n)+r(m)*w(n-m);endendendplot(f)（2）零相位子波图像：0100200300400500600-0.8-0.6-0.4-0.200.20.40.60.81零相位子波的一维地震记录（3）最小相位子波程序：t=0.002;r=3;fm=25;forn=1:51w(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:550form=1:500if(1=(n-m)&&(n-m)=51)f(n)=f(n)+r(m)*w(n-m);endendend0100200300400500600-0.8-0.6-0.4-0.200.20.40.60.81零相位子波的一维地震记录plot(f)(4)最小相位子波图像：（5）对比零相位子波和最小相位子波的一维地震记录：0100200300400500600-0.6-0.4-0.200.20.40.60.8最小相位子波的一维地震记录0100200300400500600-0.8-0.6-0.4-0.200.20.40.60.81对比最小相位和零相位子波的一维地震记录零相位最小相位250300350400450500-0.3-0.2-0.100.10.20.30.40.5对比最小相位和零相位子波的一维地震记录零相位最小相位放大如下图：局部放大可发现，最小相位子波比零相位子波的地震记录要滞后。最小相位子波的能量要稍小于零相位子波的能量。程序：t=0.002;r=3;fm=25;forn=1:51w1(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);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);forn=1:550form=1:500if(0(n-m)&&(n-m)=51)f1(n)=f1(n)+r(m)*w1(n-m);f2(n)=f2(n)+r(m)*w2(n-m);endendendplot(f1)holdplot(f2,'r')title('对比最小相位和零相位子波的一维地震记录')四、结果分析：1、离散雷克子波，取采样间隔一般是2毫秒或者1毫秒，在这里我取了2毫秒。采样点数取51个。2、代表子波宽度，随着的增大，子波能量后移。3、褶积公式：Nmmnwmrnwnrnf1)()()()()(中，r的长度为M=500，w的长度为N=51，则f的长度为N+M-1=550。所以计算过程中应该用二重循环，外层是n从1到550，内层是m从1到500.同时循环过程要满足1=n-m=51。实验二、带通滤波及频谱分析一、实验内容对复合频率信号进行频谱分析，并根据其振幅谱设计带陷滤波器，滤掉某些频率成分。二、实验步骤1、设计某一信号)(tx，包含多种频率成分（可用雷克子波tfetwmtfm2cos22/2）；2、将)(tx离散，并应用fft进行频谱分析，绘出振幅谱；3、分析振幅谱有什么特点，在频率域设计带陷滤波器（可加斜坡），以消除某频段（大于62.5Hz）的频率成分，并显示滤波后的振幅谱。（要求绘出滤波器图形）4、将滤波后的信号反变换回时间域，并绘出信号曲线，观察其与原信号的差别。5、根据所学知识对实验结果进行分析。三、实验过程以主频为25Hz的雷克子波为例。设置的子波长度为200，滤波器长度为51。（1）不加斜坡的滤波器020406080100120140160180200-0.500.51主频为25Hz的雷克子波05010015020025030000.511.522.533.54快速傅立叶变换后的振幅谱振幅05010015020025030000.10.20.30.40.50.60.70.80.91滤波器图形05010015020025030000.511.522.533.54滤波后的振幅谱振幅050100150200250300-0.3-0.2-0.100.10.20.30.40.50.60.7滤波后的信号曲线分析：关于N/2=128对称。但是吉卜斯现象波动明显。程序如下：clft=0.002;r=3;fm=25;forn=1:200x(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);endfigure(1)plot(x)title('主频为25Hz的雷克子波')forn=1:256forn=1:200x(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);endforn=201:256x(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/(N*t);k1=fx/ff;k2=N-k1;forn=1:256forn=1:k1r(n)=1;endforn=k1-1:k2r(n)=0;endforn=k2+1:Nr(n)=1;endxx(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）加斜坡的滤波器05010015020025030000.10.20.30.40.50.60.70.80.91滤波器图形05010015020025030000.511.522.533.54滤波后的振幅谱振幅050100150200250300-0.3-0.2-0.100.10.20.30.40.50.60.7滤波后的信号曲线程序：除滤波器处的程序不同外，其余相同，不做赘述。下面附上加斜坡的滤波器的程序如下：forn=1:256forn=1:k1-5r(n)=1;endforn=k1-4:k1-1r(n)=8-k1/4;endforn=k1:k2r(n)=0;endforn=k2:k2+4r(n)=k2/4-56;endforn=k2+5:Nr(n)=1;endxx(n)=r(n)*p(n);end分析，在0和1过渡时没有直接跳跃，而是设计了一个线性函数，即加了一个斜坡过渡。吉卜斯效应：当频率特性曲线是不连续函数而对滤波因子去有限项时，导致对应的频率特征曲线是一波动的曲线，频率域滤波算子反生畸变。（3）对比两种滤波器05010015020025030000.511.522.533.54加斜坡不加斜坡（4）对比两种滤波后的信号曲线050100150200250300-0.3-0.2-0.100.10.20.30.40.50.60.7对比带陷滤波器和非带陷滤波器作用加斜坡不加斜坡050100150200250300-0.500.51对比带斜坡不带斜坡原始雷克子波10203040506070-0.3-0.25-0.2-0.15-0.1-0.0500.05局部放大带斜坡不带斜坡原始雷克子波（5）对比不同斜率的斜坡的滤波作用：050100150200250300-0.3-0.2-0.100.10.20.30.40.50.60.7对比斜率不同的斜坡的滤波作用斜坡的斜率为1/7斜坡的斜率为1/41020304050607080-0.15-0.1-0.0500.05对比斜率不同的斜坡的滤波作用（放大）斜坡的斜率为1/7斜坡的斜率为1/4