1 /*
2
3 Copyright (C) 2011,2013,2015 Alois Schloegl <alois.schloegl@ist.ac.at>
4 This file is part of the "BioSig for C/C++" repository
5 (biosig4c++) at http://biosig.sf.net/
6
7 This program 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 */
13
14 #include "mex.h"
15 #include <ctype.h>
16 #include <math.h>
17 #include <stdlib.h>
18 #include <string.h>
19 #include <time.h>
20
21 #include <biosig-dev.h>
22 #include <biosig.h>
23
24
25 #ifdef NDEBUG
26 #define VERBOSE_LEVEL 0 // turn off debugging information, but its only used without NDEBUG
27 #else
28 extern int VERBOSE_LEVEL; // used for debugging
29 #endif
30
31 #ifdef tmwtypes_h
32 #if (MX_API_VER<=0x07020000)
33 typedef int mwSize;
34 #endif
35 #endif
36
getDouble(const mxArray * pm,size_t idx)37 double getDouble(const mxArray *pm, size_t idx) {
38 size_t n = mxGetNumberOfElements(pm);
39 if (n == 0) return(NAN);
40 if (n <= idx) idx = n-1;
41
42 switch (mxGetClassID(pm)) {
43 case mxCHAR_CLASS:
44 case mxLOGICAL_CLASS:
45 case mxINT8_CLASS:
46 return(*((int8_t*)mxGetData(pm) + idx));
47 case mxUINT8_CLASS:
48 return(*((uint8_t*)mxGetData(pm) + idx));
49 case mxDOUBLE_CLASS:
50 return(*((double*)mxGetData(pm) + idx));
51 case mxSINGLE_CLASS:
52 return(*((float*)mxGetData(pm) + idx));
53 case mxINT16_CLASS:
54 return(*((int16_t*)mxGetData(pm) + idx));
55 case mxUINT16_CLASS:
56 return(*((uint16_t*)mxGetData(pm) + idx));
57 case mxINT32_CLASS:
58 return(*((int32_t*)mxGetData(pm) + idx));
59 case mxUINT32_CLASS:
60 return(*((uint32_t*)mxGetData(pm) + idx));
61 case mxINT64_CLASS:
62 return(*((int64_t*)mxGetData(pm) + idx));
63 case mxUINT64_CLASS:
64 return(*((uint64_t*)mxGetData(pm) + idx));
65 /*
66 case mxFUNCTION_CLASS:
67 case mxUNKNOWN_CLASS:
68 case mxCELL_CLASS:
69 case mxSTRUCT_CLASS:
70 */
71 default:
72 return(NAN);
73 }
74 return(NAN);
75 }
76
mexFunction(int nlhs,mxArray * plhs[],int nrhs,const mxArray * prhs[])77 void mexFunction(
78 int nlhs, /* number of expected outputs */
79 mxArray *plhs[], /* array of pointers to output arguments */
80 int nrhs, /* number of inputs */
81 const mxArray *prhs[] /* array of pointers to input arguments */
82 )
83
84 {
85 int k;
86 const mxArray *arg;
87 HDRTYPE *hdr;
88 size_t count;
89 time_t T0;
90 char *FileName;
91 char tmpstr[128];
92 int status;
93 int CHAN = 0;
94 double *ChanList=NULL;
95 int NS = -1;
96 char FlagOverflowDetection = 1, FlagUCAL = 0;
97 void *data = NULL;
98 mxArray *p = NULL, *p1 = NULL, *p2 = NULL;
99
100 #ifdef CHOLMOD_H
101 cholmod_sparse RR,*rr=NULL;
102 double dummy;
103 #endif
104
105 // ToDO: output single data
106 // mxClassId FlagMXclass=mxDOUBLE_CLASS;
107
108
109 if (nrhs<1) {
110 mexPrintf(" Usage of mexSSAVE:\n");
111 mexPrintf("\tstatus=mexSSAVE(HDR,data)\n");
112 mexPrintf(" Input:\n\tHDR\tHeader structure \n");
113 mexPrintf("\tdata\tdata matrix, one channel per column\n");
114 mexPrintf("\tHDR\theader structure\n\n");
115 mexPrintf("\tstatus 0: file saves successfully\n\n");
116 mexPrintf("\tstatus <>0: file could not saved\n\n");
117 return;
118 }
119
120 /*
121 improve checks for input arguments
122 */
123 /* process input arguments */
124 if (mxIsNumeric(prhs[1]) &&
125 mxIsStruct( prhs[0])) {
126 data = (void*) mxGetData(prhs[1]);
127 // get number of channels
128 size_t NS = mxGetN (prhs[1]);
129 // get number of events
130 size_t NEvt=0;
131 if ( (p = mxGetField(prhs[0], 0, "EVENT") ) != NULL ) {
132 if ( (p1 = mxGetField(p, 0, "POS") ) != NULL ) {
133 NEvt = mxGetNumberOfElements(p1);
134 }
135 if ( (p1 = mxGetField(p, 0, "TYP") ) != NULL ) {
136 size_t n = mxGetNumberOfElements(p1);
137 if (n>NEvt) NEvt = n;
138 }
139 }
140
141 // allocate memory for header structure
142 hdr = constructHDR (NS, NEvt);
143 data = (biosig_data_type*) mxGetData (prhs[1]);
144 hdr->NRec = mxGetM (prhs[1]);
145 hdr->SPR = 1;
146 }
147 else {
148 mexErrMsgTxt("mexSSAVE(HDR,data) failed because HDR and data, are not a struct and numeric, resp.\n");
149 return;
150 }
151
152
153 for (k = 2; k < nrhs; k++) {
154 arg = prhs[k];
155
156 if (mxIsChar(arg)) {
157 #ifdef DEBUG
158 mexPrintf("arg[%i]=%s \n",k,mxArrayToString(prhs[k]));
159 #endif
160 if (!strcmp(mxArrayToString(prhs[k]), "OVERFLOWDETECTION:ON"))
161 FlagOverflowDetection = 1;
162 else if (!strcmp(mxArrayToString(prhs[k]), "OVERFLOWDETECTION:OFF"))
163 FlagOverflowDetection = 0;
164 else if (!strcmp(mxArrayToString(prhs[k]), "UCAL:ON"))
165 FlagUCAL = 1;
166 else if (!strcmp(mxArrayToString(prhs[k]), "UCAL:OFF"))
167 FlagUCAL = 0;
168 }
169 else {
170 #ifndef mexSOPEN
171 mexPrintf("mexSSAVE: argument #%i is invalid.",k+1);
172 mexErrMsgTxt("mexSSAVE failes because of unknown parameter\n");
173 #else
174 mexPrintf("mexSOPEN: argument #%i is invalid.",k+1);
175 mexErrMsgTxt("mexSOPEN fails because of unknown parameter\n");
176 #endif
177 }
178 }
179
180 /***** SET INPUT ARGUMENTS *****/
181 hdr->FLAG.OVERFLOWDETECTION = FlagOverflowDetection;
182 hdr->FLAG.UCAL = FlagUCAL;
183 #ifdef CHOLMOD_H
184 hdr->FLAG.ROW_BASED_CHANNELS = (rr!=NULL);
185 #else
186 hdr->FLAG.ROW_BASED_CHANNELS = 0;
187 #endif
188
189
190 /***** CHECK INPUT HDR STUCTURE CONVERT TO libbiosig hdr *****/
191 if (VERBOSE_LEVEL>7) mexPrintf("110: input arguments checked\n");
192
193 if ( (p = mxGetField(prhs[0], 0, "TYPE") ) != NULL ) {
194 mxGetString(p,tmpstr,sizeof(tmpstr));
195 hdr->TYPE = GetFileTypeFromString(tmpstr);
196 }
197 if ( (p = mxGetField(prhs[0], 0, "VERSION") ) != NULL ) {
198 mxGetString(p, tmpstr, sizeof(tmpstr));
199 hdr->VERSION = atof(tmpstr);
200 }
201 if ( (p = mxGetField(prhs[0], 0, "T0") ) != NULL ) hdr->T0 = (gdf_time)getDouble(p, 0);
202 if ( (p = mxGetField(prhs[0], 0, "tzmin") ) != NULL )
203 hdr->tzmin = (int16_t)getDouble(p, 0);
204 else
205 hdr->tzmin = -timezone/60;
206 if ( (p = mxGetField(prhs[0], 0, "FileName") ) != NULL ) FileName = mxArrayToString(p);
207 if ( (p = mxGetField(prhs[0], 0, "SampleRate") ) != NULL ) hdr->SampleRate = getDouble(p, 0);
208 if ( (p = mxGetField(prhs[0], 0, "NS") ) != NULL ) hdr->NS = getDouble(p, 0);
209
210 #ifdef DEBUG
211 mexPrintf("mexSSAVE [400] TYPE=<%s><%s> VERSION=%f\n",tmpstr,GetFileTypeString(hdr->TYPE),hdr->VERSION);
212 #endif
213
214 p1 = mxGetField(prhs[0], 0, "SPR");
215 p2 = mxGetField(prhs[0], 0, "NRec");
216 if ( p1 && p2) {
217 hdr->SPR = (size_t)getDouble(p1, 0);
218 hdr->NRec = (size_t)getDouble(p2, 0);
219 }
220 else if (!p1 && p2) {
221 hdr->SPR = hdr->NRec;
222 hdr->NRec = (size_t)getDouble(p2, 0);
223 hdr->SPR /= hdr->NRec;
224 }
225 else if ( p1 && !p2) {
226 hdr->SPR = (size_t)getDouble(p1, 0);
227 hdr->NRec /= hdr->SPR;
228 }
229 else if (!p1 && !p2) {
230 ; /* use default values SPR=1, NREC = size(data,1) */
231 }
232
233 if (hdr->NRec * hdr->SPR != mxGetM (prhs[1]) )
234 mexPrintf("mexSSAVE: warning HDR.NRec * HDR.SPR (%i*%i = %i) does not match number of rows (%i) in data.", hdr->NRec, hdr->SPR, hdr->NRec*hdr->SPR, mxGetM(prhs[1]) );
235
236
237 if ( (p = mxGetField(prhs[0], 0, "Label") ) != NULL ) {
238 if ( mxIsCell(p) ) {
239 for (k = 0; k < hdr->NS; k++)
240 mxGetString(mxGetCell(p,k), hdr->CHANNEL[k].Label, MAX_LENGTH_LABEL+1);
241 }
242 }
243 if ( (p = mxGetField(prhs[0], 0, "Transducer") ) != NULL ) {
244 if ( mxIsCell(p) ) {
245 for (k = 0; k < hdr->NS; k++)
246 mxGetString(mxGetCell(p,k), hdr->CHANNEL[k].Transducer, MAX_LENGTH_LABEL+1);
247 }
248 }
249
250
251 if ( (p = mxGetField(prhs[0], 0, "LowPass") ) != NULL ) {
252 for (k = 0; k < hdr->NS; k++)
253 hdr->CHANNEL[k].LowPass = (float)getDouble(p,k);
254 }
255 if ( (p = mxGetField(prhs[0], 0, "HighPass") ) != NULL ) {
256 for (k = 0; k < hdr->NS; k++)
257 hdr->CHANNEL[k].HighPass = (float)getDouble(p,k);
258 }
259 if ( (p = mxGetField(prhs[0], 0, "Notch") ) != NULL ) {
260 for (k = 0; k < hdr->NS; k++)
261 hdr->CHANNEL[k].Notch = (float)getDouble(p,k);
262 }
263 if ( (p = mxGetField(prhs[0], 0, "PhysMax") ) != NULL ) {
264 for (k = 0; k < hdr->NS; k++)
265 hdr->CHANNEL[k].PhysMax = (double)getDouble(p,k);
266 }
267 if ( (p = mxGetField(prhs[0], 0, "PhysMin") ) != NULL ) {
268 for (k = 0; k < hdr->NS; k++)
269 hdr->CHANNEL[k].PhysMin = (double)getDouble(p,k);
270 }
271 if ( (p = mxGetField(prhs[0], 0, "DigMax") ) != NULL ) {
272 for (k = 0; k < hdr->NS; k++)
273 hdr->CHANNEL[k].DigMax = (double)getDouble(p,k);
274 }
275 if ( (p = mxGetField(prhs[0], 0, "DigMin") ) != NULL ) {
276 for (k = 0; k < hdr->NS; k++)
277 hdr->CHANNEL[k].DigMin = (double)getDouble(p,k);
278 }
279
280 if ( (p = mxGetField(prhs[0], 0, "PhysDimCode") ) != NULL ) {
281 for (k = 0; k < hdr->NS; k++)
282 hdr->CHANNEL[k].PhysDimCode = (uint16_t)getDouble(p,k);
283 }
284 else if ( (p = mxGetField(prhs[0], 0, "PhysDim") ) != NULL ) {
285 if ( mxIsCell(p) ) {
286 for (k = 0; k < hdr->NS; k++)
287 mxGetString(mxGetCell(p,k), tmpstr, sizeof(tmpstr));
288 hdr->CHANNEL[k].PhysDimCode = PhysDimCode(tmpstr);
289 }
290 }
291
292 if ( (p = mxGetField(prhs[0], 0, "GDFTYP") ) != NULL ) {
293 for (k = 0; k < hdr->NS; k++)
294 hdr->CHANNEL[k].GDFTYP = (uint16_t)getDouble(p,k);
295 }
296 if ( (p = mxGetField(prhs[0], 0, "TOffset") ) != NULL ) {
297 for (k = 0; k < hdr->NS; k++)
298 hdr->CHANNEL[k].TOffset = (float)getDouble(p,k);
299 }
300 if ( (p = mxGetField(prhs[0], 0, "Impedance") ) != NULL ) {
301 for (k = 0; k < hdr->NS; k++)
302 hdr->CHANNEL[k].Impedance = (float)getDouble(p,k);
303 }
304 if ( (p = mxGetField(prhs[0], 0, "fZ") ) != NULL ) {
305 for (k = 0; k < hdr->NS; k++)
306 hdr->CHANNEL[k].fZ = (float)getDouble(p,k);
307 }
308 if ( (p = mxGetField(prhs[0], 0, "AS") ) != NULL ) {
309 if ( (p1 = mxGetField(p, 0, "SPR") ) != NULL ) {
310 // define channel-based samplingRate, HDR.SampleRate*HDR.AS.SPR(channel)/HDR.SPR;
311 for (k = 0; k < hdr->NS; k++)
312 hdr->CHANNEL[k].SPR = (double)getDouble(p1,k);
313 }
314 }
315
316 if ( (p = mxGetField(prhs[0], 0, "Patient") ) != NULL ) {
317 if ( (p1 = mxGetField(p, 0, "Id") ) != NULL )
318 if (mxIsChar(p1)) mxGetString(p1, hdr->Patient.Id, MAX_LENGTH_PID+1);
319 if ( (p1 = mxGetField(p, 0, "Name") ) != NULL )
320 if (mxIsChar(p1)) mxGetString(p1, hdr->Patient.Name, MAX_LENGTH_PID+1);
321 if ( (p1 = mxGetField(p, 0, "Sex") ) != NULL ) {
322 if (mxIsChar(p1)) {
323 char sex = toupper(*mxGetChars(p1));
324 hdr->Patient.Sex = (sex=='M') + 2*(sex=='F');
325 }
326 else
327 hdr->Patient.Sex = (int8_t)getDouble(p1,0);
328 }
329
330 if ( (p1 = mxGetField(p, 0, "Handedness") ) != NULL )
331 hdr->Patient.Handedness = (int8_t)getDouble(p1,0);
332 if ( (p1 = mxGetField(p, 0, "Smoking") ) != NULL )
333 hdr->Patient.Smoking = (int8_t)getDouble(p1,0);
334 if ( (p1 = mxGetField(p, 0, "AlcoholAbuse") ) != NULL )
335 hdr->Patient.AlcoholAbuse = (int8_t)getDouble(p1,0);
336 if ( (p1 = mxGetField(p, 0, "DrugAbuse") ) != NULL )
337 hdr->Patient.DrugAbuse = (int8_t)getDouble(p1,0);
338 if ( (p1 = mxGetField(p, 0, "Medication") ) != NULL )
339 hdr->Patient.Medication = (int8_t)getDouble(p1,0);
340 if ( (p1 = mxGetField(p, 0, "Impairment") ) != NULL ) {
341 if ( (p2 = mxGetField(p1, 0, "Visual") ) != NULL )
342 hdr->Patient.Impairment.Visual = (int8_t)getDouble(p2,0);
343 if ( (p2 = mxGetField(p1, 0, "Heart") ) != NULL )
344 hdr->Patient.Impairment.Heart = (int8_t)getDouble(p2,0);
345 }
346
347 if ( (p1 = mxGetField(p, 0, "Weight") ) != NULL )
348 hdr->Patient.Weight = (uint8_t)getDouble(p1,0);
349 if ( (p1 = mxGetField(p, 0, "Height") ) != NULL )
350 hdr->Patient.Height = (uint8_t)getDouble(p1,0);
351 if ( (p1 = mxGetField(p, 0, "Birthday") ) != NULL )
352 hdr->Patient.Birthday = (gdf_time)getDouble(p1,0);
353 }
354
355 if ( (p = mxGetField(prhs[0], 0, "ID") ) != NULL ) {
356 if ( (p1 = mxGetField(p, 0, "Recording") ) != NULL )
357 if (mxIsChar(p1)) mxGetString(p1, hdr->ID.Recording, MAX_LENGTH_RID+1);
358 if ( (p1 = mxGetField(p, 0, "Technician") ) != NULL )
359 if (mxIsChar(p1)) {
360 char* str = mxArrayToString(p1);
361 hdr->ID.Technician = strdup(str);
362 }
363 if ( (p1 = mxGetField(p, 0, "Hospital") ) != NULL && mxIsChar(p1) ) {
364 size_t len = mxGetN(p1)*mxGetN(p1);
365 hdr->ID.Hospital = (char*)realloc(hdr->ID.Hospital, len+1);
366 mxGetString(p1, hdr->ID.Hospital, len);
367 }
368 if ( (p1 = mxGetField(p, 0, "Equipment") ) != NULL )
369 memcpy(&hdr->ID.Equipment,mxGetData(p1), 8);
370 if ( (p1 = mxGetField(p, 0, "Manufacturer") ) != NULL ) {
371 uint8_t pos = 0;
372 if ( ( (p2 = mxGetField(p1, 0, "Name") ) != NULL ) && mxIsChar(p2)) {
373 //hdr->ID.Manufacturer.Name=mxGetChars(p2);
374 mxGetString(p1, hdr->ID.Manufacturer._field,MAX_LENGTH_MANUF);
375 pos = strlen(hdr->ID.Manufacturer._field)+1;
376 }
377 else {
378 hdr->ID.Manufacturer._field[pos++] = 0;
379 }
380
381 if ( ( (p2 = mxGetField(p1, 0, "Model") ) != NULL ) && mxIsChar(p2)) {
382 //hdr->ID.Manufacturer.Model=mxGetChars(p2);
383 mxGetString(p1, hdr->ID.Manufacturer._field + pos, MAX_LENGTH_MANUF);
384 pos += strlen(hdr->ID.Manufacturer._field + pos)+1;
385 }
386 else {
387 hdr->ID.Manufacturer._field[pos++] = 0;
388 }
389
390 if ( ( (p2 = mxGetField(p1, 0, "Version") ) != NULL ) && mxIsChar(p2)) {
391 //hdr->ID.Manufacturer.Version=mxGetChars(p2);
392 mxGetString(p1, hdr->ID.Manufacturer._field + pos, MAX_LENGTH_MANUF);
393 pos += strlen(hdr->ID.Manufacturer._field+pos)+1;
394 }
395 else {
396 hdr->ID.Manufacturer._field[pos++] = 0;
397 }
398
399 if ( ( (p2 = mxGetField(p1, 0, "SerialNumber") ) != NULL ) && mxIsChar(p2)) {
400 //hdr->ID.Manufacturer.SerialNumber=mxGetChars(p2);
401 mxGetString(p1, hdr->ID.Manufacturer._field + pos, MAX_LENGTH_MANUF);
402 pos += strlen(hdr->ID.Manufacturer._field+pos)+1;
403 }
404 else {
405 hdr->ID.Manufacturer._field[pos++] = 0;
406 }
407 }
408 }
409
410 if ( (p = mxGetField(prhs[0], 0, "FLAG") ) != NULL ) {
411 if ( (p1 = mxGetField(p, 0, "OVERFLOWDETECTION") ) != NULL )
412 hdr->FLAG.OVERFLOWDETECTION = (char)getDouble(p1,0);
413 if ( (p1 = mxGetField(p, 0, "UCAL") ) != NULL )
414 hdr->FLAG.UCAL = (char)getDouble(p1,0);
415 if ( (p1 = mxGetField(p, 0, "ANONYMOUS") ) != NULL )
416 hdr->FLAG.ANONYMOUS = (char)getDouble(p1,0);
417 if ( (p1 = mxGetField(p, 0, "ROW_BASED_CHANNELS") ) != NULL )
418 hdr->FLAG.ROW_BASED_CHANNELS = (char)getDouble(p1,0);
419 }
420
421 if ( (p = mxGetField(prhs[0], 0, "EVENT") ) != NULL ) {
422 if ( (p1 = mxGetField(p, 0, "SampleRate") ) != NULL ) {
423 hdr->EVENT.SampleRate = (double)getDouble(p1,0);
424 }
425 if ( (p1 = mxGetField(p, 0, "POS") ) != NULL ) {
426 size_t n = mxGetNumberOfElements(p1);
427 for (k = 0; k < n; k++)
428 hdr->EVENT.POS[k] = (uint32_t)getDouble(p1,k);
429 }
430 if ( (p1 = mxGetField(p, 0, "TYP") ) != NULL ) {
431 size_t n = mxGetNumberOfElements(p1);
432 for (k = 0; k < n; k++)
433 hdr->EVENT.TYP[k] = (uint16_t)getDouble(p1,k);
434 }
435 if ( (p1 = mxGetField(p, 0, "DUR") ) != NULL ) {
436 size_t n = mxGetNumberOfElements(p1);
437 for (k = 0; k < n; k++)
438 hdr->EVENT.DUR[k] = (uint32_t)getDouble(p1,k);
439 }
440 if ( (p1 = mxGetField(p, 0, "CHN") ) != NULL ) {
441 size_t n = mxGetNumberOfElements(p1);
442 for (k = 0; k < n; k++)
443 hdr->EVENT.CHN[k] = (uint16_t)getDouble(p1,k);
444 }
445
446 if ( (p1 = mxGetField(p, 0, "CodeDesc") ) != NULL ) {
447 size_t n = mxGetNumberOfElements(p1);
448 hdr->EVENT.LenCodeDesc = n+1;
449 // get total size for storing all CodeDesc strings
450 size_t memsiz = 1;
451 for (k = 0; k < n; k++) {
452 mxArray *p2 = mxGetCell(p1,k);
453 memsiz += mxGetN(p2)+1;
454 }
455 /* allocate memory for
456 hdr->.AS.auxBUF contains the \0-terminated CodeDesc strings,
457 hdr->EVENT.CodeDesc contains the pointer to the strings
458 */
459 hdr->EVENT.CodeDesc = (const char**) realloc(hdr->EVENT.CodeDesc, hdr->EVENT.LenCodeDesc*sizeof(char*));
460 hdr->AS.auxBUF = (uint8_t*)realloc(hdr->AS.auxBUF, memsiz*sizeof(char));
461
462 // first element is always the empty string.
463 hdr->AS.auxBUF[0] = 0;
464 hdr->EVENT.CodeDesc[0] = (char*)hdr->AS.auxBUF;
465 size_t pos = 1;
466 for (k = 0; k < n; k++) {
467 mxArray *p2 = mxGetCell(p1,k);
468 size_t buflen = mxGetN(p2)+1; //*mxGetM(p2)
469 mxGetString(p2, (char*)hdr->AS.auxBUF+pos, buflen);
470 hdr->EVENT.CodeDesc[k+1] = (char*)hdr->AS.auxBUF+pos;
471 pos += buflen;
472 }
473 }
474
475 }
476
477 hdr = sopen(FileName, "w", hdr);
478 if (serror2(hdr)) mexErrMsgTxt("mexSSAVE: sopen failed \n");
479
480 swrite((biosig_data_type*)data, hdr->NRec, hdr);
481 if (serror2(hdr)) mexErrMsgTxt("mexSSAVE: swrite failed \n");
482
483 destructHDR(hdr);
484 if (serror2(hdr)) mexErrMsgTxt("mexSSAVE: sclose failed \n");
485
486 };
487
488