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