biosig.c

c++编写的用于生物信号处理的软件库

	/*
    $Id: biosig.c,v 1.139 2008/03/14 08:30:41 schloegl Exp $
    Copyright (C) 2005,2006,2007,2008 Alois Schloegl <[e-mail]>
    This file is part of the "BioSig for C/C++" repository 
    (biosig4c++) at http://biosig.sf.net/ 


    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License
    as published by the Free Software Foundation; either version 3
    of the License, or (at your option) any later version.

 */

/*

	reading and writing of GDF files is demonstrated 
	
	Features: 
	- reading and writing of EDF, BDF, GDF1, GDF2, CWFB, HL7aECG, SCP files 
	- reading of ACQ, AINF, BKR, BrainVision, CNT, DEMG, EGI, MFER files 
	
	implemented functions: 
	- SOPEN, SREAD, SWRITE, SCLOSE, SEOF, SSEEK, STELL, SREWIND 
	
*/

#include <ctype.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

//#include <libxml/xmlreader.h>

#include "biosig-dev.h"

const int16_t GDFTYP_BYTE[] = {
	1, 1, 1, 2, 2, 4, 4, 8, 8, 4, 8, 0, 0, 0, 0, 0,   /* 0  */ 
	4, 8,16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   /* 16 */ 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   /* 32 */ 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   /* 48 */ 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   /* 64 */ 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    /* 128 */ 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2,    /* 256 - 271*/  
	0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 4,    /* 255+24 = bit24, 3 byte */ 
	0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 6, 
	0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 8, 
	0, 0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0, 0,10, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,    /* 384 - 399*/ 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
	0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2,    /* 512 - 527*/ 
	0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 4, 
	0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 6, 
	0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 8, 
	0, 0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0, 0,10};


const char *LEAD_ID_TABLE[] = { "unspecified",
	"I","II","V1","V2","V3","V4","V5","V6",
	"V7","V2R","V1","V2","V3","V4","V5","V6",
	"V7","V2R","V3R","V4R","V5R","V6R","V7R","X",
	"Y","Z","CC5","CM5","LA","RA","LL","fI",
	"fE","fC","fA","fM","fF","fH","dI",
	"dII","dV1","dV2","dV3","dV4","dV5",
	"dV6","dV7","dV2R","dV3R","dV4R","dV5R",
	"dV6R","dV7R","dX","dY","dZ","dCC5","dCM5",
	"dLA","dRA","dLL","dfI","dfE","dfC","dfA",
	"dfM","dfF","dfH","III","aVR","aVL","aVF",
	"aVRneg","V8","V9","V8R","V9R","D","A","J",
	"Defib","Extern","A1","A2","A3","A4","dV8",
	"dV9","dV8R","dV9R","dD","dA","dJ","Chest",
	"V","VR","VL","VF","MCL","MCL1","MCL2","MCL3",
	"MCL4","MCL5","MCL6","CC","CC1","CC2","CC3",
	"CC4","CC6","CC7","CM","CM1","CM2","CM3","CM4",
	"CM6","dIII","daVR","daVL","daVF","daVRneg","dChest",
	"dV","dVR","dVL","dVF","CM7","CH5","CS5","CB5","CR5",
	"ML","AB1","AB2","AB3","AB4","ES","AS","AI","S",
	"dDefib","dExtern","dA1","dA2","dA3","dA4","dMCL1",
	"dMCL2","dMCL3","dMCL4","dMCL5","dMCL6","RL","CV5RL",
	"CV6LL","CV6LU","V10","dMCL","dCC","dCC1","dCC2",
	"dCC3","dCC4","dCC6","dCC7","dCM","dCM1","dCM2",
	"dCM3","dCM4","dCM6","dCM7","dCH5","dCS5","dCB5",
	"dCR5","dML","dAB1","dAB2","dAB3","dAB4","dES",
	"dAS","dAI","dS","dRL","dCV5RL","dCV6LL","dCV6LU","dV10"
/*   EEG Leads - non consecutive index 
	,"NZ","FPZ","AFZ","FZ","FCZ","CZ","CPZ","PZ",
	"POZ","OZ","IZ","FP1","FP2","F1","F2","F3","F4",
	"F5","F6","F7","F8","F9","F10","FC1","FC2","FC3",
	"FC4","FC5","FC6","FT7","FT8","FT9","FT10","C1",
	"C2","C3","C4","C5","C6","CP1","CP2","CP3","CP4",
	"CP5","CP6","P1","P2","P3","P4","P5","P6","P9",
	"P10","O1","O2","AF3","AF4","AF7","AF8","PO3",
	"PO4","PO7","PO8","T3","T7","T4","T8","T5","P7",
	"T6","P8","T9","T10","TP7","TP8","TP9","TP10",
	"A1","A2","T1","T2","PG1","PG2","SP1","SP2",
	"E0","EL1","EL2","EL3","EL4","EL5","EL6","EL7",
	"ER1","ER2","ER3","ER4","ER5","ER6","ER7","ELL",
	"ERL","ELA","ELB","ERA","ERB"
*/
	, "\0\0" };  // stop marker 



/****************************************************************************/
/**                                                                        **/
/**                      INTERNAL FUNCTIONS                                **/
/**                                                                        **/
/****************************************************************************/
 
// greatest common divisor 
uint32_t gcd(uint32_t A, uint32_t B)
{	size_t t; 
	if (A<B) {t=B; B=A; A=t;}; 
	while (B) {
		t = B; 
		B = A%B;
		A = t; 
	}
	return(A);
};

// least common multiple - used for obtaining the common HDR.SPR
uint32_t lcm(uint32_t A, uint32_t B)
{	
	// return(A*(B/gcd(A,B)) with overflow detection 
	uint64_t A64 = A;
	A64 *= B/gcd(A,B);
	if (A64 > 0x00000000ffffffffllu) {
		fprintf(stderr,"Error: HDR.SPR=LCM(%i,%i) overflows and does not fit into uint32.\n",A,B);
		exit(-4);
	}
	return((uint32_t)A64);
};

#if __BYTE_ORDER == __BIG_ENDIAN
float l_endian_f32(float x) 
#elif __BYTE_ORDER == __LITTLE_ENDIAN
float b_endian_f32(float x) 
#endif
{
	union {
		float f32;
		uint32_t u32;
	} b1,b2; 
	b1.f32 = x; 
	b2 = b1; 
	b2.u32 = bswap_32(b1.u32);
	return(b2.f32);
}

#if __BYTE_ORDER == __BIG_ENDIAN
double l_endian_f64(double x) 
#elif __BYTE_ORDER == __LITTLE_ENDIAN
double b_endian_f64(double x) 
#endif
{
	union {
		double f64;
		uint32_t u32[2];
	} b1,b2; 
	b1.f64 = x; 
	b2 = b1; 
	b2.u32[0] = bswap_32(b1.u32[1]);
	b2.u32[1] = bswap_32(b1.u32[0]);
	return(b2.f64);
}

#ifdef __sparc__
/*    SPARC: missing alignment must be explicitly handled     */ 
uint16_t leu16p(uint8_t* i) {
	// decode little endian uint16 pointer 
	return ((*i) + (*(i+1) << 8));
}
int16_t lei16p(uint8_t* i) {
	// decode little endian int16 pointer 
	uint16_t o = ((*i) + (*(i+1) << 8));
	return(*(int16_t*)(&o)); 
}
uint32_t leu32p(uint8_t* i) {
	// decode little endian uint32 pointer 
	uint32_t o=0;
	for (char k=0; k<4; k++)
		o += (*(i+k))<<(k*8);
	return(o); 
}
int32_t lei32p(uint8_t* i) {
	// decode little endian int32 pointer 
	uint32_t o=0;
	for (char k=0; k<4; k++)
		o += (*(i+k))<<(k*8);
	return(*(int32_t*)(&o)); 
}
uint64_t leu64p(uint8_t* i) {
	// decode little endian uint64 pointer 
	uint64_t o=0;
	for (char k=0; k<8; k++)
		o += (*(i+k))<<(k*8);
	return(o); 
}
int64_t lei64p(uint8_t* i) {
	// decode little endian int64 pointer 
	int64_t o=0;
	for (char k=0; k<8; k++)
		o += (*(i+k))<<(k*8);
	return(*(int64_t*)(&o)); 
}
float lef32p(uint8_t* i) {
	// decode little endian float pointer 
	uint32_t o;
	for (char k=0, o=0; k<4; k++)
		o += (*(i+k))<<(k*8);
	return(*(float*)(&o)); 
}
double lef64p(uint8_t* i) {
	// decode little endian double pointer 
	uint64_t o=0;
	for (char k=0; k<8; k++)
		o += (*(i+k))<<(k*8);
	return(*(double*)(&o)); 
}

uint16_t beu16p(uint8_t* i) {
	// decode big endian uint16 pointer 
	return ((*i<<8) + (*(i+1)));
}
int16_t bei16p(uint8_t* i) {
	// decode big endian int16 pointer 
	uint16_t o = ((*i << 8) + (*(i+1)));
	return(*(int16_t*)(&o)); 
}
uint32_t beu32p(uint8_t* i) {
	// decode big endian uint32 pointer 
	uint32_t o=0;
	for (char k=0; k<4; k++) {
		o<<=8;
		o += *(i+k);
	}	
	return(o); 
}
int32_t bei32p(uint8_t* i) {
	// decode big endian int32 pointer 
	uint32_t o=0;
	for (char k=0; k<4; k++) {
		o<<=8;
		o += *(i+k);
	}	
	return(*(int32_t*)(&o)); 
}
uint64_t beu64p(uint8_t* i) {
	// decode big endian uint64 pointer 
	uint64_t o=0;
	for (char k=0; k<8; k++) {
		o<<=8;
		o += *(i+k);
	}	
	return(o); 
}
int64_t bei64p(uint8_t* i) {
	// decode big endian int64 pointer 
	uint64_t o=0;
	for (char k=0; k<8; k++){
		o<<=8;
		o += *(i+k);
	}	
	return(*(int64_t*)(&o)); 
}
float bef32p(uint8_t* i) {
	// decode big endian float pointer 
	uint32_t o=0;
	for (char k=0; k<4; k++){
		o<<=8;
		o += *(i+k);
	}	
	return(*(float*)(&o)); 
}
double bef64p(uint8_t* i) {
	// decode big endian double pointer 
	uint64_t o=0;
	for (char k=0; k<8; k++){
		o<<=8;
		o += *(i+k);
	}	
	return(*(double*)(&o)); 
}

void leu16a(uint16_t i, uint8_t* r) {
	i = l_endian_u16(i);
	memcpy(&i,r,sizeof(i));
}
void lei16a( int16_t i, uint8_t* r) {
	i = l_endian_i16(i);
	memcpy(&i,r,sizeof(i));
}
void leu32a(uint32_t i, uint8_t* r) {
	i = l_endian_u32(i);
	memcpy(&i,r,sizeof(i));
}
void lei32a( int32_t i, uint8_t* r) {
	i = l_endian_i32(i);
	memcpy(&i,r,sizeof(i));