用C语言读把SGY地震数据读成txt

sgy地震数据格式：
文件头总长度为3600字节，分两部分。
1.文件头第一部分
长度： 3200bytes ；组成： 80bytes*40；特性：EBCDIC字符集，参数卡，需要转换为ASCII码后才能显示。
2.文件头第二部分
长度：400bytes；数据类型：32位、16位的整型；特性：二进制头，记录数据体信息。
3.道头 240bytes。
读取：
因为我这里只需要数据，所以直接读头文件有用信息，然后跳过头文件和道头字节，把数据写入txt。

#include "stdio.h" 
#include "stdlib.h" 
#include "string.h" 
#include "math.h" 

float ibm2float(int x);
int swap4byte(int value);

void main()
{
	FILE *filein,*fileout;

	char *inputfile,*outputfile;  //输入输出文件指针
	inputfile = (char*)malloc(sizeof(char)*200);
	outputfile = (char*)malloc(sizeof(char)*200);

	int trace_length = 0;	//道长度（即每道的时间采样点数）
	long int i,len;  //len表示文件的总字节数
	int traces = 0;  //道数 
	unsigned char f3200[3200];
	short int f400[200];

	inputfile = "";  //指定文件，这里需要改成读入绝对路径
	outputfile = "";

	filein = fopen(inputfile,"rb");  //打开文件
	if(filein == NULL)   
	{   
		printf("Can not open File %s\n",inputfile);
		exit(0);
	}
	fileout = fopen(outputfile,"w");
	if(fileout == NULL)   
	{   
		printf("Can not open File %s \n",outputfile);
		exit(0);   
	}

	fread(f3200,3200,1,filein);  //读入3200字节文件头部分  
	fread(f400,200,2,filein);    //读入400字节文件头部分 

	short svalue;
	svalue=((f400[10] & 0xff00)>> 8)|((f400[10] & 0x00ff)<< 8);  //读道长度
	trace_length=svalue; 
	
	fseek(filein,0,2);  
	len=ftell(filein);	//文件长度
	traces=(len-3600)/(240+4L*trace_length);  //计算道数
	printf("len = %d\t traces = %d,trace_length = %d\n",len,traces,trace_length);

	fseek(filein,3600,SEEK_SET);
	char *header;
	header = (char*)malloc(sizeof(char)*240);
	int *data;
	data = (int*)malloc(sizeof(int)*trace_length);
	float *d1;
	d1 = (float*)malloc(sizeof(float)*trace_length);
	int sizeofTraceData = 4L *trace_length ;
	int sizeofTraceHeader = 240;
	int qw1=0;

	for(i=0;i<traces;i++)
	{

		if(i%10000==0)
		{
			printf("%d\n",i);
		}
		fread(header,sizeofTraceHeader,1,filein);  //读道头和数据体
		fread(data,sizeofTraceData,1,filein);
		int j;
		for(j = 0;j < trace_length;j++)
		{
			qw1 = swap4byte(data[j]);
			d1[j] = ibm2float(qw1) ;
			fprintf(fileout,"%f\n",d1[j]);  //把浮点型数据写入文件
		}
	}
	free(data);  //释放内存
	free(header);
	free(d1);
	fclose(filein);  //关闭文件
	fclose(fileout);
}

float ibm2float(int x)
{
	int sign,expo,mantissa,res;
	float result,mant;
	if(x){
		sign = x>>31 & 0x01; //计算IBM浮点数的符号s

		expo = x>>24 & 0x7f; //计算IBM浮点数的指数e

		mantissa = x & 0x00ffffff; //计算IBM浮点数的小数部分f

		if(expo == 127 &&(mantissa & 0x00800000)&& !(mantissa & 0x007fffff)){//判断非数字

			res=(sign<<31)| 0x7f800001;

			memcpy(&result,&res,sizeof(float));

			return result;

		}

		if(expo > 96){    //判断无穷大值

			res =(sign<<31)| 0x7f800000;

			memcpy(&result,&res,sizeof(float));

		}

		else if(expo < 32)result = 0; //判断0值

		else {

			expo = expo-64; //得到指数偏移后的指数e-64

			mant =(float)(mantissa)/ pow(2,24); //转化为十进制小数0.f

			result =(1-2*sign)* mant * pow(16,expo); //带入(2)式中得到十进制数值

		}

		return result; //返回结果
	}
	result = x; //当x=0时,返回0值
	return result;
}

int swap4byte(int value)
{
	int svalue;
	svalue =((value & 0xff000000)>> 24)|((value & 0x00ff0000)>> 8)|

		((value & 0x0000ff00)<< 8)|((value & 0x000000ff)<< 24);

	return svalue;
}