v2.0.0
Loading...
Searching...
No Matches
fiff_raw_data.cpp
Go to the documentation of this file.
1//=============================================================================================================
21
22//=============================================================================================================
23// INCLUDES
24//=============================================================================================================
25
26#include "fiff_raw_data.h"
27#include "fiff_events.h"
28#include "fiff_tag.h"
29#include "fiff_stream.h"
30#include "cstdlib"
31
32#include <stdexcept>
33//=============================================================================================================
34// USED NAMESPACES
35//=============================================================================================================
36
37using namespace FIFFLIB;
38using namespace Eigen;
39
40//=============================================================================================================
41// DEFINE MEMBER METHODS
42//=============================================================================================================
43
49
50//=============================================================================================================
51
52FiffRawData::FiffRawData(QIODevice& p_IODevice)
53: first_samp(-1)
54, last_samp(-1)
55{
56 //setup FiffRawData object
57 if (!FiffStream::setup_read_raw(p_IODevice, *this)) {
58 throw std::runtime_error("Error during fiff setup raw read");
59 }
60}
61
62//=============================================================================================================
63
64FiffRawData::FiffRawData(QIODevice& p_IODevice, bool b_littleEndian)
65: first_samp(-1)
66, last_samp(-1)
67{
68 //setup FiffRawData object
69 if (!FiffStream::setup_read_raw(p_IODevice, *this, false, b_littleEndian)) {
70 throw std::runtime_error("Error during fiff setup raw read");
71 }
72}
73
74//=============================================================================================================
75
77: file(p_FiffRawData.file)
78, info(p_FiffRawData.info)
79, first_samp(p_FiffRawData.first_samp)
80, last_samp(p_FiffRawData.last_samp)
81, cals(p_FiffRawData.cals)
82, rawdir(p_FiffRawData.rawdir)
83, proj(p_FiffRawData.proj)
84, comp(p_FiffRawData.comp)
85{
86}
87
88//=============================================================================================================
89
93
94//=============================================================================================================
95
97{
98 info.clear();
99 first_samp = -1;
100 last_samp = -1;
101 cals = RowVectorXd();
102 rawdir.clear();
103 proj = MatrixXd();
104 comp.clear();
105}
106
107//=============================================================================================================
108#include <QElapsedTimer>
109#include <QDebug>
110bool FiffRawData::read_raw_segment(MatrixXd& data,
111 MatrixXd& times,
112 fiff_int_t from,
113 fiff_int_t to,
114 const RowVectorXi& sel,
115 bool do_debug) const
116{
117 bool projAvailable = true;
118
119 if (this->proj.size() == 0) {
120 //qDebug() << "FiffRawData::read_raw_segment - No projectors setup. Consider calling MNE::setup_compensators.";
121 projAvailable = false;
122 }
123
124 if (from == -1)
125 from = this->first_samp;
126 if (to == -1)
127 to = this->last_samp;
128 //
129 // Initial checks
130 //
131 if (from < this->first_samp)
132 from = this->first_samp;
133 if (to > this->last_samp)
134 to = this->last_samp;
135 //
136 if (from > to) {
137 qWarning("No data in this range %d ... %d = %9.3f ... %9.3f secs...", from, to, (static_cast<float>(from)) / this->info.sfreq, (static_cast<float>(to)) / this->info.sfreq);
138 return false;
139 }
140 //
141 // Initialize the data and calibration vector
142 //
143 qint32 nchan = this->info.nchan;
144 qint32 dest = 0; //1;
145 qint32 i, k, r;
146
147 using T = Eigen::Triplet<double>;
148 std::vector<T> tripletList;
149 tripletList.reserve(nchan);
150 for (i = 0; i < nchan; ++i)
151 tripletList.push_back(T(i, i, this->cals[i]));
152
153 SparseMatrix<double> cal(nchan, nchan);
154 cal.setFromTriplets(tripletList.begin(), tripletList.end());
155 // cal.makeCompressed();
156
157 MatrixXd mult_full;
158 //
159 if (sel.size() == 0) {
160 data = MatrixXd(nchan, to - from + 1);
161 // data->setZero();
162 if (projAvailable || this->comp.kind != -1) {
163 if (!projAvailable)
164 mult_full = this->comp.data->data * cal;
165 else if (this->comp.kind == -1)
166 mult_full = this->proj * cal;
167 else
168 mult_full = this->proj * this->comp.data->data * cal;
169 }
170 } else {
171 data = MatrixXd(sel.size(), to - from + 1);
172 // data->setZero();
173
174 MatrixXd selVect(sel.size(), nchan);
175
176 selVect.setZero();
177
178 if (!projAvailable && this->comp.kind == -1) {
179 tripletList.clear();
180 tripletList.reserve(sel.size());
181 for (i = 0; i < sel.size(); ++i)
182 tripletList.push_back(T(i, i, this->cals[sel[i]]));
183 cal = SparseMatrix<double>(sel.size(), sel.size());
184 cal.setFromTriplets(tripletList.begin(), tripletList.end());
185 } else {
186 if (!projAvailable) {
187 for (i = 0; i < sel.size(); ++i)
188 selVect.row(i) = this->comp.data->data.block(sel[i], 0, 1, nchan);
189 mult_full = selVect * cal;
190 } else if (this->comp.kind == -1) {
191 for (i = 0; i < sel.size(); ++i)
192 selVect.row(i) = this->proj.block(sel[i], 0, 1, nchan);
193
194 mult_full = selVect * cal;
195 } else {
196 for (i = 0; i < sel.size(); ++i)
197 selVect.row(i) = this->proj.block(sel[i], 0, 1, nchan);
198
199 mult_full = selVect * this->comp.data->data * cal;
200 }
201 }
202 }
203
204 //
205 // Make mult sparse
206 //
207 tripletList.clear();
208 tripletList.reserve(mult_full.rows() * mult_full.cols());
209 for (i = 0; i < mult_full.rows(); ++i)
210 for (k = 0; k < mult_full.cols(); ++k)
211 if (mult_full(i, k) != 0)
212 tripletList.push_back(T(i, k, mult_full(i, k)));
213
214 SparseMatrix<double> mult(mult_full.rows(), mult_full.cols());
215 if (tripletList.size() > 0)
216 mult.setFromTriplets(tripletList.begin(), tripletList.end());
217 // mult.makeCompressed();
218
220 if (!this->file->device()->isOpen()) {
221 if (!this->file->device()->open(QIODevice::ReadOnly)) {
222 qWarning("Cannot open file %s", this->info.filename.toUtf8().constData());
223 }
224 fid = this->file;
225 } else {
226 fid = this->file;
227 }
228
229 MatrixXd one, newData, tmp_data;
230 FiffRawDir thisRawDir;
231 FiffTag::UPtr t_pTag;
232 fiff_int_t first_pick, last_pick, picksamp;
233 for (k = 0; k < this->rawdir.size(); ++k) {
234 thisRawDir = this->rawdir[k];
235 //
236 // Do we need this buffer
237 //
238 if (thisRawDir.last > from) {
239 if (thisRawDir.ent->kind == -1) {
240 //
241 // Take the easy route: skip is translated to zeros
242 //
243 if (do_debug)
244 qDebug("S");
245 if (sel.cols() <= 0)
246 one.resize(nchan, thisRawDir.nsamp);
247 else
248 one.resize(sel.cols(), thisRawDir.nsamp);
249
250 one.setZero();
251 } else {
252 fid->read_tag(t_pTag, thisRawDir.ent->pos);
253 //
254 // Depending on the state of the projection and selection
255 // we proceed a little bit differently
256 //
257 if (mult.cols() == 0) {
258 if (sel.cols() == 0) {
259 if (t_pTag->type == FIFFT_DAU_PACK16)
260 one = cal * (Map<const MatrixDau16>(t_pTag->toDauPack16(), nchan, thisRawDir.nsamp)).cast<double>();
261 else if (t_pTag->type == FIFFT_INT)
262 one = cal * (Map<const MatrixXi>(t_pTag->toInt(), nchan, thisRawDir.nsamp)).cast<double>();
263 else if (t_pTag->type == FIFFT_FLOAT)
264 one = cal * (Map<const MatrixXf>(t_pTag->toFloat(), nchan, thisRawDir.nsamp)).cast<double>();
265 else if (t_pTag->type == FIFFT_SHORT)
266 one = cal * (Map<const MatrixShort>(t_pTag->toShort(), nchan, thisRawDir.nsamp)).cast<double>();
267 else if (t_pTag->type == FIFFT_DOUBLE)
268 one = cal * Map<const MatrixXd>(t_pTag->toDouble(), nchan, thisRawDir.nsamp);
269 else {
270 qWarning("Data Storage Format not known yet [1]!! Type: %d\n", t_pTag->type);
271 this->file->device()->close();
272 return false;
273 }
274 } else {
275 //ToDo find a faster solution for this!! --> make cal and mul sparse like in MATLAB
276 newData.resize(sel.cols(), thisRawDir.nsamp); //ToDo this can be done much faster, without newData
277
278 if (t_pTag->type == FIFFT_DAU_PACK16) {
279 tmp_data = (Map<const MatrixDau16>(t_pTag->toDauPack16(), nchan, thisRawDir.nsamp)).cast<double>();
280
281 for (r = 0; r < sel.size(); ++r)
282 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
283 } else if (t_pTag->type == FIFFT_INT) {
284 tmp_data = (Map<const MatrixXi>(t_pTag->toInt(), nchan, thisRawDir.nsamp)).cast<double>();
285
286 for (r = 0; r < sel.size(); ++r)
287 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
288 } else if (t_pTag->type == FIFFT_FLOAT) {
289 tmp_data = (Map<const MatrixXf>(t_pTag->toFloat(), nchan, thisRawDir.nsamp)).cast<double>();
290
291 for (r = 0; r < sel.size(); ++r)
292 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
293 } else if (t_pTag->type == FIFFT_SHORT) {
294 tmp_data = (Map<const MatrixShort>(t_pTag->toShort(), nchan, thisRawDir.nsamp)).cast<double>();
295
296 for (r = 0; r < sel.size(); ++r)
297 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
298 } else if (t_pTag->type == FIFFT_DOUBLE) {
299 tmp_data = Map<const MatrixXd>(t_pTag->toDouble(), nchan, thisRawDir.nsamp);
300
301 for (r = 0; r < sel.size(); ++r)
302 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
303 } else {
304 qWarning("Data Storage Format not known yet [2]!! Type: %d\n", t_pTag->type);
305 this->file->device()->close();
306 return false;
307 }
308
309 one = cal * newData;
310 }
311 } else {
312 if (t_pTag->type == FIFFT_DAU_PACK16)
313 one = mult * (Map<const MatrixDau16>(t_pTag->toDauPack16(), nchan, thisRawDir.nsamp)).cast<double>();
314 else if (t_pTag->type == FIFFT_INT)
315 one = mult * (Map<const MatrixXi>(t_pTag->toInt(), nchan, thisRawDir.nsamp)).cast<double>();
316 else if (t_pTag->type == FIFFT_FLOAT)
317 one = mult * (Map<const MatrixXf>(t_pTag->toFloat(), nchan, thisRawDir.nsamp)).cast<double>();
318 else if (t_pTag->type == FIFFT_SHORT)
319 one = mult * (Map<const MatrixShort>(t_pTag->toShort(), nchan, thisRawDir.nsamp)).cast<double>();
320 else if (t_pTag->type == FIFFT_DOUBLE)
321 one = mult * Map<const MatrixXd>(t_pTag->toDouble(), nchan, thisRawDir.nsamp);
322 else {
323 qWarning("Data Storage Format not known yet [3]!! Type: %d\n", t_pTag->type);
324 this->file->device()->close();
325 return false;
326 }
327 }
328 }
329 //
330 // The picking logic is a bit complicated
331 //
332 if (to >= thisRawDir.last && from <= thisRawDir.first) {
333 //
334 // We need the whole buffer
335 //
336 first_pick = 0; //1;
337 last_pick = thisRawDir.nsamp - 1;
338 if (do_debug)
339 qDebug("W");
340 } else if (from > thisRawDir.first) {
341 first_pick = from - thisRawDir.first; // + 1;
342 if (to < thisRawDir.last) {
343 //
344 // Something from the middle
345 //
346 // qDebug() << "This needs to be debugged!";
347 last_pick = thisRawDir.nsamp + to - thisRawDir.last - 1; //is this alright?
348 if (do_debug)
349 qDebug("M");
350 } else {
351 //
352 // From the middle to the end
353 //
354 last_pick = thisRawDir.nsamp - 1;
355 if (do_debug)
356 qDebug("E");
357 }
358 } else {
359 //
360 // From the beginning to the middle
361 //
362 first_pick = 0; //1;
363 last_pick = to - thisRawDir.first; // + 1;
364 if (do_debug)
365 qDebug("B");
366 }
367 //
368 // Now we are ready to pick
369 //
370 picksamp = last_pick - first_pick + 1;
371
372 if (do_debug) {
373 qDebug() << "first_pick: " << first_pick;
374 qDebug() << "last_pick: " << last_pick;
375 qDebug() << "picksamp: " << picksamp;
376 }
377
378 if (picksamp > 0) {
379 // for(r = 0; r < data->rows(); ++r)
380 // for(c = 0; c < picksamp; ++c)
381 // (*data)(r,dest + c) = one(r,first_pick + c);
382 data.block(0, dest, data.rows(), picksamp) = one.block(0, first_pick, data.rows(), picksamp);
383
384 dest += picksamp;
385 }
386 }
387 //
388 // Done?
389 //
390 if (thisRawDir.last >= to) {
391 //qInfo(" [done]");
392 break;
393 }
394 }
395
396 if (!this->file->device()->isOpen()) {
397 this->file->device()->close();
398 }
399
400 times = MatrixXd(1, to - from + 1);
401
402 for (i = 0; i < times.cols(); ++i)
403 times(0, i) = static_cast<float>(from + i) / this->info.sfreq;
404
405 return true;
406}
407
408//=============================================================================================================
409
410bool FiffRawData::read_raw_segment(MatrixXd& data,
411 MatrixXd& times,
412 SparseMatrix<double>& multSegment,
413 fiff_int_t from,
414 fiff_int_t to,
415 const RowVectorXi& sel,
416 bool do_debug) const
417{
418 bool projAvailable = true;
419
420 if (this->proj.size() == 0) {
421 //qInfo() << "FiffRawData::read_raw_segment - No projectors setup. Consider calling MNE::setup_compensators.";
422 projAvailable = false;
423 }
424
425 if (from == -1)
426 from = this->first_samp;
427 if (to == -1)
428 to = this->last_samp;
429 //
430 // Initial checks
431 //
432 if (from < this->first_samp)
433 from = this->first_samp;
434 if (to > this->last_samp)
435 to = this->last_samp;
436 //
437 if (from > to) {
438 qWarning("No data in this range");
439 return false;
440 }
441 //
442 // Initialize the data and calibration vector
443 //
444 qint32 nchan = this->info.nchan;
445 qint32 dest = 0; //1;
446 qint32 i, k, r;
447
448 using T = Eigen::Triplet<double>;
449 std::vector<T> tripletList;
450 tripletList.reserve(nchan);
451 for (i = 0; i < nchan; ++i)
452 tripletList.push_back(T(i, i, this->cals[i]));
453
454 SparseMatrix<double> cal(nchan, nchan);
455 cal.setFromTriplets(tripletList.begin(), tripletList.end());
456 // cal.makeCompressed();
457
458 MatrixXd mult_full;
459 //
460 if (sel.size() == 0) {
461 data = MatrixXd(nchan, to - from + 1);
462 // data->setZero();
463 if (projAvailable || this->comp.kind != -1) {
464 if (!projAvailable)
465 mult_full = this->comp.data->data * cal;
466 else if (this->comp.kind == -1)
467 mult_full = this->proj * cal;
468 else
469 mult_full = this->proj * this->comp.data->data * cal;
470 }
471 } else {
472 data = MatrixXd(sel.size(), to - from + 1);
473 // data->setZero();
474
475 MatrixXd selVect(sel.size(), nchan);
476
477 selVect.setZero();
478
479 if (!projAvailable && this->comp.kind == -1) {
480 tripletList.clear();
481 tripletList.reserve(sel.size());
482 for (i = 0; i < sel.size(); ++i)
483 tripletList.push_back(T(i, i, this->cals[sel[i]]));
484 cal = SparseMatrix<double>(sel.size(), sel.size());
485 cal.setFromTriplets(tripletList.begin(), tripletList.end());
486 } else {
487 if (!projAvailable) {
488 for (i = 0; i < sel.size(); ++i)
489 selVect.row(i) = this->comp.data->data.block(sel[i], 0, 1, nchan);
490 mult_full = selVect * cal;
491 } else if (this->comp.kind == -1) {
492 for (i = 0; i < sel.size(); ++i)
493 selVect.row(i) = this->proj.block(sel[i], 0, 1, nchan);
494
495 mult_full = selVect * cal;
496 } else {
497 for (i = 0; i < sel.size(); ++i)
498 selVect.row(i) = this->proj.block(sel[i], 0, 1, nchan);
499
500 mult_full = selVect * this->comp.data->data * cal;
501 }
502 }
503 }
504
505 //
506 // Make mult sparse
507 //
508 tripletList.clear();
509 tripletList.reserve(mult_full.rows() * mult_full.cols());
510 for (i = 0; i < mult_full.rows(); ++i)
511 for (k = 0; k < mult_full.cols(); ++k)
512 if (mult_full(i, k) != 0)
513 tripletList.push_back(T(i, k, mult_full(i, k)));
514
515 SparseMatrix<double> mult(mult_full.rows(), mult_full.cols());
516 if (tripletList.size() > 0)
517 mult.setFromTriplets(tripletList.begin(), tripletList.end());
518 // mult.makeCompressed();
519
520 //
521
523 if (!this->file->device()->isOpen()) {
524 if (!this->file->device()->open(QIODevice::ReadOnly)) {
525 qWarning("Cannot open file %s", this->info.filename.toUtf8().constData());
526 }
527 fid = this->file;
528 } else {
529 fid = this->file;
530 }
531
532 MatrixXd one;
533 fiff_int_t first_pick, last_pick, picksamp;
534 for (k = 0; k < this->rawdir.size(); ++k) {
535 FiffRawDir thisRawDir = this->rawdir[k];
536 //
537 // Do we need this buffer
538 //
539 if (thisRawDir.last > from) {
540 if (thisRawDir.ent->kind == -1) {
541 //
542 // Take the easy route: skip is translated to zeros
543 //
544 if (do_debug)
545 qDebug("S");
546 if (sel.cols() <= 0)
547 one.resize(nchan, thisRawDir.nsamp);
548 else
549 one.resize(sel.cols(), thisRawDir.nsamp);
550
551 one.setZero();
552 } else {
553 FiffTag::UPtr t_pTag;
554 fid->read_tag(t_pTag, thisRawDir.ent->pos);
555 //
556 // Depending on the state of the projection and selection
557 // we proceed a little bit differently
558 //
559 if (mult.cols() == 0) {
560 if (sel.cols() == 0) {
561 if (t_pTag->type == FIFFT_DAU_PACK16)
562 one = cal * (Map<const MatrixDau16>(t_pTag->toDauPack16(), nchan, thisRawDir.nsamp)).cast<double>();
563 else if (t_pTag->type == FIFFT_INT)
564 one = cal * (Map<const MatrixXi>(t_pTag->toInt(), nchan, thisRawDir.nsamp)).cast<double>();
565 else if (t_pTag->type == FIFFT_FLOAT)
566 one = cal * (Map<const MatrixXf>(t_pTag->toFloat(), nchan, thisRawDir.nsamp)).cast<double>();
567 else if (t_pTag->type == FIFFT_SHORT)
568 one = cal * (Map<const MatrixShort>(t_pTag->toShort(), nchan, thisRawDir.nsamp)).cast<double>();
569 else if (t_pTag->type == FIFFT_DOUBLE)
570 one = cal * Map<const MatrixXd>(t_pTag->toDouble(), nchan, thisRawDir.nsamp);
571 else {
572 qWarning("Data Storage Format not known yet [1]!! Type: %d\n", t_pTag->type);
573 this->file->device()->close();
574 return false;
575 }
576 } else {
577 //ToDo find a faster solution for this!! --> make cal and mul sparse like in MATLAB
578 MatrixXd newData(sel.cols(), thisRawDir.nsamp); //ToDo this can be done much faster, without newData
579
580 if (t_pTag->type == FIFFT_DAU_PACK16) {
581 MatrixXd tmp_data = (Map<const MatrixDau16>(t_pTag->toDauPack16(), nchan, thisRawDir.nsamp)).cast<double>();
582
583 for (r = 0; r < sel.size(); ++r)
584 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
585 } else if (t_pTag->type == FIFFT_INT) {
586 MatrixXd tmp_data = (Map<const MatrixXi>(t_pTag->toInt(), nchan, thisRawDir.nsamp)).cast<double>();
587
588 for (r = 0; r < sel.size(); ++r)
589 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
590 } else if (t_pTag->type == FIFFT_FLOAT) {
591 MatrixXd tmp_data = (Map<const MatrixXf>(t_pTag->toFloat(), nchan, thisRawDir.nsamp)).cast<double>();
592
593 for (r = 0; r < sel.size(); ++r)
594 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
595 } else if (t_pTag->type == FIFFT_SHORT) {
596 MatrixXd tmp_data = (Map<const MatrixShort>(t_pTag->toShort(), nchan, thisRawDir.nsamp)).cast<double>();
597
598 for (r = 0; r < sel.size(); ++r)
599 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
600 } else if (t_pTag->type == FIFFT_DOUBLE) {
601 MatrixXd tmp_data = Map<const MatrixXd>(t_pTag->toDouble(), nchan, thisRawDir.nsamp);
602
603 for (r = 0; r < sel.size(); ++r)
604 newData.block(r, 0, 1, thisRawDir.nsamp) = tmp_data.block(sel[r], 0, 1, thisRawDir.nsamp);
605 } else {
606 qWarning("Data Storage Format not known yet [2]!! Type: %d\n", t_pTag->type);
607 this->file->device()->close();
608 return false;
609 }
610
611 one = cal * newData;
612 }
613 } else {
614 if (t_pTag->type == FIFFT_DAU_PACK16)
615 one = mult * (Map<const MatrixDau16>(t_pTag->toDauPack16(), nchan, thisRawDir.nsamp)).cast<double>();
616 else if (t_pTag->type == FIFFT_INT)
617 one = mult * (Map<const MatrixXi>(t_pTag->toInt(), nchan, thisRawDir.nsamp)).cast<double>();
618 else if (t_pTag->type == FIFFT_FLOAT)
619 one = mult * (Map<const MatrixXf>(t_pTag->toFloat(), nchan, thisRawDir.nsamp)).cast<double>();
620 else if (t_pTag->type == FIFFT_SHORT)
621 one = mult * (Map<const MatrixShort>(t_pTag->toShort(), nchan, thisRawDir.nsamp)).cast<double>();
622 else if (t_pTag->type == FIFFT_DOUBLE)
623 one = mult * Map<const MatrixXd>(t_pTag->toDouble(), nchan, thisRawDir.nsamp);
624 else {
625 qWarning("Data Storage Format not known yet [3]!! Type: %d\n", t_pTag->type);
626 this->file->device()->close();
627 return false;
628 }
629 }
630 }
631 //
632 // The picking logic is a bit complicated
633 //
634 if (to >= thisRawDir.last && from <= thisRawDir.first) {
635 //
636 // We need the whole buffer
637 //
638 first_pick = 0; //1;
639 last_pick = thisRawDir.nsamp - 1;
640 if (do_debug)
641 qDebug("W");
642 } else if (from > thisRawDir.first) {
643 first_pick = from - thisRawDir.first; // + 1;
644 if (to < thisRawDir.last) {
645 //
646 // Something from the middle
647 //
648 // qDebug() << "This needs to be debugged!";
649 last_pick = thisRawDir.nsamp + to - thisRawDir.last - 1; //is this alright?
650 if (do_debug)
651 qDebug("M");
652 } else {
653 //
654 // From the middle to the end
655 //
656 last_pick = thisRawDir.nsamp - 1;
657 if (do_debug)
658 qDebug("E");
659 }
660 } else {
661 //
662 // From the beginning to the middle
663 //
664 first_pick = 0; //1;
665 last_pick = to - thisRawDir.first; // + 1;
666 if (do_debug)
667 qDebug("B");
668 }
669 //
670 // Now we are ready to pick
671 //
672 picksamp = last_pick - first_pick + 1;
673
674 if (do_debug) {
675 qDebug() << "first_pick: " << first_pick;
676 qDebug() << "last_pick: " << last_pick;
677 qDebug() << "picksamp: " << picksamp;
678 }
679
680 if (picksamp > 0) {
681 // for(r = 0; r < data->rows(); ++r)
682 // for(c = 0; c < picksamp; ++c)
683 // (*data)(r,dest + c) = one(r,first_pick + c);
684 data.block(0, dest, data.rows(), picksamp) = one.block(0, first_pick, data.rows(), picksamp);
685
686 dest += picksamp;
687 }
688 }
689 //
690 // Done?
691 //
692 if (thisRawDir.last >= to) {
693 //qInfo(" [done]");
694 break;
695 }
696 }
697
698 if (mult.cols() == 0)
699 multSegment = cal;
700 else
701 multSegment = mult;
702
703 if (!this->file->device()->isOpen()) {
704 this->file->device()->close();
705 }
706
707 times = MatrixXd(1, to - from + 1);
708
709 for (i = 0; i < times.cols(); ++i)
710 times(0, i) = static_cast<float>(from + i) / this->info.sfreq;
711
712 return true;
713}
714
715//=============================================================================================================
716
718 MatrixXd& times,
719 float from,
720 float to,
721 const RowVectorXi& sel) const
722{
723 //
724 // Convert to samples
725 //
726 from = floor(static_cast<double>(from) * this->info.sfreq);
727 to = ceil(static_cast<double>(to) * this->info.sfreq);
728 //
729 // Read it
730 //
731 return this->read_raw_segment(data, times, static_cast<qint32>(from), static_cast<qint32>(to), sel);
732}
733
734//=============================================================================================================
735
736bool FiffRawData::save(QIODevice& p_IODevice,
737 const RowVectorXi& picks,
738 int decim,
739 int from,
740 int to) const
741{
742 if (decim < 1)
743 decim = 1;
744
745 int firstSamp = (from >= 0) ? from : first_samp;
746 int lastSamp = (to >= 0) ? to : last_samp;
747
748 if (firstSamp > lastSamp) {
749 qWarning() << "[FiffRawData::save] Invalid sample range.";
750 return false;
751 }
752
753 // start_writing_raw picks the channels itself; picking info here as well applied picks twice.
754 FiffInfo outInfo = info;
755
756 // Adjust sampling frequency for decimation
757 if (decim > 1) {
758 outInfo.sfreq = info.sfreq / static_cast<float>(decim);
759 }
760
761 // Use the standard start_writing_raw pipeline
762 RowVectorXd calsOut;
763 FiffStream::SPtr pStream = FiffStream::start_writing_raw(p_IODevice, outInfo, calsOut, picks);
764 if (!pStream) {
765 qWarning() << "[FiffRawData::save] Cannot start writing raw file.";
766 return false;
767 }
768
769 // Without this the copy restarts at sample 0 instead of at its position in the recording.
770 int firstOut = firstSamp / decim;
771 pStream->write_int(FIFF_FIRST_SAMPLE, &firstOut);
772
773 // Write data in blocks
774 const int blockSize = 2000;
775 int blockSamples = decim * blockSize;
776
777 for (int samp = firstSamp; samp <= lastSamp; samp += blockSamples) {
778 int nsamp = qMin(blockSamples, lastSamp - samp + 1);
779
780 MatrixXd segData;
781 MatrixXd segTimes;
782 if (!read_raw_segment(segData, segTimes, samp, samp + nsamp - 1, picks)) {
783 qWarning() << "[FiffRawData::save] Error reading data at sample" << samp;
784 pStream->finish_writing_raw();
785 return false;
786 }
787
788 // Decimate if needed
789 if (decim > 1) {
790 int nOut = (nsamp + decim - 1) / decim;
791 MatrixXd decimData(segData.rows(), nOut);
792 for (int s = 0, idx = 0; s < nsamp && idx < nOut; s += decim, ++idx) {
793 decimData.col(idx) = segData.col(s);
794 }
795 segData = decimData;
796 }
797
798 pStream->write_raw_buffer(segData, calsOut);
799 }
800
801 pStream->finish_writing_raw();
802
803 qInfo() << "[FiffRawData::save] Saved raw data from sample" << firstSamp
804 << "to" << lastSamp << "(decim=" << decim << ")";
805 return true;
806}
FIFF tag: the 16-byte tag header (kind, type, size, next) plus its decoded payload.
FIFF binary tag-stream layer: wraps a QIODevice to read and write FIFF tags, directories,...
constexpr auto FIFFT_DOUBLE
Definition fiff_file.h:224
constexpr auto FIFFT_DAU_PACK16
Definition fiff_file.h:234
constexpr auto FIFF_FIRST_SAMPLE
Definition fiff_file.h:452
constexpr auto FIFFT_SHORT
Definition fiff_file.h:221
constexpr auto FIFFT_FLOAT
Definition fiff_file.h:223
constexpr auto FIFFT_INT
Definition fiff_file.h:222
Stim-channel event list (sample, previous value, new value triples) with FIFF read/write helpers.
FIFF continuous raw recording: FiffInfo plus a directory of FIFF_DATA_BUFFER tags for random-access s...
FIFF file I/O, in-memory data structures and high-level readers/writers.
qint32 fiff_int_t
Definition fiff_types.h:86
Full FIFF measurement info: per-channel descriptors, sampling and filter setup, projectors,...
Definition fiff_info.h:90
Eigen::RowVectorXd cals
bool read_raw_segment(Eigen::MatrixXd &data, Eigen::MatrixXd &times, fiff_int_t from=-1, fiff_int_t to=-1, const Eigen::RowVectorXi &sel=defaultRowVectorXi, bool do_debug=false) const
FiffStream::SPtr file
bool save(QIODevice &p_IODevice, const Eigen::RowVectorXi &picks=Eigen::RowVectorXi(), int decim=1, int from=-1, int to=-1) const
Eigen::MatrixXd proj
bool read_raw_segment_times(Eigen::MatrixXd &data, Eigen::MatrixXd &times, float from, float to, const Eigen::RowVectorXi &sel=defaultRowVectorXi) const
QList< FiffRawDir > rawdir
FiffDirEntry::SPtr ent
QSharedPointer< FiffStream > SPtr
static bool setup_read_raw(QIODevice &p_IODevice, FiffRawData &data, bool allow_maxshield=true, bool is_littleEndian=false)
static FiffStream::SPtr start_writing_raw(QIODevice &p_IODevice, const FiffInfo &info, Eigen::RowVectorXd &cals, Eigen::MatrixXi sel=defaultMatrixXi, bool bResetRange=true)
std::unique_ptr< FiffTag > UPtr
Definition fiff_tag.h:165