MNE-CPP  beta 1.0
bci.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //*************************************************************************************************************
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include "bci.h"
43 
44 
45 //*************************************************************************************************************
46 //=============================================================================================================
47 // QT INCLUDES
48 //=============================================================================================================
49 
50 #include <QtCore/QtPlugin>
51 #include <QtCore/QTextStream>
52 #include <QDebug>
53 
54 
55 //*************************************************************************************************************
56 //=============================================================================================================
57 // USED NAMESPACES
58 //=============================================================================================================
59 
60 using namespace BCIPlugin;
61 using namespace std;
62 using namespace UTILSLIB;
63 
64 //*************************************************************************************************************
65 //=============================================================================================================
66 // DEFINE MEMBER METHODS
67 //=============================================================================================================
68 
70 : m_qStringResourcePath(qApp->applicationDirPath()+"/mne_x_plugins/resources/bci/")
71 , m_bProcessData(false)
72 {
73 }
74 
75 
76 //*************************************************************************************************************
77 
79 {
80  //std::cout << "BCI::~BCI() " << std::endl;
81 
82  //If the program is closed while the sampling is in process
83  if(this->isRunning())
84  this->stop();
85 }
86 
87 
88 //*************************************************************************************************************
89 
90 QSharedPointer<IPlugin> BCI::clone() const
91 {
92  QSharedPointer<BCI> pTMSIClone(new BCI());
93  return pTMSIClone;
94 }
95 
96 
97 //*************************************************************************************************************
98 
99 void BCI::init()
100 {
101  m_bIsRunning = false;
102  m_bTriggerActivated = false;
103 
104  // Inputs - Source estimates and sensor level
105  m_pRTSEInput = PluginInputData<RealTimeSourceEstimate>::create(this, "BCIInSource", "BCI source input data");
106  connect(m_pRTSEInput.data(), &PluginInputConnector::notify, this, &BCI::updateSource, Qt::DirectConnection);
107  m_inputConnectors.append(m_pRTSEInput);
108 
109  m_pRTMSAInput = PluginInputData<NewRealTimeMultiSampleArray>::create(this, "BCIInSensor", "SourceLab sensor input data");
110  connect(m_pRTMSAInput.data(), &PluginInputConnector::notify, this, &BCI::updateSensor, Qt::DirectConnection);
111  m_inputConnectors.append(m_pRTMSAInput);
112 
113  // Output streams
114  m_pBCIOutputOne = PluginOutputData<NewRealTimeSampleArray>::create(this, "ControlSignal", "BCI output data One");
115  m_pBCIOutputOne->data()->setArraySize(1);
116  m_pBCIOutputOne->data()->setName("Boundary");
117  m_outputConnectors.append(m_pBCIOutputOne);
118 
119  m_pBCIOutputTwo = PluginOutputData<NewRealTimeSampleArray>::create(this, "ControlSignal", "BCI output data Two");
120  m_pBCIOutputTwo->data()->setArraySize(1);
121  m_pBCIOutputTwo->data()->setName("Left electrode var");
122  m_outputConnectors.append(m_pBCIOutputTwo);
123 
124  m_pBCIOutputThree = PluginOutputData<NewRealTimeSampleArray>::create(this, "ControlSignal", "BCI output data Three");
125  m_pBCIOutputThree->data()->setArraySize(1);
126  m_pBCIOutputThree->data()->setName("Right electrode var");
127  m_outputConnectors.append(m_pBCIOutputThree);
128 
129  m_pBCIOutputFour = PluginOutputData<NewRealTimeSampleArray>::create(this, "ControlSignal", "BCI output data Four");
130  m_pBCIOutputFour->data()->setArraySize(1);
131  m_pBCIOutputFour->data()->setName("Left electrode");
132  m_outputConnectors.append(m_pBCIOutputFour);
133 
134  m_pBCIOutputFive = PluginOutputData<NewRealTimeSampleArray>::create(this, "ControlSignal", "BCI output data Five");
135  m_pBCIOutputFive->data()->setArraySize(1);
136  m_pBCIOutputFive->data()->setName("Right electrode");
137  m_outputConnectors.append(m_pBCIOutputFive);
138 
139  // Delete Buffer - will be initailzed with first incoming data
140  m_pBCIBuffer_Sensor = CircularMatrixBuffer<double>::SPtr();
141  m_pBCIBuffer_Source = CircularMatrixBuffer<double>::SPtr();
142 
143  // Delete fiff info because the initialisation of the fiff info is seen as the first data acquisition from the input stream
144  m_pFiffInfo_Sensor = FiffInfo::SPtr();
145 
146  // Intitalise GUI stuff
147  m_bSubtractMean = true;
148  m_bUseFilter = true;
149  m_bUseSensorData = true;
150  m_bUseSourceData = false;
151  m_bDisplayFeatures = true;
152  m_bUseArtefactThresholdReduction = false;
153  m_dSlidingWindowSize = 0.5;
154  m_dTimeBetweenWindows = 0.5;
155  m_dDisplayRangeBoundary = 15;
156  m_dDisplayRangeVariances = 5;
157  m_dDisplayRangeElectrodes = 7;
158  m_iNumberFeatures = 1;
159  m_dThresholdValue = 30;
160  m_iNumberFeaturesToDisplay = 30;
161  m_iFeatureCalculationType = 0;
162 
163  // Intitalise feature selection
164  m_slChosenFeatureSensor << "LA4" << "RA4"; //<< "TEST";
165 
166  // Initalise sliding window stuff
167  m_bFillSensorWindowFirstTime = true;
168 
169  // Initialise filter stuff
170  m_filterOperator = QSharedPointer<FilterData>(new FilterData());
171 
172  m_dFilterLowerBound = 7.0;
173  m_dFilterUpperBound = 14.0;
174  m_dParcksWidth = m_dFilterLowerBound-1; // (m_dFilterUpperBound-m_dFilterLowerBound)/2;
175  m_iFilterOrder = 256;
176 
177  // Init BCIFeatureWindow for visualization
178  m_BCIFeatureWindow = QSharedPointer<BCIFeatureWindow>(new BCIFeatureWindow(this));
179 }
180 
181 
182 //*************************************************************************************************************
183 
185 {
186 
187 }
188 
189 
190 //*************************************************************************************************************
191 
193 {
194  // Do not start BCI if the number of chosen electrodes/features is uneven
195  if(m_slChosenFeatureSensor.size()%2 != 0)
196  {
197  if(!m_slChosenFeatureSensor.contains("TEST"))
198  {
199  QMessageBox msgBox;
200  msgBox.setText("The number of selected electrodes needs to be even.");
201  msgBox.setInformativeText("Please select feautres again");
202  msgBox.setStandardButtons(QMessageBox::Ok);
203  msgBox.exec();
204  return false;
205  }
206  }
207 
208  // Initialise index
209  m_iTBWIndexSensor = 0;
210  m_iNumberOfCalculatedFeatures = 0;
211 
212  // BCIFeatureWindow show and init
213  if(m_bDisplayFeatures)
214  {
215  m_BCIFeatureWindow->initGui();
216  m_BCIFeatureWindow->show();
217  }
218 
219  // Init debug output stream
220  QString path("BCIDebugFile.txt");
221  path.prepend(m_qStringResourcePath);
222  m_outStreamDebug.open(path.toStdString(), ios::trunc);
223 
224  m_pFiffInfo_Sensor = FiffInfo::SPtr();
225 
226  m_bFillSensorWindowFirstTime = true;
227 
228 // // Set display ranges for output channels
229 // m_pBCIOutputOne->data()->setMaxValue(m_dDisplayRangeBoundary);
230 // m_pBCIOutputOne->data()->setMinValue(-m_dDisplayRangeBoundary);
231 
232 // m_pBCIOutputTwo->data()->setMaxValue(m_dDisplayRangeVariances);
233 // m_pBCIOutputTwo->data()->setMinValue(0);
234 
235 // m_pBCIOutputThree->data()->setMaxValue(m_dDisplayRangeVariances);
236 // m_pBCIOutputThree->data()->setMinValue(0);
237 
238 // m_pBCIOutputFour->data()->setMaxValue(m_dDisplayRangeElectrodes);
239 // m_pBCIOutputFour->data()->setMinValue(-m_dDisplayRangeElectrodes);
240 
241 // m_pBCIOutputFive->data()->setMaxValue(m_dDisplayRangeElectrodes);
242 // m_pBCIOutputFive->data()->setMinValue(-m_dDisplayRangeElectrodes);
243 
244  // variance
245  m_pBCIOutputOne->data()->setMaxValue(5);
246  m_pBCIOutputOne->data()->setMinValue(-5);
247 
248  m_pBCIOutputTwo->data()->setMaxValue(10);
249  m_pBCIOutputTwo->data()->setMinValue(0);
250 
251  m_pBCIOutputThree->data()->setMaxValue(10);
252  m_pBCIOutputThree->data()->setMinValue(0);
253 
254  m_pBCIOutputFour->data()->setMaxValue(1e-05);
255  m_pBCIOutputFour->data()->setMinValue(-1e-05);
256 
257  m_pBCIOutputFive->data()->setMaxValue(1e-05);
258  m_pBCIOutputFive->data()->setMinValue(-1e-05);
259 
260 // // log variance
261 // m_pBCIOutputOne->data()->setMaxValue(10);
262 // m_pBCIOutputOne->data()->setMinValue(-10);
263 
264 // m_pBCIOutputTwo->data()->setMaxValue(15);
265 // m_pBCIOutputTwo->data()->setMinValue(0);
266 
267 // m_pBCIOutputThree->data()->setMaxValue(15);
268 // m_pBCIOutputThree->data()->setMinValue(0);
269 
270 // m_pBCIOutputFour->data()->setMaxValue(10e-04);
271 // m_pBCIOutputFour->data()->setMinValue(-10e-04);
272 
273 // m_pBCIOutputFive->data()->setMaxValue(10e-04);
274 // m_pBCIOutputFive->data()->setMinValue(-10e-04);
275 
276  m_bIsRunning = true;
277 
278  QThread::start();
279 
280  return true;
281 }
282 
283 
284 //*************************************************************************************************************
285 
286 bool BCI::stop()
287 {
288  m_bIsRunning = false;
289 
290  // Get data buffers out of idle state if they froze in the acquire or release function
291  //In case the semaphore blocks the thread -> Release the QSemaphore and let it exit from the pop function (acquire statement)
292 
293  if(m_bProcessData) // Only clear if buffers have been initialised
294  {
295  m_pBCIBuffer_Sensor->releaseFromPop();
296  m_pBCIBuffer_Sensor->releaseFromPush();
297  // m_pBCIBuffer_Source->releaseFromPop();
298  // m_pBCIBuffer_Source->releaseFromPush();
299  }
300 
301  // Stop filling buffers with data from the inputs
302  m_bProcessData = false;
303 
304  // Reset trigger
305  m_bTriggerActivated = false;
306 
307  // Clear stream
308  m_outStreamDebug.close();
309  m_outStreamDebug.clear();
310 
311  // Hide feature visualization window
312  if(m_bDisplayFeatures)
313  m_BCIFeatureWindow->hide();
314 
315  // Delete all features and classification results
316  clearFeatures();
318 
319  return true;
320 }
321 
322 
323 //*************************************************************************************************************
324 
326 {
327  return _IAlgorithm;
328 }
329 
330 
331 //*************************************************************************************************************
332 
333 QString BCI::getName() const
334 {
335  return "BCI EEG";
336 }
337 
338 
339 //*************************************************************************************************************
340 
342 {
343  BCISetupWidget* setupWidget = new BCISetupWidget(this);//widget is later destroyed by CentralWidget - so it has to be created everytime new
344 
345  //init properties dialog
346  setupWidget->initGui();
347 
348  return setupWidget;
349 }
350 
351 
352 //*************************************************************************************************************
353 
355 {
356  QSharedPointer<NewRealTimeMultiSampleArray> pRTMSA = pMeasurement.dynamicCast<NewRealTimeMultiSampleArray>();
357  if(pRTMSA)
358  {
359  //Check if buffer initialized
360  if(!m_pBCIBuffer_Sensor)
361  m_pBCIBuffer_Sensor = CircularMatrixBuffer<double>::SPtr(new CircularMatrixBuffer<double>(64, pRTMSA->getNumChannels(), pRTMSA->getMultiArraySize()));
362 
363  // Load Fiff information on sensor level
364  if(!m_pFiffInfo_Sensor)
365  {
366  m_pFiffInfo_Sensor = pRTMSA->getFiffInfo();
367 
368  // Adjust working matrixes (sliding window and time between windows matrix) size so that the samples from the tmsi plugin stream fit in the matrix perfectly
369  int arraySize = pRTMSA->getMultiArraySize();
370  int modulo = int(m_pFiffInfo_Sensor->sfreq*m_dSlidingWindowSize) % arraySize;
371  int rows = m_slChosenFeatureSensor.size();
372  int cols = m_pFiffInfo_Sensor->sfreq*m_dSlidingWindowSize-modulo;
373  m_matSlidingWindowSensor.resize(rows, cols);
374 
375  m_matStimChannelSensor.resize(1,cols);
376 
377  modulo = int(m_pFiffInfo_Sensor->sfreq*m_dTimeBetweenWindows) % arraySize;
378  rows = m_slChosenFeatureSensor.size();
379  cols = m_pFiffInfo_Sensor->sfreq*m_dTimeBetweenWindows-modulo;
380  m_matTimeBetweenWindowsSensor.resize(rows, cols);
381 
382  m_matTimeBetweenWindowsStimSensor.resize(1,cols);
383 
384  // Build filter operator
385  double dCenterFreqNyq = (m_dFilterLowerBound+((m_dFilterUpperBound - m_dFilterLowerBound)/2))/(m_pFiffInfo_Sensor->sfreq/2);
386  double dBandwidthNyq = (m_dFilterUpperBound - m_dFilterLowerBound)/(m_pFiffInfo_Sensor->sfreq/2);
387  double dParksWidth = m_dParcksWidth/(m_pFiffInfo_Sensor->sfreq/2);
388 
389 // // Calculate needed fft length
390 // int exponent = ceil(log10(m_matSlidingWindowSensor.cols())/log10(2));
391 // int fftLength = pow(2,exponent+1);
392 
393  // Initialise filter operator
394  m_filterOperator = QSharedPointer<FilterData>(new FilterData(QString("BPF"),FilterData::BPF,m_iFilterOrder,dCenterFreqNyq,dBandwidthNyq,dParksWidth,m_matSlidingWindowSensor.cols()+m_iFilterOrder)); // letztes Argument muss 2er potenz sein - fft länge
395 
396  // Write filter coefficients to debug file
397  for(int i = 0; i<m_filterOperator->m_dCoeffA.cols(); i++)
398  m_outStreamDebug << m_filterOperator->m_dFFTCoeffA(0,i).real() <<"+" << m_filterOperator->m_dFFTCoeffA(0,i).imag() << "i " << endl;
399 
400  m_outStreamDebug << endl << endl;
401 
402  for(int i = 0; i<m_filterOperator->m_dCoeffA.cols(); i++)
403  m_outStreamDebug << m_filterOperator->m_dCoeffA(0,i) << endl;
404 
405  m_outStreamDebug << "---------------------------------------------------------------------" << endl;
406  }
407 
408  // Only process data when fiff info has been initialised in run() method
409  if(m_bProcessData)
410  {
411  MatrixXd t_mat(pRTMSA->getNumChannels(), pRTMSA->getMultiArraySize());
412 
413  for(unsigned char i = 0; i < pRTMSA->getMultiArraySize(); ++i)
414  t_mat.col(i) = pRTMSA->getMultiSampleArray()[i];
415 
416  m_pBCIBuffer_Sensor->push(&t_mat);
417  }
418  }
419 }
420 
421 
422 //*************************************************************************************************************
423 
425 {
426  cout<<"update source"<<endl;
427  QSharedPointer<RealTimeSourceEstimate> pRTSE = pMeasurement.dynamicCast<RealTimeSourceEstimate>();
428  if(pRTSE)
429  {
430  //Check if buffer initialized
431  if(!m_pBCIBuffer_Source)
432  m_pBCIBuffer_Source = CircularMatrixBuffer<double>::SPtr(new CircularMatrixBuffer<double>(64, pRTSE->getValue().size(), pRTSE->getArraySize()));
433 
434  if(m_bProcessData)
435  {
436  MatrixXd t_mat(pRTSE->getValue().size(), pRTSE->getArraySize());
437 
438  for(unsigned char i = 0; i < pRTSE->getArraySize(); ++i)
439  t_mat.col(i) = pRTSE->getStc().data.col(i);
440 
441  m_pBCIBuffer_Source->push(&t_mat);
442  }
443  }
444 }
445 
446 
447 //*************************************************************************************************************
448 
449 void BCI::applyMeanCorrectionConcurrently(QPair<int, RowVectorXd> &chdata)
450 {
451  chdata.second = chdata.second.array() - chdata.second.mean(); // Row vector needs to be transformed to a array in order to subract a scalar
452 }
453 
454 
455 //*************************************************************************************************************
456 
457 void BCI::applyFilterOperatorConcurrently(QPair<int, RowVectorXd> &chdata)
458 {
459  chdata.second = m_filterOperator->applyFFTFilter(chdata.second);
460 }
461 
462 
463 //*************************************************************************************************************
464 
465 QPair< int,QList<double> > BCI::applyFeatureCalcConcurrentlyOnSensorLevel(const QPair<int, RowVectorXd> &chdata)
466 {
467  RowVectorXd data = chdata.second;
468  QList<double> features;
469 
470  // TODO: Divide into subsignals
471  switch(m_iFeatureCalculationType)
472  {
473  case 0:
474  features << data.squaredNorm(); // Compute variance
475  break;
476  case 1:
477  features << abs(log10(data.squaredNorm())); // Compute log of variance
478  break;
479  default:
480  features << data.squaredNorm(); // Compute variance
481  break;
482  }
483 
484  return QPair< int,QList<double> >(chdata.first, features);
485 }
486 
487 
488 //*************************************************************************************************************
489 
491 {
492  return classificationBoundaryValue(featData);
493 }
494 
495 
496 //*************************************************************************************************************
497 
498 double BCI::classificationBoundaryValue(const QList<double> &featData)
499 {
500  double return_val = 0;
501 
502  if(featData.size() == m_vLoadedSensorBoundary[1].size())
503  {
504  VectorXd feat_temp(featData.size());
505  for(int i = 0; i<featData.size(); i++)
506  feat_temp(i) = featData.at(i);
507 
508  return_val = m_vLoadedSensorBoundary[0](0) + m_vLoadedSensorBoundary[1].dot(feat_temp);
509  }
510 
511  return return_val;
512 }
513 
514 
515 //*************************************************************************************************************
516 
518 {
519  m_qMutex.lock();
520  m_lFeaturesSensor.clear();
521  m_qMutex.unlock();
522 }
523 
524 
525 //*************************************************************************************************************
526 
528 {
529  m_qMutex.lock();
530  m_lClassResultsSensor.clear();
531  m_qMutex.unlock();
532 }
533 
534 
535 //*************************************************************************************************************
536 
537 bool BCI::hasThresholdArtefact(const QList< QPair<int,RowVectorXd> > &data)
538 {
539  // Perform simple threshold artefact reduction
540  double max = 0;
541  double min = 0;
542 
543  if(m_bUseArtefactThresholdReduction)
544  {
545  // find min max in current m_matSlidingWindowSensor/qlMatrixRows after mean was subtracted
546  for(int i = 0; i<data.size(); i++)
547  {
548  if(data.at(i).second.maxCoeff() > max)
549  max = data.at(i).second.maxCoeff();
550 
551  if(data.at(i).second.minCoeff() < min)
552  min = data.at(i).second.minCoeff();
553  }
554  }
555 
556 // cout<<"max: "<<max<<endl;
557 // cout<<"min: "<<min<<endl;
558 
559  if(max<m_dThresholdValue*1e-06 && min>m_dThresholdValue*-1e-06) // If max is outside the threshold -> completley discard m_matSlidingWindowSensor
560  return false;
561  else
562  {
563  //cout<<"Rejected artefact"<<endl;
564  return true;
565  }
566 }
567 
568 
569 //*************************************************************************************************************
570 
571 bool BCI::lookForTrigger(const MatrixXd &data)
572 {
573  // Check if capacitive touch trigger signal was received - Note that there can also be "beep" triggers in the received data, which are only 1 sample wide -> therefore look for 2 samples with a value of 254 each
574  for(int i = 0; i<data.cols()-1; i++)
575  {
576  if(data(0,i) == 254 && data(0,i+1) == 254) // data - corresponds with channel 136 which is the trigger channel
577  return true;
578  }
579 
580  return false;
581 }
582 
583 
584 //*************************************************************************************************************
585 
586 void BCI::run()
587 {
588  while(m_bIsRunning)
589  {
590  // Decide which data to use - sensor or source level data
591  if(m_bUseSensorData)
593  else
594  if(m_bUseSourceData)
596  }
597 }
598 
599 
600 //*************************************************************************************************************
601 
603 {
604  // Wait for fiff Info if not yet received - this is needed because we have to wait until the buffers are firstly initiated in the update functions
605  while(!m_pFiffInfo_Sensor)
606  msleep(10);
607 
608  // Start filling buffers with data from the inputs
609  m_bProcessData = true;
610 
611  // Sensor level: Fill matrices with data
612  if(m_bFillSensorWindowFirstTime) // Sensor level: Fill m_matSlidingWindowSensor with data for the first time
613  {
614  if(m_iTBWIndexSensor < m_matSlidingWindowSensor.cols())
615  {
616  //cout<<"About to pop matrix"<<endl;
617  MatrixXd t_mat = m_pBCIBuffer_Sensor->pop();
618  //cout<<"poped matrix"<<endl;
619 
620  // Get only the rows from the matrix which correspond with the selected features, namely electrodes on sensor level and destrieux clustered regions on source level
621  for(int i = 0; i < m_matSlidingWindowSensor.rows(); i++)
622  m_matSlidingWindowSensor.block(i, m_iTBWIndexSensor, 1, t_mat.cols()) = t_mat.block(m_mapElectrodePinningScheme[m_slChosenFeatureSensor.at(i)], 0, 1, t_mat.cols());
623 
624  m_matStimChannelSensor.block(0, m_iTBWIndexSensor, 1, t_mat.cols()) = t_mat.block(136, 0, 1, t_mat.cols());
625 
626  m_iTBWIndexSensor = m_iTBWIndexSensor + t_mat.cols();
627  }
628  else // m_matSlidingWindowSensor is full for the first time
629  {
630  m_iTBWIndexSensor = 0;
631  m_bFillSensorWindowFirstTime = false;
632  }
633  }
634  else // Sensor level: Fill m_matTimeBetweenWindowsSensor matrix until full -> Then -> recalculate m_matSlidingWindowSensor, calculate features and classify
635  {
636  if(m_iTBWIndexSensor < m_matTimeBetweenWindowsSensor.cols())
637  {
638  //cout<<"About to pop matrix"<<endl;
639  MatrixXd t_mat = m_pBCIBuffer_Sensor->pop();
640  //cout<<"poped matrix"<<endl;
641 
642  // Get only the rows from the matrix which correspond with the selected features, namely electrodes on sensor level and destrieux clustered regions on source level
643  for(int i = 0; i < m_matTimeBetweenWindowsSensor.rows(); i++)
644  m_matTimeBetweenWindowsSensor.block(i, m_iTBWIndexSensor, 1, t_mat.cols()) = t_mat.block(m_mapElectrodePinningScheme[m_slChosenFeatureSensor.at(i)], 0, 1, t_mat.cols());
645 
646  m_matTimeBetweenWindowsStimSensor.block(0, m_iTBWIndexSensor, 1, t_mat.cols()) = t_mat.block(136, 0, 1, t_mat.cols());
647 
648  m_iTBWIndexSensor = m_iTBWIndexSensor + t_mat.cols();
649  }
650  else // Recalculate m_matSlidingWindowSensor -> Calculate features, classify and store results
651  {
652  // ----1---- Recalculate m_matSlidingWindowSensor
653  //cout<<"----1----"<<endl;
654  //Move block from the right to the left -> use eval() to resolve Eigens aliasing problem
655  m_matSlidingWindowSensor.block(0, 0, m_matSlidingWindowSensor.rows(), m_matSlidingWindowSensor.cols()-m_matTimeBetweenWindowsSensor.cols()) = m_matSlidingWindowSensor.block(0, m_matTimeBetweenWindowsSensor.cols(), m_matSlidingWindowSensor.rows(), m_matSlidingWindowSensor.cols()-m_matTimeBetweenWindowsSensor.cols()).eval();
656  m_matStimChannelSensor.block(0, 0, m_matStimChannelSensor.rows(), m_matStimChannelSensor.cols()-m_matTimeBetweenWindowsStimSensor.cols()) = m_matStimChannelSensor.block(0, m_matTimeBetweenWindowsStimSensor.cols(), m_matStimChannelSensor.rows(), m_matStimChannelSensor.cols()-m_matTimeBetweenWindowsStimSensor.cols()).eval();
657 
658  // push m_matTimeBetweenWindowsSensor from the right
659  m_matSlidingWindowSensor.block(0, m_matSlidingWindowSensor.cols()-m_matTimeBetweenWindowsSensor.cols(), m_matTimeBetweenWindowsSensor.rows(), m_matTimeBetweenWindowsSensor.cols()) = m_matTimeBetweenWindowsSensor;
660  m_matStimChannelSensor.block(0, m_matStimChannelSensor.cols()-m_matTimeBetweenWindowsStimSensor.cols(), m_matTimeBetweenWindowsStimSensor.rows(), m_matTimeBetweenWindowsStimSensor.cols()) = m_matTimeBetweenWindowsStimSensor;
661 
662  //cout<<m_matStimChannelSensor;
663 
664  // Test if data is correctly streamed to this plugin
665  if(m_slChosenFeatureSensor.contains("TEST"))
666  {
667  cout<<"Recalculate matrix"<<endl;
668 
669  for(int i = 0; i<m_matSlidingWindowSensor.cols() ; i++)
670  cout << m_matSlidingWindowSensor(m_matSlidingWindowSensor.rows()-1,i) <<endl;
671 
672 // for(int i = 1; i<m_matSlidingWindowSensor.cols() ; i++)
673 // {
674 // cout << m_matSlidingWindowSensor(m_matSlidingWindowSensor.rows()-1,i-1) <<endl;
675 // if(m_matSlidingWindowSensor(m_matSlidingWindowSensor.rows()-1,i) - m_matSlidingWindowSensor(m_matSlidingWindowSensor.rows()-1,i-1) != 1)
676 // cout<<"Sequence error while streaming from tmsi plugin at position: "<<i<<endl;
677 // }
678  }
679 
680  // ----2---- Transform matrix into QList structure, so that QTConcurrent can handle it properly
681  //cout<<"----2----"<<endl;
682  QList< QPair<int,RowVectorXd> > qlMatrixRows;
683  for(int i = 0; i< m_matSlidingWindowSensor.rows(); i++)
684  qlMatrixRows << QPair<int,RowVectorXd>(i, m_matSlidingWindowSensor.row(i));
685 
686  int iNumberOfFeatures = qlMatrixRows.size();
687 
688  // ----3---- Subtract mean in m_matSlidingWindowSensor concurrently using map()
689  //cout<<"----3----"<<endl;
690  if(m_bSubtractMean)
691  {
692  QFuture<void> futureMean = QtConcurrent::map(qlMatrixRows,[this](QPair<int,RowVectorXd>& rowdata) {
694  });
695 
696  futureMean.waitForFinished();
697  }
698 
699  // ----4---- Do simple threshold artefact reduction
700  //cout<<"----4----"<<endl;
701  if(hasThresholdArtefact(qlMatrixRows) == false)
702  {
703  // Look for trigger flag
704  if(lookForTrigger(m_matStimChannelSensor) && !m_bTriggerActivated)
705  {
706  // cout << "Trigger activated" << endl;
707  //QFuture<void> future = QtConcurrent::run(Beep, 450, 700);
708  m_bTriggerActivated = true;
709  }
710 
711  // ----5---- Filter data in m_matSlidingWindowSensor concurrently using map()
712  //cout<<"----5----"<<endl;
713  // TODO: work only on qlMatrixRows -> filteredRows doesnt need to be created -> more efficient
714  QList< QPair<int,RowVectorXd> > filteredRows = qlMatrixRows;
715 
716  if(m_bUseFilter)
717  {
718  QFuture<void> futureFilter = QtConcurrent::map(filteredRows,[this](QPair<int,RowVectorXd>& chdata) {
720  });
721 
722  futureFilter.waitForFinished();
723  }
724 
725 // // Write data before and after filtering to debug file
726 // for(int i=0; i<qlMatrixRows.size(); i++)
727 // m_outStreamDebug << "qlMatrixRows at row " << i << ": " << qlMatrixRows.at(i).second << "\n";
728 
729 // m_outStreamDebug<<endl;
730 
731 // for(int i=0; i<filteredRows.size(); i++)
732 // m_outStreamDebug << "filteredRows at row " << i << ": " << filteredRows.at(i).second << "\n";
733 
734 // m_outStreamDebug << endl << "---------------------------------------------------------" << endl << endl;
735 
736 // // Write filtered data continously to file
737 // for(int i = 0; i<filteredRows.at(0).second.size() ;i++)
738 // m_outStreamDebug << filteredRows.at(0).second(i)<<endl;
739 
740  // ----6---- Calculate features concurrently using mapped()
741  //cout<<"----6----"<<endl;
742  std::function< QPair<int,QList<double> > (QPair<int,RowVectorXd>&)> applyOpsFeatures = [this](QPair<int,RowVectorXd>& chdata) -> QPair< int,QList<double> > {
744  };
745 
746  QFuture< QPair< int,QList<double> > > futureCalculatedFeatures = QtConcurrent::mapped(filteredRows.begin(), filteredRows.end(), applyOpsFeatures);
747 
748  futureCalculatedFeatures.waitForFinished();
749 
750  m_iNumberOfCalculatedFeatures++;
751 
752  // ----7---- Store features
753  //cout<<"----7----"<<endl;
754  m_lFeaturesSensor.append(futureCalculatedFeatures.results());
755 
756  // ----8---- If enough features (windows) have been calculated (processed) -> classify all features and average results
757  //cout<<"----8----"<<endl;
758  if(m_iNumberOfCalculatedFeatures == m_iNumberFeatures)
759  {
760  // Transform m_lFeaturesSensor into an easier file structure -> create feature points
761  QList< QList<double> > lFeaturesSensor_new;
762 
763  for(int i = 0; i<m_lFeaturesSensor.size()-iNumberOfFeatures+1; i = i + iNumberOfFeatures) // iterate over QPair feature List
764  for(int z = 0; z<m_lFeaturesSensor.at(0).second.size(); z++) // iterate over number of sub signals
765  {
766  QList<double> temp;
767  for(int t = 0; t<iNumberOfFeatures; t++) // iterate over chosen features (electrodes)
768  temp.append(m_lFeaturesSensor.at(i+t).second.at(z));
769  lFeaturesSensor_new.append(temp);
770  }
771 
772 // //Check sizes
773 // cout<<"lFeaturesSensor_new.size()"<<lFeaturesSensor_new.size()<<endl;
774 // cout<<"lFeaturesSensor_new.at(0).size()"<<lFeaturesSensor_new.at(0).size()<<endl;
775 // cout<<"m_lFeaturesSensor.size()"<<m_lFeaturesSensor.size()<<endl;
776 
777  // Display features
778  if(m_bDisplayFeatures)
779  emit paintFeatures((MyQList)lFeaturesSensor_new, m_bTriggerActivated);
780 
781  // Reset trigger
782  m_bTriggerActivated = false;
783 
784  // ----9---- Classify features concurrently using mapped() ----------
785  //cout<<"----9----"<<endl;
786 
787  std::function<double (QList<double>&)> applyOpsClassification = [this](QList<double>& featData){
789  };
790 
791  QFuture<double> futureClassificationResults = QtConcurrent::mapped(lFeaturesSensor_new.begin(), lFeaturesSensor_new.end(), applyOpsClassification);
792 
793  futureClassificationResults.waitForFinished();
794 
795  // ----10---- Generate final classification result -> average all classification results
796  //cout<<"----10----"<<endl;
797  double dfinalResult = 0;
798 
799  for(int i = 0; i<futureClassificationResults.resultCount() ;i++)
800  dfinalResult += futureClassificationResults.resultAt(i);
801 
802  dfinalResult = dfinalResult/futureClassificationResults.resultCount();
803  cout << "dfinalResult: " << dfinalResult << endl << endl;
804 
805  // ----11---- Store final result
806  //cout<<"----11----"<<endl;
807  m_lClassResultsSensor.append(dfinalResult);
808 
809  // ----12---- Send result to the output stream, i.e. which is connected to the triggerbox
810  //cout<<"----12----"<<endl;
811  VectorXd variances(iNumberOfFeatures);
812  variances.setZero();
813 
814  for(int i = 0; i<lFeaturesSensor_new.size(); i++)
815  for(int t = 0; t<iNumberOfFeatures; t++)
816  variances(t) = variances(t) + lFeaturesSensor_new.at(i).at(t);
817 
818  variances = variances/lFeaturesSensor_new.size();
819 
820  m_pBCIOutputOne->data()->setValue(dfinalResult);
821  m_pBCIOutputTwo->data()->setValue(variances(0));
822  m_pBCIOutputThree->data()->setValue(variances(1));
823 
824  for(int i = 0; i<filteredRows.at(0).second.cols() ; i++)
825  {
826  m_pBCIOutputFour->data()->setValue(filteredRows.at(0).second(0,i));
827  m_pBCIOutputFive->data()->setValue(filteredRows.at(1).second(0,i));
828  }
829 
830  // Clear classifications
831  clearFeatures();
832 
833  // Reset counter
834  m_iNumberOfCalculatedFeatures = 0;
835  } // End if enough features (windows) have been calculated (processed)
836  } // End if artefact reduction
837  else
838  {
839  // If trial has been rejected -> plot zeros as result and the filtered electrode channel
840  m_pBCIOutputOne->data()->setValue(0);
841  m_pBCIOutputTwo->data()->setValue(0);
842  m_pBCIOutputThree->data()->setValue(0);
843 
844  QList< QPair<int,RowVectorXd> > filteredRows = qlMatrixRows;
845 
846  if(m_bUseFilter)
847  {
848  QFuture<void> futureFilter = QtConcurrent::map(filteredRows,[this](QPair<int,RowVectorXd>& chdata) {
850  });
851 
852  futureFilter.waitForFinished();
853  }
854 
855  for(int i = 0; i<filteredRows.at(0).second.cols() ; i++)
856  {
857  m_pBCIOutputFour->data()->setValue(filteredRows.at(0).second(0,i));
858  m_pBCIOutputFive->data()->setValue(filteredRows.at(1).second(0,i));
859  }
860  }
861 
862  m_iTBWIndexSensor = 0;
863  }
864  }
865 }
866 
867 
868 //*************************************************************************************************************
869 
871 {
872 }
double classificationBoundaryValue(const QList< double > &featData)
Definition: bci.cpp:498
virtual IPlugin::PluginType getType() const
Definition: bci.cpp:325
QPair< int, QList< double > > applyFeatureCalcConcurrentlyOnSensorLevel(const QPair< int, RowVectorXd > &chdata)
Definition: bci.cpp:465
void applyFilterOperatorConcurrently(QPair< int, RowVectorXd > &chdata)
Definition: bci.cpp:457
OutputConnectorList m_outputConnectors
Definition: IPlugin.h:222
virtual void run()
Definition: bci.cpp:586
The BCIFeatureWindow class provides a visualization tool for calculated features. ...
The TMSISetupWidget class provides the TMSI configuration window.
void BCIOnSourceLevel()
Definition: bci.cpp:870
QSharedPointer< CircularMatrixBuffer > SPtr
void updateSource(XMEASLIB::NewMeasurement::SPtr pMeasurement)
Definition: bci.cpp:424
QSharedPointer< FiffInfo > SPtr
Definition: fiff_info.h:99
virtual void unload()
Definition: bci.cpp:184
Definition: bci.h:75
void BCIOnSensorLevel()
Definition: bci.cpp:602
double applyClassificationCalcConcurrentlyOnSensorLevel(QList< double > &featData)
Definition: bci.cpp:490
void push(const Matrix< _Tp, Dynamic, Dynamic > *pMatrix)
virtual bool stop()
Definition: bci.cpp:286
bool lookForTrigger(const MatrixXd &data)
Definition: bci.cpp:571
void applyMeanCorrectionConcurrently(QPair< int, RowVectorXd > &rowdata)
Definition: bci.cpp:449
virtual void init()
Definition: bci.cpp:99
Matrix< _Tp, Dynamic, Dynamic > pop()
void clearFeatures()
Definition: bci.cpp:517
bool hasThresholdArtefact(const QList< QPair< int, RowVectorXd > > &data)
Definition: bci.cpp:537
virtual ~BCI()
Definition: bci.cpp:78
virtual QWidget * setupWidget()
Definition: bci.cpp:341
Real-time source estimate measurement.
QSharedPointer< NewMeasurement > SPtr
virtual QString getName() const
Definition: bci.cpp:333
void clearClassifications()
Definition: bci.cpp:527
virtual QSharedPointer< IPlugin > clone() const
Definition: bci.cpp:90
static QSharedPointer< PluginOutputData< T > > create(IPlugin *parent, const QString &name, const QString &descr)
virtual bool start()
Definition: bci.cpp:192
Contains the declaration of the BCI class.
InputConnectorList m_inputConnectors
Definition: IPlugin.h:221
void updateSensor(XMEASLIB::NewMeasurement::SPtr pMeasurement)
Definition: bci.cpp:354
static QSharedPointer< PluginInputData< T > > create(IPlugin *parent, const QString &name, const QString &descr)
The RealTimeMultiSampleArrayNew class is the base class of every RealTimeMultiSampleArrayNew Measurem...