长江大学matlab程序设计报告

MATLAB程序设计报告学院：地球物理与石油资源学院班级：物探姓名：学号：班内编号：指导教师：完成日期：2015年6月30日2一、课程设计题目和要求题目：地震数据（segy格式数据）读写日期：2014年3月27日至2014年4月27日主要内容：1、理解segy格式2、二维动态数组3、将读取的数据显示出来要求：1、源代码格式规范2、运行结果以SURFER的图片形式提交3、独立完成，严禁抄袭，逐一验收二、编程思路分析：1、segy格式的了解：文件前3600字节为格式属性等的描述，跳过3216字节可读取采样率，跳过3220字节可读取采样数等。3600字节之后为每一道检波器检测的信息，每道前240字节为每一道的描述，之后数据才是实际测得数据！所以运用fopen,freed,fseek,ftell,fclose可以获得采样率，采样数等信息。2、读取检波器测的数据并保存在一个矩阵[SISp]中，其中SI为每道采样数，Sp为道数。因为读到的数据是ibm格式，所以我们还要将ibm格式的数据转换成数值格式（调用ibm2num函数）。这里就有两种选择：第一，将数据读完后再进行格式转换，第二，边读数据边转换格式。经过测试第二种方案速度更快所以选择第二种方案。33、数据的显示，这里也有两种方案：第一，如果直接用读到的矩阵[SISp]进行画图会出现一个问题，那就是各道的波形就会叠加在一起，根本得不到地震剖面。所以我们要让各道波形分开，通过n*100*ones[SISp],其中n为道号，由于没有读出道间距我就假设道间距为100，然后n*100*ones[SISp]+[SISp]就相当于每道的波形错开一个道间距，再显示。第二，直接调用wiggle或imagesc显示。三、源代码主程序：segyid=fopen('101.seg','rb');if~segyid{disp('can''topenfile!');exit;};endfseek(segyid,3220,'bof');%读取采样点数SI=fread(segyid,1,'int16','b');fseek(segyid,3216,'bof');%读取采样率Sp=fread(segyid,1,'int16','b');fseek(segyid,0,'eof');%计算总文件字节数BytesInFile=ftell(segyid);Tn=(BytesInFile-3600)/(SI*4+240);%计算道数fclose(segyid);seismic_data=zeros(SI,Tn);segyid=fopen('101.seg','r','b');if~segyid{disp('can''topenfile!');exit;};end%读取地震数据forn=1:Tnfseek(segyid,3600+n*240+(n-1)*SI*4,'bof');seismic=fread(segyid,[1SI],'uint32','b');seismic_data(:,n)=ibm2num(uint32(seismic'));4endfclose(segyid);%h=0:Sp/1000:Sp/1000*(SI-1);%forn=1:Tn%q=100*n*ones(SI,Tn);%c=q(:,n);%seismic_data11=seismic_data1(:,n);%seismic_data2=seismic_data11+c;%end;%plot(seismic_data2,h);%holdon;wiggle(seismic_data);Xlabel('炮间距');Ylabel('采样点时间/s');Title('地震剖面图');%选择第50道得到它的波形和频谱%零相位滤波filter一般滤波data=seismic_data(:,50)';t=0:1/250:4.496;[a,df]=getFrequencySpectrum(data,1/250);[b,a]=butter(4,20/(250/2),'low');%低通yy=filtfilt(b,a,data);%零相位滤波filter一般滤波figure;plot(t,data-500,'r',t,yy+500,'g');xlabel('时间/s');ylabel('振幅/m');legend('干绕波波形','有效波波形',1);title('第50道地震波干绕波和有效波波形图');getFrequencySpectrum(yy,1/250);Ibm2num函数代码：functionx=ibm2num(b)x=repmat(NaN,size(b));sign=bitget(b,32);%getsignfromfirstbitsign=double(sign);%formathexexp=bitand(b,uint32(hex2dec('7f000000')));%getexponentfromfirstbyte,last7bitsexp=bitshift(exp,-24);%formatlongexp=double(exp)-64;%removebiasfromexponent%formathex5frac=bitand(b,uint32(hex2dec('00ffffff')));%getmantissafromlast3bytes%formatlongfrac=double(frac);frac=frac/2^24;x=(1-2*sign).*16.^exp.*frac;err=frac==0&(exp~=-64|sign~=0);%biasremovalisincorrectforzeroifany(err)%TMH19/06/2003disp(['WARNING:',mfilename,'Invalidzeroinput--SuredataareIBMFLOATformatted?'])return;enderr=frac~=0&(frac1/16|frac=1);ifany(err)%TMH19/06/2003disp(['WARNING:',mfilename,'Invalidmantissainput--SuredataareIBMFLOATformatted?'])return;endWiggle函数代码：Functionwiggle(x,t,Data,style,dmax,showmax,plImage,imageax,example_plot);if(nargin==9);np=3;subplot(np,np,1);wiggle(Data);subplot(np,np,2);wiggle(Data,dmax);subplot(np,np,3);wiggle(x,t,Data);subplot(np,np,4);wiggle(x,t,Data,style,dmax);subplot(np,np,5);wiggle(x,t,Data,style,dmax,showmax);subplot(np,np,6);wiggle(x,t,Data,style,dmax,showmax,plImage);ifisempty(dmax),dmax=max(abs(Data(:)));endsubplot(np,np,7);wiggle(x,t,Data,style,dmax,showmax,plImage,dmax./10);returnendshowmax_def=100;style_def='wiggle';ifnargin==1,6Data=x;t=[1:1:size(Data,1)];x=[1:1:size(Data,2)];dmax=max(Data(:));style=style_def;showmax=showmax_def;endifnargin==2,Data=x;dmax=t;t=[1:1:size(Data,1)];x=[1:1:size(Data,2)];style=style_def;showmax=showmax_def;endifnargin==3,style=style_def;dmax=max(abs(Data(:)));showmax=showmax_def;endifnargin==4,dmax=max(abs(Data(:)));showmax=showmax_def;endifnargin==5,showmax=showmax_def;endifnargin7plImage=0;endifisempty(dmax),%Setscalingfactordmaxifemptydmax=max(abs(Data(:)));endifisempty(showmax),showmax=100;endifnargin==7,imageax=[-11].*dmax;endifplImage==1,imagesc(x,t,Data);if(length(imageax)==1)imageax=[-11].*abs(imageax);7endcaxis(imageax);holdonendif(showmax0)iflength(x)1,dx=x(2)-x(1);endntraces=length(x);d=ntraces/showmax;ifd=1;d=1;endd=round(d);dmax=dmax/d;LineWidth=0.0001;EdgeColor=[000];fori=1:d:length(x)xt=dx*Data(:,i)'./dmax;if(strmatch('VA',style)==1)xt1=xt;xt1(find(xt10))=0;f1=fill(x(i)+[xt,fliplr(xt1)],[t,fliplr(t)],[000]);set(f1,'LineWidth',LineWidth)set(f1,'EdgeColor',EdgeColor)%set(f1,'EdgeAlpha',[0]);%GIVESROCKYIMAGESholdonif(strmatch('VA2',style,'exact')==1)xt2=xt;xt2(find(xt20))=0;f2=fill(x(i)+[xt,fliplr(xt2)],[t,fliplr(t)],[100]);set(f2,'LineWidth',LineWidth)set(f2,'EdgeColor',EdgeColor)%set(f2,'EdgeAlpha',[0])endelseplot(xt+x(i),t,'k-','linewidth',.05);endifi==1,holdon;endendendholdoff;tryaxis([min(x)-(x(2)-x(1))max(x)+(x(2)-x(1))min(t)max(t)])catchaxis([min(x)max(x)min(t)max(t)])endset(gca,'Ydir','revers')%getFrequencySpectrum函数得到频谱：8function[a,df]=getFrequencySpectrum(data,dt)N=length(data);t=[1:N]*dt;nfft=2^nextpow2(N);f=1/dt*(0:(nfft/2))/nfft;W=fft(data,nfft);a=abs(W);a=a(1:(nfft/2+1));a=a/max(a);p=angle(W);p=p(1:(nfft/2+1));NN=length(a);df=1/dt/(NN-1)/2;figure;subplot(3,1,1);plot(dt*[0:N-1],data);set(gca,'ygrid','on');xlabel('时间/s');ylabel('振幅/m');title('第50道地震波波形图');subplot(3