gusucode.com > DSISoft是由加拿大地质调查局发布的用于垂直地震剖面(VSP)数据处理的免费软件包 > dsisoftv3/dsisoftv3/dsisoftv3/main/segy2mat.m
function [dataout]=segy2mat(filename,cross_wire_file,platform) %[dataout]=segy2mat(filename,cross_wire_file,platform) % %This function reads data stored in segy format, the oil industry standard %for seismic data exchange, and converts it into a MATLAB variable in DSI %format. This module has been written to use a crosswire file to deliniate %trace header locations as there are many variations on the segy format. %The crosswire file is and ASCII file consisting of numbers only separted %into columns using either spaces or Tabs. MATLAB-style comments can be %included in the file provided that they begin with a percent character '%'. %The file should consist of 5 columns: segy byte to start at (byte 1 is %beginning of a trace header), the format of this element (1 for 'char', %2 for 'int16', and 3 for 'int32'), a multiplier, a value to add to the %element, and the trace header word to store this value in the DSI variable. % %See the following paper for a description of the standard segy format: %Barry, K.M., Cavers, D.A., and Kneale, C.W., Recomended Standards for %Digital Tape Formats, Geophysics, Vol. 40, No. 2 (April 1975), P. 344-352. % %INPUT %filename - name of segy file. Must be enclosed in single quotes. %cross_wire_file - name of crosswire file. Must be enclosed in single quotes. %platform = 'l' for little-endian (PC) files % 'b' for big-endian (Unix) files % %OUTPUT %dataout - data in DSI format % %DSISoft Customized VSP Processing software %written by K.S. Beaty June, 1998 %$Id: segy2mat.m,v 2.0 1999/05/21 18:46:30 mah Exp $ %$Log: segy2mat.m,v $ %Revision 2.0 1999/05/21 18:46:30 mah %Release 2 % %Revision 1.4 1999/01/25 18:58:31 kay %Ignore EBCDIC header %more portable loading of crosswire file % %Revision 1.3 1999/01/12 16:17:34 kay %Fixed a couple of typos... % %Revision 1.2 1999/01/11 19:14:55 kay %added support of ibm->ieee formats. % %Revision 1.1 1999/01/06 19:09:08 kay %Initial revision % % %Copyright (C) 1998 Seismology and Electromagnetic Section/ %Continental Geosciences Division/Geological Survey of Canada % %This library is free software; you can redistribute it and/or %modify it under the terms of the GNU Library General Public %License as published by the Free Software Foundation; either %version 2 of the License, or (at your option) any later version. % %This library is distributed in the hope that it will be useful, %but WITHOUT ANY WARRANTY; without even the implied warranty of %MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU %Library General Public License for more details. % %You should have received a copy of the GNU Library General Public %License along with this library; if not, write to the %Free Software Foundation, Inc., 59 Temple Place - Suite 330, %Boston, MA 02111-1307, USA. % %DSI Consortium %Continental Geosciences Division %Geological Survey of Canada %615 Booth St. %Ottawa, Ontario %K1A 0E9 % %email: dsi@cg.nrcan.gc.ca disp('[dataout]=segy2mat(filename,cross_wire_file,platform)') platform=lower(platform); EBCDICSIZE=3200; %EBCDIC/ASCII header length (bytes) RHLEN=400; %reel header length (bytes) THLEN=240; %trace header length (bytes) % open the files fid=fopen(filename,'r',platform); if fid== -1 error('Unable to open data file'); return; end %if fseek(fid,0,1); endoffile=ftell(fid); %crs=load(cross_wire_file); eval (['load ',cross_wire_file]); eval (['crs=',cross_wire_file(1:length(cross_wire_file)-4),';']); eval (['clear ', cross_wire_file(1:length(cross_wire_file)-4)]); for c=1:size(crs,1) if crs(c,2)==1 crsformat{c}='char'; elseif crs(c,2)==2 crsformat{c}='int16'; else crsformat{c}='int32'; end %if end %for fseek(fid,0,-1); %read reel id header 1 reel_id_header1=fread(fid,EBCDICSIZE,'char'); %ind=find(reel_id_header1~=32); %line_name=char(reel_id_header1(1:max(ind)))'; %disp('1st reel header read'); %disp(line_name); %read reel header 2 (file header) %initialize dataout file header dataout.fh{13}=0; %max record fold fseek(fid,3204,-1); dataout.fh{2}=fread(fid,1,'int32'); %line number fseek(fid,3208,-1); dataout.fh{3}=fread(fid,1,'int32'); %reel number fseek(fid,3212,-1); ntr=fread(fid,1,'int16'); %number of traces per record fseek(fid,3216,-1); smp=fread(fid,1,'int16')/1000000; %sampling rate in seconds dataout.fh{8}=smp; dataout.fh{9}=0; %start time (s) fseek(fid,3220,-1); npts=fread(fid,1,'int16'); %number of points per trace dataout.fh{7}=npts; %number of points per trace dataout.fh{10}=dataout.fh{9}+(npts-1)*smp; %end time (s) fseek(fid,3224,-1); format=fread(fid,1,'int16'); if format==1 dformat='float32'; dformatlen=4; elseif format ==2 dformat='int32'; dformatlen=4; %4 bytes elseif format==3 dformat='int16'; dformatlen=2; %2 bytes else error(['Don"t know how to deal with data format ',num2str(format)]); end %if/else fseek(fid,3254,-1); units=fread(fid,1,'int16'); if units==2 %measurements in feet instead of metres units=0.3048; %conversion factor (m/ft) else units=1; %metres end %if TRLEN=npts*dformatlen; %number of bytes taken up by each trace (data only) %disp('2nd reel header read'); if ntr==0 ntr=1; end %if rec=1; tr=0; while ~(ftell(fid)==endoffile) dataout.th{rec}=zeros(64,ntr); %initialize trace headers dataout.dat{rec}=zeros(npts,ntr); %initialize data matrix for k=1:ntr %read trace header buff=EBCDICSIZE+RHLEN+THLEN*(k-1+tr)+TRLEN*(k-1+tr); %use crosswire file to read in trace headers for c=1:size(crs,1) fseek(fid,buff+crs(c,1)-1,-1); dataout.th{rec}(crs(c,5),k)=fread(fid,1,crsformat{c})*crs(c,3)+crs(c,4); end %for %read trace data buff=buff+THLEN; fseek(fid,buff,-1); if (format==1) ibm1=fread(fid,npts,'uint'); data=ibm2ieee(ibm1); count=npts; else [data,count]=fread(fid,npts,dformat); if count~=npts error(['Trouble reading data for trace ',num2str(k),' of record ',num2str(rec),'.']) end %if end %if dataout.dat{rec}(:,k)=data.*units; end %loop over traces dataout.th{rec}(12,1)=ntr; %number of traces this record dataout.th{rec}(13,:)=1:ntr; %trace sequential number deadtr=find(dataout.th{rec}(6,:)==2); dataout.th{rec}(6,deadtr)=-1; %convert segy kill flags to DSI kill flags tr=tr+ntr; %number of traces read so far rec=rec+1; end %while dataout.fh{12}=rec-1; %number of records read dataout.fh{1}=tr; %total number of traces in file dataout.fh{13}=ntr; %max record fold fclose(fid); disp('done reading file')