读取32位IBM浮点数的Qt C++程序

SEGY是地震勘探存储数据的标准格式，最早在1975年发布revision 0。由于那时候地震数据的处理一般是在服务器上进行的，所以SEGY格式打上了时代印记，例如SEGY文件采用大字节序，数据点是IBM 32位浮点数。2002年发布的版本revision 1增加了IEEE 32位浮点数，但还是采用大字节序。

PC的处理器都是小字节序的，32位浮点数是IEEE浮点数，因此在PC上读取SEGY格式的文件时就比较麻烦，主要体现在：

• 读取SEGY文件里的整数或IEEE浮点数之后，必须进行字节序交换，转换为小字节序，才能正确地表示为PC上的数。

• 32位的IBM浮点数与32位IEEE浮点数格式不一样，不是简单的交换字节序的问题，需要经过计算才能将读取的IBM浮点数转换为PC上的浮点数。

本人在编写的一个数据处理软件中，需要读取其他采集设备转换出的标准SEGY文件，有的SEGY文件采用IBM浮点数，就涉及到将IBM浮点数转换为IEEE浮点数值。

关于32位IBM和IEEE浮点数结构及其转换方法可参考下面的一个文献：

符茂松《32位IEEE和IBM浮点数结构及其转换方法》，工程地球物理学报，2011年12月。

图A.  32位IBM浮点数结构

32位IBM浮点数的结构如图A所示，其具体解释见参考文献。

下面是读取一个4字节IBM浮点数，然后将其转换为IEEE浮点数的示意代码（只是从文件读取IBM浮点数、字节序交换、转换为IEEE浮点数数值的示意代码，不是读取一个文件内容的完整代码）。

   QFile aFile(aFileName);  //以文件方式读出
   if (!(aFile.open(QIODevice::ReadOnly)))
      return false;

    unsigned long  IBMFloatBytes; //即quint32
    aFile.read((char *)&IBMFloatBytes, 4);//用无符号4字节整数读出IBM 浮点数的4个字节
// 交换字节，变成如图A的字节序
    IBMFloatBytes=(((IBMFloatBytes>>24)&0xff) | ((IBMFloatBytes&0xff)<<24) |
                  ((IBMFloatBytes>>8)&0xff00) | ((IBMFloatBytes&0xff00)<<8));

   float  aSampValue=IBM2IEEE_float(IBMFloatBytes);//4字节IBM浮点数的实际数值

• 从文件中读取4字节IBM浮点数时，采用是读取一个4字节无符号整数的方法，这样32位数据原始不动的保存到变量IBMFloatBytes里。

• 变量IBMFloatBytes的4个字节实际是大字节序的IBM浮点数，将其进行字节序交换，变成小字节序，也就是图A所示的顺序。

• 调用函数IBM2IEEE_float()将变量IBMFloatBytes存储的小字节序的4个字节转换为IEEE浮点数。PC上C语言里的float就是4字节IEEE浮点数。
  

下面是函数IBM2IEEE_float( ) 的完整代码，它的功能是返回如图A存储的4个字节表示的浮点数的数值。

float QsgyFileReader::IBM2IEEE_float(const unsigned long IBMFloatBytes)
{ //IBM 浮点数转换为IEEE浮点数， IBMFloatBytes是已经转换为小字节序的 IEEE浮点数32位编码
   unsigned long ibmCode=IBMFloatBytes;

   float ieeeFloat=0;

//IBM浮点数： SEEEEEEE MMMMMMMM MMMMMMMM MMMMMMMM
//浮点数值 Value = (-1)^s * M * 16^(E-64)

   unsigned long   signCode=(ibmCode>>31);//获取符号位, S=0或1， sign=00 00 00 0000000S
   int sign=1; //正数
   if (signCode==1)
      sign=-1;//负数

   ibmCode= (ibmCode<<1); // 左移移出符号位, ibmCode= EEEEEEEM MMMMMMMM MMMMMMMM MMMMMMM0

   int exponent = (int)(ibmCode >> 25);  // 获取阶数, exponent=00 00 00 0EEEEEEE

   unsigned long fraction= (ibmCode << 7);  //移出符号位和阶数剩余的部分：尾数部分fraction=MMMMMMMM MMMMMMMM MMMMMMM 00000000
   fraction= (fraction >> 8);  //fraction=00000000 MMMMMMMM MMMMMMMM MMMMMMM

   static double ratio=qPow(2,24); //使用静态局部变量，避免重复计算
   if ((exponent==0) && (fraction==0))  //00000000 00 00 00 或10000000 00 00 00 都表示0
      ieeeFloat=0;
   else //IBMfloat= (sign)*Fraction*16^(exponent-64);
   {
      double   P2=qPow(16, exponent-64);
      double   P1=fraction/ratio;
      ieeeFloat=sign*P1*P2;
   }
    return   ieeeFloat;
}

IBM2IEEE_float( )函数的代码通过多次的移位操作，得到浮点数的符号、指数部分和尾数部分，然后计算出其表示的浮点数。代码不详细解释了，看注释即懂。

这里只实现了读取IBM浮点数，转换为的IEEE浮点数是什么样的存储结构由CPU自动处理。在导出标准格式SEGY文件时，使用IEEE浮点数格式，避免进行IEEE到IBM格式的转换。

图B是本人编写的数据处理软件中的数据波形显示，读取的数据文件是数据点格式为32位IBM浮点数的标准SEGY文件。

图B 读取数据点格式为IBM浮点数的标准SEGY文件并显示其波形