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.