sgyread.m

这是matlab在地球物理数据处理方面的源码

function [traces,t,headasc,headbin,trahead] = sgyread(sgyfile)
% Program to automatically read a SEGY file produced for the radar data sets.
% Program returns matrix traces with amplitudes
% and time base t.
% Can easily add more variables if required
% Refer to SEG digital tape standards booklet.
% Written by D. Schmitt, April 30, 1996
% The input file must have the full filename including
% the standard .sgy designation at the end.
% To run [tracematrix,timebase,headasc,headbin,trahead] = sgyread('filename');
% tracematrix is the set of all traces in matrix form with each column
% equal to a trace.
% timebase is the corresponding times in seconds for the traces.
% headasc is the ascii header(3200 bytes)
% headbin is the 400 byte binary header which contains important information.
% Note that VISTA does not properly write out everything necessary for
% this file to work (the number of traces in the file) and the program
% must be modified when a Vista SEGY file is examined.
% trahead is the 101 by numtraces matrix which contains the trace header
% information. See the SEG-Y technical manual for additional details.

fclose('all');
% fid = fopen(sgyfile,'r','b')  % for workstation mode
%fid = fopen(sgyfile,'r','l');  % for standard PC Vista mode
fid = fopen(sgyfile,'r','b')

headasc = reshape(fread(fid,3200,'char'),80,40);  
headasc = setstr(headasc');

headbin(1:3) = fread(fid,3,'int32');
headbin(4:197) = fread(fid,194,'int16'); 

numsamps = headbin(8);  % Number of samples per trace
samprate = headbin(6);  % Sample rate in microseconds per sample

if headbin(10) == 1
   dataformat = ['float32'];
   datasize = 4;
elseif headbin(10) == 2
      dataformat = ['int32'];
      datasize = 4;
elseif headbin(10) == 3
         dataformat = ['int16'];
         datasize = 2;
else
         error('Error - Data format not specified in binary header, end program')
end
      
if numsamps == 0
      error('Error - number of samples per trace not specified, cannot find numtrace, end program')
end 
   presentpos = ftell(fid);
   fseek(fid,0,1) % Move to end of file
   numbytes = ftell(fid) % How many bytes in total in file
   fseek(fid,presentpos,-1); % Move back to point before trace info.
   numtrace = (numbytes - 3200 - 400)/(240 + datasize*numsamps);
   
trahead = zeros(71,numtrace);  % There are 71 values assigned
				% in the standard SEG-Y format
traces = zeros(numsamps,numtrace);
            
for i = 1:numtrace
   
	trahead(1:7,i) = fread(fid,7,'int32');
	trahead(8:11,i) = fread(fid,4,'int16');
	trahead(12:19,i) = fread(fid,8,'int32');
	trahead(20:21,i) = fread(fid,2,'int16');
	trahead(22:25,i) = fread(fid,4,'int32');
	trahead(26:101,i) = fread(fid,76,'int16');
	[traces(:,i),count2] = fread(fid,numsamps,dataformat); % Read trace data
end % i

t = 1e-6*[0:samprate:(numsamps-1)*samprate]; % Make timebase in seconds

fclose(fid);