50 using namespace FIFFLIB;
51 using namespace Eigen;
72 printf(
"\tError during fiff setup raw read.\n");
87 printf(
"\tError during fiff setup raw read.\n");
96 : file(p_FiffRawData.file)
97 , info(p_FiffRawData.info)
98 , first_samp(p_FiffRawData.first_samp)
99 , last_samp(p_FiffRawData.last_samp)
100 , cals(p_FiffRawData.cals)
101 , rawdir(p_FiffRawData.rawdir)
102 , proj(p_FiffRawData.proj)
103 , comp(p_FiffRawData.comp)
120 cals = RowVectorXd();
127 #include <QElapsedTimer>
133 const RowVectorXi& sel,
136 bool projAvailable =
true;
138 if (this->
proj.size() == 0) {
140 projAvailable =
false;
157 printf(
"No data in this range %d ... %d = %9.3f ... %9.3f secs...", from, to, ((
float)from)/this->
info.
sfreq, ((
float)to)/this->info.sfreq);
168 typedef Eigen::Triplet<double> T;
169 std::vector<T> tripletList;
170 tripletList.reserve(nchan);
171 for(i = 0; i < nchan; ++i)
172 tripletList.push_back(T(i, i, this->
cals[i]));
174 SparseMatrix<double> cal(nchan, nchan);
175 cal.setFromTriplets(tripletList.begin(), tripletList.end());
182 data = MatrixXd(nchan, to-from+1);
184 if (projAvailable || this->
comp.
kind != -1)
187 mult_full = this->
comp.
data->data*cal;
189 mult_full = this->
proj*cal;
196 data = MatrixXd(sel.size(),to-from+1);
199 MatrixXd selVect(sel.size(), nchan);
203 if (!projAvailable && this->
comp.
kind == -1)
206 tripletList.reserve(sel.size());
207 for(i = 0; i < sel.size(); ++i)
208 tripletList.push_back(T(i, i, this->
cals[sel[i]]));
209 cal = SparseMatrix<double>(sel.size(), sel.size());
210 cal.setFromTriplets(tripletList.begin(), tripletList.end());
216 qDebug() <<
"This has to be debugged! #1";
217 for( i = 0; i < sel.size(); ++i)
218 selVect.row(i) = this->
comp.
data->data.block(sel[i],0,1,nchan);
219 mult_full = selVect*cal;
223 for( i = 0; i < sel.size(); ++i)
224 selVect.row(i) = this->
proj.block(sel[i],0,1,nchan);
226 mult_full = selVect*cal;
230 qDebug() <<
"This has to be debugged! #3";
231 for( i = 0; i < sel.size(); ++i)
232 selVect.row(i) = this->
proj.block(sel[i],0,1,nchan);
234 mult_full = selVect*this->
comp.
data->data*cal;
243 tripletList.reserve(mult_full.rows()*mult_full.cols());
244 for(i = 0; i < mult_full.rows(); ++i)
245 for(
k = 0;
k < mult_full.cols(); ++
k)
246 if(mult_full(i,
k) != 0)
247 tripletList.push_back(T(i,
k, mult_full(i,
k)));
249 SparseMatrix<double> mult(mult_full.rows(),mult_full.cols());
250 if(tripletList.size() > 0)
251 mult.setFromTriplets(tripletList.begin(), tripletList.end());
255 if (!this->
file->device()->isOpen())
257 if (!this->
file->device()->open(QIODevice::ReadOnly))
259 printf(
"Cannot open file %s",this->
info.
filename.toUtf8().constData());
268 MatrixXd one, newData, tmp_data;
271 fiff_int_t first_pick, last_pick, picksamp;
278 if (thisRawDir.
last > from)
280 if (thisRawDir.
ent->kind == -1)
288 one.resize(nchan,thisRawDir.
nsamp);
290 one.resize(sel.cols(),thisRawDir.
nsamp);
296 fid->read_tag(t_pTag, thisRawDir.
ent->pos);
301 if (mult.cols() == 0)
305 if (t_pTag->type == FIFFT_DAU_PACK16)
306 one = cal*(Map< MatrixDau16 >( t_pTag->toDauPack16(),nchan, thisRawDir.
nsamp)).cast<double>();
307 else if(t_pTag->type == FIFFT_INT)
308 one = cal*(Map< MatrixXi >( t_pTag->toInt(),nchan, thisRawDir.
nsamp)).cast<double>();
309 else if(t_pTag->type == FIFFT_FLOAT)
310 one = cal*(Map< MatrixXf >( t_pTag->toFloat(),nchan, thisRawDir.
nsamp)).cast<double>();
311 else if(t_pTag->type == FIFFT_SHORT)
312 one = cal*(Map< MatrixShort >( t_pTag->toShort(),nchan, thisRawDir.
nsamp)).cast<double>();
314 printf(
"Data Storage Format not known yet [1]!! Type: %d\n", t_pTag->type);
319 newData.resize(sel.cols(), thisRawDir.
nsamp);
321 if (t_pTag->type == FIFFT_DAU_PACK16)
323 tmp_data = (Map< MatrixDau16 > ( t_pTag->toDauPack16(),nchan, thisRawDir.
nsamp)).cast<double>();
325 for(r = 0; r < sel.size(); ++r)
326 newData.block(r,0,1,thisRawDir.
nsamp) = tmp_data.block(sel[r],0,1,thisRawDir.
nsamp);
328 else if(t_pTag->type == FIFFT_INT)
330 tmp_data = (Map< MatrixXi >( t_pTag->toInt(),nchan, thisRawDir.
nsamp)).cast<double>();
332 for(r = 0; r < sel.size(); ++r)
333 newData.block(r,0,1,thisRawDir.
nsamp) = tmp_data.block(sel[r],0,1,thisRawDir.
nsamp);
335 else if(t_pTag->type == FIFFT_FLOAT)
337 tmp_data = (Map< MatrixXf > ( t_pTag->toFloat(),nchan, thisRawDir.
nsamp)).cast<double>();
339 for(r = 0; r < sel.size(); ++r)
340 newData.block(r,0,1,thisRawDir.
nsamp) = tmp_data.block(sel[r],0,1,thisRawDir.
nsamp);
342 else if(t_pTag->type == FIFFT_SHORT)
344 tmp_data = (Map< MatrixShort > ( t_pTag->toShort(),nchan, thisRawDir.
nsamp)).cast<double>();
346 for(r = 0; r < sel.size(); ++r)
347 newData.block(r,0,1,thisRawDir.
nsamp) = tmp_data.block(sel[r],0,1,thisRawDir.
nsamp);
351 printf(
"Data Storage Format not known yet [2]!! Type: %d\n", t_pTag->type);
359 if (t_pTag->type == FIFFT_DAU_PACK16)
360 one = mult*(Map< MatrixDau16 >( t_pTag->toDauPack16(),nchan, thisRawDir.
nsamp)).cast<double>();
361 else if(t_pTag->type == FIFFT_INT)
362 one = mult*(Map< MatrixXi >( t_pTag->toInt(),nchan, thisRawDir.
nsamp)).cast<double>();
363 else if(t_pTag->type == FIFFT_FLOAT)
364 one = mult*(Map< MatrixXf >( t_pTag->toFloat(),nchan, thisRawDir.
nsamp)).cast<double>();
366 printf(
"Data Storage Format not known yet [3]!! Type: %d\n", t_pTag->type);
372 if (to >= thisRawDir.
last && from <= thisRawDir.
first)
378 last_pick = thisRawDir.
nsamp - 1;
382 else if (from > thisRawDir.
first)
384 first_pick = from - thisRawDir.
first;
385 if(to < thisRawDir.
last)
391 last_pick = thisRawDir.
nsamp + to - thisRawDir.
last - 1;
400 last_pick = thisRawDir.
nsamp - 1;
411 last_pick = to - thisRawDir.
first;
418 picksamp = last_pick - first_pick + 1;
422 qDebug() <<
"first_pick: " << first_pick;
423 qDebug() <<
"last_pick: " << last_pick;
424 qDebug() <<
"picksamp: " << picksamp;
432 data.block(0,dest,data.rows(),picksamp) = one.block(0, first_pick, data.rows(), picksamp);
440 if (thisRawDir.
last >= to)
447 if (!this->
file->device()->isOpen()) {
448 this->
file->device()->close();
451 times = MatrixXd(1, to-from+1);
453 for (i = 0; i < times.cols(); ++i)
454 times(0, i) = ((float)(from+i)) / this->
info.
sfreq;
463 SparseMatrix<double>& multSegment,
466 const RowVectorXi& sel,
469 bool projAvailable =
true;
471 if (this->
proj.size() == 0) {
473 projAvailable =
false;
490 printf(
"No data in this range\n");
501 typedef Eigen::Triplet<double> T;
502 std::vector<T> tripletList;
503 tripletList.reserve(nchan);
504 for(i = 0; i < nchan; ++i)
505 tripletList.push_back(T(i, i, this->
cals[i]));
507 SparseMatrix<double> cal(nchan, nchan);
508 cal.setFromTriplets(tripletList.begin(), tripletList.end());
515 data = MatrixXd(nchan, to-from+1);
517 if (projAvailable || this->
comp.
kind != -1)
520 mult_full = this->
comp.
data->data*cal;
522 mult_full = this->
proj*cal;
529 data = MatrixXd(sel.size(),to-from+1);
532 MatrixXd selVect(sel.size(), nchan);
536 if (!projAvailable && this->
comp.
kind == -1)
539 tripletList.reserve(sel.size());
540 for(i = 0; i < sel.size(); ++i)
541 tripletList.push_back(T(i, i, this->
cals[sel[i]]));
542 cal = SparseMatrix<double>(sel.size(), sel.size());
543 cal.setFromTriplets(tripletList.begin(), tripletList.end());
549 qDebug() <<
"This has to be debugged! #1";
550 for( i = 0; i < sel.size(); ++i)
551 selVect.row(i) = this->
comp.
data->data.block(sel[i],0,1,nchan);
552 mult_full = selVect*cal;
556 for( i = 0; i < sel.size(); ++i)
557 selVect.row(i) = this->
proj.block(sel[i],0,1,nchan);
559 mult_full = selVect*cal;
563 qDebug() <<
"This has to be debugged! #3";
564 for( i = 0; i < sel.size(); ++i)
565 selVect.row(i) = this->
proj.block(sel[i],0,1,nchan);
567 mult_full = selVect*this->
comp.
data->data*cal;
576 tripletList.reserve(mult_full.rows()*mult_full.cols());
577 for(i = 0; i < mult_full.rows(); ++i)
578 for(
k = 0;
k < mult_full.cols(); ++
k)
579 if(mult_full(i,
k) != 0)
580 tripletList.push_back(T(i,
k, mult_full(i,
k)));
582 SparseMatrix<double> mult(mult_full.rows(),mult_full.cols());
583 if(tripletList.size() > 0)
584 mult.setFromTriplets(tripletList.begin(), tripletList.end());
590 if (!this->
file->device()->isOpen())
592 if (!this->
file->device()->open(QIODevice::ReadOnly))
594 printf(
"Cannot open file %s",this->
info.
filename.toUtf8().constData());
604 fiff_int_t first_pick, last_pick, picksamp;
611 if (thisRawDir.
last > from)
613 if (thisRawDir.
ent->kind == -1)
621 one.resize(nchan,thisRawDir.
nsamp);
623 one.resize(sel.cols(),thisRawDir.
nsamp);
630 fid->read_tag(t_pTag, thisRawDir.
ent->pos);
635 if (mult.cols() == 0)
639 if (t_pTag->type == FIFFT_DAU_PACK16)
640 one = cal*(Map< MatrixDau16 >( t_pTag->toDauPack16(),nchan, thisRawDir.
nsamp)).cast<double>();
641 else if(t_pTag->type == FIFFT_INT)
642 one = cal*(Map< MatrixXi >( t_pTag->toInt(),nchan, thisRawDir.
nsamp)).cast<double>();
643 else if(t_pTag->type == FIFFT_FLOAT)
644 one = cal*(Map< MatrixXf >( t_pTag->toFloat(),nchan, thisRawDir.
nsamp)).cast<double>();
645 else if(t_pTag->type == FIFFT_SHORT)
646 one = cal*(Map< MatrixShort >( t_pTag->toShort(),nchan, thisRawDir.
nsamp)).cast<double>();
648 printf(
"Data Storage Format not known yet [1]!! Type: %d\n", t_pTag->type);
654 MatrixXd newData(sel.cols(), thisRawDir.
nsamp);
656 if (t_pTag->type == FIFFT_DAU_PACK16)
658 MatrixXd tmp_data = (Map< MatrixDau16 > ( t_pTag->toDauPack16(),nchan, thisRawDir.
nsamp)).cast<double>();
660 for(r = 0; r < sel.size(); ++r)
661 newData.block(r,0,1,thisRawDir.
nsamp) = tmp_data.block(sel[r],0,1,thisRawDir.
nsamp);
663 else if(t_pTag->type == FIFFT_INT)
665 MatrixXd tmp_data = (Map< MatrixXi >( t_pTag->toInt(),nchan, thisRawDir.
nsamp)).cast<double>();
667 for(r = 0; r < sel.size(); ++r)
668 newData.block(r,0,1,thisRawDir.
nsamp) = tmp_data.block(sel[r],0,1,thisRawDir.
nsamp);
670 else if(t_pTag->type == FIFFT_FLOAT)
672 MatrixXd tmp_data = (Map< MatrixXf > ( t_pTag->toFloat(),nchan, thisRawDir.
nsamp)).cast<double>();
674 for(r = 0; r < sel.size(); ++r)
675 newData.block(r,0,1,thisRawDir.
nsamp) = tmp_data.block(sel[r],0,1,thisRawDir.
nsamp);
677 else if(t_pTag->type == FIFFT_SHORT)
679 MatrixXd tmp_data = (Map< MatrixShort > ( t_pTag->toShort(),nchan, thisRawDir.
nsamp)).cast<double>();
681 for(r = 0; r < sel.size(); ++r)
682 newData.block(r,0,1,thisRawDir.
nsamp) = tmp_data.block(sel[r],0,1,thisRawDir.
nsamp);
686 printf(
"Data Storage Format not known yet [2]!! Type: %d\n", t_pTag->type);
694 if (t_pTag->type == FIFFT_DAU_PACK16)
695 one = mult*(Map< MatrixDau16 >( t_pTag->toDauPack16(),nchan, thisRawDir.
nsamp)).cast<double>();
696 else if(t_pTag->type == FIFFT_INT)
697 one = mult*(Map< MatrixXi >( t_pTag->toInt(),nchan, thisRawDir.
nsamp)).cast<double>();
698 else if(t_pTag->type == FIFFT_FLOAT)
699 one = mult*(Map< MatrixXf >( t_pTag->toFloat(),nchan, thisRawDir.
nsamp)).cast<double>();
701 printf(
"Data Storage Format not known yet [3]!! Type: %d\n", t_pTag->type);
707 if (to >= thisRawDir.
last && from <= thisRawDir.
first)
713 last_pick = thisRawDir.
nsamp - 1;
717 else if (from > thisRawDir.
first)
719 first_pick = from - thisRawDir.
first;
720 if(to < thisRawDir.
last)
726 last_pick = thisRawDir.
nsamp + to - thisRawDir.
last - 1;
735 last_pick = thisRawDir.
nsamp - 1;
746 last_pick = to - thisRawDir.
first;
753 picksamp = last_pick - first_pick + 1;
757 qDebug() <<
"first_pick: " << first_pick;
758 qDebug() <<
"last_pick: " << last_pick;
759 qDebug() <<
"picksamp: " << picksamp;
767 data.block(0,dest,data.rows(),picksamp) = one.block(0, first_pick, data.rows(), picksamp);
775 if (thisRawDir.
last >= to)
787 if (!this->
file->device()->isOpen()) {
788 this->
file->device()->close();
791 times = MatrixXd(1, to-from+1);
793 for (i = 0; i < times.cols(); ++i)
794 times(0, i) = ((float)(from+i)) / this->
info.
sfreq;
805 const RowVectorXi& sel)
const