sopen_scp_read.c

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

			en1064.Section4.len_ms	= leu16p(PtrCurSect+curSectPos);
			en1064.Section4.fiducial_sample	= leu16p(PtrCurSect+curSectPos+2);
			en1064.Section4.N	= leu16p(PtrCurSect+curSectPos+4);
			en1064.Section4.SPR	= hdr->SPR/4;

			en1064.Section4.beat	= (typeof(en1064.Section4.beat))malloc(en1064.Section4.N*sizeof(*en1064.Section4.beat));

			curSectPos += 6;
			for (i=0; i < en1064.Section4.N; i++) {
				en1064.Section4.beat[i].btyp = leu16p(PtrCurSect+curSectPos);
				en1064.Section4.beat[i].SB   = leu32p(PtrCurSect+curSectPos+2);
				en1064.Section4.beat[i].fcM  = leu32p(PtrCurSect+curSectPos+6);
				en1064.Section4.beat[i].SE   = leu32p(PtrCurSect+curSectPos+10);
				curSectPos += 14;
			}	
			for (i=0; i < en1064.Section4.N; i++) {
				en1064.Section4.beat[i].QB   = leu32p(PtrCurSect+curSectPos);
				en1064.Section4.beat[i].QE   = leu32p(PtrCurSect+curSectPos+4);
				curSectPos += 8;
				en1064.Section4.SPR += en1064.Section4.beat[i].QE-en1064.Section4.beat[i].QB-1;
			}	
			if (en1064.Section4.len_ms==0) {
				en1064.FLAG.REF_BEAT = 0;
			}	
		}

		/**** SECTION 5 ****/
		else if (curSect==5)  {
			Cal5 			= leu16p(PtrCurSect+curSectPos);
			en1064.Section5.AVM	= leu16p(PtrCurSect+curSectPos);
			en1064.Section5.dT_us	= leu16p(PtrCurSect+curSectPos+2);
			en1064.Section5.DIFF 	= *(PtrCurSect+curSectPos+4);
			en1064.Section5.Length  = (1000L * en1064.Section4.len_ms) / en1064.Section5.dT_us; // hdr->SPR;
			en1064.Section5.inlen	= (typeof(en1064.Section5.inlen))malloc(hdr->NS*sizeof(*en1064.Section5.inlen));
			for (i=0; i < hdr->NS; i++) {
				en1064.Section5.inlen[i] = leu16p(PtrCurSect+curSectPos+6+2*i);
				if (!section[4].length && (en1064.Section5.Length<en1064.Section5.inlen[i]))
					en1064.Section5.Length = en1064.Section5.inlen[i];
			}
			if (!section[4].length && en1064.FLAG.HUFFMAN) {
				 en1064.Section5.Length *= 5; // decompressed data might need more space 
				 fprintf(stderr,"Warning SCPOPEN: Section 4 not defined - size of Sec5 can be only guessed (%i allocated)\n",en1064.Section5.Length);
			}

			if (en1064.FLAG.REF_BEAT) {
				en1064.Section5.datablock = (int32_t*)malloc(4 * hdr->NS * en1064.Section5.Length); 

				Ptr2datablock           = (PtrCurSect+curSectPos+6+2*hdr->NS);
				for (i=0; i < hdr->NS; i++) {
					en1064.Section5.inlen[i] = leu16p(PtrCurSect+curSectPos+6+2*i);
					if (en1064.FLAG.HUFFMAN) {
						if (B4C_ERRNUM) {
							deallocEN1064(en1064);
							return(-1);
						}
						DecodeHuffman(HTrees, Huffman, Ptr2datablock, en1064.Section5.inlen[i], en1064.Section5.datablock + en1064.Section5.Length*i, en1064.Section5.Length);
					}
					else {
						for (k1=0; k1<en1064.Section5.Length; k1++)
							en1064.Section5.datablock[i*en1064.Section5.Length+k1] = lei16p(Ptr2datablock + 2*(i*en1064.Section5.Length + k1));
					}
					Ptr2datablock += en1064.Section5.inlen[i];
				}	
				size_t ix;
				data = en1064.Section5.datablock;
				if (en1064.Section5.DIFF==1)
					for (k1 = 0; k1 < hdr->NS; k1++)
					for (ix = k1*en1064.Section5.Length+1; ix < (k1+1)*en1064.Section5.Length; ix++)
						data[ix] += data[ix-1];

				else if (en1064.Section5.DIFF==2)
					for (k1 = 0; k1 < hdr->NS; k1++)
					for (ix = k1*en1064.Section5.Length+2; ix < (k1+1)*en1064.Section5.Length; ix++)
						data[ix] += 2*data[ix-1] - data[ix-2];
			}
		}

		/**** SECTION 6 ****/
		else if (curSect==6)  {
			hdr->NRec = 1;

			en1064.Section6.inlen	= (typeof(en1064.Section6.inlen))malloc(hdr->NS*sizeof(*en1064.Section6.inlen));
			uint16_t gdftyp 	= 5;	// int32: internal raw data type   
			hdr->AS.rawdata = (uint8_t*)malloc(4 * hdr->NS * hdr->SPR * hdr->NRec); 
			data = (int32_t*)hdr->AS.rawdata;

			en1064.Section6.AVM	= leu16p(PtrCurSect+curSectPos);
			en1064.Section6.dT_us	= leu16p(PtrCurSect+curSectPos+2);
			hdr->SampleRate 	= 1e6/en1064.Section6.dT_us;
			en1064.Section6.DIFF 	= *(PtrCurSect+curSectPos+4);
			en1064.FLAG.DIFF 	= *(PtrCurSect+curSectPos+4);
			en1064.Section6.BIMODAL = *(PtrCurSect+curSectPos+5);
			en1064.FLAG.BIMODAL     = *(PtrCurSect+curSectPos+5);

			Cal6 			= leu16p(PtrCurSect+curSectPos);
			en1064.Section6.dT_us	= leu16p(PtrCurSect+curSectPos+2);
			aECG->FLAG.DIFF 	= *(PtrCurSect+curSectPos+4);
			aECG->FLAG.BIMODAL = *(PtrCurSect+curSectPos+5);

			if ((section[5].length>4) &&  en1064.Section5.dT_us) 
				dT_us = en1064.Section5.dT_us;
			else 	
				dT_us = en1064.Section6.dT_us;
			hdr->SampleRate 	= 1e6/dT_us; 

			typeof(hdr->SPR) SPR  = ( en1064.FLAG.BIMODAL ? en1064.Section4.SPR : hdr->SPR);  

			if      (Cal5==0 && Cal6 >0) Cal0 = Cal6;
			else if (Cal5 >0 && Cal6==0) Cal0 = Cal5;
			else if (Cal5 >0 && Cal6 >0) Cal0 = gcd(Cal5,Cal6);
			uint16_t cal5 = Cal5/Cal0; 
			uint16_t cal6 = Cal6/Cal0; 

			Ptr2datablock = (PtrCurSect+curSectPos + 6 + hdr->NS*2);   // pointer for huffman decoder
			len = 0;  
			size_t ix; 			
			for (i=0; i < hdr->NS; i++) {
				hdr->CHANNEL[i].SPR 	    = hdr->SPR;
				hdr->CHANNEL[i].PhysDimCode = 4275; // PhysDimCode("uV") physical unit "uV"
				hdr->CHANNEL[i].Cal 	    = Cal0 * 1e-3;
				hdr->CHANNEL[i].Off         = 0;
				hdr->CHANNEL[i].OnOff       = 1;    // 1: ON 0:OFF
				hdr->CHANNEL[i].Transducer[0] = '\0';
				hdr->CHANNEL[i].GDFTYP      = gdftyp;  

				// ### these values should represent the true saturation values ### //
				hdr->CHANNEL[i].DigMax      = ldexp(1.0,20)-1;
				hdr->CHANNEL[i].DigMin      = ldexp(-1.0,20);
				hdr->CHANNEL[i].PhysMax     = hdr->CHANNEL[i].DigMax * hdr->CHANNEL[i].Cal;
				hdr->CHANNEL[i].PhysMin     = hdr->CHANNEL[i].DigMin * hdr->CHANNEL[i].Cal;
				
#ifndef WITHOUT_SCP_DECODE
//				if (AS_DECODE > 0) continue; 
#endif

				en1064.Section6.inlen[i]    = leu16p(PtrCurSect+curSectPos+6+2*i);
				if (en1064.FLAG.HUFFMAN) {
					if (B4C_ERRNUM) {
						deallocEN1064(en1064);
						return(-1);
					}
					DecodeHuffman(HTrees, Huffman, Ptr2datablock, en1064.Section6.inlen[i], data + i*hdr->SPR, hdr->SPR);
				}
				else {
					for (k1=0, ix = i*hdr->SPR; k1 < SPR; k1++)
						data[ix+k1] = lei16p(Ptr2datablock + 2*k1);
				}
				len += en1064.Section6.inlen[i];
				Ptr2datablock += en1064.Section6.inlen[i]; 

				if (aECG->FLAG.DIFF==1) {
					for (ix = i*hdr->SPR+1; ix < i*hdr->SPR + SPR; ix++)
						data[ix] += data[ix-1];
				}		
				else if (aECG->FLAG.DIFF==2) {
					for (ix = i*hdr->SPR+2; ix < i*hdr->SPR + SPR; ix++)
						data[ix] += 2*data[ix-1] - data[ix-2];
				}

#ifndef WITHOUT_SCP_DECODE
				if (aECG->FLAG.BIMODAL || en1064.FLAG.REF_BEAT) { 	
//				if (aECG->FLAG.BIMODAL) {
//				if (aECG->FLAG.REF_BEAT {
					/*	this is experimental work
						Bimodal and RefBeat decompression are under development. 
						"continue" ignores code below 
						AS_DECODE=1 will call later SCP-DECODE instead  
					*/ 	
					AS_DECODE = 1; continue; 
				}
#endif

				if (aECG->FLAG.BIMODAL) {
					// ### FIXME ### 
					ix = i*hdr->SPR;		// memory offset
					k1 = en1064.Section4.SPR-1; 	// SPR of decimated data 
					k2 = hdr->SPR;			// SPR of sample data
					uint32_t k3 = en1064.Section4.N-1; // # of protected zones 
					uint8_t  k4 = 4;		// decimation factor 
					do {
						--k2;
						data[ix + k2] = data[ix + k1];
						if (k2 > en1064.Section4.beat[k3].QE) { // outside protected zone 
							if (--k4==0) {k4=4; --k1; };
						}
						else {	// inside protected zone 
							--k1;
							if (k2<en1064.Section4.beat[k3].QB) {--k3; k4=4;};
						}
					} while (k2 && (k1>=0));
				}
			
				if (en1064.FLAG.REF_BEAT) {
					/* Add reference beats */
					// ### FIXME ### 
					for (k1 = 0; k1 < en1064.Section4.N; k1++) {
						if (en1064.Section4.beat[k1].btyp == 0)
						for (ix = 0; ix < en1064.Section5.Length; ix++) {
							int32_t ix1 = en1064.Section4.beat[k1].SB - en1064.Section4.beat[k1].fcM + ix;
							int32_t ix2 = i*hdr->SPR + ix1;
							if ((en1064.Section4.beat[k1].btyp==0) && (ix1 >= 0) && (ix1 < hdr->SPR))
								data[ix2] = data[ix2] * cal6 + en1064.Section5.datablock[i*en1064.Section5.Length+ix] * cal5;
						}
					}
				}
			}

			en1064.Section6.datablock = data; 
			hdr->FLAG.SWAP  = 0;

			curSectPos += 6 + 2*hdr->NS + len;
		}

		/**** SECTION 7 ****/
		else if (curSect==7)  {
			uint16_t N_QRS = *(uint8_t*)(PtrCurSect+curSectPos)-1;
			uint8_t  N_PaceMaker = *(uint8_t*)(PtrCurSect+curSectPos+1);
			uint16_t RRI = leu16p(PtrCurSect+curSectPos+2);
			uint16_t PPI = leu16p(PtrCurSect+curSectPos+4);
			curSectPos += 6;
			size_t curSectPos0 = curSectPos; // backup of pointer 
			
			// skip data on QRS measurements 
			/*
			// ### FIXME ### 
			It seems that the P,QRS, and T wave events can not be reconstructed 
			because they refer to the reference beat and not to the overall signal data.  	
			Maybe Section 4 information need to be used. However, EN1064 does not mention this.
			
			hdr->EVENT.POS = (uint32_t*)realloc(hdr->EVENT.POS, (hdr->EVENT.N+5*N_QRS+N_PaceMaker)*sizeof(*hdr->EVENT.POS));
			hdr->EVENT.TYP = (uint16_t*)realloc(hdr->EVENT.TYP, (hdr->EVENT.N+5*N_QRS+N_PaceMaker)*sizeof(*hdr->EVENT.TYP));
			hdr->EVENT.DUR = (uint32_t*)realloc(hdr->EVENT.DUR, (hdr->EVENT.N+5*N_QRS+N_PaceMaker)*sizeof(*hdr->EVENT.DUR));
			hdr->EVENT.CHN = (uint16_t*)realloc(hdr->EVENT.CHN, (hdr->EVENT.N+5*N_QRS+N_PaceMaker)*sizeof(*hdr->EVENT.CHN));
			for (i=0; i < 5*N_QRS; i++) {
				hdr->EVENT.DUR[hdr->EVENT.N+i] = 0;
				hdr->EVENT.CHN[hdr->EVENT.N+i] = 0;
			}
			for (i=0; i < 5*N_QRS; i+=5) {
				uint8_t typ = *(PtrCurSect+curSectPos+i);
				hdr->EVENT.TYP[hdr->EVENT.N]   = 0x0502;
				hdr->EVENT.TYP[hdr->EVENT.N+1] = 0x8502;
				hdr->EVENT.TYP[hdr->EVENT.N+2] = 0x0503;
				hdr->EVENT.TYP[hdr->EVENT.N+3] = 0x8503;
				hdr->EVENT.TYP[hdr->EVENT.N+4] = 0x8506;
				hdr->EVENT.POS[hdr->EVENT.N]   = leu16p(PtrCurSect+curSectPos0);
				hdr->EVENT.POS[hdr->EVENT.N+1] = leu16p(PtrCurSect+curSectPos0+2);
				hdr->EVENT.POS[hdr->EVENT.N+2] = leu16p(PtrCurSect+curSectPos0+4);
				hdr->EVENT.POS[hdr->EVENT.N+3] = leu16p(PtrCurSect+curSectPos0+6);
				hdr->EVENT.POS[hdr->EVENT.N+4] = leu16p(PtrCurSect+curSectPos0+8);
				hdr->EVENT.N+= 5;
				curSectPos0 += 16;
			}
			*/
			curSectPos += N_QRS*16;
				// pace maker information is stored in sparse sampling channel
			if (N_PaceMaker>0) {
				hdr->EVENT.POS = (uint32_t*)realloc(hdr->EVENT.POS, (hdr->EVENT.N+N_PaceMaker)*sizeof(*hdr->EVENT.POS));
				hdr->EVENT.TYP = (uint16_t*)realloc(hdr->EVENT.TYP, (hdr->EVENT.N+N_PaceMaker)*sizeof(*hdr->EVENT.TYP));
				hdr->EVENT.DUR = (uint32_t*)realloc(hdr->EVENT.DUR, (hdr->EVENT.N+N_PaceMaker)*sizeof(*hdr->EVENT.DUR));
				hdr->EVENT.CHN = (uint16_t*)realloc(hdr->EVENT.CHN, (hdr->EVENT.N+N_PaceMaker)*sizeof(*hdr->EVENT.CHN));
				/* add pacemaker channel */
				hdr->CHANNEL = (CHANNEL_TYPE *) realloc(hdr->CHANNEL,(++hdr->NS)*sizeof(CHANNEL_TYPE));
				i = hdr->NS;
				hdr->CHANNEL[i].SPR 	    = 0;    // sparse event channel 
				hdr->CHANNEL[i].PhysDimCode = 4275; // PhysDimCode("uV") physical unit "uV"
				hdr->CHANNEL[i].Cal 	    = 1;
				hdr->CHANNEL[i].Off         = 0;
				hdr->CHANNEL[i].OnOff       = 1;    // 1: ON 0:OFF
				strcpy(hdr->CHANNEL[i].Transducer,"Pacemaker");
				hdr->CHANNEL[i].GDFTYP      = 3;  

				// ### these values should represent the true saturation values ###//
				hdr->CHANNEL[i].DigMax      = ldexp(1.0,15)-1;
				hdr->CHANNEL[i].DigMin      = ldexp(-1.0,15);
				hdr->CHANNEL[i].PhysMax     = hdr->CHANNEL[i].DigMax * hdr->CHANNEL[i].Cal;
				hdr->CHANNEL[i].PhysMin     = hdr->CHANNEL[i].DigMin * hdr->CHANNEL[i].Cal;
			}
			// skip pacemaker spike measurements  
			for (i=0; i < N_PaceMaker; i++) {
				++hdr->EVENT.N;
				hdr->EVENT.TYP[hdr->EVENT.N] = 0x7fff;
				hdr->EVENT.CHN[hdr->EVENT.N] = hdr->NS;
				hdr->EVENT.POS[hdr->EVENT.N] = (uint32_t)(leu16p(PtrCurSect+curSectPos)*hdr->SampleRate*1e-3);
				hdr->EVENT.DUR[hdr->EVENT.N] = leu16p(PtrCurSect+curSectPos+2);
				curSectPos += 4;
			}
			// skip pacemaker spike information section  
			curSectPos += N_PaceMaker*6;
			
			// QRS type information 
			N_QRS = leu16p(PtrCurSect+curSectPos);
			curSectPos += 2;

		}

		/**** SECTION 8 ****/
		else if (curSect==8)  {
		}

		/**** SECTION 9 ****/
		else if (curSect==9)  {
		}

		/**** SECTION 10 ****/
		else if (curSect==10)  {
		}

		/**** SECTION 11 ****/
		else if (curSect==11)  {
		}
		else {
		}
	}	
	hdr->Dur[0] = hdr->SPR * dT_us; 
	hdr->Dur[1] = 1000000L; 

	/* free allocated memory */ 
	deallocEN1064(en1064);	


#ifndef WITHOUT_SCP_DECODE
   	if (AS_DECODE==0) return(0); 
    		
/*
---------------------------------------------------------------------------
Copyright (C) 2006  Eugenio Cervesato.
Developed at the Associazione per la Ricerca in Cardiologia - Pordenone - Italy,

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 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
---------------------------------------------------------------------------
*/

	/* Fall back method: 





		+ implements Huffman, reference beat and Bimodal compression. 
		- uses piece-wise file access
		- defines intermediate data structure
	*/	

	fprintf(stdout, "\nUse SCP_DECODE (Huffman=%i RefBeat=%i Bimodal=%i)\n", aECG->FLAG.HUFFMAN, aECG->FLAG.REF_BEAT, aECG->FLAG.BIMODAL);

	textual.des.acquiring.protocol_revision_number = aECG->Section1.Tag14.VERSION;
	textual.des.analyzing.protocol_revision_number = aECG->Section1.Tag15.VERSION; 

	decode.flag_Res.bimodal = (aECG->Section1.Tag14.VERSION > 10 ? aECG->FLAG.BIMODAL : 0);  
	decode.Reconstructed    = (int32_t*) hdr->AS.rawdata; 

	if (scp_decode(hdr, section, decode, record, textual, add_filter)) {
		if (Cal0>1)
			for (i=0; i < hdr->NS * hdr->SPR * hdr->NRec; ++i)
				data[i] /= Cal0;
		hdr->FLAG.SWAP  = 0;
	}
	else { 
		B4C_ERRNUM = B4C_CANNOT_OPEN_FILE;
		B4C_ERRMSG = "SCP-DECODE can not read file"; 
		return(0);
	}
	
	 // end of fall back method 
	return(1);
#endif 
};