读取和写入地震segy文件的matlab代码

function writesegyfile(sgyfile, datain, sampint, numsamps, numtraces, segfmt, byteorder) if nargin < 3 error('Invalid number of arguments: at least 3 required (segyfile, datain, sampint)'); end if nargin < 4 numsamps = size(datain,1); end if nargin < 5 numtraces = size(datain,2); end if nargin < 6 segfmt = 5; end if nargin < 7 byteorder = 'be'; end switch lower(byteorder) case {'be', 'bigendian', 'ieee-be'} platform='ieee-be'; case {'le', 'littleendian', 'ieee-le'} platform = 'ieee-le'; otherwise error(['Invalid byte order parameter: ' byteorder]); end fid = fopen(sgyfile, 'w', platform); textLine{1} = char(['File generated ' datestr(now)]); switch segfmt case 1 textLine{1} = 'IBM Floating point fmt (not implemented)'; dformat = 'float32'; % 4 bytes, should be IBM floating point case 2 textLine{2} = '4-byte integer format'; dformat = 'int32'; % 4 bytes, signed. case 3 textLine{3} = '2-byte integer format'; dformat = 'int16'; % 2 bytes, signed. case 4 textLine{4} = 'WARNING BAD DATA: Fixed point w/gain code format (not implemented)' error('Can not write this format. (Fixed point with gain code.)'); case 5 textLine{5} = 'IEEE Floating point format'; dformat = 'float32'; % 4 bytes presumably IEEE floating point case 8 textLine{8} = '1-byte integer format'; dformat = 'uchar8'; % 4 bytes presumably IEEE floating point otherwise error(['invalid format specified: ', num2str(segfmt)]); end textLine{39} = 'SEG Y rev1'; % ooh, bleeding edge textLine{40} = 'END EBCDIC'; % % Next line should contain 80 spaces textLine{41} = ' '; for i = 1:40 if isempty(textLine{i}) textLine{i} = ' '; end end textHdrBlock = char(textLine); textHdrBlock = [reshape(sprintf('C%2d ',1:40),4,40)' textHdrBlock(1:40,1:76)]; textHdrBlock = ascii2ebcdic(textHdrBlock); count = fwrite(fid, textHdrBlock', 'uchar'); if count ~= 3200 error(['EBCDIC header is too short. Size is only ', num2str(count), '.']); end temp = zeros(27,1); temp(6) = sampint*1000000; temp(8) = numsamps; temp(10) = segfmt; if (fwrite(fid, temp(1:3), 'ulong') ~= 3) | ... (fwrite(fid, temp(4:27), 'ushort') ~= 24) | ... (fwrite(fid, zeros(170,1), 'ushort') ~= 170) ; error('The binary header has been truncated.') end clear temp tracewrit = 1; trcount = 0; traceheader = zeros(120,1); traceheader(floor((29+1)/2)) = 1; % seismic data traceheader(floor((35+1)/2)) = 1; % production data traceheader(floor((115+1)/2)) = numsamps; % number of samples for this trace while tracewrit <= numtraces trcount = fwrite(fid, traceheader, 'ushort'); count = fwrite(fid, datain(:, tracewrit), dformat); if count ~= numsamps error(['Problem writing trace ', num2str(tracewrit), '.']) end tracewrit = tracewrit + 1; end fclose(fid);