How to read 32-bit *IBM* (not IEEE) floating point values from a binary file?

Question

I am trying to read a binary file. This file records seismic data. Basically its first 3600 bytes are for the general file information, following are for seismic traces, each begins with 240 bytes for information about this trace.

For the file I am reading, it has 3617 traces, each trace has 6240 bytes, with the first 240 bytes for trace headers. Following the trace headers are some velocity values for underground formation, it has 1500 data points. I can view this data using some commercial software... the values should be around 4999.99 feet/s(velocity of sea water) to 14800 feet/s. But what I read range from 590.125 to 743.25, which is obviously not correct.

Later I found out that the binary file I am trying to read stores IBM-format floating point values.

If I read it as an IEEE binary file, for integer data, it is all right, but for floating point data, what I read is totally wrong.

The trace head I read is about the geographic location of the traces, and it is correctly read by Mathematica, giving {0, 20, 40, 60, ..., 72320} feet. I don't need to read all the trace headers.

What I care about most is reading the correct velocity values.

Here's the code I have so far (the input files are from timodel_vp.segy.gz (22 MB) from here: http://software.seg.org/datasets/2D/Hess_VTI/)

nByteFile = 
  FileByteCount[
   "timodel_c11.segy"];

nByteTrace = (nByteFile - 3600)/3617;

nTraceSample = (nByteTrace - 240)/4;


str = OpenRead[
   "timodel_vp.segy", 
   BinaryFormat -> True];
pos = Range[3600, nByteFile - 1, nByteTrace];

traces = {SetStreamPosition[str, # + 76]; 
  BinaryRead[str, "Integer32", ByteOrdering -> +1], 
  SetStreamPosition[str, # + 240]; 
  Table[BinaryRead[str, "Real32", 
    ByteOrdering -> +1], {nTraceSample}]} & /@ pos; 

Answer 1

I imported the data this way:

data = Import["ftp://software.seg.org/pub/datasets/2D/Hess_VTI/timodel_vp.segy.gz", 
              "Binary"];

nByteFile = Length@data
nByteTrace = (nByteFile - 3600)/3617
nTraceSample = (nByteTrace - 240)/4

22573680
6240
1500

You have to convert bytes to IBM 32-bit floating-point yourself. This takes a list of bytes, divides into 4-byte groups and converts them.

ibmFloat32[bytes_List] := 
  With[{sign = BitGet[#[[1]], 7],
        exponent = BitClear[#[[1]], 7], 
        significand = FromDigits[#[[2 ;; 4]], 2^8]},
     N @ (-1)^sign * significand/2^24 * 16^(exponent - 64)] & /@ 
   Partition[bytes, 4];

Then we can find the traces. It takes about a minute. (The code could be made more efficient).

pos = Range[3600, nByteFile - 1, nByteTrace];

traces = {FromDigits[data[[# + 76 + 1 ;; # + 80]], 2^8],
          ibmFloat32[data[[# + 240 + 1 ;; # + nByteTrace]]]} & /@ pos;

Example

traces[[2]] // Short

{20, {4999.99, <<1498>>, 11878.}}

This is in the range specified, so I hope it's working.

---

This is almost 100 times faster, but less readable:

newibmFloat32 = Compile[{{bytes, _Integer, 1}},
   (2. UnitStep[127.5 - #[[1]]] - 1) *
   16^(#[[1]] - UnitStep[#[[1]] - 127.5] 128. - 64.) *
   {2.^-8, 2.^-16, 2.^-24}.#[[2 ;; 4]] &@
      Transpose@Partition[N@bytes, 4]
   ];

traces2 = {FromDigits[data[[# + 76 + 1 ;; # + 80]], 2^8],
           newibmFloat32[data[[# + 240 + 1 ;; # + nByteTrace]]]} & /@ pos;

traces2 == traces

True