1 /*
2
3 Copyright (C) 2008-2013,2018 Alois Schloegl <alois.schloegl@gmail.com>
4 This file is part of the "BioSig for C/C++" repository
5 (biosig4c++) at http://biosig.sf.net/
6
7 BioSig is free software; you can redistribute it and/or
8 modify it under the terms of the GNU General Public License
9 as published by the Free Software Foundation; either version 3
10 of the License, or (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program. If not, see <http://www.gnu.org/licenses/>.
19
20
21 */
22
23 #include <assert.h>
24 #include <ctype.h>
25 #include <stdlib.h>
26 #include <string.h>
27 #ifdef _WIN32
28 // Can't include sys/stat.h or sopen is declared twice.
29 #include <sys/types.h>
30 struct stat {
31 _dev_t st_dev;
32 _ino_t st_ino;
33 unsigned short st_mode;
34 short st_nlink;
35 short st_uid;
36 short st_gid;
37 _dev_t st_rdev;
38 _off_t st_size;
39 time_t st_atime;
40 time_t st_mtime;
41 time_t st_ctime;
42 };
43 int __cdecl stat(const char *_Filename,struct stat *_Stat);
44 #else
45 #include <sys/stat.h>
46 #endif
47
48 #include "../biosig.h"
49
50 /* TODO:
51 - need to separate sopen_heka() and sread_heka()
52 - data swapping
53 */
54
55 #define min(a,b) (((a) < (b)) ? (a) : (b))
56 #define max(a,b) (((a) > (b)) ? (a) : (b))
57
58 /****************************************************************************
59 rational :
60 computes the rational approximation of a floating point number
61 such that n/d is an approximation for r with an relative
62 error smaller than tol
63
64 see Octave's rat.m
65 ****************************************************************************/
rational(double x,double tol,long * n,long * d)66 void rational (double x, double tol, long *n, long *d) {
67
68 if (x != x) { // i.e. isnan(x)
69 *n = 0;
70 *d = 0;
71 return;
72 }
73
74 if (!finite(x)) {
75 *n = x>0; // i.e. sign(x)
76 *d = 0;
77 return;
78 }
79
80 tol *= fabs(x);
81 *n = lround(x);
82 *d = 1;
83 double frac = x - *n;
84 long lastn = 1, lastd = 0;
85
86 while (fabs((*d) * x - (*n) ) >= fabs((*d) * tol)) {
87 double flip = 1.0/frac;
88 long step = lround(flip);
89 frac = flip - step;
90
91 long nextn = *n, nextd = *d;
92 *n = *n * step + lastn;
93 *d = *d * step + lastd;
94 lastn = nextn;
95 lastd = nextd;
96 }
97
98 if (*d < 0) {
99 *n = - *n;
100 *d = - *d;
101 }
102 }
103
104
105 /****************************************************************************
106 heka2gdftime
107 converts heka time format into gdftime
108 ****************************************************************************/
heka2gdftime(double t)109 gdf_time heka2gdftime(double t) {
110 t -= 1580970496; if (t<0) t += 4294967296; t += 9561652096;
111 return (uint64_t)ldexp(t/(24.0*60*60) + 584755, 32); // +datenum(1601,1,1));
112 }
113
114 /****************************************************************************
115 sopen_heka
116 reads heka format
117
118 if itx is not null, the file is converted into an ITX formated file
119 and streamed to itx, too.
120 ****************************************************************************/
sopen_heka(HDRTYPE * hdr,FILE * itx)121 void sopen_heka(HDRTYPE* hdr, FILE *itx) {
122 size_t count = hdr->HeadLen;
123
124 if (hdr->TYPE==HEKA && hdr->VERSION > 1) {
125
126 int32_t Levels=0;
127 uint16_t k;
128 //int32_t *Sizes=NULL;
129 int32_t Counts[5], counts[5]; //, Sizes[5];
130 memset(Counts,0,20);
131 memset(counts,0,20);
132 //memset(Sizes,0,20);
133 uint32_t StartOfData=0,StartOfPulse=0;
134
135 union {
136 struct {
137 int32_t Root;
138 int32_t Group;
139 int32_t Series;
140 int32_t Sweep;
141 int32_t Trace;
142 } Rec;
143 int32_t all[5];
144 } Sizes;
145
146
147 // HEKA PatchMaster file format
148
149 count = hdr->HeadLen;
150 struct stat FileBuf;
151 stat(hdr->FileName,&FileBuf);
152 hdr->AS.Header = (uint8_t*)realloc(hdr->AS.Header, FileBuf.st_size);
153 if (count < 1024)
154 count += ifread(hdr->AS.Header+count, 1, 1024-count, hdr);
155 hdr->HeadLen = count;
156
157 hdr->FILE.LittleEndian = *(uint8_t*)(hdr->AS.Header+52) > 0;
158 char SWAP = ( hdr->FILE.LittleEndian && (__BYTE_ORDER == __BIG_ENDIAN)) \
159 || (!hdr->FILE.LittleEndian && (__BYTE_ORDER == __LITTLE_ENDIAN));
160
161 if (SWAP) {
162 biosigERROR(hdr, B4C_FORMAT_UNSUPPORTED, "Heka/Patchmaster format requires data swapping - this is not supported yet.");
163 return;
164 }
165 SWAP = 0; // might be useful for compile time optimization
166
167 /* get file size and read whole file */
168 count += ifread(hdr->AS.Header+count, 1, FileBuf.st_size - count, hdr);
169
170 if (VERBOSE_LEVEL>7) fprintf(stdout,"%s (line %i) %s(...): %i bytes read\n",__FILE__,__LINE__,__func__, (int)count);
171
172 // double oTime;
173 uint32_t nItems;
174 if (hdr->FILE.LittleEndian) {
175 // oTime = lef64p(hdr->AS.Header+40); // not used
176 nItems = leu32p(hdr->AS.Header+48);
177 }
178 else {
179 // oTime = bef64p(hdr->AS.Header+40); // not used
180 nItems = beu32p(hdr->AS.Header+48);
181 }
182
183 if (VERBOSE_LEVEL>7) fprintf(stdout,"%s (line %i) %s(...): nItems=%i\n",__FILE__,__LINE__,__func__, nItems);
184
185 if (hdr->VERSION == 1) {
186 Sizes.Rec.Root = 544;
187 Sizes.Rec.Group = 128;
188 Sizes.Rec.Series = 1120;
189 Sizes.Rec.Sweep = 160;
190 Sizes.Rec.Trace = 296;
191 }
192 else if (hdr->VERSION == 2)
193 for (k=0; k < min(12,nItems); k++) {
194
195 if (VERBOSE_LEVEL>7) fprintf(stdout,"%s (line %i): HEKA nItems=%i\n",__func__,__LINE__, k);
196
197 uint32_t start = *(uint32_t*)(hdr->AS.Header+k*16+64);
198 uint32_t length = *(uint32_t*)(hdr->AS.Header+k*16+64+4);
199 if (SWAP) {
200 start = bswap_32(start);
201 length = bswap_32(length);
202 }
203 uint8_t *ext = hdr->AS.Header + k*16 + 64 + 8;
204
205 if (VERBOSE_LEVEL>7) fprintf(stdout,"%s (line %i): HEKA #%i: <%s> [%i:+%i]\n",__func__,__LINE__,k,ext,start,length);
206
207 if (!start) break;
208
209 if ((start+8) > count) {
210 biosigERROR(hdr, B4C_INCOMPLETE_FILE, "Heka/Patchmaster: file is corrupted - segment with pulse data is not available!");
211 return;
212 }
213
214 if (!memcmp(ext,".pul\0\0\0\0",8)) {
215 // find pulse data
216 ifseek(hdr, start, SEEK_SET);
217
218 //magic = *(int32_t*)(hdr->AS.Header+start);
219 Levels = *(int32_t*)(hdr->AS.Header+start+4);
220 if (SWAP) Levels = bswap_32(Levels);
221 if (Levels>5) {
222 biosigERROR(hdr, B4C_FORMAT_UNSUPPORTED, "Heka/Patchmaster format with more than 5 levels not supported");
223 return;
224 }
225
226 if (VERBOSE_LEVEL>7) fprintf(stdout,"%s (line %i): HEKA #%i Levels=%i\n",__func__,__LINE__,k,Levels);
227
228 memcpy(Sizes.all,hdr->AS.Header+start+8,sizeof(int32_t)*Levels);
229 if (SWAP) {
230 int l;
231 for (l=0; l < Levels; l++) Sizes.all[l] = bswap_32(Sizes.all[l]);
232 }
233
234 if (VERBOSE_LEVEL>7) {int l; for (l=0; l < Levels; l++) fprintf(stdout,"%s (line %i): HEKA #%i %i\n",__func__,__LINE__,l, Sizes.all[l]); }
235
236 StartOfPulse = start + 8 + 4 * Levels;
237 }
238
239 else if (!memcmp(ext,".dat\0\0\0\0",8)) {
240 StartOfData = start;
241 }
242 }
243
244 if (VERBOSE_LEVEL>7) fprintf(stdout,"%s (line %i) %s(...): \n",__FILE__,__LINE__,__func__);
245
246 // if (!Sizes) free(Sizes); Sizes=NULL;
247
248 /* DONE: HEKA, check channel number and label
249 pass 1:
250 + get number of sweeps
251 + get number of channels
252 + check whether all traces of a single sweep have the same SPR, and Fs
253 + check whether channelnumber (TrAdcChannel), scaling (DataScaler) and Label fit among all sweeps
254 + extract the total number of samples
255 + physical units
256 + level 4 may have no children
257 + count event descriptions Level2/SeLabel
258 pass 2:
259 + initialize data to NAN
260 + skip sweeps if selected channel is not in it
261 + Y scale, physical scale
262 + Event.CodeDescription, Events,
263 resampling
264 */
265
266 uint32_t k1=0, k2=0, k3=0, k4=0;
267 uint32_t K1=0, K2=0, K3=0, K4=0, K5=0;
268 double t;
269 size_t pos;
270
271 // read K1
272 if (SWAP) {
273 K1 = bswap_32(*(uint32_t*)(hdr->AS.Header + StartOfPulse + Sizes.Rec.Root));
274 hdr->VERSION = bswap_32(*(uint32_t*)(hdr->AS.Header + StartOfPulse));
275 union {
276 double f64;
277 uint64_t u64;
278 } c;
279 c.u64 = bswap_64(*(uint64_t*)(hdr->AS.Header + StartOfPulse + 520));
280 t = c.f64;
281 } else {
282 K1 = (*(uint32_t*)(hdr->AS.Header + StartOfPulse + Sizes.Rec.Root));
283 hdr->VERSION = (*(uint32_t*)(hdr->AS.Header + StartOfPulse));
284 t = (*(double*)(hdr->AS.Header + StartOfPulse + 520));
285 }
286
287 if (VERBOSE_LEVEL>7) fprintf(stdout,"%s (line %i) %s(...): \n",__FILE__,__LINE__,__func__);
288
289 hdr->T0 = heka2gdftime(t); // this is when when heka was started, data is recorded later.
290 hdr->SampleRate = 0.0;
291 double *DT = NULL; // list of sampling intervals per channel
292 hdr->SPR = 0;
293
294 if (VERBOSE_LEVEL>7) fprintf(stdout,"%s (line %i) %s(...): %p\n",__FILE__,__LINE__,__func__,hdr->EVENT.CodeDesc);
295
296 /*******************************************************************************************************
297 HEKA: read structural information
298 *******************************************************************************************************/
299
300 pos = StartOfPulse + Sizes.Rec.Root + 4;
301 size_t EventN=0;
302
303 for (k1=0; k1<K1; k1++) {
304 // read group
305
306 if (VERBOSE_LEVEL>7) fprintf(stdout,"HEKA L1 @%i=\t%i/%i \n",(int)(pos+StartOfData),k1,K1);
307
308 pos += Sizes.Rec.Group+4;
309 // read number of children
310 K2 = (*(uint32_t*)(hdr->AS.Header+pos-4));
311
312 hdr->AS.auxBUF = (uint8_t*)realloc(hdr->AS.auxBUF,K2*33); // used to store name of series
313 for (k2=0; k2<K2; k2++) {
314 // read series
315 union {
316 double f64;
317 uint64_t u64;
318 } Delay;
319 char *SeLabel = (char*)(hdr->AS.Header+pos+4); // max 32 bytes
320 strncpy((char*)hdr->AS.auxBUF + 33*k2, (char*)hdr->AS.Header+pos+4, 32); hdr->AS.auxBUF[33*k2+32] = 0;
321 SeLabel = (char*)hdr->AS.auxBUF + 33*k2;
322 double tt = *(double*)(hdr->AS.Header+pos+136); // time of series. TODO: this time should be taken into account
323 Delay.u64 = bswap_64(*(uint64_t*)(hdr->AS.Header+pos+472+176));
324
325 gdf_time t = heka2gdftime(tt);
326
327 if (VERBOSE_LEVEL>7) {
328 char tmp[60];
329 snprintf_gdfdate(tmp, sizeof(tmp), t);
330 fprintf(stdout,"HEKA L2 @%i=%s %f\t%i/%i %i/%i t=%.17g %s\n",(int)(pos+StartOfData),SeLabel,Delay.f64,k1,K1,k2,K2,ldexp(t,-32),tmp);
331 }
332
333 pos += Sizes.Rec.Series + 4;
334 // read number of children
335 K3 = (*(uint32_t*)(hdr->AS.Header+pos-4));
336 if (EventN <= hdr->EVENT.N + K3 + 2) {
337 EventN = max(max(16,EventN),hdr->EVENT.N+K3+2) * 2;
338 if (reallocEventTable(hdr, EventN) == SIZE_MAX) {
339 biosigERROR(hdr, B4C_MEMORY_ALLOCATION_FAILED, "Allocating memory for event table failed.");
340 return;
341 };
342 }
343
344 if (!hdr->AS.SegSel[0] && !hdr->AS.SegSel[1] && !hdr->AS.SegSel[2]) {
345 // in case of reading the whole file (no sweep selection), include marker for start of series
346 FreeTextEvent(hdr, hdr->EVENT.N, SeLabel);
347 hdr->EVENT.POS[hdr->EVENT.N] = hdr->SPR; // within reading the structure, hdr->SPR is used as a intermediate variable counting the number of samples
348 #if (BIOSIG_VERSION >= 10500)
349 hdr->EVENT.TimeStamp[hdr->EVENT.N] = t;
350 #endif
351 hdr->EVENT.N++;
352 }
353
354 for (k3=0; k3<K3; k3++) {
355 // read sweep
356 hdr->NRec++; // increase number of sweeps
357 size_t SPR = 0, spr = 0;
358 gdf_time t = heka2gdftime(*(double*)(hdr->AS.Header+pos+48)); // time of sweep. TODO: this should be taken into account
359
360 if (VERBOSE_LEVEL>7) {
361 char tmp[60];
362 snprintf_gdfdate(tmp, sizeof(tmp), t);
363 fprintf(stdout,"HEKA L3 @%i= %fHz\t%i/%i %i/%i %i/%i %s\n",(int)(pos+StartOfData),hdr->SampleRate,k1,K1,k2,K2,k3,K3,tmp);
364 }
365
366 char flagSweepSelected = (hdr->AS.SegSel[0]==0 || k1+1==hdr->AS.SegSel[0])
367 && (hdr->AS.SegSel[1]==0 || k2+1==hdr->AS.SegSel[1])
368 && (hdr->AS.SegSel[2]==0 || k3+1==hdr->AS.SegSel[2]);
369
370 // hdr->SPR
371 if (hdr->SPR==0)
372 hdr->T0 = t; // start time of first recording determines the start time of the recording
373 else if (flagSweepSelected && hdr->SPR > 0) {
374 // marker for start of sweep
375 hdr->EVENT.POS[hdr->EVENT.N] = hdr->SPR; // within reading the structure, hdr->SPR is used as a intermediate variable counting the number of samples
376 hdr->EVENT.TYP[hdr->EVENT.N] = 0x7ffe;
377 #if (BIOSIG_VERSION >= 10500)
378 hdr->EVENT.TimeStamp[hdr->EVENT.N] = t;
379 #endif
380 hdr->EVENT.N++;
381 }
382
383 pos += Sizes.Rec.Sweep + 4;
384 // read number of children
385 K4 = (*(uint32_t*)(hdr->AS.Header+pos-4));
386 for (k4=0; k4<K4; k4++) {
387 // read trace
388 double DigMin, DigMax;
389 uint16_t gdftyp = 0;
390 uint32_t ns = (*(uint32_t*)(hdr->AS.Header+pos+36));
391 uint32_t DataPos = (*(uint32_t*)(hdr->AS.Header+pos+40));
392 spr = (*(uint32_t*)(hdr->AS.Header+pos+44));
393 double DataScaler= (*(double*)(hdr->AS.Header+pos+72));
394 double Toffset = (*(double*)(hdr->AS.Header+pos+80)); // time offset of
395 uint16_t pdc = PhysDimCode((char*)(hdr->AS.Header + pos + 96));
396 double dT = (*(double*)(hdr->AS.Header+pos+104));
397 //double XStart = (*(double*)(hdr->AS.Header+pos+112));
398 uint16_t XUnits = PhysDimCode((char*)(hdr->AS.Header+pos+120));
399 double YRange = (*(double*)(hdr->AS.Header+pos+128));
400 double YOffset = (*(double*)(hdr->AS.Header+pos+136));
401 double Bandwidth = (*(double*)(hdr->AS.Header+pos+144));
402 //double PipetteResistance = (*(double*)(hdr->AS.Header+pos+152));
403 double RsValue = (*(double*)(hdr->AS.Header+pos+192));
404
405 uint8_t ValidYRange = hdr->AS.Header[pos+220];
406 uint16_t AdcChan = (*(uint16_t*)(hdr->AS.Header+pos+222));
407 /* obsolete: range is defined by DigMin/DigMax * DataScaler + YOffset
408 double PhysMin = (*(double*)(hdr->AS.Header+pos+224));
409 double PhysMax = (*(double*)(hdr->AS.Header+pos+232));
410 */
411
412 if (VERBOSE_LEVEL>7) fprintf(stdout, "%s (line %i): %i %i %i %i %i %g %g 0x%x xUnits=%i %g %g %g %g %i %i\n", __FILE__,__LINE__, k1, k2, k3, k4, ns, DataScaler, Toffset, pdc, XUnits, YRange, YOffset, Bandwidth, RsValue, ValidYRange, AdcChan);
413
414 switch (hdr->AS.Header[pos+70]) {
415 case 0: gdftyp = 3; //int16
416 /*
417 It seems that the range is 1.024*(2^15-1)/2^15 nA or V
418 and symetric around zero. i.e. YOffset is zero
419 */
420 DigMax = ldexp(1.0,15) - 1.0;
421 DigMin = -DigMax;
422 break;
423 case 1: gdftyp = 5; //int32
424 DigMax = ldexp(1.0, 31) - 1.0;
425 DigMin = -DigMax;
426 break;
427 case 2: gdftyp = 16; //float32
428 DigMax = 1e9;
429 DigMin = -1e9;
430 break;
431 case 3: gdftyp = 17; //float64
432 DigMax = 1e9;
433 DigMin = -1e9;
434 break;
435 default:
436 DigMax = NAN;
437 DigMin = NAN;
438 biosigERROR(hdr, B4C_FORMAT_UNSUPPORTED, "Heka/Patchmaster: data type not supported");
439 };
440
441 if (SWAP) {
442 AdcChan = bswap_16(AdcChan);
443 ns = bswap_32(ns);
444 DataPos = bswap_32(DataPos);
445 spr = bswap_32(spr);
446 // avoid breaking strict-aliasing rules
447 union {
448 double f64;
449 uint64_t u64;
450 } c;
451 c.f64 = dT; c.u64 = bswap_64(c.u64); dT = c.f64;
452 c.f64 = YRange; c.u64 = bswap_64(c.u64); YRange = c.f64;
453 c.f64 = YOffset; c.u64 = bswap_64(c.u64); YOffset = c.f64;
454 //c.f64 = PhysMax; c.u64 = bswap_64(c.u64); PhysMax = c.f64;
455 //c.f64 = PhysMin; c.u64 = bswap_64(c.u64); PhysMin = c.f64;
456 c.f64 = Toffset; c.u64 = bswap_64(c.u64); Toffset = c.f64;
457 }
458
459 if (YOffset != 0.0)
460 fprintf(stderr,"!!! WARNING !!! HEKA: the offset is not zero - "
461 "this case is not tested and might result in incorrect scaling of "
462 "the data,\n!!! YOU ARE WARNED !!!\n");
463
464 // scale to standard units - no prefix
465 double Cal = DataScaler * PhysDimScale(pdc);
466 double Off = YOffset * PhysDimScale(pdc);
467 pdc &= 0xffe0;
468 float Fs = 1.0 / ( dT * PhysDimScale(XUnits) ) ; // float is used to avoid spurios accuracy, round to single precision accuracy
469
470 if (flagSweepSelected) {
471
472 if (hdr->SampleRate <= 0.0) hdr->SampleRate = Fs;
473 if (fabs(hdr->SampleRate - Fs) > 1e-9*Fs) {
474 long DIV1 = 1, DIV2 = 1;
475 rational(hdr->SampleRate*dT*PhysDimScale(XUnits), 1e-6, &DIV2, &DIV1);
476
477 if (DIV1 > 1) {
478 if ( ((size_t)DIV1 * hdr->SPR) > 0xffffffffffffffff) {
479 fprintf(stderr,"!!! WARNING sopen_heka(%s) !!! due to resampling, the data will have more then 2^31 samples !!!\n", hdr->FileName);
480 biosigERROR(hdr,B4C_FORMAT_UNSUPPORTED,"HEKA file has more than 2^32 samples - this is not supported yet");
481 }
482 hdr->SPR *= DIV1;
483 hdr->SampleRate *= DIV1;
484 hdr->EVENT.SampleRate = hdr->SampleRate;
485 size_t n = 0;
486 while (n < hdr->EVENT.N)
487 hdr->EVENT.POS[n++] *= DIV1;
488 }
489 if (DIV2 > 1) spr *= DIV2;
490 }
491
492 // samples per sweep
493 if (k4==0) SPR = spr;
494 else if (SPR != spr) {
495 biosigERROR(hdr, B4C_FORMAT_UNSUPPORTED, "Heka/Patchmaster: number of samples among channels within a single sweep do not match.");
496 return;
497 }
498 }
499
500 char *Label = (char*)hdr->AS.Header+pos+4;
501 for (ns=0; ns < hdr->NS; ns++) {
502 if (!strcmp(hdr->CHANNEL[ns].Label,Label)) break;
503 }
504
505 if (VERBOSE_LEVEL>7) fprintf(stdout,"HEKA L4 @%i= #%i,%i, %s %f %fHz\t%i/%i %i/%i %i/%i %i/%i \n",(int)(pos+StartOfData),ns,AdcChan,Label,hdr->SampleRate,Fs,k1,K1,k2,K2,k3,K3,k4,K4);
506
507 CHANNEL_TYPE *hc;
508 if (ns >= hdr->NS) {
509 hdr->NS = ns + 1;
510 #ifdef WITH_TIMESTAMPCHANNEL
511 // allocate memory for an extra time stamp channel, which is define only after the end of the channel loop - see below
512 hdr->CHANNEL = (CHANNEL_TYPE*) realloc(hdr->CHANNEL, (hdr->NS + 1) * sizeof(CHANNEL_TYPE));
513 #else
514 hdr->CHANNEL = (CHANNEL_TYPE*) realloc(hdr->CHANNEL, hdr->NS * sizeof(CHANNEL_TYPE));
515 #endif
516 hc = hdr->CHANNEL + ns;
517 strncpy(hc->Label, Label, max(32, MAX_LENGTH_LABEL));
518 hc->Label[max(32,MAX_LENGTH_LABEL)] = 0;
519 hc->Transducer[0] = 0;
520 hc->SPR = 1;
521 hc->PhysDimCode = pdc;
522 hc->OnOff = 1;
523 hc->GDFTYP = gdftyp;
524 hc->LeadIdCode = 0;
525 hc->DigMin = DigMin;
526 hc->DigMax = DigMax;
527
528 // TODO: case of non-zero YOffset is not tested //
529 hc->PhysMax = DigMax * Cal + Off;
530 hc->PhysMin = DigMin * Cal + Off;
531
532 hc->Cal = Cal;
533 hc->Off = Off;
534 hc->TOffset = Toffset;
535
536 #ifndef NDEBUG
537 double Cal2 = (hc->PhysMax - hc->PhysMin) / (hc->DigMax - hc->DigMin);
538 double Off2 = hc->PhysMin - Cal2 * hc->DigMin;
539 double Off3 = hc->PhysMax - Cal2 * hc->DigMax;
540
541 assert(fabs(Cal-Cal2) < 1e-8 * Cal);
542 assert(fabs(Off-Off2) < 1e-8 * Cal);
543 assert(fabs(Off-Off3) < 1e-8 * Cal);
544 #endif
545
546 /* TODO: fix remaining channel header */
547 /* LowPass, HighPass, Notch, Impedance, */
548 hc->HighPass = NAN;
549 hc->LowPass = (Bandwidth > 0) ? Bandwidth : NAN;
550 hc->Notch = NAN;
551 hc->Impedance = (RsValue > 0) ? RsValue : NAN;
552
553 DT = (double*) realloc(DT, hdr->NS*sizeof(double));
554 DT[ns] = dT;
555 }
556 else {
557 /*
558 channel has been already defined in earlier sweep.
559 check compatibility and adapt internal format when needed
560 */
561 hc = hdr->CHANNEL + ns;
562 double PhysMax = DigMax * Cal + Off;
563 double PhysMin = DigMin * Cal + Off;
564 // get max value to avoid false positive saturation detection when scaling changes
565 if (hc->PhysMax < PhysMax) hc->PhysMax = PhysMax;
566 if (hc->PhysMin > PhysMin) hc->PhysMin = PhysMin;
567
568 if (hc->GDFTYP < gdftyp) {
569 /* when data type changes, use the largest data type */
570 if (4 < hc->GDFTYP && hc->GDFTYP < 9 && gdftyp==16)
571 /* (U)INT32, (U)INT64 + FLOAT32 -> DOUBLE */
572 hc->GDFTYP = 17;
573 else
574 hc->GDFTYP = gdftyp;
575 }
576 else if (hc->GDFTYP > gdftyp) {
577 /* when data type changes, use the largest data type */
578 if (4 < gdftyp && gdftyp < 9 && hc->GDFTYP==16)
579 /* (U)INT32, (U)INT64 + FLOAT32 -> DOUBLE */
580 hc->GDFTYP = 17;
581 }
582
583 if (fabs(hc->Cal - Cal) > 1e-9*Cal) {
584 /* when scaling changes from sweep to sweep, use floating point numbers internally. */
585 if (hc->GDFTYP < 5) // int16 or smaller
586 hc->GDFTYP = 16;
587 else if (hc->GDFTYP < 9) // int32, int64 -> double
588 hc->GDFTYP = 17;
589 }
590
591 if ((pdc & 0xFFE0) != (hc->PhysDimCode & 0xFFE0)) {
592 fprintf(stdout, "ERROR: [%i,%i,%i,%i] Yunits in %s do not match %04x(%s) ! %04x(%s)\n",k1,k2,k3,k4, Label, pdc, PhysDim3(pdc), hc->PhysDimCode, PhysDim3(hc->PhysDimCode));
593 biosigERROR(hdr, B4C_FORMAT_UNSUPPORTED, "Heka/Patchmaster: Yunits do not match");
594 }
595 if ( ( VERBOSE_LEVEL > 7 ) && ( fabs( DT[ns] - dT) > 1e-9 * dT) ) {
596 fprintf(stdout, "%s (line %i) different sampling rates [%i,%i,%i,%i]#%i,%f/%f \n",__FILE__,__LINE__,(int)k1,(int)k2,(int)k3,(int)k4,(int)ns, 1.0/DT[ns],1.0/dT);
597 }
598 }
599
600 if (YOffset) {
601 biosigERROR(hdr, B4C_FORMAT_UNSUPPORTED, "Heka/Patchmaster: YOffset is not zero");
602 }
603 if (hdr->AS.Header[pos+220] != 1) {
604 fprintf(stderr,"WARNING Heka/Patchmaster: ValidYRange not set to 1 but %i in sweep [%i,%i,%i,%i]\n", hdr->AS.Header[pos+220],k1+1,k2+1,k3+1,k4+1);
605 }
606
607 if (VERBOSE_LEVEL>7) fprintf(stdout,"HEKA L6 @%i= #%i,%i, %s %f-%fHz\t%i/%i %i/%i %i/%i %i/%i \n",(int)(pos+StartOfData),ns,AdcChan,Label,hdr->SampleRate,Fs,k1,K1,k2,K2,k3,K3,k4,K4);
608
609 pos += Sizes.Rec.Trace+4;
610 // read number of children -- this should be 0 - ALWAYS;
611 K5 = (*(uint32_t*)(hdr->AS.Header+pos-4));
612 if (K5) {
613 biosigERROR(hdr, B4C_FORMAT_UNSUPPORTED, "Heka/Patchmaster: Level 4 has some children");
614 }
615 } // end loop k4
616
617 // if sweep is selected, add number of samples to counter
618 if (flagSweepSelected) {
619 if ( hdr->SPR > 0xffffffffffffffffu-SPR) {
620 biosigERROR(hdr,B4C_FORMAT_UNSUPPORTED,"HEKA file has more than 2^32 samples - this is not supported yet");
621 }
622 hdr->SPR += SPR;
623 }
624 } // end loop k3
625 } // end loop k2
626 } // end loop k1
627
628 #ifndef NO_BI
629 if (DT) free(DT);
630 #else
631 size_t *BI = (size_t*) DT; // DT is not used anymore, use space for BI
632 #endif
633 DT = NULL;
634
635 #ifdef WITH_TIMESTAMPCHANNEL
636 {
637 /*
638 define time stamp channel, memory is already allocated above
639 */
640 CHANNEL_TYPE *hc = hdr->CHANNEL + hdr->NS;
641 hc->GDFTYP = 7; // corresponds to int64_t, gdf_time
642 strcpy(hc->Label,"Timestamp");
643 hc->Transducer[0]=0;
644 hc->PhysDimCode = 2272; // units: days [d]
645 hc->LeadIdCode = 0;
646 hc->SPR = 1;
647 hc->Cal = ldexp(1.0, -32);
648 hc->Off = 0.0;
649 hc->OnOff = 1;
650 hc->DigMax = ldexp( 1.0, 61);
651 hc->DigMin = 0;
652 hc->PhysMax = hc->DigMax * hc->Cal;
653 hc->PhysMin = hc->DigMin * hc->Cal;
654 hc->TOffset = 0.0;
655 hc->Impedance = NAN;
656 hc->HighPass = NAN;
657 hc->LowPass = NAN;
658 hc->Notch = NAN;
659 hc->XYZ[0] = 0.0;
660 hc->XYZ[1] = 0.0;
661 hc->XYZ[2] = 0.0;
662
663 hdr->NS++;
664 }
665 #endif
666
667 hdr->NRec = 1;
668 hdr->AS.bpb = 0;
669 for (k = 0; k < hdr->NS; k++) {
670 CHANNEL_TYPE *hc = hdr->CHANNEL + k;
671
672 hc->Cal = (hc->PhysMax - hc->PhysMin) / (hc->DigMax - hc->DigMin);
673 hc->Off = hc->PhysMin - hc->DigMin * hc->Cal;
674 #ifndef NO_BI
675 hc->bi = hdr->AS.bpb;
676 #else
677 BI[k] = hdr->AS.bpb;
678 #endif
679 hc->SPR = hdr->SPR;
680 hdr->AS.bpb += hc->SPR * (GDFTYP_BITS[hc->GDFTYP]>>3); // multiplation must not exceed 32 bit limit
681 }
682
683 if (hdr->AS.B4C_ERRNUM) {
684 #ifdef NO_BI
685 if (BI) free(BI);
686 #endif
687 return;
688 }
689 hdr->ID.Manufacturer.Name = "HEKA/Patchmaster";
690
691
692
693 /******************************************************************************
694 SREAD_HEKA
695
696 void sread_heka(HDRTYPE* hdr, FILE *itx, ... ) {
697
698 ******************************************************************************/
699
700 if (VERBOSE_LEVEL > 7) fprintf(stdout,"HEKA: 400: %"PRIi64" %"PRIi32" %"PRIi64"\n",hdr->NRec, hdr->AS.bpb, hdr->NRec * (size_t)hdr->AS.bpb);
701
702 size_t sz = hdr->NRec * (size_t)hdr->AS.bpb;
703 if (sz/hdr->NRec < hdr->AS.bpb) {
704 biosigERROR(hdr, B4C_MEMORY_ALLOCATION_FAILED, "memory allocation failed - more than 2GB required but platform supports only 32 bit!");
705 return;
706 }
707
708 void* tmpptr = realloc(hdr->AS.rawdata, sz);
709 if (tmpptr!=NULL)
710 hdr->AS.rawdata = (uint8_t*) tmpptr;
711 else {
712 biosigERROR(hdr, B4C_MEMORY_ALLOCATION_FAILED, "memory allocation failed - not enough memory!");
713 return;
714 }
715 assert(hdr->NRec >= 0);
716 memset(hdr->AS.rawdata, 0xff, hdr->NRec * (size_t)hdr->AS.bpb); // initialize with NAN's
717
718
719 #ifdef NO_BI
720 #define _BI (BI[k])
721 #else
722 #define _BI (hc->bi)
723 #endif
724 /* initialize with NAN's */
725 for (k=0; k<hdr->NS; k++) {
726 size_t k1;
727 CHANNEL_TYPE *hc = hdr->CHANNEL+k;
728 switch (hc->GDFTYP) {
729 case 3:
730 for (k1=0; k1<hc->SPR; k1++) {
731 *(uint16_t*)(hdr->AS.rawdata + _BI + k1 * 2) = 0x8000;
732 }
733 break;
734 case 5:
735 for (k1=0; k1<hc->SPR; k1++) *(uint32_t*)(hdr->AS.rawdata + _BI + k1 * 4) = 0x80000000;
736 break;
737 case 7:
738 for (k1=0; k1<hc->SPR; k1++) *(int64_t*)(hdr->AS.rawdata + _BI + k1 * 4) = 0x8000000000000000LL;
739 break;
740 case 16:
741 for (k1=0; k1<hc->SPR; k1++) *(float*)(hdr->AS.rawdata + _BI + k1 * 4) = NAN;
742 break;
743 case 17:
744 for (k1=0; k1<hc->SPR; k1++) *(double*)(hdr->AS.rawdata + _BI + k1 * 8) = NAN;
745 break;
746 }
747 }
748 #undef _BI
749
750 char *WAVENAME = NULL;
751 if (itx) {
752 fprintf(itx, "IGOR\r\nX Silent 1\r\n");
753 const char *fn = strrchr(hdr->FileName,'\\');
754 if (fn) fn++;
755 else fn = strrchr(hdr->FileName,'/');
756 if (fn) fn++;
757 else fn = hdr->FileName;
758
759 size_t len = strspn(fn,".");
760 WAVENAME = (char*)malloc(strlen(hdr->FileName)+7);
761 if (len)
762 strncpy(WAVENAME, fn, len);
763 else
764 strcpy(WAVENAME, fn); // Flawfinder: ignore
765 }
766
767 if (VERBOSE_LEVEL>7) hdr2ascii(hdr,stdout,4);
768
769 /*******************************************************************************************************
770 HEKA: read data blocks
771 *******************************************************************************************************/
772 uint32_t SPR = 0;
773 pos = StartOfPulse + Sizes.Rec.Root + 4;
774 for (k1=0; k1<K1; k1++) {
775 // read group
776
777 if (VERBOSE_LEVEL>7) fprintf(stdout,"HEKA+L1 @%i=\t%i/%i \n",(int)(pos+StartOfData),k1,K1);
778
779 pos += Sizes.Rec.Group+4;
780 // read number of children
781 K2 = (*(uint32_t*)(hdr->AS.Header+pos-4));
782
783 for (k2=0; k2<K2; k2++) {
784 // read series
785 union {
786 double f64;
787 uint64_t u64;
788 } Delay;
789 uint32_t spr = 0;
790 char *SeLabel = (char*)(hdr->AS.Header+pos+4); // max 32 bytes
791 Delay.u64 = bswap_64(*(uint64_t*)(hdr->AS.Header+pos+472+176));
792
793 if (VERBOSE_LEVEL>7) fprintf(stdout,"HEKA+L2 @%i=%s %f\t%i/%i %i/%i \n",(int)(pos+StartOfData),SeLabel,Delay.f64,k1,K1,k2,K2);
794
795 /* move to reading of data */
796 pos += Sizes.Rec.Series+4;
797 // read number of children
798 K3 = (*(uint32_t*)(hdr->AS.Header+pos-4));
799 for (k3=0; k3<K3; k3++) {
800 #if defined(WITH_TIMESTAMPCHANNEL)
801
802 #ifdef NO_BI
803 #define _BI (BI[hdr->NS-1])
804 #else
805 #define _BI (hdr->CHANNEL[hdr->NS-1].bi)
806 #endif
807 gdf_time t = heka2gdftime(*(double*)(hdr->AS.Header+pos+48)); // time of sweep. TODO: this should be taken into account
808 *(int64_t*)(hdr->AS.rawdata + _BI + SPR * 8) = t;
809 #undef _BI
810
811 #endif // WITH_TIMESTAMPCHANNEL
812 // read sweep
813 char flagSweepSelected = (hdr->AS.SegSel[0]==0 || k1+1==hdr->AS.SegSel[0])
814 && (hdr->AS.SegSel[1]==0 || k2+1==hdr->AS.SegSel[1])
815 && (hdr->AS.SegSel[2]==0 || k3+1==hdr->AS.SegSel[2]);
816
817 if (VERBOSE_LEVEL>7) fprintf(stdout,"HEKA+L3 @%i=\t%i/%i %i/%i %i/%i sel=%i\n",(int)(pos+StartOfData),k1,K1,k2,K2,k3,K3,flagSweepSelected);
818
819 pos += Sizes.Rec.Sweep + 4;
820 // read number of children
821 K4 = (*(uint32_t*)(hdr->AS.Header+pos-4));
822 size_t DIV=1;
823 for (k4=0; k4<K4; k4++) {
824 if (!flagSweepSelected) {
825 pos += Sizes.Rec.Trace+4;
826 continue;
827 }
828
829 // read trace
830 uint16_t gdftyp = 0;
831 uint32_t ns = (*(uint32_t*)(hdr->AS.Header+pos+36));
832 uint32_t DataPos = (*(uint32_t*)(hdr->AS.Header+pos+40));
833 spr = (*(uint32_t*)(hdr->AS.Header+pos+44));
834 double DataScaler= (*(double*)(hdr->AS.Header+pos+72));
835 double Toffset = (*(double*)(hdr->AS.Header+pos+80)); // time offset of
836 uint16_t pdc = PhysDimCode((char*)(hdr->AS.Header + pos + 96));
837 char *physdim = (char*)(hdr->AS.Header + pos + 96);
838 double dT = (*(double*)(hdr->AS.Header+pos+104));
839 // double XStart = (*(double*)(hdr->AS.Header+pos+112));
840 // uint16_t XUnits = PhysDimCode((char*)(hdr->AS.Header+pos+120));
841 double YRange = (*(double*)(hdr->AS.Header+pos+128));
842 double YOffset = (*(double*)(hdr->AS.Header+pos+136));
843 // double Bandwidth = (*(double*)(hdr->AS.Header+pos+144));
844 uint16_t AdcChan = (*(uint16_t*)(hdr->AS.Header+pos+222));
845 /*
846 double PhysMin = (*(double*)(hdr->AS.Header+pos+224));
847 double PhysMax = (*(double*)(hdr->AS.Header+pos+232));
848 */
849 switch (hdr->AS.Header[pos+70]) {
850 case 0: gdftyp = 3; break; // int16
851 case 1: gdftyp = 5; break; // int32
852 case 2: gdftyp = 16; break; // float32
853 case 3: gdftyp = 17; break; // float64
854 default:
855 biosigERROR(hdr, B4C_FORMAT_UNSUPPORTED, "Heka/Patchmaster unknown data type is used");
856 };
857
858 if (SWAP) {
859 AdcChan = bswap_16(AdcChan);
860 ns = bswap_32(ns);
861 DataPos = bswap_32(DataPos);
862 spr = bswap_32(spr);
863 // avoid breaking strict-aliasing rules
864 union {
865 double f64;
866 uint64_t u64;
867 } c;
868 c.f64 = dT; c.u64 = bswap_64(c.u64); dT = c.f64;
869 c.f64 = YRange; c.u64 = bswap_64(c.u64); YRange = c.f64;
870 c.f64 = YOffset; c.u64 = bswap_64(c.u64); YOffset = c.f64;
871 /*
872 c.f64 = PhysMax; c.u64 = bswap_64(c.u64); PhysMax = c.f64;
873 c.f64 = PhysMin; c.u64 = bswap_64(c.u64); PhysMin = c.f64;
874 */
875 c.f64 = Toffset; c.u64 = bswap_64(c.u64); Toffset = c.f64;
876 }
877 double Fs = round(1.0 / dT);
878 DIV = round(hdr->SampleRate / Fs);
879
880 char *Label = (char*)(hdr->AS.Header+pos+4);
881 for (ns=0; ns < hdr->NS; ns++) {
882 if (!strcmp(hdr->CHANNEL[ns].Label, Label)) break;
883 }
884 CHANNEL_TYPE *hc = hdr->CHANNEL+ns;
885
886 if (VERBOSE_LEVEL>7) fprintf(stdout,"HEKA+L4 @%i= #%i,%i,%i/%i %s\t%i/%i %i/%i %i/%i %i/%i DIV=%i,%i,%i\n",(int)(pos+StartOfData),ns,AdcChan,spr,SPR,Label,k1,K1,k2,K2,k3,K3,k4,K4,(int)DIV,gdftyp,hc->GDFTYP);
887
888 if (itx) {
889 uint32_t k5;
890 double Cal = DataScaler;
891
892 assert(hdr->CHANNEL[ns].Off==0.0);
893 double Off = 0.0;
894
895 fprintf(itx, "\r\nWAVES %s_%i_%i_%i_%i\r\nBEGIN\r\n", WAVENAME,k1+1,k2+1,k3+1,k4+1);
896 switch (hc->GDFTYP) {
897 case 3:
898 for (k5 = 0; k5 < spr; ++k5)
899 fprintf(itx,"% e\n", (double)*(int16_t*)(hdr->AS.Header + DataPos + k5 * 2) * Cal + Off);
900 break;
901 case 5:
902 for (k5 = 0; k5 < spr; ++k5)
903 fprintf(itx,"% e\n", (double)*(int32_t*)(hdr->AS.Header + DataPos + k5 * 4) * Cal + Off);
904 break;
905 case 16:
906 for (k5 = 0; k5 < spr; ++k5)
907 fprintf(itx,"% e\n", (double)*(float*)(hdr->AS.Header + DataPos + k5 * 4) * Cal + Off);
908 break;
909 case 17:
910 for (k5 = 0; k5 < spr; ++k5)
911 fprintf(itx,"% e\n", *(double*)(hdr->AS.Header + DataPos + k5 * 8) * Cal + Off);
912 break;
913 }
914 fprintf(itx, "END\r\nX SetScale/P x, %g, %g, \"s\", %s_%i_%i_%i_%i\r\n", Toffset, dT, WAVENAME, k1+1,k2+1,k3+1,k4+1);
915 fprintf(itx, "X SetScale y,0,0,\"%s\", %s_%i_%i_%i_%i\n", physdim, WAVENAME, k1+1,k2+1,k3+1,k4+1);
916 }
917
918 #ifdef NO_BI
919 #define _BI (BI[ns])
920 #else
921 #define _BI (hc->bi)
922 #endif
923 // no need to check byte order because File.Endian is set and endian conversion is done in sread
924 if ((DIV==1) && (gdftyp == hc->GDFTYP)) {
925 uint16_t sz = GDFTYP_BITS[hc->GDFTYP]>>3;
926 memcpy(hdr->AS.rawdata + _BI + SPR * sz, hdr->AS.Header + DataPos, spr * sz);
927 }
928 else if (1) {
929 double Cal = DataScaler * PhysDimScale(pdc) / hdr->CHANNEL[ns].Cal;
930 assert(Cal==1.0 || hc->GDFTYP > 15); // when scaling changes, target data type is always float/double -> see above
931 uint32_t k5,k6;
932 switch (gdftyp) {
933 case 3:
934 switch (hc->GDFTYP) {
935 case 3:
936 for (k5 = 0; k5 < spr; ++k5) {
937 int16_t ival = *(int16_t*)(hdr->AS.Header + DataPos + k5 * 2);
938 for (k6 = 0; k6 < DIV; ++k6)
939 *(int16_t*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 2) = ival;
940 }
941 break;
942 case 5:
943 for (k5 = 0; k5 < spr; ++k5) {
944 int16_t ival = *(int16_t*)(hdr->AS.Header + DataPos + k5 * 2);
945 for (k6 = 0; k6 < DIV; ++k6)
946 *(int32_t*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 4) = (int32_t)ival;
947 }
948 break;
949 case 16:
950 for (k5 = 0; k5 < spr; ++k5) {
951 int16_t ival = *(int16_t*)(hdr->AS.Header + DataPos + k5 * 2);
952 for (k6 = 0; k6 < DIV; ++k6)
953 *(float*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 4) = (float)ival * Cal;
954 }
955 break;
956 case 17:
957 for (k5 = 0; k5 < spr; ++k5) {
958 int16_t ival = *(int16_t*)(hdr->AS.Header + DataPos + k5 * 2);
959 for (k6 = 0; k6 < DIV; ++k6)
960 *(double*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 8) = (double)ival * Cal;
961 }
962 break;
963 }
964 break;
965 case 5:
966 switch (hc->GDFTYP) {
967 case 5:
968 for (k5 = 0; k5 < spr; ++k5) {
969 int32_t ival = *(int32_t*)(hdr->AS.Header + DataPos + k5 * 4);
970 for (k6 = 0; k6 < DIV; ++k6)
971 *(int32_t*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 4) = ival;
972 }
973 break;
974 case 16:
975 for (k5 = 0; k5 < spr; ++k5) {
976 int32_t ival = *(int32_t*)(hdr->AS.Header + DataPos + k5 * 4);
977 for (k6 = 0; k6 < DIV; ++k6)
978 *(float*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 4) = (float)ival * Cal;
979 }
980 break;
981 case 17:
982 for (k5 = 0; k5 < spr; ++k5) {
983 int32_t ival = *(int32_t*)(hdr->AS.Header + DataPos + k5 * 4);
984 for (k6 = 0; k6 < DIV; ++k6)
985 *(double*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 8) = (double)ival * Cal;
986 }
987 break;
988 }
989 break;
990 case 16:
991 switch (hc->GDFTYP) {
992 case 16:
993 for (k5 = 0; k5 < spr; ++k5) {
994 float ival = *(float*)(hdr->AS.Header + DataPos + k5 * 4);
995 for (k6 = 0; k6 < DIV; ++k6)
996 *(float*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 4) = ival * Cal;
997 }
998 break;
999 case 17:
1000 for (k5 = 0; k5 < spr; ++k5) {
1001 float ival = *(float*)(hdr->AS.Header + DataPos + k5 * 4);
1002 for (k6 = 0; k6 < DIV; ++k6)
1003 *(double*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 8) = (double)ival * Cal;
1004 }
1005 break;
1006 }
1007 break;
1008 case 17:
1009 switch (hc->GDFTYP) {
1010 case 17:
1011 for (k5 = 0; k5 < spr; ++k5) {
1012 double ival = *(double*)(hdr->AS.Header + DataPos + k5 * 8);
1013 for (k6 = 0; k6 < DIV; ++k6)
1014 *(double*)(hdr->AS.rawdata + _BI + (SPR + k5*DIV + k6) * 8) = ival * Cal;
1015 }
1016 break;
1017 }
1018 break;
1019 }
1020 }
1021 #undef _BI
1022 pos += Sizes.Rec.Trace+4;
1023
1024 }
1025 if (flagSweepSelected) SPR += spr * DIV;
1026 }
1027 }
1028 }
1029 #ifdef NO_BI
1030 if (BI) free(BI);
1031 #endif
1032 hdr->AS.first = 0;
1033 hdr->AS.length = hdr->NRec;
1034 free(hdr->AS.Header);
1035 hdr->AS.Header = NULL;
1036
1037 if (VERBOSE_LEVEL>7) fprintf(stdout,"End of SOPEN_HEKA\n");
1038 }
1039
1040 else {
1041 biosigERROR(hdr, B4C_FORMAT_UNSUPPORTED, "Heka/Patchmaster format has unsupported version number");
1042 }
1043 }
1044
1045
1046 #ifdef __cplusplus
1047 }
1048 #endif
1049
1050