#!/usr/bin/perl



if ($#ARGV != 2) { 

   print STDERR "Usage: despike cutoff input.sgy output.sgy \n\n";

   print STDERR "If the cutoff is positive, values higher than the\n";

   print STDERR "cutoff will be clipped. If the cutoff is negative,\n";

   print STDERR "values lower than the cutoff will be clipped. The\n";

   print STDERR "output file cannot be the same as the input file.\n";

}



if ($ARGV[1] eq $ARGV[2]) { die "Output file can't be same as input file" }



$ebcdicLength = 3200; 

$binaryLength = 400; 

$traceHeaderLength = 240;



open (IN,  "<$ARGV[1]") || die "Could not open $ARGV[1] $! \n"; 

open (OUT, ">$ARGV[2]") || die "Could not open $ARGV[2] $! \n";



$cutoff = $ARGV[0];



# This next line makes $cutoff into a packed binary string 

# containing a single floating point number.



$cutoffp = pack("f", $cutoff);



# These lines just copy the EBCDIC header from the input file to the

# output file. Why are we using read() and syswrite() rather than <> 

# and print()? Because we might be reading from tape instead of disk,

# in which case the block structure must be preserved.



sysread  (IN,  $ebcdic, $ebcdicLength); 

syswrite (OUT, $ebcdic, $ebcdicLength); 

sysread  (IN,  $binary, $binaryLength); 

syswrite (OUT, $binary, $binaryLength);



# These two variables are set to the values of the short integers 

# at byte offsets 20 and 24 in the $binary packed string.



($numberSamples, $sampleFormat) = unpack("@20s@24s", $binary);



if ($sampleFormat == 3) { $traceLength = 2 * $numberSamples }

else                    { $traceLength = 4 * $numberSamples }

$traceBufferLength = $traceLength + $traceHeaderLength; 



while (!eof(IN)) {



    # We read with sysread() to get the exact number of bytes in a 

    # tape block. If the block is short, we either have a bad tape or

    # we're at the end of the tape. Either way, we're done. 

    $bytesRead = sysread(IN, $traceBuffer, $traceBufferLength); 

    last if $bytesRead < $traceBufferLength;



    $tracesRead++;



    # Announce every hundredth trace 

    if ($tracesRead % 100 == 0) { 

        print STDERR "Processed $tracesRead traces \n"; 

    }



    for $i (1 .. $numberSamples) { 

        $currentSample = $traceHeaderLength + ( ($i - 1) * 4);



        # Take the four bytes beginning at $currentSample; 

        $sample = substr($traceBuffer, $currentSample, 4);



        # ...and unpack() them as a single float 

        $amp = unpack("f", $sample);



        if ($cutoff > 0) { 

            if ($amp > $cutoff) {

                # If it's too high, replace it with the cutoff 

                substr($traceBuffer, $currentSample, 4) = $cutoffp; 

            } 

        } else { 

            if ($amp < $cutoff) {

                # If it's too low, replace it with the cutoff 

                substr($traceBuffer, $currentSample, 4) = $cutoffp; 

            } 

        } 

    } 

    syswrite(OUT, $traceBuffer, $traceBufferLength); 

}