This MATLAB script was used to do frequency analysis on 1 mil
"latest.shot.segy". The script doesn't care about the shot number.
It simply reads the trace specified.
#!/bin/csh -f
if( $#argv != 2 ) then
echo "Usage: fanal filename fno ftr"
exit 1
endif
set FILE = $1
set FTR = $2
matlab << eof1
a = (0:.9765625:500)';
rdtrsegy $FILE -1 $FTR;
x = fft(ans,1024);
b = abs(x(1:513,1));
plot(a,b,'b');
xlabel ('frequency (Hz.)');
title ('file $FILE trace $FTR');
pause(120)
print -dpsc plot.ps;
eof1
A similar script for 2 mil data is:
#!/bin/csh -f
if( $#argv != 2 ) then
echo "Usage: fanal filename fno ftr"
exit 1
endif
set FILE = $1
set FTR = $2
matlab << eof1
a = (0:.4882812:250)';
rdtrsegy $FILE -1 $FTR;
x = fft(ans,1024);
b = abs(x(1:513,1));
plot(a,b,'b');
xlabel ('frequency (Hz.)');
title ('file $FILE trace $FTR');
pause(20)
print -dpsc plot.ps;
eof1
The MATLAB m file rdtrcsegy.m is:
function [seis,ntr,fmt,dsrt,dom,x,t0,nsamps,dt,nx,sh, ...
tr,rp,trp,trcid]=readsegy(filenm,a1,a2);
%
% Reads an SEG-Y trace into "ans"
% arguments:
% filename - The name of the SEG-Y file to be read
% a1 - The shot/rp number of the trace to be read
% a2 - The trace number to be read.
%
% Short program-driver for reading SEGY files
% The SEGY binary tape header parameters passed here are
% (all of them of 1x1 size, because it is a tape header,
% so you have to write it out just once):
% seis(nsamps by nx) Seismogram data
% ntr : Number of traces per shot (or record)
% fmt : SEGY format type, >4, "host" floating point
% dsrt : How is data sorted. 0 or 1 by shots; 2 by CDP gathers
% dom : The domain of the data (0 or 1 - time, 7 - tau-p, etc.)
% Variables are in the two temporary arrays -
% lbuf holds long integers and ibuf holds short integers
% The header parameters passed here are:
% x (nx by 1) : Range in km
% t0(nx by 1) : Start time in sec
%%% id(nx by 1) : Trace index integer (ie 1 through nx)
% nsamps(nx by 1) : Number of samples per trace
% dt(1 by 1) : Sample interval in sec/sample
% nx : Number of traces written in
% sh : Shot number
% tr : Trace within the shot
% rp : RP or CDP number
% trp : Trace number within RP or CDP
% trcfid : Trace ID; live = 1
fno = str2num(a1);
ftr = str2num(a2);
% Open file and read the "tape" header
fid = fopen(filenm, 'r');
if fid == -1
error('error opening SEG-Y file')
end
% Read the EBCDIC header, putting into temp arrays
% Useful information for subsequent writing is stored in appropriate variables
thbuf1 = fread(fid,3200)';
%disp( 'EBCDIC tape header ')
thbuf1(1,1);
%fseek(fid,3200,0);
thbuf2 = fread(fid,200,'short')';
ntr = thbuf2(7);
fmt = thbuf2(13);
dsrt = thbuf2(15);
dom = thbuf2(31);
tst = thbuf2(200);
if (fmt == 1)
exit('Cannot read IBM floating point format');
end
% Next read the trace headers and data, putting into temp arrays
% Reads each header twice, first putting into ibuf in short int format
% then into lbuf in long int format
% Array ibuf is (nx by 120) and array lbuf is (nx by 60)
buftr = 40; % Size of buffer increments
nx = 0; % Number of traces read
nread = 0;
doit = 1;
while (doit)
% T1 = fread(fid,120,'short')';
T1 = fread(fid,120,'int16')';
if (feof(fid) == 1) break; end
nread = nread + 1;
if (nread == 1 ) % Preassign Buffer space for trace headers/data
nsamps = T1(58);
bibuf = zeros(buftr,120);
blbuf = zeros(buftr, 60);
bseis = zeros(buftr,nsamps);
end
nsamps = T1(58);
% fseek(fid,-240,0)
fseek(fid,-240,'cof');
% T2 = fread(fid,60,'long')';
T2 = fread(fid,60,'int32')';
if (feof(fid) == 1)
fprintf(1,'reread of header failed.\n');
break;
end
if( fmt == 2 ) T3 = fread(fid,nsamps,'long')'; end;
if( fmt == 3 ) T3 = fread(fid,nsamps,'short')'; end;
if( fmt > 4 ) T3 = fread(fid,nsamps,'float')'; end;
if (feof(fid) == 1)
fprintf(1,'read of data failed.\n');
break;
end
nx = nx + 1;
bibuf(nx,:) = T1;
blbuf(nx,:) = T2;
if( blbuf(nx,7) == 0 )
no = blbuf(nx,3);
tr = blbuf(nx,4);
else
no = blbuf(nx,6);
tr = blbuf(nx,7);
end
bseis(nx,:) = T3;
if( tr == ftr )
doit = 0;
else
nx = nx - 1;
end
% Allocate additional storage space for data arrays
if ( nx > 0 )
if( rem(nx,buftr) == 0)
seis = [seis;bseis];
ibuf = [ibuf;bibuf];
lbuf = [lbuf;blbuf];
end
end
end
if ( nx == 0 )
quit
end
seis = bseis(1:nx,:);
ibuf = bibuf(1:nx,:);
lbuf = blbuf(1:nx,:);
fclose(fid);
% identify the variables inside the ibuf - the short integers
% or lbuf - the long integers, arrays to pass on
dt = bibuf(:,59) * 0.000001;
nsamps = bibuf(:,58);
x = blbuf(:,10) * 0.001;
t0 = bibuf(:,55) * 0.001;
sh = blbuf(:,3);
tr = blbuf(:,4);
rp = blbuf(:,6);
trp = blbuf(:,7);
trcid = bibuf(:,15);
%id = [:]'
seis = bseis(:,:)';
Return to EW0008 example.
Return to SIOSEIS examples.
Go to the list of seismic processes.
Go to SIOSEIS introduction.