MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_raw_data.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "mne_raw_data.h"
42 
43 #include <QFile>
44 
45 #include <Eigen/Core>
46 
47 #define _USE_MATH_DEFINES
48 #include <math.h>
49 
50 //=============================================================================================================
51 // USED NAMESPACES
52 //=============================================================================================================
53 
54 using namespace Eigen;
55 using namespace FIFFLIB;
56 using namespace MNELIB;
57 
58 #ifndef TRUE
59 #define TRUE 1
60 #endif
61 
62 #ifndef FALSE
63 #define FALSE 0
64 #endif
65 
66 #ifndef FAIL
67 #define FAIL -1
68 #endif
69 
70 #ifndef OK
71 #define OK 0
72 #endif
73 
74 #define MALLOC_36(x,t) (t *)malloc((x)*sizeof(t))
75 #define REALLOC_36(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
76 
77 #define ALLOC_CMATRIX_36(x,y) mne_cmatrix_36((x),(y))
78 
79 #if defined(_WIN32) || defined(_WIN64)
80 #define snprintf _snprintf
81 #define vsnprintf _vsnprintf
82 #define strcasecmp _stricmp
83 #define strncasecmp _strnicmp
84 #endif
85 
86 static void matrix_error(int kind, int nr, int nc)
87 
88 {
89  if (kind == 1)
90  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
91  else if (kind == 2)
92  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
93  else
94  printf("Allocation error for a %d x %d matrix\n",nr,nc);
95  if (sizeof(void *) == 4) {
96  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
97  printf("Please consider moving to a 64-bit platform.");
98  }
99  printf("Cannot continue. Sorry.\n");
100  exit(1);
101 }
102 
103 float **mne_cmatrix_36(int nr,int nc)
104 
105 {
106  int i;
107  float **m;
108  float *whole;
109 
110  m = MALLOC_36(nr,float *);
111  if (!m) matrix_error(1,nr,nc);
112  whole = MALLOC_36(nr*nc,float);
113  if (!whole) matrix_error(2,nr,nc);
114 
115  for(i=0;i<nr;i++)
116  m[i] = whole + i*nc;
117  return m;
118 }
119 
120 #define FREE_36(x) if ((char *)(x) != NULL) free((char *)(x))
121 
122 #define FREE_CMATRIX_36(m) mne_free_cmatrix_36((m))
123 
124 void mne_free_cmatrix_36 (float **m)
125 {
126  if (m) {
127  FREE_36(*m);
128  FREE_36(m);
129  }
130 }
131 
132 namespace MNELIB
133 {
134 
135 typedef struct {
136  int size; /* Size of this buffer in floats */
137  float *data; /* The allocated buffer */
138  float ***datap; /* The matrix which uses this */
140 
141 typedef struct {
142  ringBufBuf_36 *bufs;
143  int nbuf;
144  int next;
146 
147 }
148 
149 void mne_free_ring_buffer_36(void *thisp)
150 
151 {
152  int k;
153  ringBuf_36 this_buf = (ringBuf_36)thisp;
154 
155  if (!this_buf)
156  return;
157 
158  for (k = 0; k < this_buf->nbuf; k++)
159  FREE_36(this_buf->bufs[k]->data);
160  FREE_36(this_buf->bufs);
161  FREE_36(this_buf);
162  return;
163 }
164 
165 void mne_free_event(mneEvent e)
166 {
167  if (!e)
168  return;
169  FREE_36(e->comment);
170  FREE_36(e);
171  return;
172 }
173 
174 void mne_free_event_list_36(mneEventList list)
175 
176 {
177  int k;
178  if (!list)
179  return;
180  for (k = 0; k < list->nevent; k++)
181  mne_free_event(list->events[k]);
182  FREE_36(list->events);
183  FREE_36(list);
184  return;
185 }
186 
187 void mne_free_name_list_36(char **list, int nlist)
188 /*
189  * Free a name list array
190  */
191 {
192  int k;
193  if (list == NULL || nlist == 0)
194  return;
195  for (k = 0; k < nlist; k++) {
196 #ifdef FOO
197  printf("%d %s\n",k,list[k]);
198 #endif
199  FREE_36(list[k]);
200  }
201  FREE_36(list);
202  return;
203 }
204 
205 //============================= misc_util.c =============================
206 
207 void mne_string_to_name_list(const QString& s, QStringList& listp,int &nlistp)
208 /*
209  * Convert a colon-separated list into a string array
210  */
211 {
212  QStringList list;
213 
214  if (!s.isEmpty() && s.size() > 0) {
216  //list = s.split(":");
217  }
218  listp = list;
219  nlistp = list.size();
220  return;
221 }
222 
223 QString mne_name_list_to_string(const QStringList& list)
224 /*
225  * Convert a string array to a colon-separated string
226  */
227 {
228  int nlist = list.size();
229  QString res;
230  if (nlist == 0 || list.isEmpty())
231  return res;
232 // res[0] = '\0';
233  for (int k = 0; k < nlist-1; k++) {
234  res += list[k];
235  res += ":";
236  }
237  res += list[nlist-1];
238  return res;
239 }
240 
241 QString mne_channel_names_to_string(const QList<FIFFLIB::FiffChInfo>& chs,
242  int nch)
243 /*
244  * Make a colon-separated string out of channel names
245  */
246 {
247  QStringList names;
248  QString res;
249  if (nch <= 0)
250  return res;
251  for (int k = 0; k < nch; k++)
252  names.append(chs[k].ch_name);
253  res = mne_name_list_to_string(names);
254  return res;
255 }
256 
257 void mne_channel_names_to_name_list(const QList<FIFFLIB::FiffChInfo>& chs,
258  int nch,
259  QStringList& listp,
260  int &nlistp)
261 
262 {
263  QString s = mne_channel_names_to_string(chs,nch);
264  mne_string_to_name_list(s,listp,nlistp);
265  return;
266 }
267 
268 //============================= mne_apply_filter.c =============================
269 
270 namespace MNELIB
271 {
272 
273 typedef struct {
274  float *freq_resp; /* Frequency response */
275  float *eog_freq_resp; /* Frequency response (EOG) */
276  float *precalc; /* Precalculated data for FFT */
277  int np; /* Length */
278  float nprec;
280 
281 }
282 
283 static void filter_data_free(void *datap)
284 
285 {
286  filterData data = (filterData)datap;
287  if (!data)
288  return;
289  FREE_36(data->freq_resp);
290  FREE_36(data->eog_freq_resp);
291  FREE_36(data->precalc);
292  FREE_36(data);
293  return;
294 }
295 
296 static filterData new_filter_data()
297 
298 {
299  filterData data = MALLOC_36(1,filterDataRec);
300 
301  data->freq_resp = NULL;
302  data->eog_freq_resp = NULL;
303  data->precalc = NULL;
304  data->np = 0;
305  return data;
306 }
307 
308 int mne_compare_filters(mneFilterDef f1,
309  mneFilterDef f2)
310 /*
311  * Return 0 if the two filter definitions are same, 1 otherwise
312  */
313 {
314  if (f1->filter_on != f2->filter_on ||
315  std::fabs(f1->lowpass - f2->lowpass) > 0.1 ||
316  std::fabs(f1->lowpass_width - f2->lowpass_width) > 0.1 ||
317  std::fabs(f1->highpass - f2->highpass) > 0.1 ||
318  std::fabs(f1->highpass_width - f2->highpass_width) > 0.1 ||
319  std::fabs(f1->eog_lowpass - f2->eog_lowpass) > 0.1 ||
320  std::fabs(f1->eog_lowpass_width - f2->eog_lowpass_width) > 0.1 ||
321  std::fabs(f1->eog_highpass - f2->eog_highpass) > 0.1 ||
322  std::fabs(f1->eog_highpass_width - f2->eog_highpass_width) > 0.1)
323  return 1;
324  else
325  return 0;
326 }
327 
328 //============================= mne_fft.c =============================
329 
330 void mne_fft_ana(float *data,int np, float **precalcp)
331 /*
332  * FFT analysis for real data
333  */
334 {
335  float *precalc;
336 
337  printf("##################### DEBUG Error: FFT analysis needs to be implemented");
338 
339  // if (precalcp && *precalcp)
340  // precalc = *precalcp;
341  // else {
342  // precalc = MALLOC(2*np+15,float);
343  // rffti(&np,precalc);
344  // if (precalcp)
345  // *precalcp = precalc;
346  // }
347  // rfftf(&np,data,precalc);
348 // if (!precalcp)
349 // FREE_36(precalc);
350  return;
351 }
352 
353 void mne_fft_syn(float *data,int np, float **precalcp)
354 /*
355  * FFT synthesis for real data
356  */
357 {
358  float *precalc;
359  float mult;
360 
361  printf("##################### DEBUG Error: FFT synthesis needs to be implemented");
362 
363  // if (precalcp && *precalcp)
364  // precalc = *precalcp;
365  // else {
366  // precalc = MALLOC(2*np+15,float);
367  // rffti(&np,precalc);
368  // if (precalcp)
369  // *precalcp = precalc;
370  // }
371  // rfftb(&np,data,precalc);
372  // /*
373  // * Normalization
374  // */
375  // mult = 1.0/np;
376  // mne_scale_vector(mult,data,np);
377 
378 // if (!precalcp)
379 // FREE_36(precalc);
380  return;
381 }
382 
383 int mne_apply_filter(mneFilterDef filter, void *datap, float *data, int ns, int zero_pad, float dc_offset, int kind)
384 /*
385  * Do the magick trick
386  */
387 {
388  int k,p,n;
389  filterData d = (filterData)datap;
390  float *freq_resp;
391 
392  if (ns != filter->size + 2*filter->taper_size) {
393  printf("Incorrect data length in apply_filter");
394  return FAIL;
395  }
396  /*
397  * Zero padding
398  */
399  if (zero_pad) {
400  for (k = 0; k < filter->taper_size; k++)
401  data[k] = 0.0;
402  for (k = filter->taper_size + filter->size; k < ns; k++)
403  data[k] = 0.0;
404  }
405  if (!filter->filter_on) /* Nothing else to do */
406  return OK;
407  /*
408  * Make things nice by compensating for the dc offset
409  */
410  if (dc_offset != 0.0) {
411  for (k = filter->taper_size; k < filter->taper_size + filter->size; k++)
412  data[k] = data[k] - dc_offset;
413  }
414  if (!d)
415  return OK;
416  if (!d->freq_resp)
417  return OK;
418  /*
419  * Next comes the FFT
420  */
421  mne_fft_ana(data,ns,&d->precalc);
422  /*
423  * Multiply with the frequency response
424  * See FFTpack doc for details of the arrangement
425  */
426  n = ns % 2 == 0 ? ns/2 : (ns+1)/2;
427  p = 0;
428  /*
429  * No imaginary part for the DC component
430  */
431  if (kind == FIFFV_EOG_CH)
432  freq_resp = d->eog_freq_resp;
433  else
434  freq_resp = d->freq_resp;
435  data[p] = data[p]*freq_resp[0]; p++;
436  /*
437  * The other components
438  */
439  for (k = 1 ; k < n ; k++) {
440  data[p] = data[p]*freq_resp[k]; p++;
441  data[p] = data[p]*freq_resp[k]; p++;
442  }
443  /*
444  * Then the last value
445  */
446  if (ns % 2 == 0)
447  data[p] = data[p]*freq_resp[k];
448 
449  mne_fft_syn(data,ns,&d->precalc);
450 
451  return OK;
452 }
453 
454 void mne_create_filter_response(mneFilterDef filter,
455  float sfreq,
456  void **filter_datap,
457  mneUserFreeFunc *filter_data_freep,
458  int *highpass_effective)
459 /*
460  * Create a frequency response and return also the function to free it
461  */
462 {
463  int resp_size;
464  int k,s,w,f;
465  int highpasss,lowpasss;
466  int highpass_widths,lowpass_widths;
467  float lowpass,highpass,lowpass_width,highpass_width;
468  float *freq_resp;
469  float pi4 = M_PI/4.0;
470  float mult,add,c;
471  filterData filter_data;
472 
473  resp_size = (filter->size + 2*filter->taper_size)/2 + 1;
474 
475  filter_data = new_filter_data();
476  filter_data->freq_resp = MALLOC_36(resp_size,float);
477  filter_data->eog_freq_resp = MALLOC_36(resp_size,float);
478  filter_data->np = resp_size;
479 
480  for (k = 0; k < resp_size; k++) {
481  filter_data->freq_resp[k] = 1.0;
482  filter_data->eog_freq_resp[k] = 1.0;
483  }
484  *highpass_effective = FALSE;
485 
486  for (f = 0; f < 2; f++) {
487  highpass = f == 0 ? filter->highpass : filter->eog_highpass;
488  highpass_width = f == 0 ? filter->highpass_width : filter->eog_highpass_width;
489  lowpass = f == 0 ? filter->lowpass : filter->eog_lowpass;
490  lowpass_width = f == 0 ? filter->lowpass_width : filter->eog_lowpass_width;
491  freq_resp = f == 0 ? filter_data->freq_resp : filter_data->eog_freq_resp;
492  /*
493  * Start simple first
494  */
495  highpasss = ((resp_size-1)*highpass)/(0.5*sfreq);
496  lowpasss = ((resp_size-1)*lowpass)/(0.5*sfreq);
497 
498  lowpass_widths = ((resp_size-1)*lowpass_width)/(0.5*sfreq);
499  lowpass_widths = (lowpass_widths+1)/2; /* What user specified */
500 
501  if (filter->highpass_width > 0.0) {
502  highpass_widths = ((resp_size-1)*highpass_width)/(0.5*sfreq);
503  highpass_widths = (highpass_widths+1)/2; /* What user specified */
504  }
505  else
506  highpass_widths = 3; /* Minimal */
507 
508  if (filter->filter_on) {
509  printf("filter : %7.3f ... %6.1f Hz bins : %d ... %d of %d hpw : %d lpw : %d\n",
510  highpass,
511  lowpass,
512  highpasss,
513  lowpasss,
514  resp_size,
515  highpass_widths,
516  lowpass_widths);
517  }
518  if (highpasss > highpass_widths + 1) {
519  w = highpass_widths;
520  mult = 1.0/w;
521  add = 3.0;
522  for (k = 0; k < highpasss-w+1; k++)
523  freq_resp[k] = 0.0;
524  for (k = -w+1, s = highpasss-w+1; k < w; k++, s++) {
525  if (s >= 0 && s < resp_size) {
526  c = cos(pi4*(k*mult+add));
527  freq_resp[s] = freq_resp[s]*c*c;
528  }
529  }
530  *highpass_effective = TRUE;
531  }
532  else
533  *highpass_effective = *highpass_effective || (filter->highpass == 0.0);
534 
535  if (lowpass_widths > 0) {
536  w = lowpass_widths;
537  mult = 1.0/w;
538  add = 1.0;
539  for (k = -w+1, s = lowpasss-w+1; k < w; k++, s++) {
540  if (s >= 0 && s < resp_size) {
541  c = cos(pi4*(k*mult+add));
542  freq_resp[s] = freq_resp[s]*c*c;
543  }
544  }
545  for (k = s; k < resp_size; k++)
546  freq_resp[k] = 0.0;
547  }
548  else {
549  for (k = lowpasss; k < resp_size; k++)
550  freq_resp[k] = 0.0;
551  }
552  if (filter->filter_on) {
553  if (*highpass_effective)
554  printf("Highpass filter will work as specified.\n");
555  else
556  printf("NOTE: Highpass filter omitted due to a too low corner frequency.\n");
557  }
558  else
559  printf("NOTE: Filter is presently switched off.\n");
560  }
561  *filter_datap = filter_data;
562  *filter_data_freep = filter_data_free;
563  return;
564 }
565 
566 //============================= mne_ringbuffer.c =============================
567 
568 namespace MNELIB
569 {
570 
571 typedef struct {
572  int size; /* Size of this buffer in floats */
573  float *data; /* The allocated buffer */
574  float ***datap; /* The matrix which uses this */
576 
577 typedef struct {
578  ringBufBuf *bufs;
579  int nbuf;
580  int next;
582 
583 }
584 
585 void mne_free_ring_buffer(void *thisp)
586 
587 {
588  int k;
589  ringBuf this_buf = (ringBuf)thisp;
590 
591  if (!this_buf)
592  return;
593 
594  for (k = 0; k < this_buf->nbuf; k++)
595  FREE_36(this_buf->bufs[k]->data);
596  FREE_36(this_buf->bufs);
597  FREE_36(this_buf);
598  return;
599 }
600 
601 void *mne_initialize_ring(int nbuf)
602 
603 {
604  int k;
605  ringBuf ring;
606 
607  ring = MALLOC_36(1,ringBufRec);
608  ring->bufs = MALLOC_36(nbuf,ringBufBuf);
609  ring->nbuf = nbuf;
610 
611  for (k = 0; k < nbuf; k++) {
612  ring->bufs[k] = MALLOC_36(1,ringBufBufRec);
613  ring->bufs[k]->size = 0;
614  ring->bufs[k]->data = NULL;
615  ring->bufs[k]->datap = NULL;
616  }
617  ring->next = 0;
618 
619 #ifdef DEBUG
620  printf("Ring buffer structure with %d entries initialized\n",ring->nbuf);
621 #endif
622 
623  return ring;
624 }
625 
626 void mne_allocate_from_ring(void *ringp, int nrow, int ncol, float ***res)
627 /*
628  * Get a new buffer
629  */
630 {
631  float **mat;
632  int j;
633  ringBufBuf buf;
634  ringBuf ring = (ringBuf)ringp;
635 
636  if (ring->next > ring->nbuf-1)
637  ring->next = 0;
638 
639 #ifdef DEBUG
640  printf("Allocating buf # %d\n",ring->next);
641 #endif
642 
643  buf = ring->bufs[ring->next++];
644 
645  if (buf->datap) { /* Clear the reference */
646  FREE_36(*buf->datap);
647  *buf->datap = NULL;
648  }
649  *res = mat = MALLOC_36(nrow,float *);
650  if (buf->size < nrow*ncol)
651  buf->data = REALLOC_36(buf->data,nrow*ncol,float);
652 
653  for (j = 0; j < nrow; j++)
654  mat[j] = buf->data + j*ncol;
655 
656  buf->datap = res;
657 
658  return;
659 }
660 
661 //============================= mne_raw_routines.c =============================
662 
663 int mne_read_raw_buffer_t(//fiffFile in, /* Input file */
664  FiffStream::SPtr& stream,
665  const FiffDirEntry::SPtr& ent, /* The directory entry to read */
666  float **data, /* Allocated for npick x nsamp samples */
667  int nchan, /* Number of channels in the data */
668  int nsamp, /* Expected number of samples */
669  const QList<FIFFLIB::FiffChInfo>& chs, /* Channel info for ALL channels */
670  int *pickno, /* Which channels to pick */
671  int npick) /* How many */
672 
673 {
674  FiffTag::SPtr t_pTag;
675 // fiffTagRec tag;
676  fiff_short_t *this_samples;
677  fiff_float_t *this_samplef;
678  fiff_int_t *this_sample;
679 
680  int s,c;
681  int do_all;
682  float *mult;
683 
684 // tag.data = NULL;
685 
686  if (npick == 0) {
687  pickno = MALLOC_36(nchan, int);
688  for (c = 0; c < nchan; c++)
689  pickno[c] = c;
690  do_all = TRUE;
691  npick = nchan;
692  }
693  else
694  do_all = FALSE;
695 
696  mult = MALLOC_36(npick,float);
697  for (c = 0; c < npick; c++)
698  mult[c] = chs[pickno[c]].cal*chs[pickno[c]].range;
699 
700 // if (fiff_read_this_tag(in->fd,ent->pos,&tag) == FIFF_FAIL)
701 // goto bad;
702  if (!stream->read_tag(t_pTag,ent->pos))
703  goto bad;
704 
705  if (ent->type == FIFFT_FLOAT) {
706  if ((int)(t_pTag->size()/(sizeof(fiff_float_t)*nchan)) != nsamp) {
707  printf("Incorrect number of samples in buffer.");
708  goto bad;
709  }
710  qDebug() << "ToDo: Check whether this_samplef contains the right stuff!!! - use VectorXf instead";
711  this_samplef = t_pTag->toFloat();
712  for (s = 0; s < nsamp; s++, this_samplef += nchan) {
713  for (c = 0; c < npick; c++)
714  data[c][s] = mult[c]*this_samplef[pickno[c]];
715  }
716  }
717  else if (ent->type == FIFFT_SHORT || ent->type == FIFFT_DAU_PACK16) {
718  if ((int)(t_pTag->size()/(sizeof(fiff_short_t)*nchan)) != nsamp) {
719  printf("Incorrect number of samples in buffer.");
720  goto bad;
721  }
722  qDebug() << "ToDo: Check whether this_samples contains the right stuff!!! - use VectorXi instead";
723  this_samples = (fiff_short_t *)t_pTag->data();
724  for (s = 0; s < nsamp; s++, this_samples += nchan) {
725  for (c = 0; c < npick; c++)
726  data[c][s] = mult[c]*this_samples[pickno[c]];
727  }
728  }
729  else if (ent->type == FIFFT_INT) {
730  if ((int)(t_pTag->size()/(sizeof(fiff_int_t)*nchan)) != nsamp) {
731  printf("Incorrect number of samples in buffer.");
732  goto bad;
733  }
734  qDebug() << "ToDo: Check whether this_sample contains the right stuff!!! - use VectorXi instead";
735  this_sample = t_pTag->toInt();
736  for (s = 0; s < nsamp; s++, this_sample += nchan) {
737  for (c = 0; c < npick; c++)
738  data[c][s] = mult[c]*this_sample[pickno[c]];
739  }
740  }
741  else {
742  printf("We are not prepared to handle raw data type: %d",ent->type);
743  goto bad;
744  }
745  if (do_all)
746  FREE_36(pickno);
747  FREE_36(mult);
748 // FREE_36(tag.data);
749  return OK;
750 
751 bad : {
752  if (do_all)
753  FREE_36(pickno);
754 // FREE_36(tag.data);
755  return FAIL;
756  }
757 }
758 
759 //============================= mne_process_bads.c =============================
760 
761 int mne_read_bad_channel_list_from_node(FiffStream::SPtr& stream,
762  const FiffDirNode::SPtr& pNode, QStringList& listp, int& nlistp)
763 {
764  FiffDirNode::SPtr node,bad;
765  QList<FiffDirNode::SPtr> temp;
766  QStringList list;
767  int nlist = 0;
768  FiffTag::SPtr t_pTag;
769  QString names;
770 
771  if (pNode->isEmpty())
772  node = stream->dirtree();
773  else
774  node = pNode;
775 
776  temp = node->dir_tree_find(FIFFB_MNE_BAD_CHANNELS);
777  if (temp.size() > 0) {
778  bad = temp[0];
779 
780  bad->find_tag(stream, FIFF_MNE_CH_NAME_LIST, t_pTag);
781  if (t_pTag) {
782  names = t_pTag->toString();
783  mne_string_to_name_list(names,list,nlist);
784  }
785  }
786  listp = list;
787  nlistp = nlist;
788  return OK;
789 }
790 
791 int mne_read_bad_channel_list(const QString& name, QStringList& listp, int& nlistp)
792 
793 {
794  QFile file(name);
795  FiffStream::SPtr stream(new FiffStream(&file));
796 
797  int res;
798 
799  if(!stream->open())
800  return FAIL;
801 
802  res = mne_read_bad_channel_list_from_node(stream,stream->dirtree(),listp,nlistp);
803 
804  stream->close();
805 
806  return res;
807 }
808 
809 int mne_sparse_vec_mult2(FiffSparseMatrix* mat, /* The sparse matrix */
810  float *vector, /* Vector to be multiplied */
811  float *res) /* Result of the multiplication */
812 /*
813  * Multiply a vector by a sparse matrix.
814  */
815 {
816  int i,j;
817 
818  if (mat->coding == FIFFTS_MC_RCS) {
819  for (i = 0; i < mat->m; i++) {
820  res[i] = 0.0;
821  for (j = mat->ptrs[i]; j < mat->ptrs[i+1]; j++)
822  res[i] += mat->data[j]*vector[mat->inds[j]];
823  }
824  return 0;
825  }
826  else if (mat->coding == FIFFTS_MC_CCS) {
827  for (i = 0; i < mat->m; i++)
828  res[i] = 0.0;
829  for (i = 0; i < mat->n; i++)
830  for (j = mat->ptrs[i]; j < mat->ptrs[i+1]; j++)
831  res[mat->inds[j]] += mat->data[j]*vector[i];
832  return 0;
833  }
834  else {
835  printf("mne_sparse_vec_mult2: unknown sparse matrix storage type: %d",mat->coding);
836  return -1;
837  }
838 }
839 
840 int mne_sparse_mat_mult2(FiffSparseMatrix* mat, /* The sparse matrix */
841  float **mult, /* Matrix to be multiplied */
842  int ncol, /* How many columns in the above */
843  float **res) /* Result of the multiplication */
844 /*
845  * Multiply a dense matrix by a sparse matrix.
846  */
847 {
848  int i,j,k;
849  float val;
850 
851  if (mat->coding == FIFFTS_MC_RCS) {
852  for (i = 0; i < mat->m; i++) {
853  for (k = 0; k < ncol; k++) {
854  val = 0.0;
855  for (j = mat->ptrs[i]; j < mat->ptrs[i+1]; j++)
856  val += mat->data[j]*mult[mat->inds[j]][k];
857  res[i][k] = val;
858  }
859  }
860  }
861  else if (mat->coding == FIFFTS_MC_CCS) {
862  for (k = 0; k < ncol; k++) {
863  for (i = 0; i < mat->m; i++)
864  res[i][k] = 0.0;
865  for (i = 0; i < mat->n; i++)
866  for (j = mat->ptrs[i]; j < mat->ptrs[i+1]; j++)
867  res[mat->inds[j]][k] += mat->data[j]*mult[i][k];
868  }
869  }
870  else {
871  printf("mne_sparse_mat_mult2: unknown sparse matrix storage type: %d",mat->coding);
872  return -1;
873  }
874  return 0;
875 }
876 
877 #define APPROX_RING_BUF_SIZE (600*1024*1024)
878 
879 static int approx_ring_buf_size = APPROX_RING_BUF_SIZE;
880 
881 //=============================================================================================================
882 // DEFINE MEMBER METHODS
883 //=============================================================================================================
884 
885 MneRawData::MneRawData()
886 :info(NULL)
887 ,bufs(NULL)
888 ,nbuf(0)
889 ,proj(NULL)
890 ,bad(NULL)
891 ,nbad(0)
892 ,first_samp(0)
893 ,omit_samp(0)
894 ,omit_samp_old(0)
895 ,event_list(NULL)
896 ,max_event(0)
897 ,dig_trigger_mask(0)
898 ,ring(NULL)
899 ,filt_ring(NULL)
900 ,filt_bufs(NULL)
901 ,nfilt_buf(0)
902 ,first_sample_val(NULL)
903 ,filter(NULL)
904 ,filter_data(NULL)
905 ,filter_data_free(NULL)
906 ,offsets(NULL)
907 ,deriv(NULL)
908 ,deriv_matched(NULL)
909 ,deriv_offsets(NULL)
910 ,user(NULL)
911 ,user_free(NULL)
912 ,comp(NULL)
913 ,comp_file(MNE_CTFV_NOGRAD)
914 ,comp_now(MNE_CTFV_NOGRAD)
915 ,sss(NULL)
916 {
917 }
918 
919 //=============================================================================================================
920 
922 {
923 // fiff_close(this->file);
924  this->stream->close();
925  this->filename.clear();
926  this->ch_names.clear();
927 
928  MneRawBufDef::free_bufs(this->bufs,this->nbuf);
929  mne_free_ring_buffer_36(this->ring);
930 
931  MneRawBufDef::free_bufs(this->filt_bufs,this->nfilt_buf);
932  mne_free_ring_buffer_36(this->filt_ring);
933 
934  if(this->proj)
935  delete this->proj;
936  this->badlist.clear();
937  FREE_36(this->first_sample_val);
938  FREE_36(this->bad);
939  FREE_36(this->offsets);
940  if(this->comp)
941  delete this->comp;
942  if(this->sss)
943  delete this->sss;
944 
945  if (this->filter_data_free)
946  this->filter_data_free(this->filter_data);
947  if (this->user_free)
948  this->user_free(this->user);
949  this->dig_trigger.clear();
950  mne_free_event_list_36(this->event_list);
951 
952  if(this->info)
953  delete this->info;
954 
955  delete this->deriv;
956  delete this->deriv_matched;
957  FREE_36(this->deriv_offsets);
958 }
959 
960 //=============================================================================================================
961 
962 void MneRawData::mne_raw_add_filter_response(MneRawData *data, int *highpass_effective)
963 /*
964  * Add the standard filter frequency response function
965  */
966 {
967  if (!data)
968  return;
969  /*
970  * Free the previous filter definition
971  */
972  if (data->filter_data_free)
973  data->filter_data_free(data->filter_data);
974  data->filter_data = NULL;
975  data->filter_data_free = NULL;
976  /*
977  * Nothing more to do if there is no filter
978  */
979  if (!data->filter)
980  return;
981  /*
982  * Create a new one
983  */
984  mne_create_filter_response(data->filter,
985  data->info->sfreq,
986  &data->filter_data,
987  &data->filter_data_free,
988  highpass_effective);
989 }
990 
991 //=============================================================================================================
992 
993 void MneRawData::setup_filter_bufs(MneRawData *data)
994 /*
995  * These will hold the filtered data
996  */
997 {
998  mneFilterDef filter;
999  int nfilt_buf;
1000  MneRawBufDef* bufs;
1001  int j,k;
1002  int firstsamp;
1003  int nring_buf;
1004  int highpass_effective;
1005 
1006  MneRawBufDef::free_bufs(data->filt_bufs,data->nfilt_buf);
1007  data->filt_bufs = NULL;
1008  data->nfilt_buf = 0;
1009  mne_free_ring_buffer(data->filt_ring);
1010  data->filt_ring = NULL;
1011 
1012  if (!data || !data->filter)
1013  return;
1014  filter = data->filter;
1015 
1016  for (nfilt_buf = 0, firstsamp = data->first_samp-filter->taper_size;
1017  firstsamp < data->nsamp + data->first_samp;
1018  firstsamp = firstsamp + filter->size)
1019  nfilt_buf++;
1020 #ifdef DEBUG
1021  printf("%d filter buffers needed\n",nfilt_buf);
1022 #endif
1023  bufs = MALLOC_36(nfilt_buf,MneRawBufDef);
1024  for (k = 0, firstsamp = data->first_samp-filter->taper_size; k < nfilt_buf; k++,
1025  firstsamp = firstsamp + filter->size) {
1026  bufs[k].ns = filter->size + 2*filter->taper_size;
1027  bufs[k].firsts = firstsamp;
1028  bufs[k].lasts = firstsamp + bufs[k].ns - 1;
1029  // bufs[k].ent = NULL;
1030  bufs[k].nchan = data->info->nchan;
1031  bufs[k].is_skip = FALSE;
1032  bufs[k].vals = NULL;
1033  bufs[k].valid = FALSE;
1034  bufs[k].ch_filtered = MALLOC_36(data->info->nchan,int);
1035  bufs[k].comp_status = MNE_CTFV_NOGRAD;
1036 
1037  for (j = 0; j < data->info->nchan; j++)
1038  bufs[k].ch_filtered[j] = FALSE;
1039  }
1040  data->filt_bufs = bufs;
1041  data->nfilt_buf = nfilt_buf;
1042  nring_buf = approx_ring_buf_size/((2*filter->taper_size+filter->size)*
1043  data->info->nchan*sizeof(float));
1044  data->filt_ring = mne_initialize_ring(nring_buf);
1045  mne_raw_add_filter_response(data,&highpass_effective);
1046 
1047  return;
1048 }
1049 
1050 //=============================================================================================================
1051 
1052 int MneRawData::load_one_buffer(MneRawData *data, MneRawBufDef *buf)
1053 /*
1054  * load just one
1055  */
1056 {
1057  if (buf->ent->kind == FIFF_DATA_SKIP) {
1058  printf("Cannot load a skip");
1059  return FAIL;
1060  }
1061  if (!buf->vals) { /* The data space may have been reused */
1062  buf->valid = FALSE;
1063  mne_allocate_from_ring(data->ring, buf->nchan,buf->ns,&buf->vals);
1064  }
1065  if (buf->valid)
1066  return OK;
1067 
1068 #ifdef DEBUG
1069  printf("Read buffer %d .. %d\n",buf->firsts,buf->lasts);
1070 #endif
1071 
1072  if (mne_read_raw_buffer_t(data->stream,
1073  buf->ent,
1074  buf->vals,
1075  buf->nchan,
1076  buf->ns,
1077  data->info->chInfo,
1078  NULL,0) != OK) {
1079  buf->valid = FALSE;
1080  return FAIL;
1081  }
1082  buf->valid = TRUE;
1083  buf->comp_status = data->comp_file;
1084  return OK;
1085 }
1086 
1087 //=============================================================================================================
1088 
1089 int MneRawData::compensate_buffer(MneRawData *data, MneRawBufDef *buf)
1090 /*
1091  * Apply compensation channels
1092  */
1093 {
1094  MneCTFCompData* temp;
1095 
1096  if (!data->comp)
1097  return OK;
1098  if (!data->comp->undo && !data->comp->current)
1099  return OK;
1100  if (buf->comp_status == data->comp_now)
1101  return OK;
1102  if (!buf->vals)
1103  return OK;
1104  /*
1105  * Have to do the hard job
1106  */
1107  if (data->comp->undo) {
1108  temp = data->comp->current;
1109  data->comp->current = data->comp->undo;
1110  data->comp->undo = temp;
1111  /*
1112  * Undo the previous compensation
1113  */
1114  if (MneCTFCompDataSet::mne_apply_ctf_comp_t(data->comp,FALSE,buf->vals,data->info->nchan,buf->ns) != OK) {
1115  temp = data->comp->undo;
1116  data->comp->undo = data->comp->current;
1117  data->comp->current = temp;
1118  goto bad;
1119  }
1120  temp = data->comp->undo;
1121  data->comp->undo = data->comp->current;
1122  data->comp->current = temp;
1123  }
1124  if (data->comp->current) {
1125  /*
1126  * Apply new compensation
1127  */
1128  if (MneCTFCompDataSet::mne_apply_ctf_comp_t(data->comp,TRUE,buf->vals,data->info->nchan,buf->ns) != OK)
1129  goto bad;
1130  }
1131  buf->comp_status = data->comp_now;
1132  return OK;
1133 
1134 bad :
1135  return FAIL;
1136 }
1137 
1138 //=============================================================================================================
1139 
1140 int MneRawData::mne_raw_pick_data(MneRawData *data, mneChSelection sel, int firsts, int ns, float **picked)
1141 /*
1142  * Data from a selection of channels
1143  */
1144 {
1145  int k,s,p,start,c,fills;
1146  int ns2,s2;
1147  MneRawBufDef* this_buf;
1148  float *values;
1149  int need_some;
1150 
1151  float **deriv_vals = NULL;
1152  int deriv_ns = 0;
1153  int nderiv = 0;
1154 
1155  if (firsts < data->first_samp) {
1156  for (s = 0, p = firsts; p < data->first_samp; s++, p++) {
1157  if (sel)
1158  for (c = 0; c < sel->nchan; c++)
1159  picked[c][s] = 0.0;
1160  else
1161  for (c = 0; c < data->info->nchan; c++)
1162  picked[c][s] = 0.0;
1163  }
1164  ns = ns - s;
1165  firsts = data->first_samp;
1166  }
1167  else
1168  s = 0;
1169  /*
1170  * There is possibly nothing to do
1171  */
1172  if (sel) {
1173  for (c = 0, need_some = FALSE; c < sel->nchan; c++) {
1174  if (sel->pick[c] >= 0 || sel->pick_deriv[c] >= 0) {
1175  need_some = TRUE;
1176  break;
1177  }
1178  }
1179  if (!need_some)
1180  return OK;
1181  }
1182  /*
1183  * Have to to the hard work
1184  */
1185  for (k = 0, this_buf = data->bufs, s = 0; k < data->nbuf; k++, this_buf++) {
1186  if (this_buf->lasts >= firsts) {
1187  start = firsts - this_buf->firsts;
1188  if (start < 0)
1189  start = 0;
1190  if (this_buf->is_skip) {
1191  for (p = start; p < this_buf->ns && ns > 0; p++, ns--, s++) {
1192  if (sel) {
1193  for (c = 0; c < sel->nchan; c++)
1194  if (sel->pick[c] >= 0)
1195  picked[c][s] = 0.0;
1196  }
1197  else {
1198  for (c = 0; c < data->info->nchan; c++)
1199  picked[c][s] = 0.0;
1200  }
1201  }
1202  }
1203  else {
1204  /*
1205  * Load the buffer
1206  */
1207  if (load_one_buffer(data,this_buf) != OK)
1208  return FAIL;
1209  /*
1210  * Apply compensation
1211  */
1212  if (compensate_buffer(data,this_buf) != OK)
1213  return FAIL;
1214  ns2 = s2 = 0;
1215  if (sel) {
1216  /*
1217  * Do we need the derived channels?
1218  */
1219  if (sel->nderiv > 0 && data->deriv_matched) {
1220  if (deriv_ns < this_buf->ns || nderiv != data->deriv_matched->deriv_data->nrow) {
1221  FREE_CMATRIX_36(deriv_vals);
1222  deriv_vals = ALLOC_CMATRIX_36(data->deriv_matched->deriv_data->nrow,this_buf->ns);
1223  nderiv = data->deriv_matched->deriv_data->nrow;
1224  deriv_ns = this_buf->ns;
1225  }
1226  if (mne_sparse_mat_mult2(data->deriv_matched->deriv_data->data,this_buf->vals,this_buf->ns,deriv_vals) == FAIL) {
1227  FREE_CMATRIX_36(deriv_vals);
1228  return FAIL;
1229  }
1230  }
1231  for (c = 0; c < sel->nchan; c++) {
1232  /*
1233  * First pick the ordinary channels...
1234  */
1235  if (sel->pick[c] >= 0) {
1236  values = this_buf->vals[sel->pick[c]];
1237  for (p = start, s2 = s, ns2 = ns; p < this_buf->ns && ns2 > 0; p++, ns2--, s2++)
1238  picked[c][s2] = values[p];
1239  }
1240  /*
1241  * ...then the derived ones
1242  */
1243  else if (sel->pick_deriv[c] >= 0 && data->deriv_matched) {
1244  values = deriv_vals[sel->pick_deriv[c]];
1245  for (p = start, s2 = s, ns2 = ns; p < this_buf->ns && ns2 > 0; p++, ns2--, s2++)
1246  picked[c][s2] = values[p];
1247  }
1248  }
1249  }
1250  else {
1251  for (c = 0; c < data->info->nchan; c++)
1252  for (p = start, s2 = s, ns2 = ns; p < this_buf->ns && ns2 > 0; p++, ns2--, s2++)
1253  picked[c][s2] = this_buf->vals[c][p];
1254  }
1255  s = s2;
1256  ns = ns2;
1257  }
1258  if (ns == 0)
1259  break;
1260  }
1261  }
1262  /*
1263  * Extend with the last available sample or zero if the request is beyond the data
1264  */
1265  if (s > 0) {
1266  fills = s-1;
1267  for (; ns > 0; ns--, s++) {
1268  if (sel)
1269  for (c = 0; c < sel->nchan; c++)
1270  picked[c][s] = picked[c][fills];
1271  else
1272  for (c = 0; c < data->info->nchan; c++)
1273  picked[c][s] = picked[c][fills];
1274  }
1275  }
1276  else {
1277  for (; ns > 0; ns--, s++) {
1278  if (sel)
1279  for (c = 0; c < sel->nchan; c++)
1280  picked[c][s] = 0;
1281  else
1282  for (c = 0; c < data->info->nchan; c++)
1283  picked[c][s] = 0;
1284  }
1285  }
1286  FREE_CMATRIX_36(deriv_vals);
1287  return OK;
1288 }
1289 
1290 //=============================================================================================================
1291 
1292 int MneRawData::mne_raw_pick_data_proj(MneRawData *data, mneChSelection sel, int firsts, int ns, float **picked)
1293 /*
1294  * Data from a set of channels, apply projection
1295  */
1296 {
1297  int k,s,p,start,c,fills;
1298  MneRawBufDef* this_buf;
1299  float **values;
1300  float *pvalues;
1301  float *deriv_pvalues = NULL;
1302 
1303  if (!data->proj || (sel && !MneProjOp::mne_proj_op_affect(data->proj,sel->chspick,sel->nchan) && !MneProjOp::mne_proj_op_affect(data->proj,sel->chspick_nospace,sel->nchan)))
1304  return mne_raw_pick_data(data,sel,firsts,ns,picked);
1305 
1306  if (firsts < data->first_samp) {
1307  for (s = 0, p = firsts; p < data->first_samp; s++, p++) {
1308  if (sel)
1309  for (c = 0; c < sel->nchan; c++)
1310  picked[c][s] = 0.0;
1311  else
1312  for (c = 0; c < data->info->nchan; c++)
1313  picked[c][s] = 0.0;
1314  }
1315  ns = ns - s;
1316  firsts = data->first_samp;
1317  }
1318  else
1319  s = 0;
1320  pvalues = MALLOC_36(data->info->nchan,float);
1321  for (k = 0, this_buf = data->bufs; k < data->nbuf; k++, this_buf++) {
1322  if (this_buf->lasts >= firsts) {
1323  start = firsts - this_buf->firsts;
1324  if (start < 0)
1325  start = 0;
1326  if (this_buf->is_skip) {
1327  for (p = start; p < this_buf->ns && ns > 0; p++, ns--, s++) {
1328  if (sel) {
1329  for (c = 0; c < sel->nchan; c++)
1330  if (sel->pick[c] >= 0)
1331  picked[c][s] = 0.0;
1332  }
1333  else {
1334  for (c = 0; c < data->info->nchan; c++)
1335  picked[c][s] = 0.0;
1336  }
1337  }
1338  }
1339  else {
1340  /*
1341  * Load the buffer
1342  */
1343  if (load_one_buffer(data,this_buf) != OK)
1344  return FAIL;
1345  /*
1346  * Apply compensation
1347  */
1348  if (compensate_buffer(data,this_buf) != OK)
1349  return FAIL;
1350  /*
1351  * Apply projection
1352  */
1353  values = this_buf->vals;
1354  if (sel && sel->nderiv > 0 && data->deriv_matched)
1355  deriv_pvalues = MALLOC_36(data->deriv_matched->deriv_data->nrow,float);
1356  for (p = start; p < this_buf->ns && ns > 0; p++, ns--, s++) {
1357  for (c = 0; c < data->info->nchan; c++)
1358  pvalues[c] = values[c][p];
1359  if (MneProjOp::mne_proj_op_proj_vector(data->proj,pvalues,data->info->nchan,TRUE) != OK)
1360  qWarning()<<"Error";
1361  if (sel) {
1362  if (sel->nderiv > 0 && data->deriv_matched) {
1363  if (mne_sparse_vec_mult2(data->deriv_matched->deriv_data->data,pvalues,deriv_pvalues) == FAIL)
1364  return FAIL;
1365  }
1366  for (c = 0; c < sel->nchan; c++) {
1367  /*
1368  * First try the ordinary channels...
1369  */
1370  if (sel->pick[c] >= 0)
1371  picked[c][s] = pvalues[sel->pick[c]];
1372  /*
1373  * ...then the derived ones
1374  */
1375  else if (sel->pick_deriv[c] >= 0 && data->deriv_matched)
1376  picked[c][s] = deriv_pvalues[sel->pick_deriv[c]];
1377  }
1378  }
1379  else {
1380  for (c = 0; c < data->info->nchan; c++) {
1381  picked[c][s] = pvalues[c];
1382  }
1383  }
1384  }
1385  }
1386  if (ns == 0)
1387  break;
1388  }
1389  }
1390  FREE_36(deriv_pvalues);
1391  FREE_36(pvalues);
1392  /*
1393  * Extend with the last available sample or zero if the request is beyond the data
1394  */
1395  if (s > 0) {
1396  fills = s-1;
1397  for (; ns > 0; ns--, s++) {
1398  if (sel)
1399  for (c = 0; c < sel->nchan; c++)
1400  picked[c][s] = picked[c][fills];
1401  else
1402  for (c = 0; c < data->info->nchan; c++)
1403  picked[c][s] = picked[c][fills];
1404  }
1405  }
1406  else {
1407  for (; ns > 0; ns--, s++) {
1408  if (sel)
1409  for (c = 0; c < sel->nchan; c++)
1410  picked[c][s] = 0;
1411  else
1412  for (c = 0; c < data->info->nchan; c++)
1413  picked[c][s] = 0;
1414  }
1415  }
1416  return OK;
1417 }
1418 
1419 //=============================================================================================================
1420 
1421 int MneRawData::load_one_filt_buf(MneRawData *data, MneRawBufDef *buf)
1422 /*
1423  * Load and filter one buffer
1424  */
1425 {
1426  int k;
1427  int res;
1428 
1429  float **vals;
1430 
1431  if (!buf->vals) {
1432  buf->valid = FALSE;
1433  mne_allocate_from_ring(data->filt_ring, buf->nchan, buf->ns,&buf->vals);
1434  }
1435  if (buf->valid)
1436  return OK;
1437 
1438  vals = MALLOC_36(buf->nchan,float *);
1439  for (k = 0; k < buf->nchan; k++) {
1440  buf->ch_filtered[k] = FALSE;
1441  vals[k] = buf->vals[k] + data->filter->taper_size;
1442  }
1443 
1444  res = mne_raw_pick_data_proj(data,NULL,buf->firsts + data->filter->taper_size,buf->ns - 2*data->filter->taper_size,vals);
1445 
1446  FREE_36(vals);
1447 
1448 #ifdef DEBUG
1449  if (res == OK)
1450  printf("Loaded filtered buffer %d...%d %d %d last = %d\n",
1451  buf->firsts,buf->lasts,buf->lasts-buf->firsts+1,buf->ns,data->first_samp + data->nsamp);
1452 #endif
1453  buf->valid = res == OK;
1454  return res;
1455 }
1456 
1457 //=============================================================================================================
1458 
1459 int MneRawData::mne_raw_pick_data_filt(MneRawData *data, mneChSelection sel, int firsts, int ns, float **picked)
1460 /*
1461  * Data for a selection (filtered and picked)
1462  */
1463 {
1464  int k,s,bs,c;
1465  int bs1,bs2,s1,s2,lasts;
1466  MneRawBufDef* this_buf;
1467  float *values;
1468  float **deriv_vals = NULL;
1469  float *dc = NULL;
1470  float dc_offset;
1471  int deriv_ns = 0;
1472  int nderiv = 0;
1473  int filter_was;
1474 
1475  if (!data->filter || !data->filter->filter_on)
1476  return mne_raw_pick_data_proj(data,sel,firsts,ns,picked);
1477 
1478  if (sel) {
1479  for (s = 0; s < ns; s++)
1480  for (c = 0; c < sel->nchan; c++)
1481  picked[c][s] = 0.0;
1482  }
1483  else {
1484  for (s = 0; s < ns; s++)
1485  for (c = 0; c < data->info->nchan; c++)
1486  picked[c][s] = 0.0;
1487  }
1488  lasts = firsts + ns - 1;
1489  /*
1490  * Take into account the initial dc offset (compensate and project)
1491  */
1492  if (data->first_sample_val) {
1493  dc = MALLOC_36(data->info->nchan,float);
1494 
1495  for (k = 0; k < data->info->nchan; k++)
1496  dc[k] = data->first_sample_val[k];
1497  /*
1498  * Is this correct??
1499  */
1500  if (data->comp && data->comp->current)
1501  if (MneCTFCompDataSet::mne_apply_ctf_comp(data->comp,TRUE,dc,data->info->nchan,NULL,0) != OK)
1502  goto bad;
1503  if (data->proj)
1504  if (MneProjOp::mne_proj_op_proj_vector(data->proj,dc,data->info->nchan,TRUE) != OK)
1505  goto bad;
1506  }
1507  filter_was = data->filter->filter_on;
1508  /*
1509  * Find the first buffer to consider
1510  */
1511  for (k = 0, this_buf = data->filt_bufs; k < data->nfilt_buf; k++, this_buf++) {
1512  if (this_buf->lasts >= firsts)
1513  break;
1514  }
1515  for (; k < data->nfilt_buf && this_buf->firsts <= lasts; k++, this_buf++) {
1516 #ifdef DEBUG
1517  printf("this_buf (%d): %d..%d\n",k,this_buf->firsts,this_buf->lasts);
1518 #endif
1519  /*
1520  * Load the buffer first and apply projection
1521  */
1522  if (load_one_filt_buf(data,this_buf) != OK)
1523  goto bad;
1524  /*
1525  * Then filter all relevant channels (not stimuli)
1526  */
1527  if (sel) {
1528  for (c = 0; c < sel->nchan; c++) {
1529  if (sel->pick[c] >= 0) {
1530  if (!this_buf->ch_filtered[sel->pick[c]]) {
1531  /*
1532  * Do not filter stimulus channels
1533  */
1534  dc_offset = 0.0;
1535  if (data->info->chInfo[sel->pick[c]].kind == FIFFV_STIM_CH)
1536  data->filter->filter_on = FALSE;
1537  else if (dc)
1538  dc_offset = dc[sel->pick[c]];
1539  if (mne_apply_filter(data->filter,data->filter_data,this_buf->vals[sel->pick[c]],this_buf->ns,TRUE,
1540  dc_offset,data->info->chInfo[sel->pick[c]].kind) != OK) {
1541  data->filter->filter_on = filter_was;
1542  goto bad;
1543  }
1544  this_buf->ch_filtered[sel->pick[c]] = TRUE;
1545  data->filter->filter_on = filter_was;
1546  }
1547  }
1548  }
1549  /*
1550  * Also check channels included in derivations if they are used
1551  */
1552  if (sel->nderiv > 0 && data->deriv_matched) {
1553  MneDeriv* der = data->deriv_matched;
1554  for (c = 0; c < der->deriv_data->ncol; c++) {
1555  if (der->in_use[c] > 0 &&
1556  !this_buf->ch_filtered[c]) {
1557  /*
1558  * Do not filter stimulus channels
1559  */
1560  dc_offset = 0.0;
1561  if (data->info->chInfo[c].kind == FIFFV_STIM_CH)
1562  data->filter->filter_on = FALSE;
1563  else if (dc)
1564  dc_offset = dc[c];
1565  if (mne_apply_filter(data->filter,data->filter_data,this_buf->vals[c],this_buf->ns,TRUE,
1566  dc_offset,data->info->chInfo[c].kind) != OK) {
1567  data->filter->filter_on = filter_was;
1568  goto bad;
1569  }
1570  this_buf->ch_filtered[c] = TRUE;
1571  data->filter->filter_on = filter_was;
1572  }
1573  }
1574  }
1575  }
1576  else {
1577  /*
1578  * Simply filter all channels if there is no selection
1579  */
1580  for (c = 0; c < data->info->nchan; c++) {
1581  if (!this_buf->ch_filtered[c]) {
1582  /*
1583  * Do not filter stimulus channels
1584  */
1585  dc_offset = 0.0;
1586  if (data->info->chInfo[c].kind == FIFFV_STIM_CH)
1587  data->filter->filter_on = FALSE;
1588  else if (dc)
1589  dc_offset = dc[c];
1590  if (mne_apply_filter(data->filter,data->filter_data,this_buf->vals[c],this_buf->ns,TRUE,
1591  dc_offset,data->info->chInfo[c].kind) != OK) {
1592  data->filter->filter_on = filter_was;
1593  goto bad;
1594  }
1595  this_buf->ch_filtered[c] = TRUE;
1596  data->filter->filter_on = filter_was;
1597  }
1598  }
1599  }
1600  /*
1601  * Decide the picking limits
1602  */
1603  if (firsts >= this_buf->firsts) {
1604  bs1 = firsts - this_buf->firsts;
1605  s1 = 0;
1606  }
1607  else {
1608  bs1 = 0;
1609  s1 = this_buf->firsts - firsts;
1610  }
1611  if (lasts >= this_buf->lasts) {
1612  bs2 = this_buf->ns;
1613  s2 = this_buf->lasts - lasts + ns;
1614  }
1615  else {
1616  bs2 = lasts - this_buf->lasts + this_buf->ns;
1617  s2 = ns;
1618  }
1619 #ifdef DEBUG
1620  printf("buf : %d..%d %d\n",bs1,bs2,bs2-bs1);
1621  printf("dest : %d..%d %d\n",s1,s2,s2-s1);
1622 #endif
1623  /*
1624  * Then pick data from all relevant channels
1625  */
1626  if (sel) {
1627  if (sel->nderiv > 0 && data->deriv_matched) {
1628  /*
1629  * Compute derived data if we need it
1630  */
1631  if (deriv_ns < this_buf->ns || nderiv != data->deriv_matched->deriv_data->nrow) {
1632  FREE_CMATRIX_36(deriv_vals);
1633  deriv_vals = ALLOC_CMATRIX_36(data->deriv_matched->deriv_data->nrow,this_buf->ns);
1634  nderiv = data->deriv_matched->deriv_data->nrow;
1635  deriv_ns = this_buf->ns;
1636  }
1637  if (mne_sparse_mat_mult2(data->deriv_matched->deriv_data->data,this_buf->vals,this_buf->ns,deriv_vals) == FAIL)
1638  goto bad;
1639  }
1640  for (c = 0; c < sel->nchan; c++) {
1641  /*
1642  * First the ordinary channels
1643  */
1644  if (sel->pick[c] >= 0) {
1645  values = this_buf->vals[sel->pick[c]];
1646  for (s = s1, bs = bs1; s < s2; s++, bs++)
1647  picked[c][s] += values[bs];
1648  }
1649  else if (sel->pick_deriv[c] >= 0 && data->deriv_matched) {
1650  values = deriv_vals[sel->pick_deriv[c]];
1651  for (s = s1, bs = bs1; s < s2; s++, bs++)
1652  picked[c][s] += values[bs];
1653  }
1654  }
1655  }
1656  else {
1657  for (c = 0; c < data->info->nchan; c++) {
1658  values = this_buf->vals[c];
1659  for (s = s1, bs = bs1; s < s2; s++, bs++)
1660  picked[c][s] += values[bs];
1661  }
1662  }
1663  }
1664  FREE_CMATRIX_36(deriv_vals);
1665  FREE_36(dc);
1666  return OK;
1667 
1668 bad : {
1669  FREE_CMATRIX_36(deriv_vals);
1670  FREE_36(dc);
1671  return FAIL;
1672  }
1673 }
1674 
1675 //=============================================================================================================
1676 
1677 MneRawData *MneRawData::mne_raw_open_file_comp(const QString& name,
1678  int omit_skip,
1679  int allow_maxshield,
1680  mneFilterDef filter,
1681  int comp_set)
1682 /*
1683  * Open a raw data file
1684  */
1685 {
1686  MneRawInfo* info = NULL;
1687  MneRawData* data = NULL;
1688 
1689  QFile file(name);
1690  FiffStream::SPtr stream(new FiffStream(&file));
1691  // fiffFile in = NULL;
1692 
1693  FiffDirEntry::SPtr dir;
1694  QList<FiffDirEntry::SPtr> dir0;
1695  // fiffTagRec tag;
1696  FiffTag::SPtr t_pTag;
1697  FiffChInfo ch;
1698  MneRawBufDef* bufs;
1699  int k, b, nbuf, ndir, nnames;
1700  int current_dir0 = 0;
1701 
1702  // tag.data = NULL;
1703 
1704  if (MneRawInfo::mne_load_raw_info(name,allow_maxshield,&info) == FAIL)
1705  goto bad;
1706 
1707  for (k = 0; k < info->nchan; k++) {
1708  ch = info->chInfo.at(k);
1709  if (QString::compare(ch.ch_name,MNE_DEFAULT_TRIGGER_CH) == 0) {
1710  if (std::fabs(1.0 - ch.range) > 1e-5) {
1711  ch.range = 1.0;
1712  printf("%s range set to %f\n",MNE_DEFAULT_TRIGGER_CH,ch.range);
1713  }
1714  }
1715  /*
1716  * Take care of the nonzero unit multiplier
1717  */
1718  if (ch.unit_mul != 0) {
1719  ch.cal = pow(10.0,(double)(ch.unit_mul))*ch.cal;
1720  printf("Ch %s unit multiplier %d -> 0\n",ch.ch_name.toLatin1().data(),ch.unit_mul);
1721  ch.unit_mul = 0;
1722  }
1723  }
1724  // if ((in = fiff_open(name)) == NULL)
1725  // goto bad;
1726  if(!stream->open())
1727  goto bad;
1728 
1729  data = new MneRawData;
1730  data->filename = name;
1731  data->stream = stream;
1732  data->info = info;
1733  /*
1734  * Add the channel name list
1735  */
1736  mne_channel_names_to_name_list(info->chInfo,
1737  info->nchan,
1738  data->ch_names,
1739  nnames);
1740  if (nnames != info->nchan) {
1741  printf("Channel names were not translated correctly into a name list");
1742  goto bad;
1743  }
1744  /*
1745  * Compensation data
1746  */
1747  data->comp = MneCTFCompDataSet::mne_read_ctf_comp_data(data->filename);
1748  if (data->comp) {
1749  if (data->comp->ncomp > 0)
1750  printf("Read %d compensation data sets from %s\n",data->comp->ncomp,data->filename.toUtf8().constData());
1751  else
1752  printf("No compensation data in %s\n",data->filename.toUtf8().constData());
1753  }
1754  else
1755  qWarning() << "err_print_error()";
1756  if ((data->comp_file = MneCTFCompDataSet::mne_get_ctf_comp(data->info->chInfo,data->info->nchan)) == FAIL)
1757  goto bad;
1758  printf("Compensation in file : %s\n",MneCTFCompDataSet::mne_explain_ctf_comp(MneCTFCompDataSet::mne_map_ctf_comp_kind(data->comp_file)));
1759  if (comp_set < 0)
1760  data->comp_now = data->comp_file;
1761  else
1762  data->comp_now = comp_set;
1763 
1765  data->comp_now,
1766  data->info->chInfo,
1767  data->info->nchan,
1768  QList<FIFFLIB::FiffChInfo>(),
1769  0) == FAIL)
1770  goto bad;
1771  /*
1772  * SSS data
1773  */
1774  data->sss = MneSssData::read_sss_data(data->filename);
1775  if (data->sss && data->sss->job != FIFFV_SSS_JOB_NOTHING && data->sss->ncomp > 0) {
1776  printf("SSS data read from %s :\n",data->filename.toUtf8().constData());
1777  data->sss->print(stderr);
1778  }
1779  else {
1780  printf("No SSS data in %s\n",data->filename.toUtf8().constData());
1781  if(data->sss)
1782  delete data->sss;
1783  data->sss = NULL;
1784  }
1785  /*
1786  * Buffers
1787  */
1788  dir0 = data->info->rawDir;
1789  ndir = data->info->ndir;
1790  /*
1791  * Take into account the first sample
1792  */
1793  if (dir0[current_dir0]->kind == FIFF_FIRST_SAMPLE) {
1794  // if (fiff_read_this_tag(in->fd,dir0->pos,&tag) == FIFF_FAIL)
1795  // goto bad;
1796  if (!stream->read_tag(t_pTag,dir0[current_dir0]->pos))
1797  goto bad;
1798  data->first_samp = *t_pTag->toInt();
1799  current_dir0++;
1800  ndir--;
1801  }
1802  if (dir0[current_dir0]->kind == FIFF_DATA_SKIP) {
1803  int nsamp_skip;
1804  // if (fiff_read_this_tag(in->fd,dir0->pos,&tag) == FIFF_FAIL)
1805  // goto bad;
1806  if (!stream->read_tag(t_pTag,dir0[current_dir0]->pos))
1807  goto bad;
1808  nsamp_skip = data->info->buf_size*(*t_pTag->toInt());
1809  printf("Data skip of %d samples in the beginning\n",nsamp_skip);
1810  current_dir0++;
1811  ndir--;
1812  if (dir0[current_dir0]->kind == FIFF_FIRST_SAMPLE) {
1813  // if (fiff_read_this_tag(in->fd,dir0->pos,&tag) == FIFF_FAIL)
1814  // goto bad;
1815  if (!stream->read_tag(t_pTag,dir0[current_dir0]->pos))
1816  goto bad;
1817  data->first_samp += *t_pTag->toInt();
1818  current_dir0++;
1819  ndir--;
1820  }
1821  if (omit_skip) {
1822  data->omit_samp = data->first_samp + nsamp_skip;
1823  data->omit_samp_old = nsamp_skip;
1824  data->first_samp = 0;
1825  }
1826  else {
1827  data->first_samp = data->first_samp + nsamp_skip;
1828  }
1829  }
1830  else if (omit_skip) {
1831  data->omit_samp = data->first_samp;
1832  data->first_samp = 0;
1833  }
1834 #ifdef DEBUG
1835  printf("data->first_samp = %d\n",data->first_samp);
1836 #endif
1837  /*
1838  * Figure out the buffers
1839  */
1840  // for (k = 0, dir = dir0, nbuf = 0; k < ndir; k++, dir++)
1841  for (k = 0, nbuf = 0; k < ndir; k++)
1842  if (dir0[k]->kind == FIFF_DATA_BUFFER ||
1843  dir0[k]->kind == FIFF_DATA_SKIP)
1844  nbuf++;
1845  bufs = MALLOC_36(nbuf,MneRawBufDef);
1846 
1847  // for (k = 0, nbuf = 0, dir = dir0; k < ndir; k++, dir++)
1848  for (k = 0, nbuf = 0; k < ndir; k++)
1849  if (dir0[k]->kind == FIFF_DATA_BUFFER ||
1850  dir0[k]->kind == FIFF_DATA_SKIP) {
1851  bufs[nbuf].ns = 0;
1852  bufs[nbuf].ent = dir0[k];
1853  bufs[nbuf].nchan = data->info->nchan;
1854  bufs[nbuf].is_skip = dir0[k]->kind == FIFF_DATA_SKIP;
1855  bufs[nbuf].vals = NULL;
1856  bufs[nbuf].valid = FALSE;
1857  bufs[nbuf].ch_filtered = NULL;
1858  bufs[nbuf].comp_status = data->comp_file;
1859  nbuf++;
1860  }
1861  data->bufs = bufs;
1862  data->nbuf = nbuf;
1863  data->nsamp = 0;
1864  for (k = 0; k < nbuf; k++) {
1865  dir = bufs[k].ent;
1866  if (dir->kind == FIFF_DATA_BUFFER) {
1867  if (dir->type == FIFFT_DAU_PACK16 || dir->type == FIFFT_SHORT)
1868  bufs[k].ns = dir->size/(data->info->nchan*sizeof(fiff_dau_pack16_t));
1869  else if (dir->type == FIFFT_FLOAT)
1870  bufs[k].ns = dir->size/(data->info->nchan*sizeof(fiff_float_t));
1871  else if (dir->type == FIFFT_INT)
1872  bufs[k].ns = dir->size/(data->info->nchan*sizeof(fiff_int_t));
1873  else {
1874  printf("We are not prepared to handle raw data type: %d",dir->type);
1875  goto bad;
1876  }
1877  }
1878  else if (dir->kind == FIFF_DATA_SKIP) {
1879  // if (fiff_read_this_tag(in->fd,dir->pos,&tag) == FIFF_FAIL)
1880  // goto bad;
1881  if (!stream->read_tag(t_pTag,dir->pos))
1882  goto bad;
1883  bufs[k].ns = data->info->buf_size*(*t_pTag->toInt());
1884  }
1885  bufs[k].firsts = k == 0 ? data->first_samp : bufs[k-1].lasts + 1;
1886  bufs[k].lasts = bufs[k].firsts + bufs[k].ns - 1;
1887  data->nsamp += bufs[k].ns;
1888  }
1889  // FREE_36(tag.data);
1890  /*
1891  * Set up the first sample values
1892  */
1893  data->bad = MALLOC_36(data->info->nchan,int);
1894  data->offsets = MALLOC_36(data->info->nchan,float);
1895  for (k = 0; k < data->info->nchan; k++) {
1896  data->bad[k] = FALSE;
1897  data->offsets[k] = 0.0;
1898  }
1899  /*
1900  * Th bad channel stuff
1901  */
1902  {
1903  if (mne_read_bad_channel_list(name,data->badlist,data->nbad) == OK) {
1904  for (b = 0; b < data->nbad; b++) {
1905  for (k = 0; k < data->info->nchan; k++) {
1906  if (QString::compare(data->info->chInfo[k].ch_name,data->badlist[b],Qt::CaseInsensitive) == 0) {
1907  data->bad[k] = TRUE;
1908  break;
1909  }
1910  }
1911  }
1912  printf("%d bad channels read from %s%s",data->nbad,name.toUtf8().constData(),data->nbad > 0 ? ":\n" : "\n");
1913  if (data->nbad > 0) {
1914  printf("\t");
1915  for (k = 0; k < data->nbad; k++)
1916  printf("%s%c",data->badlist[k].toUtf8().constData(),k < data->nbad-1 ? ' ' : '\n');
1917  }
1918  }
1919  }
1920  /*
1921  * Initialize the raw data buffers
1922  */
1923  nbuf = approx_ring_buf_size/(data->info->buf_size*data->info->nchan*sizeof(float));
1924  data->ring = mne_initialize_ring(nbuf);
1925  /*
1926  * Initialize the filter buffers
1927  */
1928  data->filter = MALLOC_36(1,mneFilterDefRec);
1929  *data->filter = *filter;
1930  setup_filter_bufs(data);
1931 
1932  {
1933  float **vals = ALLOC_CMATRIX_36(data->info->nchan,1);
1934 
1935  if (mne_raw_pick_data(data,NULL,data->first_samp,1,vals) == FAIL)
1936  goto bad;
1937  data->first_sample_val = MALLOC_36(data->info->nchan,float);
1938  for (k = 0; k < data->info->nchan; k++)
1939  data->first_sample_val[k] = vals[k][0];
1940  FREE_CMATRIX_36(vals);
1941  printf("Initial dc offsets determined\n");
1942  }
1943  printf("Raw data file %s:\n",name.toUtf8().constData());
1944  printf("\tnchan = %d\n",data->info->nchan);
1945  printf("\tnsamp = %d\n",data->nsamp);
1946  printf("\tsfreq = %-8.3f Hz\n",data->info->sfreq);
1947  printf("\tlength = %-8.3f sec\n",data->nsamp/data->info->sfreq);
1948 
1949  return data;
1950 
1951 bad : {
1952  if (data)
1953  delete(data);
1954  else
1955  if(info)
1956  delete info;
1957 
1958  return NULL;
1959  }
1960 }
1961 
1962 //=============================================================================================================
1963 
1964 MneRawData *MneRawData::mne_raw_open_file(const QString& name, int omit_skip, int allow_maxshield, mneFilterDef filter)
1965 /*
1966  * Wrapper for mne_raw_open_file to work as before
1967  */
1968 {
1969  return mne_raw_open_file_comp(name,omit_skip,allow_maxshield,filter,-1);
1970 }
#define FIFF_FIRST_SAMPLE
Definition: fiff_file.h:461
Information about raw data in fiff file.
Definition: mne_raw_info.h:80
fiff_float_t range
Definition: fiff_ch_info.h:122
QSharedPointer< FiffDirEntry > SPtr
Channel info descriptor.
Definition: fiff_ch_info.h:74
static QStringList split_name_list(QString p_sNameList)
QList< FIFFLIB::FiffDirEntry::SPtr > rawDir
Definition: mne_raw_info.h:146
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
static int mne_ctf_set_compensation(MneCTFCompDataSet *set, int compensate_to, QList< FIFFLIB::FiffChInfo > &chs, int nchan, QList< FIFFLIB::FiffChInfo > comp_chs, int ncomp_chan)
FIFFLIB::fiff_int_t coding
One item in a derivation data set.
Definition: mne_deriv.h:73
A comprehensive raw data structure.
Definition: mne_raw_data.h:86
One MNE CTF compensation description.
MneRawData class declaration.
Information about raw data in fiff file.
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
void print(FILE *f) const
FIFF File I/O routines.
Definition: fiff_stream.h:104
static MneSssData * read_sss_data(const QString &name)
#define FIFF_DATA_SKIP
Definition: fiff_file.h:559
#define FIFF_DATA_BUFFER
Definition: fiff_file.h:558
FIFFLIB::fiff_int_t * ptrs
FIFFLIB::fiff_int_t * inds
#define FIFFV_SSS_JOB_NOTHING
Definition: fiff_file.h:540
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
fiff_float_t cal
Definition: fiff_ch_info.h:123
FIFFLIB::fiff_float_t * data
QList< FIFFLIB::FiffChInfo > chInfo
Definition: mne_raw_info.h:130
Data associated with MNE computations for each mneMeasDataSet.