sgyread2.m

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

function [traces,t,headasc,headbin,trahead] = sgyread(sgyfile)
% Program to automatically read a SEGY file.
% 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
% 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

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')  


numtrace = headbin(4)+headbin(5);  % Number of traces in file
if numtrace == 0
	numtrace = 480; % Default OYO channels, data + aux.
	disp(['Warning **** Binary Header did not contain trace num info - default to 48 data + 2 aux files'])
end % if

numsamps = headbin(8);
samprate = headbin(6);

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,'float'); % Read trace data
end % i

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

fclose(fid);