1 /*
2
3 Copyright (C) 2007-2019 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 <math.h>
16 #include <stdlib.h>
17 #include <string.h>
18 #include <strings.h>
19 #ifdef HAVE_CHOLMOD
20 # if !defined(__APPLE__)
21 # include <suitesparse/cholmod.h>
22 # else
23 # include <cholmod.h>
24 # endif
25 #endif
26 #include <biosig-dev.h>
27 #include <biosig.h>
28
29 #ifdef NDEBUG
30 #define VERBOSE_LEVEL 0 // turn off debugging information, but its only used without NDEBUG
31 #else
32 extern int VERBOSE_LEVEL; // used for debugging
33 #endif
34
35
36 #ifdef tmwtypes_h
37 #if (MX_API_VER<=0x07020000)
38 typedef int mwSize;
39 #endif
40 #endif
41
42 #ifndef TRUE
43 #define TRUE (1)
44 #endif
45
46 #ifdef CHOLMOD_H
47 //#include "cholmod/matlab/cholmod_matlab.h"
48 /*
49 The function sputil_get_sparse and its license was downloaded on Oct 16, 2009 from
50 http://www.cise.ufl.edu/research/sparse/cholmod/CHOLMOD/MATLAB/cholmod_matlab.c
51 http://www.cise.ufl.edu/research/sparse/cholmod/CHOLMOD/MATLAB/License.txt
52 */
53 /*
54 CHOLMOD/MATLAB Module.
55 Copyright (C) 2005-2006, Timothy A. Davis
56 CHOLMOD is also available under other licenses; contact authors for details.
57 MATLAB(tm) is a Registered Trademark of The MathWorks, Inc.
58 http://www.cise.ufl.edu/research/sparse
59
60 Note that this license is for the CHOLMOD/MATLAB module only.
61 All CHOLMOD modules are licensed separately.
62
63
64 --------------------------------------------------------------------------------
65
66
67 This Module is free software; you can redistribute it and/or
68 modify it under the terms of the GNU General Public License
69 as published by the Free Software Foundation; either version 2
70 of the License, or (at your option) any later version.
71
72 This Module is distributed in the hope that it will be useful,
73 but WITHOUT ANY WARRANTY; without even the implied warranty of
74 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
75 GNU General Public License for more details.
76
77 You should have received a copy of the GNU General Public License
78 along with this Module; if not, write to the Free Software
79 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
80 */
81
82 /* ========================================================================== */
83 /* === sputil_get_sparse ==================================================== */
84 /* ========================================================================== */
85
86 /* Create a shallow CHOLMOD copy of a MATLAB sparse matrix. No memory is
87 * allocated. The resulting matrix A must not be modified.
88 */
89
sputil_get_sparse(const mxArray * Amatlab,cholmod_sparse * A,double * dummy,mwSize stype)90 cholmod_sparse *sputil_get_sparse
91 (
92 const mxArray *Amatlab, /* MATLAB version of the matrix */
93 cholmod_sparse *A, /* CHOLMOD version of the matrix */
94 double *dummy, /* a pointer to a valid scalar double */
95 mwSize stype /* -1: lower, 0: unsymmetric, 1: upper */
96 )
97 {
98 mwSize *Ap ;
99 A->nrow = mxGetM (Amatlab) ;
100 A->ncol = mxGetN (Amatlab) ;
101 A->p = (mwSize *) mxGetJc (Amatlab) ;
102 A->i = (mwSize *) mxGetIr (Amatlab) ;
103 Ap = (mwSize*)A->p ;
104 A->nzmax = Ap [A->ncol] ;
105 A->packed = TRUE ;
106 A->sorted = TRUE ;
107 A->nz = NULL ;
108 A->itype = CHOLMOD_LONG ; /* was CHOLMOD_INT in v1.6 and earlier */
109 A->dtype = CHOLMOD_DOUBLE ;
110 A->stype = stype ;
111
112 #ifndef MATLAB6p1_OR_EARLIER
113
114 if (mxIsLogical (Amatlab))
115 {
116 A->x = NULL ;
117 A->z = NULL ;
118 A->xtype = CHOLMOD_PATTERN ;
119 }
120 else if (mxIsEmpty (Amatlab))
121 {
122 /* this is not dereferenced, but the existence (non-NULL) of these
123 * pointers is checked in CHOLMOD */
124 A->x = dummy ;
125 A->z = dummy ;
126 A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
127 }
128 else if (mxIsDouble (Amatlab))
129 {
130 A->x = mxGetPr (Amatlab) ;
131 A->z = mxGetPi (Amatlab) ;
132 A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
133 }
134 else
135 {
136 /* only logical and complex/real double matrices supported */
137 //sputil_error (ERROR_INVALID_TYPE, 0) ; // modified by AS, Oct 2009
138 }
139
140 #else
141
142 if (mxIsEmpty (Amatlab))
143 {
144 /* this is not dereferenced, but the existence (non-NULL) of these
145 * pointers is checked in CHOLMOD */
146 A->x = dummy ;
147 A->z = dummy ;
148 A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
149 }
150 else
151 {
152 /* in MATLAB 6.1, the matrix is sparse, so it must be double */
153 A->x = mxGetPr (Amatlab) ;
154 A->z = mxGetPi (Amatlab) ;
155 A->xtype = mxIsComplex (Amatlab) ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL ;
156 }
157
158 #endif
159
160 return (A) ;
161 }
162 /* ========================================================================== */
163 /* === end of sputil_get_sparse ============================================= */
164 /* ========================================================================== */
165 #endif
166
167 #ifdef WITH_PDP
168 void sopen_pdp_read(HDRTYPE *hdr);
169 #endif
170
171 //#define VERBOSE_LEVEL 9
172 //extern int VERBOSE_LEVEL;
173 //#define DEBUG
174
mexFunction(int nlhs,mxArray * plhs[],int nrhs,const mxArray * prhs[])175 void mexFunction(
176 int nlhs, /* number of expected outputs */
177 mxArray *plhs[], /* array of pointers to output arguments */
178 int nrhs, /* number of inputs */
179 const mxArray *prhs[] /* array of pointers to input arguments */
180 )
181
182 {
183 size_t k,k1;
184 const mxArray *arg;
185 mxArray *HDR;
186 HDRTYPE *hdr;
187 CHANNEL_TYPE* cp;
188 size_t count;
189 time_t T0;
190 char *FileName=NULL;
191 int status;
192 int CHAN = 0;
193 int TARGETSEGMENT = 1;
194 double *ChanList=NULL;
195 int NS = -1;
196 char FlagOverflowDetection = 1, FlagUCAL = 0;
197 int argSweepSel = -1;
198
199 #if (BIOSIG_VERSION >= 10905)
200 biosig_options_type biosig_options;
201 biosig_options.free_text_event_limiter="\0"; // default value
202 #endif
203
204 #ifdef CHOLMOD_H
205 cholmod_sparse RR,*rr=NULL;
206 double dummy;
207 #endif
208
209 // ToDO: output single data
210 // mxClassId FlagMXclass=mxDOUBLE_CLASS;
211
212
213 if (nrhs<1) {
214 #ifdef mexSOPEN
215 mexPrintf(" Usage of mexSOPEN:\n");
216 mexPrintf("\tHDR = mexSOPEN(f)\n");
217 mexPrintf(" Input:\n\tf\tfilename\n");
218 mexPrintf("\t... = mexSLOAD(f,'--free-text-event-limiter',';')\n");
219 #if (BIOSIG_VERSION >= 10905)
220 mexPrintf("\t'--free-text-event-limiter',';' : free text limited by first \";\", remainder is ignored."
221 "\n\t\tThis can help to reduce the number of distinct free text events.\n\n");
222 #endif
223 mexPrintf(" Output:\n\tHDR\theader structure\n\n");
224 #else
225 mexPrintf(" Usage of mexSLOAD:\n");
226 mexPrintf("\t[s,HDR]=mexSLOAD(f)\n");
227 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan)\n\t\tchan must be sorted in ascending order\n");
228 #ifdef CHOLMOD_H
229 mexPrintf("\t[s,HDR]=mexSLOAD(f,ReRef)\n\t\treref is a (sparse) matrix for rerefencing\n");
230 #endif
231 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan,'...')\n");
232 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan,'OVERFLOWDETECTION:ON')\n");
233 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan,'OVERFLOWDETECTION:OFF')\n");
234 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan,'UCAL:ON')\n");
235 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan,'UCAL:OFF')\n");
236 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan,'OUTPUT:SINGLE')\n");
237 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan,'TARGETSEGMENT:<N>')\n");
238 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan,'SWEEP',[NE, NG, NS])\n");
239 mexPrintf("\t[s,HDR]=mexSLOAD(f,chan,'--free-text-event-limiter',';')\n");
240 mexPrintf(" Input:\n\tf\tfilename\n");
241 mexPrintf("\tchan\tlist of selected channels; 0=all channels [default]\n");
242 mexPrintf("\tUCAL\tON: do not calibrate data; default=OFF\n");
243 // mexPrintf("\tOUTPUT\tSINGLE: single precision; default='double'\n");
244 mexPrintf("\tOVERFLOWDETECTION\tdefault = ON\n\t\tON: values outside dynamic range are not-a-number (NaN)\n");
245 mexPrintf("\tTARGETSEGMENT:<N>\n\t\tselect segment <N> in multisegment files (like Nihon-Khoden), default=1\n\t\tIt has no effect for other data formats.\n");
246 mexPrintf("\t[NE, NG, NS] are the number of the experiment, the series and the sweep, resp. for sweep selection in HEKA/PatchMaster files. (0 indicates all)\n");
247 mexPrintf("\t\t examples: [1,2,3] the 3rd sweep from the 2nd series of experiment 1; [1,3,0] selects all sweeps from experiment=1, series=3. \n\n");
248 #if (BIOSIG_VERSION >= 10905)
249 mexPrintf("\t'--free-text-event-limiter',';' : free text limited by first \";\", remainder is ignored."
250 "\n\t\tThis can help to reduce the number of distinct free text events.\n\n");
251 #endif
252 mexPrintf(" Output:\n\ts\tsignal data, each column is one channel\n");
253 mexPrintf("\tHDR\theader structure\n\n");
254 #endif
255 return;
256 }
257
258 /*
259 improve checks for input arguments
260 */
261 /* process input arguments */
262 for (k = 0; k < nrhs; k++) {
263 arg = prhs[k];
264 if (mxIsEmpty(arg) && (k>0)) {
265 #ifdef DEBUG
266 mexPrintf("arg[%i] Empty\n",k);
267 #endif
268 }
269 else if ((k==0) && mxIsCell(arg) && mxGetNumberOfElements(arg)==1 && mxGetCell(arg,0) && mxIsChar(mxGetCell(arg,0))) {
270 FileName = mxArrayToString(mxGetCell(arg,0));
271 #ifdef DEBUG
272 mexPrintf("arg[%i] IsCell\n",k);
273 #endif
274 }
275 else if ((k==0) && mxIsStruct(arg)) {
276 FileName = mxArrayToString(mxGetField(prhs[k],0,"FileName"));
277 #ifdef DEBUG
278 mexPrintf("arg[%i] IsStruct\n",k);
279 #endif
280 }
281 else if ((k==1) && mxIsSparse(arg)) {
282 #ifdef CHOLMOD_H
283 rr = sputil_get_sparse(arg,&RR,&dummy,0);
284 #else
285 mexErrMsgTxt("This version of mexSLOAD does not support re-referencing matrix - recompile with -DWITH_CHOLMOD -lcholmod \n");
286 #endif
287 }
288 else if ((k==1) && mxIsNumeric(arg)) {
289 #ifdef DEBUG
290 mexPrintf("arg[%i] IsNumeric\n",k);
291 #endif
292 ChanList = (double*)mxGetData(prhs[k]);
293 NS = mxGetNumberOfElements(prhs[k]);
294 }
295 else if (mxIsChar(arg)) {
296 #ifdef DEBUG
297 mexPrintf("arg[%i]=%s \n",k,mxArrayToString(prhs[k]));
298 #endif
299 if (k==0)
300 FileName = mxArrayToString(prhs[k]);
301 else if (!strcmp(mxArrayToString(prhs[k]), "CNT32"))
302 ; // obsolete - supported for backwards compatibility
303 else if (!strcmp(mxArrayToString(prhs[k]), "OVERFLOWDETECTION:ON"))
304 FlagOverflowDetection = 1;
305 else if (!strcmp(mxArrayToString(prhs[k]), "OVERFLOWDETECTION:OFF"))
306 FlagOverflowDetection = 0;
307 else if (!strcmp(mxArrayToString(prhs[k]), "UCAL:ON"))
308 FlagUCAL = 1;
309 else if (!strcmp(mxArrayToString(prhs[k]), "UCAL:OFF"))
310 FlagUCAL = 0;
311 // else if (!strcmp(mxArrayToString(prhs[k]),"OUTPUT:SINGLE"))
312 // FlagMXclass = mxSINGLE_CLASS;
313 else if (!strncmp(mxArrayToString(prhs[k]),"TARGETSEGMENT:",14))
314 TARGETSEGMENT = atoi(mxArrayToString(prhs[k])+14);
315 else if (!strcasecmp(mxArrayToString(prhs[k]), "SWEEP") && (prhs[k+1] != NULL) && mxIsNumeric(prhs[k+1]))
316 argSweepSel = ++k;
317 #if (BIOSIG_VERSION >= 10905)
318 else if (!strcasecmp(mxArrayToString(prhs[k]), "--free-text-event-limiter") && (prhs[k+1] != NULL) && mxIsChar(prhs[k+1]))
319 biosig_options.free_text_event_limiter = mxArrayToString(prhs[++k]);
320 #endif
321 }
322 else {
323 #ifndef mexSOPEN
324 mexPrintf("mexSLOAD: argument #%i is invalid.",k+1);
325 mexErrMsgTxt("mexSLOAD fails because of unknown parameter\n");
326 #else
327 mexPrintf("mexSOPEN: argument #%i is invalid.",k+1);
328 mexErrMsgTxt("mexSOPEN fails because of unknown parameter\n");
329 #endif
330 }
331 }
332
333 if (VERBOSE_LEVEL>7)
334 mexPrintf("110: input arguments checked\n");
335
336 hdr = constructHDR(0,0);
337
338 #ifdef __LIBBIOSIG2_H__
339
340 unsigned flags = (!!FlagOverflowDetection)*BIOSIG_FLAG_OVERFLOWDETECTION + (!!FlagUCAL)*BIOSIG_FLAG_UCAL;
341 #ifdef CHOLMOD_H
342 flags += (rr!=NULL)*BIOSIG_FLAG_ROW_BASED_CHANNELS;
343 #else
344 biosig_reset_flag(hdr, BIOSIG_FLAG_ROW_BASED_CHANNELS);
345 #endif
346 biosig_set_flag(hdr, flags);
347
348 biosig_set_targetsegment(hdr, TARGETSEGMENT);
349
350 // sweep selection for Heka format
351 if (argSweepSel>0) {
352 double *SZ = (double*) mxGetData(prhs[argSweepSel]);
353 k = 0;
354 while (k < mxGetNumberOfElements(prhs[argSweepSel]) && k < 5) {
355 biosig_set_segment_selection(hdr, k+1, (uint32_t)SZ[k]);
356 k++;
357 }
358 }
359
360 #else //__LIBBIOSIG2_H__
361
362
363 hdr->FLAG.OVERFLOWDETECTION = FlagOverflowDetection;
364 hdr->FLAG.UCAL = FlagUCAL;
365 #ifdef CHOLMOD_H
366 hdr->FLAG.ROW_BASED_CHANNELS = (rr!=NULL);
367 #else
368 hdr->FLAG.ROW_BASED_CHANNELS = 0;
369 #endif
370 hdr->FLAG.TARGETSEGMENT = TARGETSEGMENT;
371
372 // sweep selection for Heka format
373 if (argSweepSel>0) {
374 double *SZ = (double*) mxGetData(prhs[argSweepSel]);
375 k = 0;
376 while (k < mxGetNumberOfElements(prhs[argSweepSel]) && k < 5) {
377 hdr->AS.SegSel[k] = (uint32_t)SZ[k];
378 k++;
379 }
380 }
381
382
383 #endif // __LIBBIOSIG2_H__ : TODO: below, nothing is converted to level-2 interface, yet
384
385
386 if (VERBOSE_LEVEL>7)
387 mexPrintf("120: going to sopen\n");
388
389 #if (BIOSIG_VERSION >= 10905)
390 hdr = sopen_extended(FileName, "r", hdr, &biosig_options);
391 #else
392 hdr = sopen(FileName, "r", hdr);
393 #endif
394
395 /*
396 #ifdef WITH_PDP
397 if (hdr->AS.B4C_ERRNUM) {
398 hdr->AS.B4C_ERRNUM = 0;
399 sopen_pdp_read(hdr);
400 }
401 #endif
402 */
403 if (VERBOSE_LEVEL>7)
404 mexPrintf("121: sopen done\n");
405
406 if ((status=serror2(hdr))) {
407
408 const char* fields[]={"TYPE","VERSION","FileName","FLAG","ErrNum","ErrMsg"};
409 HDR = mxCreateStructMatrix(1, 1, 6, fields);
410 #ifdef __LIBBIOSIG2_H__
411 mxSetField(HDR,0,"FileName",mxCreateString(biosig_get_filename(hdr)));
412 const char *FileTypeString = GetFileTypeString(biosig_get_filetype(hdr));
413 mxSetField(HDR,0,"VERSION",mxCreateDoubleScalar(biosig_get_version(hdr)));
414 #else
415 mxSetField(HDR,0,"FileName",mxCreateString(hdr->FileName));
416 const char *FileTypeString = GetFileTypeString(hdr->TYPE);
417 mxSetField(HDR,0,"VERSION",mxCreateDoubleScalar(hdr->VERSION));
418 #endif
419 mxArray *errnum = mxCreateNumericMatrix(1,1,mxUINT8_CLASS,mxREAL);
420 *(uint8_t*)mxGetData(errnum) = (uint8_t)status;
421 mxSetField(HDR,0,"ErrNum",errnum);
422
423 #ifdef HAVE_OCTAVE
424 // handle bug in octave: mxCreateString(NULL) causes segmentation fault
425 // Octave 3.2.3 causes a seg-fault in mxCreateString(NULL)
426 if (FileTypeString) FileTypeString="\0";
427 #endif
428 mxSetField(HDR,0,"TYPE",mxCreateString(FileTypeString));
429
430
431 char *msg = (char*)malloc(72+23+strlen(FileName)); // 72: max length of constant text, 23: max length of GetFileTypeString()
432 if (msg == NULL)
433 mxSetField(HDR,0,"ErrMsg",mxCreateString("Error mexSLOAD: Cannot open file\n"));
434 else {
435 if (status==B4C_CANNOT_OPEN_FILE)
436 sprintf(msg,"Error mexSLOAD: file %s not found.\n",FileName); /* Flawfinder: ignore *** sufficient memory is allocated above */
437 else if (status==B4C_FORMAT_UNKNOWN)
438 sprintf(msg,"Error mexSLOAD: Cannot open file %s - format %s not known.\n",FileName,FileTypeString); /* Flawfinder: ignore *** sufficient memory is allocated above */
439 else if (status==B4C_FORMAT_UNSUPPORTED)
440 sprintf(msg,"Error mexSLOAD: Cannot open file %s - format %s not supported [%s].\n", FileName, FileTypeString, hdr->AS.B4C_ERRMSG); /* Flawfinder: ignore *** sufficient memory is allocated above */
441 else
442 sprintf(msg,"Error %i mexSLOAD: Cannot open file %s - format %s not supported [%s].\n", status, FileName, FileTypeString, hdr->AS.B4C_ERRMSG); /* Flawfinder: ignore *** sufficient memory is allocated above */
443
444 mxSetField(HDR,0,"ErrMsg",mxCreateString(msg));
445 free(msg);
446 }
447
448 if (VERBOSE_LEVEL>7)
449 mexPrintf("737: abort mexSLOAD - sopen failed\n");
450
451 destructHDR(hdr);
452
453 if (VERBOSE_LEVEL>7)
454 mexPrintf("757: abort mexSLOAD - sopen failed\n");
455
456 #ifdef mexSOPEN
457 plhs[0] = HDR;
458 #else
459 plhs[0] = mxCreateDoubleMatrix(0,0, mxREAL);
460 plhs[1] = HDR;
461 #endif
462 if (VERBOSE_LEVEL>7)
463 mexPrintf("777: abort mexSLOAD - sopen failed\n");
464
465 return;
466 }
467
468 #ifdef CHOLMOD_H
469 RerefCHANNEL(hdr,rr,2);
470 #endif
471
472 if (!hdr->FLAG.OVERFLOWDETECTION && FlagOverflowDetection)
473 mexPrintf("Warning %s: Overflowdetection not supported in file %s\n", __FILE__, hdr->FileName);
474 if (hdr->FLAG.UCAL != FlagUCAL)
475 mexPrintf("Warning %s: Flag UCAL is %i instead of %i (%s)\n", __FILE__, hdr->FLAG.UCAL, FlagUCAL, hdr->FileName);
476
477
478 if (VERBOSE_LEVEL>7)
479 fprintf(stderr,"[112] SOPEN-R finished NS=%i %i\n",hdr->NS,NS);
480
481 // convert2to4_eventtable(hdr);
482
483 #ifdef CHOLMOD_H
484 if (hdr->Calib != NULL) {
485 NS = hdr->Calib->ncol;
486 }
487 else
488 #endif
489 if ((NS<0) || ((NS==1) && (ChanList[0] == 0.0))) { // all channels
490 for (k=0, NS=0; k<hdr->NS; ++k) {
491 if (hdr->CHANNEL[k].OnOff) NS++;
492 }
493 }
494 else {
495 for (k=0; k<hdr->NS; ++k)
496 hdr->CHANNEL[k].OnOff = 0; // reset
497 for (k=0; k<NS; ++k) {
498 int ch = (int)ChanList[k];
499 if ((ch < 1) || (ch > hdr->NS))
500 mexPrintf("Invalid channel number CHAN(%i) = %i!\n",k+1,ch);
501 else
502 hdr->CHANNEL[ch-1].OnOff = 1; // set
503 }
504 }
505
506 if (VERBOSE_LEVEL>7)
507 fprintf(stderr,"[113] NS=%i %i\n",hdr->NS,NS);
508
509 #ifndef mexSOPEN
510 if (hdr->FLAG.ROW_BASED_CHANNELS)
511 plhs[0] = mxCreateDoubleMatrix(NS, hdr->NRec*hdr->SPR, mxREAL);
512 else
513 plhs[0] = mxCreateDoubleMatrix(hdr->NRec*hdr->SPR, NS, mxREAL);
514
515 count = sread(mxGetPr(plhs[0]), 0, hdr->NRec, hdr);
516 hdr->NRec = count;
517 #endif
518 sclose(hdr);
519 #ifdef CHOLMOD_H
520 if (hdr->Calib && hdr->rerefCHANNEL) {
521 hdr->NS = hdr->Calib->ncol;
522 free(hdr->CHANNEL);
523 hdr->CHANNEL = hdr->rerefCHANNEL;
524 hdr->rerefCHANNEL = NULL;
525 hdr->Calib = NULL;
526 }
527 #endif
528 if ((status=serror2(hdr))) return;
529
530 if (VERBOSE_LEVEL>7)
531 fprintf(stderr,"\n[129] SREAD/SCLOSE on %s successful [%i,%i] [%i,%i] %i.\n",hdr->FileName,(int)hdr->data.size[0],(int)hdr->data.size[1],(int)hdr->NRec,(int)count,(int)NS);
532
533
534 // hdr2ascii(hdr,stderr,4);
535
536 #ifndef mexSOPEN
537
538 if (nlhs>1) {
539 #endif
540
541 char* mexFileName = (char*)mxMalloc(strlen(hdr->FileName)+1);
542
543 mxArray *tmp, *tmp2, *Patient, *Manufacturer, *ID, *EVENT, *Filter, *Flag, *FileType;
544 uint16_t numfields;
545 const char *fnames[] = {"TYPE","VERSION","FileName","T0","tzmin","Patient",\
546 "HeadLen","NS","SPR","NRec","SampleRate", "FLAG", \
547 "EVENT","Label","LeadIdCode","PhysDimCode","PhysDim","Filter",\
548 "PhysMax","PhysMin","DigMax","DigMin","Transducer","Cal","Off","GDFTYP","TOffset",\
549 "ELEC","Impedance","fZ","AS","Dur","REC","Manufacturer",NULL};
550
551 for (numfields=0; fnames[numfields++] != NULL; );
552 HDR = mxCreateStructMatrix(1, 1, --numfields, fnames);
553
554 mxSetField(HDR,0,"TYPE",mxCreateString(GetFileTypeString(hdr->TYPE)));
555 mxSetField(HDR,0,"HeadLen",mxCreateDoubleScalar(hdr->HeadLen));
556 mxSetField(HDR,0,"VERSION",mxCreateDoubleScalar(hdr->VERSION));
557 mxSetField(HDR,0,"NS",mxCreateDoubleScalar(NS));
558 mxSetField(HDR,0,"SPR",mxCreateDoubleScalar(hdr->SPR));
559 mxSetField(HDR,0,"NRec",mxCreateDoubleScalar(hdr->NRec));
560 mxSetField(HDR,0,"SampleRate",mxCreateDoubleScalar(hdr->SampleRate));
561 mxSetField(HDR,0,"Dur",mxCreateDoubleScalar(hdr->SPR/hdr->SampleRate));
562 mxSetField(HDR,0,"FileName",mxCreateString(hdr->FileName));
563
564 mxSetField(HDR,0,"T0",mxCreateDoubleScalar(ldexp(hdr->T0,-32)));
565 mxSetField(HDR,0,"tzmin",mxCreateDoubleScalar(hdr->tzmin));
566
567 /* Channel information */
568 #ifdef CHOLMOD_H
569 /*
570 if (hdr->Calib == NULL) { // is refering to &RR, do not destroy
571 mxArray *Calib = mxCreateDoubleMatrix(hdr->Calib->nrow, hdr->Calib->ncol, mxREAL);
572
573 }
574 */
575 #endif
576 mxArray *LeadIdCode = mxCreateDoubleMatrix(1,NS, mxREAL);
577 mxArray *PhysDimCode = mxCreateDoubleMatrix(1,NS, mxREAL);
578 mxArray *GDFTYP = mxCreateDoubleMatrix(1,NS, mxREAL);
579 mxArray *PhysMax = mxCreateDoubleMatrix(1,NS, mxREAL);
580 mxArray *PhysMin = mxCreateDoubleMatrix(1,NS, mxREAL);
581 mxArray *DigMax = mxCreateDoubleMatrix(1,NS, mxREAL);
582 mxArray *DigMin = mxCreateDoubleMatrix(1,NS, mxREAL);
583 mxArray *Cal = mxCreateDoubleMatrix(1,NS, mxREAL);
584 mxArray *Off = mxCreateDoubleMatrix(1,NS, mxREAL);
585 mxArray *Toffset = mxCreateDoubleMatrix(1,NS, mxREAL);
586 mxArray *ELEC_POS = mxCreateDoubleMatrix(NS,3, mxREAL);
587 /*
588 mxArray *ELEC_Orient = mxCreateDoubleMatrix(NS,3, mxREAL);
589 mxArray *ELEC_Area = mxCreateDoubleMatrix(NS,1, mxREAL);
590 */
591 mxArray *LowPass = mxCreateDoubleMatrix(1,NS, mxREAL);
592 mxArray *HighPass = mxCreateDoubleMatrix(1,NS, mxREAL);
593 mxArray *Notch = mxCreateDoubleMatrix(1,NS, mxREAL);
594 mxArray *Impedance = mxCreateDoubleMatrix(1,NS, mxREAL);
595 mxArray *fZ = mxCreateDoubleMatrix(1,NS, mxREAL);
596 mxArray *SPR = mxCreateDoubleMatrix(1,NS, mxREAL);
597 mxArray *Label = mxCreateCellMatrix(NS,1);
598 mxArray *Transducer = mxCreateCellMatrix(NS,1);
599 mxArray *PhysDim1 = mxCreateCellMatrix(NS,1);
600
601 for (k=0,k1=0; k1<NS; ++k)
602 if (hdr->CHANNEL[k].OnOff) {
603 *(mxGetPr(LeadIdCode)+k1) = (double)hdr->CHANNEL[k].LeadIdCode;
604 *(mxGetPr(PhysDimCode)+k1) = (double)hdr->CHANNEL[k].PhysDimCode;
605 *(mxGetPr(GDFTYP)+k1) = (double)hdr->CHANNEL[k].GDFTYP;
606 *(mxGetPr(PhysMax)+k1) = (double)hdr->CHANNEL[k].PhysMax;
607 *(mxGetPr(PhysMin)+k1) = (double)hdr->CHANNEL[k].PhysMin;
608 *(mxGetPr(DigMax)+k1) = (double)hdr->CHANNEL[k].DigMax;
609 *(mxGetPr(DigMin)+k1) = (double)hdr->CHANNEL[k].DigMin;
610 *(mxGetPr(Toffset)+k1) = (double)hdr->CHANNEL[k].TOffset;
611 *(mxGetPr(Cal)+k1) = (double)hdr->CHANNEL[k].Cal;
612 *(mxGetPr(Off)+k1) = (double)hdr->CHANNEL[k].Off;
613 *(mxGetPr(SPR)+k1) = (double)hdr->CHANNEL[k].SPR;
614 *(mxGetPr(LowPass)+k1) = (double)hdr->CHANNEL[k].LowPass;
615 *(mxGetPr(HighPass)+k1) = (double)hdr->CHANNEL[k].HighPass;
616 *(mxGetPr(Notch)+k1) = (double)hdr->CHANNEL[k].Notch;
617 *(mxGetPr(Impedance)+k1) = (double)hdr->CHANNEL[k].Impedance;
618 *(mxGetPr(fZ)+k1) = (double)hdr->CHANNEL[k].fZ;
619 *(mxGetPr(ELEC_POS)+k1) = (double)hdr->CHANNEL[k].XYZ[0];
620 *(mxGetPr(ELEC_POS)+k1+NS) = (double)hdr->CHANNEL[k].XYZ[1];
621 *(mxGetPr(ELEC_POS)+k1+NS*2) = (double)hdr->CHANNEL[k].XYZ[2];
622 /*
623 *(mxGetPr(ELEC_Orient)+k1) = (double)hdr->CHANNEL[k].Orientation[0];
624 *(mxGetPr(ELEC_Orient)+k1+NS) = (double)hdr->CHANNEL[k].Orientation[1];
625 *(mxGetPr(ELEC_Orient)+k1+NS*2) = (double)hdr->CHANNEL[k].Orientation[2];
626 *(mxGetPr(ELEC_Area)+k1) = (double)hdr->CHANNEL[k].Area;
627 */
628 mxSetCell(Label,k1,mxCreateString(hdr->CHANNEL[k].Label ? hdr->CHANNEL[k].Label : ""));
629 mxSetCell(Transducer,k1,mxCreateString(hdr->CHANNEL[k].Transducer ? hdr->CHANNEL[k].Transducer : ""));
630
631 mxSetCell(PhysDim1,k1,mxCreateString(PhysDim3(hdr->CHANNEL[k].PhysDimCode)));
632 k1++;
633 }
634
635 mxSetField(HDR,0,"LeadIdCode",LeadIdCode);
636 mxSetField(HDR,0,"PhysDimCode",PhysDimCode);
637 mxSetField(HDR,0,"GDFTYP",GDFTYP);
638 mxSetField(HDR,0,"PhysMax",PhysMax);
639 mxSetField(HDR,0,"PhysMin",PhysMin);
640 mxSetField(HDR,0,"DigMax",DigMax);
641 mxSetField(HDR,0,"DigMin",DigMin);
642 mxSetField(HDR,0,"TOffset",Toffset);
643 mxSetField(HDR,0,"Cal",Cal);
644 mxSetField(HDR,0,"Off",Off);
645 mxSetField(HDR,0,"Impedance",Impedance);
646 mxSetField(HDR,0,"fZ",fZ);
647 mxSetField(HDR,0,"Off",Off);
648 mxSetField(HDR,0,"PhysDim",PhysDim1);
649 mxSetField(HDR,0,"Transducer",Transducer);
650 mxSetField(HDR,0,"Label",Label);
651
652 const char* field[] = {"XYZ",NULL};
653 for (numfields=0; field[numfields++] != 0; );
654 tmp = mxCreateStructMatrix(1, 1, --numfields, field);
655 mxSetField(tmp,0,"XYZ",ELEC_POS);
656 /*
657 mxSetField(tmp,0,"Orientation",ELEC_Orient);
658 mxSetField(tmp,0,"Area",ELEC_Area);
659 */
660 mxSetField(HDR,0,"ELEC",tmp);
661
662 const char* field2[] = {"SPR",NULL};
663 for (numfields=0; field2[numfields++] != 0; );
664 tmp2 = mxCreateStructMatrix(1, 1, --numfields, field2);
665 mxSetField(tmp2,0,"SPR",SPR);
666 if (hdr->AS.bci2000!=NULL) {
667 mxAddField(tmp2, "BCI2000");
668 mxSetField(tmp2,0,"BCI2000",mxCreateString(hdr->AS.bci2000));
669 }
670 if (hdr->TYPE==Sigma) {
671 mxAddField(tmp2, "H1");
672 mxSetField(tmp2,0,"H1",mxCreateString((char*)hdr->AS.Header));
673 }
674 mxSetField(HDR,0,"AS",tmp2);
675
676 /* FLAG */
677 const char* field3[] = {"UCAL","OVERFLOWDETECTION","ROW_BASED_CHANNELS",NULL};
678 for (numfields=0; field3[numfields++] != 0; );
679 Flag = mxCreateStructMatrix(1, 1, --numfields, field3);
680 #ifdef MX_API_VER
681 //#if 1
682 // Matlab, Octave 3.6.1
683 mxSetField(Flag,0,"UCAL",mxCreateLogicalScalar(hdr->FLAG.UCAL));
684 mxSetField(Flag,0,"OVERFLOWDETECTION",mxCreateLogicalScalar(hdr->FLAG.OVERFLOWDETECTION));
685 mxSetField(Flag,0,"ROW_BASED_CHANNELS",mxCreateLogicalScalar(hdr->FLAG.ROW_BASED_CHANNELS));
686 #else
687 // mxCreateLogicalScalar are not included in Octave 3.0
688 mxSetField(Flag,0,"UCAL",mxCreateDoubleScalar(hdr->FLAG.UCAL));
689 mxSetField(Flag,0,"OVERFLOWDETECTION",mxCreateDoubleScalar(hdr->FLAG.OVERFLOWDETECTION));
690 mxSetField(Flag,0,"ROW_BASED_CHANNELS",mxCreateDoubleScalar(hdr->FLAG.ROW_BASED_CHANNELS));
691 #endif
692 mxSetField(HDR,0,"FLAG",Flag);
693
694 /* Filter */
695 const char *filter_fields[] = {"HighPass","LowPass","Notch",NULL};
696 for (numfields=0; filter_fields[numfields++] != 0; );
697 Filter = mxCreateStructMatrix(1, 1, --numfields, filter_fields);
698 mxSetField(Filter,0,"LowPass",LowPass);
699 mxSetField(Filter,0,"HighPass",HighPass);
700 mxSetField(Filter,0,"Notch",Notch);
701 mxSetField(HDR,0,"Filter",Filter);
702
703 /* annotation, marker, event table */
704 const char *event_fields[] = {"SampleRate","TYP","POS","DUR","CHN","Desc",NULL};
705
706 if (hdr->EVENT.DUR == NULL)
707 EVENT = mxCreateStructMatrix(1, 1, 3, event_fields);
708 else {
709 EVENT = mxCreateStructMatrix(1, 1, 5, event_fields);
710 mxArray *DUR = mxCreateDoubleMatrix(hdr->EVENT.N,1, mxREAL);
711 mxArray *CHN = mxCreateDoubleMatrix(hdr->EVENT.N,1, mxREAL);
712 for (k=0; k<hdr->EVENT.N; ++k) {
713 *(mxGetPr(DUR)+k) = (double)hdr->EVENT.DUR[k];
714 *(mxGetPr(CHN)+k) = (double)hdr->EVENT.CHN[k]; // channels use a 1-based index, 0 indicates all channels
715 }
716 mxSetField(EVENT,0,"DUR",DUR);
717 mxSetField(EVENT,0,"CHN",CHN);
718 }
719
720 if (hdr->EVENT.CodeDesc != NULL) {
721 mxAddField(EVENT, "CodeDesc");
722 mxArray *CodeDesc = mxCreateCellMatrix(hdr->EVENT.LenCodeDesc-1,1);
723 for (k=1; k < hdr->EVENT.LenCodeDesc; ++k) {
724 mxSetCell(CodeDesc,k-1,mxCreateString(hdr->EVENT.CodeDesc[k]));
725 }
726 mxSetField(EVENT,0,"CodeDesc",CodeDesc);
727 }
728
729 mxArray *TYP = mxCreateDoubleMatrix(hdr->EVENT.N,1, mxREAL);
730 mxArray *POS = mxCreateDoubleMatrix(hdr->EVENT.N,1, mxREAL);
731
732 for (k=0; k<hdr->EVENT.N; ++k) {
733 *(mxGetPr(TYP)+k) = (double)hdr->EVENT.TYP[k];
734 *(mxGetPr(POS)+k) = (double)hdr->EVENT.POS[k]+1; // conversion from 0-based to 1-based indexing
735 }
736 mxSetField(EVENT,0,"TYP",TYP);
737 mxSetField(EVENT,0,"POS",POS);
738
739 #if (BIOSIG_VERSION >= 10500)
740 if (hdr->EVENT.TimeStamp) {
741 mxArray *TimeStamp = mxCreateDoubleMatrix(hdr->EVENT.N,1, mxREAL);
742 for (k=0; k<hdr->EVENT.N; ++k) {
743 *(mxGetPr(TimeStamp)+k) = ldexp(hdr->EVENT.TimeStamp[k],-32);
744 }
745 mxAddField(EVENT, "TimeStamp");
746 mxSetField(EVENT,0,"TimeStamp",TimeStamp);
747 }
748 #endif
749 mxSetField(EVENT,0,"SampleRate",mxCreateDoubleScalar(hdr->EVENT.SampleRate));
750 mxSetField(HDR,0,"EVENT",EVENT);
751
752 /* Record identification */
753 const char *ID_fields[] = {"Recording","Technician","Hospital","Equipment","IPaddr",NULL};
754 for (numfields=0; ID_fields[numfields++] != 0; );
755 ID = mxCreateStructMatrix(1, 1, --numfields, ID_fields);
756 mxSetField(ID,0,"Recording",mxCreateString(hdr->ID.Recording));
757 mxSetField(ID,0,"Technician",mxCreateString(hdr->ID.Technician));
758 mxSetField(ID,0,"Hospital",mxCreateString(hdr->ID.Hospital));
759 mxSetField(ID,0,"Equipment",mxCreateString((char*)&hdr->ID.Equipment));
760 int len = 4;
761 uint8_t IPv6=0;
762 for (k=4; k<16; k++) IPv6 |= hdr->IPaddr[k];
763 if (IPv6) len=16;
764 mxArray *IPaddr = mxCreateNumericMatrix(1,len,mxUINT8_CLASS,mxREAL);
765 memcpy(mxGetData(IPaddr),hdr->IPaddr,len);
766 mxSetField(ID,0,"IPaddr",IPaddr);
767 mxSetField(HDR,0,"REC",ID);
768
769 /* Patient Information */
770 const char *patient_fields[] = {"Sex","Handedness","Id","Name","Weight","Height","Birthday",NULL};
771 for (numfields=0; patient_fields[numfields++] != 0; );
772 Patient = mxCreateStructMatrix(1, 1, --numfields, patient_fields);
773 const char *strarray;
774 #ifdef __LIBBIOSIG2_H__
775 strarray = biosig_get_patient_name(hdr);
776 mxSetField(Patient,0,"Name",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
777
778 strarray = biosig_get_patient_id(hdr);
779 mxSetField(Patient,0,"Id",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
780 #else
781 strarray = hdr->Patient.Name;
782 mxSetField(Patient,0,"Name",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
783
784 strarray = hdr->Patient.Id;
785 mxSetField(Patient,0,"Id",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
786 #endif
787 mxSetField(Patient,0,"Handedness",mxCreateDoubleScalar(hdr->Patient.Handedness));
788
789 mxSetField(Patient,0,"Sex",mxCreateDoubleScalar(hdr->Patient.Sex));
790 mxSetField(Patient,0,"Weight",mxCreateDoubleScalar((double)hdr->Patient.Weight));
791 mxSetField(Patient,0,"Height",mxCreateDoubleScalar((double)hdr->Patient.Height));
792 mxSetField(Patient,0,"Birthday",mxCreateDoubleScalar(ldexp(hdr->Patient.Birthday,-32)));
793
794 double d;
795 if (hdr->Patient.Weight==0) d = NAN; // not-a-number
796 else if (hdr->Patient.Weight==255) d = INFINITY; // Overflow
797 else d = (double)hdr->Patient.Weight;
798 mxSetField(Patient,0,"Weight",mxCreateDoubleScalar(d));
799
800 if (hdr->Patient.Height==0) d = NAN; // not-a-number
801 else if (hdr->Patient.Height==255) d = INFINITY; // Overflow
802 else d = (double)hdr->Patient.Height;
803 mxSetField(Patient,0,"Height",mxCreateDoubleScalar(d));
804
805 /* Manufacturer Information */
806 const char *manufacturer_fields[] = {"Name","Model","Version","SerialNumber",NULL};
807 for (numfields=0; manufacturer_fields[numfields++] != 0; );
808 Manufacturer = mxCreateStructMatrix(1, 1, --numfields, manufacturer_fields);
809
810 #ifdef __LIBBIOSIG2_H__
811 strarray = biosig_get_manufacturer_name(hdr);
812 mxSetField(Manufacturer,0,"Name",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
813
814 strarray = biosig_get_manufacturer_model(hdr);
815 mxSetField(Manufacturer,0,"Model",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
816
817 strarray = biosig_get_manufacturer_version(hdr);
818 mxSetField(Manufacturer,0,"Version",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
819
820 strarray = biosig_get_manufacturer_serial_number(hdr);
821 mxSetField(Manufacturer,0,"SerialNumber",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
822
823 #else
824 strarray = hdr->ID.Manufacturer.Name;
825 mxSetField(Manufacturer,0,"Name",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
826
827 strarray = hdr->ID.Manufacturer.Model;
828 mxSetField(Manufacturer,0,"Model",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
829
830 strarray = hdr->ID.Manufacturer.Version;
831 mxSetField(Manufacturer,0,"Version",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
832
833 strarray = hdr->ID.Manufacturer.SerialNumber;
834 mxSetField(Manufacturer,0,"SerialNumber",mxCreateCharMatrixFromStrings(strarray!=NULL, &strarray));
835 #endif
836 mxSetField(HDR,0,"Manufacturer",Manufacturer);
837
838
839 if (VERBOSE_LEVEL>7)
840 fprintf(stdout,"[148] going for SCLOSE\n");
841
842 mxSetField(HDR,0,"Patient",Patient);
843
844 #ifndef mexSOPEN
845 plhs[1] = HDR;
846 }
847 #else
848 plhs[0] = HDR;
849 #endif
850
851 if (VERBOSE_LEVEL>7) fprintf(stdout,"[151] going for SCLOSE\n");
852 #ifdef CHOLMOD_H
853 hdr->Calib = NULL; // is refering to &RR, do not destroy
854 #endif
855 if (VERBOSE_LEVEL>7) fprintf(stdout,"[156] SCLOSE finished\n");
856 destructHDR(hdr);
857 hdr = NULL;
858 if (VERBOSE_LEVEL>7) fprintf(stdout,"[157] SCLOSE finished\n");
859 };
860
861