SEGD2ASCII and Header Dumper in Perl

Kode Perl di bawah ini berguna untuk menampilkan header dan mengkonversi data dengan format SEGD (32bits) ke ASCII.

Untuk memahami konfigurasi dari data SEGD silakan lihat di link ini.

Screendump di bawah ini menunjukkan aplikasi kode Perl tersebut.

Berikut kodenya:
#!/usr/bin/perl
$fname=$ARGV [0];
sysopen(IN,$fname,O_RDONLY) or die "Can't Open $fname: $! ";
open my $fh, '<', $fname or die "Can't open ";
read($fh,$cacing,928,0);

$size= -s IN;
$bytewidth=4;

#############################################################
#General header #1   32bytes, some BCDs some binaries
#############################################################
$F= unpack( "H*",substr($cacing,1-1,2));
print " file number: $F","\n ";

$Y= unpack( "H*",substr($cacing,3-1,2));



$K=unpack( "H*",substr($cacing,5-1,6));
print "general constant:$K","\n ";

$YR= unpack( "H*",substr($cacing,11-1,1));
print "year: $YR","\n ";

$gd= unpack( "B*",substr($cacing,12-1,2));
$GH=bin2dec(substr($gd,0,4));
print "revision: $GH","\n ";

for my $i (1..3)
{
$DY.= bin2dec(substr($gd,4*$i,4));
}

print "julian day: $DY","\n";


$H= unpack( "H*",substr($cacing,14-1,1));
print " hour of day: $H","\n ";

$MI= unpack( "H*",substr($cacing,15-1,1));
print "minute of hour: $MI","\n ";

$SE= unpack( "H*",substr($cacing,16-1,1));
print "second of minute: $SE","\n ";

$Mc= unpack( "H*",substr($cacing,17-1,1));
print "manufacturer's code: $Mc","\n ";

$Ms= unpack( "H*",substr($cacing,18-1,2));
print "manufacturer's serial no: Ms","\n ";

$B= unpack( "H*",substr($cacing,20-1,3));
print "byte per scan: $B","\n ";

$bsi= unpack( "B*",substr($cacing,23-1,1));
for my $i (0..7)
{
$sr1=substr($bsi,$i,1);
$I+= (8/(2**$i))*$sr1, "\n";
}
print "base scan interval: $I","\n ";

$P1= unpack( "B*",substr($cacing,24-1,1));
$P=substr($P1,0,4);

if($P==0000)
   {
   print "polarity: Untested","\n ";
   } elsif($P==0001)
   {
   print "polarity: Zero","\n ";
   } elsif($P==0010)
   {
   print "polarity: 45 Deg","\n ";
   } elsif($P==0011)
   {
   print "polarity: 90 Deg","\n ";
   } elsif($P==0100)
   {
   print "polarity: 135 Deg","\n ";
   } elsif($P==0101)
   {
   print "polarity: 180 Deg","\n ";
   } elsif($P==0110)
   {
   print "polarity: 225 Deg","\n ";
   } elsif($P==0111)
   {
   print "polarity: 270 Deg","\n ";
   } elsif($P==1000)
   {
   print "polarity: 315 Deg","\n ";
   } elsif($P==1100)
   {
   print "polarity: unassigned","\n ";
   } else
   {
   print "polarity: $P","\n ";
   }

$SBX=bin2dec(substr($P1,4,4));
$SB= bin2dec(unpack( "B*",substr($cacing,25-1,1)));
$SB1=$SB*(2**$SBX);
print "S/B exponent: $SB1","\n ";
$Z= substr((unpack( "B*",substr($cacing,26-1,1))),0,4);

if($Z==0010)
   {
   print "record type: Test record","\n ";
   } elsif($Z==0100)
   {
   print "record type: Parallel channel test","\n ";
   } elsif($Z==0110)
   {
   print "record type: Direct channel test","\n ";
   } elsif($Z==1001)
   {
   print "record type: Normal record","\n ";
   } elsif($Z==0001)
   {
   print "record type: Other","\n ";
   } else
   {
   print "record type: $Z","\n ";
   }


#############################################################
#General header #2   32bytes, some BCDs some binaries
#############################################################

$ECX= bin2dec(unpack( "B*",substr($cacing,32+6-1,2)));
print "extended header blocks: $ECX","\n ";
$EH= bin2dec(unpack( "B*",substr($cacing,32+8-1,2)));
print "external header blocks: $EH","\n ";

$generalheader=32+32+32*$ECX;


#############################################################
#Channel set descriptor   32bytes, some BCDs some binaries
#############################################################


$TF= bin2dec(unpack( "B*",substr($cacing,$generalheader+3-1,2)));
print "channel set starting time: $TF","\n ";

$TE= 2*(bin2dec(unpack( "B*",substr($cacing,$generalheader+5-1,2))));
print "channel set end time: $TE","\n ";

$ECS= bin2dec(unpack( "B*",substr($cacing,$generalheader+27-1,2)));
print "extended channel set number: $ECS","\n ";

$headandchan=(32+32+32*$ECX)+($ECS*32)+($EH*32);

$samp=$TE/2;
$sr=$TE/$samp;



#############################################################
# Demux trace header  20bytes, some BCDs some binaries
#############################################################

$THE=bin2dec(unpack( "B*",substr($cacing,$headandchan+10-1,1)));
print "trace header extension: $THE","\n ";
$trtot=20+(32*$THE);
print "************************************************","\n ";
print "Key Information                                  ","\n ";
print "SEGD Revision: $GH","\n ";

if($Y=='8058')
   {
   print "format code: $Y >>> 32 bit IEEE demultiplexed","\n ";
   } elsif($Y=='0058')
   {
   print "format code: $Y >>> 32 bit IEEE multiplexed","\n ";
   } else
   {
   print "format code: $Y >>> Check SEGD Convention ","\n ";
   }

print "no of sample each trace: $samp","\n ";
print "sampling rate: $sr ms","\n ";
print "total header general, channel, external (byte): $headandchan","\n ";
print "trace header length (byte): $trtot","\n ";
$tracecount=($size-$headandchan)/($trtot+($bytewidth*$samp));
$tracecount = sprintf("%.0f", $tracecount);
print "no of traces : $tracecount\n";
print "************************************************","\n ";

#############################################################
# Reading data >>>> traces
#############################################################


print " Type 1 = IBM 32-Bits floating point or 2 = IEEE 32 bits floating point\n";
$formatcode= <STDIN>;
if($formatcode==1)
   {
   print " Converting $fname to ASCII with IBM 32-Bits floating point format...\n";
   } elsif($formatcode==2)
   {
   print " Converting $fname to ASCII with IEEE 32-Bits floating point format...\n";
   } else
   {
   print " Type 1 = IBM 32-Bits floating point or 2 = IEEE 32 bits floating point\n";
   exit;
   }

seek (IN,$headandchan+$trtot,1);
my $buffer;
my $i=1;
open OUT," > out.txt" or die "$!\n";
while (sysread(IN, $buffer, $trtot+($samp*$bytewidth),0))
{
for ($j=0; $j<=$samp-1; $j++)
{
$amp[$j]=unpack("B*",substr($buffer,$j*4));

if($formatcode==1)
   {
   $trace_value[$j] = conv_bit_string_2_ibm32float($amp[$j]);
   } elsif($formatcode==2)
   {
   $trace_value[$j] = conv_bit_string_2_ieee32float($amp[$j]);
   }
print OUT $trace_value[$j],"\t";
}
print OUT "\n";
progress_bar( $i, $tracecount, 25, '=' );
$i++;
}
print "\e[?25h";
print "\n";
print "...process done!:   out.txt (rows(traces)=$tracecount cols(samples)=$samp) has been created\n";
sub progress_bar {
    my ( $got, $total, $width, $char ) = @_;
    $width ||= 25;
    $char  ||= '=';
    my $num_width = length $total;
    local $| = 1;
    print "\e[?25l";
    printf "|%-${width}s| write %${num_width}s traces of %s (%.2f%%)\r",
        $char x (($width-1)*$got/$total). '>', $got, $total, 100*$got/$total;
}
sub bin2dec {
return unpack("N", pack("B32", substr("0" x 32 . shift, -32)));
}

sub conv_bit_string_2_ibm32float {
$bit_string_ibm = shift;
$first_digit_ibm = substr($bit_string_ibm, 0, 1);
if ( $first_digit_ibm eq "0" ) {
$sign_bit_ibm = 1;
} elsif ( $first_digit_ibm eq "1" ) {
$sign_bit_ibm = -1;
}
$bin_exponent_ibm = substr($bit_string_ibm, 1, 7);

$exponent_ibm = bin2dec($bin_exponent_ibm);
$bin_fraction_ibm = substr($bit_string_ibm, 8, 31);
@bit_chars_ibm = unpack("A1" x length($bin_fraction_ibm), $bin_fraction_ibm);
$place_holder_ibm = -1;
$fraction_ibm = 0;
foreach $bit_ibm ( @bit_chars_ibm ) {
$fraction_ibm += $bit_ibm * (2 ** $place_holder_ibm);
$place_holder_ibm += -1;
}
$ibm_float_ibm = ($sign_bit_ibm ** 1) * (16 ** ($exponent_ibm - 64)) *($fraction_ibm);
return sprintf("%.4f", $ibm_float_ibm);
}

sub conv_bit_string_2_ieee32float {
$bit_string_ieee = shift;
$first_digit_ieee = substr($bit_string_ieee, 0, 1);
if ( $first_digit_ieee eq "0" ) {
$sign_bit_ieee = 1;
} elsif ( $first_digit_ieee eq "1" ) {
$sign_bit_ieee = -1;
}
$bin_exponent_ieee = substr($bit_string_ieee, 1, 8);
$exponent_ieee = bin2dec($bin_exponent_ieee);
$bin_fraction_ieee = substr($bit_string_ieee, 9, 31);
@bit_chars_ieee = unpack("A1" x length($bin_fraction_ieee), $bin_fraction_ieee);
$place_holder_ieee = -1;
$fraction_ieee = 0;
foreach $bit_ieee ( @bit_chars_ieee ) {
$fraction_ieee += $bit_ieee * (2 ** $place_holder_ieee);
$place_holder_ieee += -1;
}
$ieee_float_ieee = ($sign_bit_ieee ** 1) * (1+$fraction_ieee)*(2**($exponent_ieee - 127));
return sprintf("%.4f", $ieee_float_ieee);
}

exit;



GNUPLOT untuk transpose dan menampilkan hasil:

#!/usr/bin/perl
my $file = "out.txt";
open( FILE, $file ) or die "Can't open file '$file': $!";
while( <FILE> ) {
chomp;
my @row = split;
push @data, \@row; 
}
close( FILE );


open OUT1," > out1.txt" or die "$!\n";
for my $j (0..$#{$data[0]})
{
print OUT1 "$j","\t";
for my $i (0..$#data)
{
print OUT1 "$data[$i][$j]","\t";
}
print  OUT1 "\n";
}



########## PLOT#############################################
open my $GP, '|-', 'gnuplot -persist';
print {$GP}  <<'__GNUPLOT__';

set yrange[] reverse;
plot for [i=2:5] 'out1.txt' using i:1 with lines title ""


__GNUPLOT__
close $GP;