MNE-CPP  0.1.9
A Framework for Electrophysiology
filter.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //=============================================================================================================
37 // INCLUDES
38 //=============================================================================================================
39 
40 #include "filter.h"
41 
42 #include <utils/mnemath.h>
43 #include <fiff/fiff_raw_data.h>
44 
45 //=============================================================================================================
46 // QT INCLUDES
47 //=============================================================================================================
48 
49 #include <QDebug>
50 
51 //=============================================================================================================
52 // EIGEN INCLUDES
53 //=============================================================================================================
54 
55 #include <Eigen/Dense>
56 #include <Eigen/Core>
57 
58 //=============================================================================================================
59 // USED NAMESPACES
60 //=============================================================================================================
61 
62 using namespace RTPROCESSINGLIB;
63 using namespace Eigen;
64 using namespace FIFFLIB;
65 using namespace UTILSLIB;
66 
67 //=============================================================================================================
68 // DEFINE GLOBAL RTPROCESSINGLIB METHODS
69 //=============================================================================================================
70 
71 bool RTPROCESSINGLIB::filterFile(QIODevice &pIODevice,
72  QSharedPointer<FiffRawData> pFiffRawData,
73  int type,
74  double dCenterfreq,
75  double bandwidth,
76  double dTransition,
77  double dSFreq,
78  int iOrder,
79  int designMethod,
80  const RowVectorXi& vecPicks,
81  bool bUseThreads)
82 {
83  // Normalize cut off frequencies to nyquist
84  dCenterfreq = dCenterfreq/(dSFreq/2.0);
85  bandwidth = bandwidth/(dSFreq/2.0);
86  dTransition = dTransition/(dSFreq/2.0);
87 
88  // create filter
89  FilterKernel filter = FilterKernel("filter_kernel",
90  type,
91  iOrder,
92  dCenterfreq,
93  bandwidth,
94  dTransition,
95  dSFreq,
96  designMethod);
97 
98  return filterFile(pIODevice,
99  pFiffRawData,
100  filter,
101  vecPicks,
102  bUseThreads);
103 }
104 
105 //=============================================================================================================
106 
107 bool RTPROCESSINGLIB::filterFile(QIODevice &pIODevice,
108  QSharedPointer<FiffRawData> pFiffRawData,
109  const FilterKernel& filterKernel,
110  const RowVectorXi& vecPicks,
111  bool bUseThreads)
112 {
113  int iOrder = filterKernel.getFilterOrder();
114 
115  RowVectorXd cals;
116  SparseMatrix<double> mult;
117  RowVectorXi sel;
118  FiffStream::SPtr outfid = FiffStream::start_writing_raw(pIODevice, pFiffRawData->info, cals);
119 
120  //Setup reading parameters
121  fiff_int_t from = pFiffRawData->first_samp;
122  fiff_int_t to = pFiffRawData->last_samp;
123 
124  // slice input data into data junks with proper length so that the slices are always >= the filter order
125  float fFactor = 2.0f;
126  int iSize = fFactor * iOrder;
127  int residual = (to - from) % iSize;
128  while(residual < iOrder) {
129  fFactor = fFactor - 0.1f;
130  iSize = fFactor * iOrder;
131  residual = (to - from) % iSize;
132 
133  if((iSize < iOrder)) {
134  qInfo() << "[Filter::filterData] Sliced data block size is too small. Filtering whole block at once.";
135  iSize = to - from;
136  break;
137  }
138  }
139 
140  float quantum_sec = iSize/pFiffRawData->info.sfreq;
141  fiff_int_t quantum = ceil(quantum_sec*pFiffRawData->info.sfreq);
142 
143  // Read, filter and write the data
144  bool first_buffer = true;
145 
146  fiff_int_t first, last;
147  MatrixXd matData, matDataOverlap;
148  MatrixXd times;
149 
150  for(first = from; first < to; first+=quantum) {
151  last = first+quantum-1;
152  if (last > to) {
153  last = to;
154  }
155 
156  if (!pFiffRawData->read_raw_segment(matData, times, mult, first, last, sel)) {
157  qWarning("[Filter::filterData] Error during read_raw_segment\n");
158  return false;
159  }
160 
161  qInfo() << "Filtering and writing block" << first << "to" << last;
162 
163  if (first_buffer) {
164  if (first > 0) {
165  outfid->write_int(FIFF_FIRST_SAMPLE,&first);
166  }
167  first_buffer = false;
168  }
169 
170  matData = filterDataBlock(matData,
171  vecPicks,
172  filterKernel,
173  bUseThreads);
174 
175  if(first == from) {
176  outfid->write_raw_buffer(matData.block(0,iOrder/2,matData.rows(),matData.cols()-iOrder), cals);
177  } else if(first + quantum >= to) {
178  matData.block(0,0,matData.rows(),iOrder) += matDataOverlap;
179  outfid->write_raw_buffer(matData.block(0,0,matData.rows(),matData.cols()-iOrder), cals);
180  } else {
181  matData.block(0,0,matData.rows(),iOrder) += matDataOverlap;
182  outfid->write_raw_buffer(matData.block(0,0,matData.rows(),matData.cols()-iOrder), cals);
183  }
184 
185  matDataOverlap = matData.block(0,matData.cols()-iOrder,matData.rows(),iOrder);
186  }
187 
188  outfid->finish_writing_raw();
189 
190  return true;
191 }
192 
193 //=============================================================================================================
194 
195 MatrixXd RTPROCESSINGLIB::filterData(const MatrixXd& mataData,
196  int type,
197  double dCenterfreq,
198  double bandwidth,
199  double dTransition,
200  double dSFreq,
201  int iOrder,
202  int designMethod,
203  const RowVectorXi& vecPicks,
204  bool bUseThreads,
205  bool bKeepOverhead)
206 {
207  // Check for size of data
208  if(mataData.cols() < iOrder){
209  qWarning() << QString("[Filter::filterData] Filter length/order is bigger than data length. Returning.");
210  return mataData;
211  }
212 
213  // Normalize cut off frequencies to nyquist
214  dCenterfreq = dCenterfreq/(dSFreq/2.0);
215  bandwidth = bandwidth/(dSFreq/2.0);
216  dTransition = dTransition/(dSFreq/2.0);
217 
218  // create filter
219  FilterKernel filter = FilterKernel("filter_kernel",
220  type,
221  iOrder,
222  dCenterfreq,
223  bandwidth,
224  dTransition,
225  dSFreq,
226  designMethod);
227 
228  return filterData(mataData,
229  filter,
230  vecPicks,
231  bUseThreads,
232  bKeepOverhead);
233 }
234 
235 //=============================================================================================================
236 
237 MatrixXd RTPROCESSINGLIB::filterData(const MatrixXd& mataData,
238  const FilterKernel& filterKernel,
239  const RowVectorXi& vecPicks,
240  bool bUseThreads,
241  bool bKeepOverhead)
242 {
243  int iOrder = filterKernel.getFilterOrder();
244 
245  // Check for size of data
246  if(mataData.cols() < iOrder){
247  qWarning() << "[Filter::filterData] Filter length/order is bigger than data length. Returning.";
248  return mataData;
249  }
250 
251  // Create output matrix with size of input matrix
252  MatrixXd matDataOut(mataData.rows(), mataData.cols()+iOrder);
253  matDataOut.setZero();
254  MatrixXd sliceFiltered;
255 
256  // slice input data into data junks with proper length so that the slices are always >= the filter order
257  float fFactor = 2.0f;
258  int iSize = fFactor * iOrder;
259  int residual = mataData.cols() % iSize;
260  while(residual < iOrder) {
261  fFactor = fFactor - 0.1f;
262  iSize = fFactor * iOrder;
263  residual = mataData.cols() % iSize;
264 
265  if(iSize < iOrder) {
266  iSize = mataData.cols();
267  break;
268  }
269  }
270 
271  if(mataData.cols() > iSize) {
272  int from = 0;
273  int numSlices = ceil(float(mataData.cols())/float(iSize)); //calculate number of data slices
274 
275  for (int i = 0; i < numSlices; i++) {
276  if(i == numSlices-1) {
277  //catch the last one that might be shorter than the other blocks
278  iSize = mataData.cols() - (iSize * (numSlices -1));
279  }
280 
281  // Filter the data block. This will return data with a fitler delay of iOrder/2 in front and back
282  sliceFiltered = filterDataBlock(mataData.block(0,from,mataData.rows(),iSize),
283  vecPicks,
284  filterKernel,
285  bUseThreads);
286 
287  // Perform overlap add
288  if(i == 0) {
289  matDataOut.block(0,0,mataData.rows(),sliceFiltered.cols()) += sliceFiltered;
290  } else {
291  matDataOut.block(0,from,mataData.rows(),sliceFiltered.cols()) += sliceFiltered;
292  }
293 
294  from += iSize;
295  }
296  } else {
297  matDataOut = filterDataBlock(mataData,
298  vecPicks,
299  filterKernel,
300  bUseThreads);
301  }
302 
303  if(bKeepOverhead) {
304  return matDataOut;
305  } else {
306  return matDataOut.block(0,iOrder/2,matDataOut.rows(),mataData.cols());
307  }
308 }
309 
310 //=============================================================================================================
311 
312 MatrixXd RTPROCESSINGLIB::filterDataBlock(const MatrixXd& mataData,
313  const RowVectorXi& vecPicks,
314  const FilterKernel& filterKernel,
315  bool bUseThreads)
316 {
317  int iOrder = filterKernel.getFilterOrder();
318 
319  // Check for size of data
320  if(mataData.cols() < iOrder){
321  qWarning() << QString("[Filter::filterDataBlock] Filter length/order is bigger than data length. Returning.");
322  return mataData;
323  }
324 
325  // Setup filters to the correct length, so we do not have to do this everytime we call the FFT filter function
326  FilterKernel filterKernelSetup = filterKernel;
327  filterKernelSetup.prepareFilter(mataData.cols());
328 
329  // Do the concurrent filtering
330  RowVectorXi vecPicksNew = vecPicks;
331  if(vecPicksNew.cols() == 0) {
332  vecPicksNew = RowVectorXi::LinSpaced(mataData.rows(), 0, mataData.rows());
333  }
334 
335  // Generate QList structure which can be handled by the QConcurrent framework
336  QList<FilterObject> timeData;
337 
338  // Only select channels specified in vecPicksNew
339  FilterObject data;
340  for(qint32 i = 0; i < vecPicksNew.cols(); ++i) {
341  data.filterKernel = filterKernelSetup;
342  data.iRow = vecPicksNew[i];
343  data.vecData = mataData.row(vecPicksNew[i]);
344  timeData.append(data);
345  }
346 
347  // Copy in data from last data block. This is necessary in order to also delay channels which are not filtered
348  MatrixXd matDataOut(mataData.rows(), mataData.cols()+iOrder);
349  matDataOut.setZero();
350  matDataOut.block(0, iOrder/2, mataData.rows(), mataData.cols()) = mataData;
351 
352  if(bUseThreads) {
353  QFuture<void> future = QtConcurrent::map(timeData,
354  filterChannel);
355  future.waitForFinished();
356  } else {
357  for(int i = 0; i < timeData.size(); ++i) {
358  filterChannel(timeData[i]);
359  }
360  }
361 
362  // Do the overlap add method and store in matDataOut
363  RowVectorXd tempData;
364 
365  for(int r = 0; r < timeData.size(); r++) {
366  // Write the newly calculated filtered data to the filter data matrix. This data has a delay of iOrder/2 in front and back
367  matDataOut.row(timeData.at(r).iRow) = timeData.at(r).vecData;
368  }
369 
370  return matDataOut;
371 }
372 
373 //=============================================================================================================
374 
376 {
377  //channelDataTime.vecData = channelDataTime.first.at(i).applyConvFilter(channelDataTime.vecData, true);
378  channelDataTime.filterKernel.applyFftFilter(channelDataTime.vecData, true); //FFT Convolution for rt is not suitable. FFT make the signal filtering non causal.
379 }
380 
381 //=============================================================================================================
382 // DEFINE MEMBER METHODS
383 //=============================================================================================================
384 
385 MatrixXd FilterOverlapAdd::calculate(const MatrixXd& mataData,
386  int type,
387  double dCenterfreq,
388  double bandwidth,
389  double dTransition,
390  double dSFreq,
391  int iOrder,
392  int designMethod,
393  const RowVectorXi& vecPicks,
394  bool bFilterEnd,
395  bool bUseThreads,
396  bool bKeepOverhead)
397 {
398  // Check for size of data
399  if(mataData.cols() < iOrder){
400  qWarning() << QString("[Filter::filterData] Filter length/order is bigger than data length. Returning.");
401  return mataData;
402  }
403 
404  // Normalize cut off frequencies to nyquist
405  dCenterfreq = dCenterfreq/(dSFreq/2.0);
406  bandwidth = bandwidth/(dSFreq/2.0);
407  dTransition = dTransition/(dSFreq/2.0);
408 
409  // create filter
410  FilterKernel filter = FilterKernel("filter_kernel",
411  type,
412  iOrder,
413  dCenterfreq,
414  bandwidth,
415  dTransition,
416  dSFreq,
417  designMethod);
418 
419  return calculate(mataData,
420  filter,
421  vecPicks,
422  bFilterEnd,
423  bUseThreads,
424  bKeepOverhead);
425 }
426 
427 //=============================================================================================================
428 
429 MatrixXd FilterOverlapAdd::calculate(const MatrixXd& mataData,
430  const FilterKernel& filterKernel,
431  const RowVectorXi& vecPicks,
432  bool bFilterEnd,
433  bool bUseThreads,
434  bool bKeepOverhead)
435 {
436  int iOrder = filterKernel.getFilterOrder();
437 
438  // Check for size of data
439  if(mataData.cols() < iOrder){
440  qWarning() << "[Filter::filterData] Filter length/order is bigger than data length. Returning.";
441  return mataData;
442  }
443 
444  // Init overlaps from last block
445  if(m_matOverlapBack.cols() != iOrder || m_matOverlapBack.rows() < mataData.rows()) {
446  m_matOverlapBack.resize(mataData.rows(), iOrder);
447  m_matOverlapBack.setZero();
448  }
449 
450  if(m_matOverlapFront.cols() != iOrder || m_matOverlapFront.rows() < mataData.rows()) {
451  m_matOverlapFront.resize(mataData.rows(), iOrder);
452  m_matOverlapFront.setZero();
453  }
454 
455  // Create output matrix with size of input matrix
456  MatrixXd matDataOut(mataData.rows(), mataData.cols()+iOrder);
457  matDataOut.setZero();
458  MatrixXd sliceFiltered;
459 
460  // slice input data into data junks with proper length so that the slices are always >= the filter order
461  float fFactor = 2.0f;
462  int iSize = fFactor * iOrder;
463  int residual = mataData.cols() % iSize;
464  while(residual < iOrder) {
465  fFactor = fFactor - 0.1f;
466  iSize = fFactor * iOrder;
467  residual = mataData.cols() % iSize;
468 
469  if(iSize < iOrder) {
470  iSize = mataData.cols();
471  break;
472  }
473  }
474 
475  if(mataData.cols() > iSize) {
476  int from = 0;
477  int numSlices = ceil(float(mataData.cols())/float(iSize)); //calculate number of data slices
478 
479  for (int i = 0; i < numSlices; i++) {
480  if(i == numSlices-1) {
481  //catch the last one that might be shorter than the other blocks
482  iSize = mataData.cols() - (iSize * (numSlices -1));
483  }
484 
485  // Filter the data block. This will return data with a fitler delay of iOrder/2 in front and back
486  sliceFiltered = filterDataBlock(mataData.block(0,from,mataData.rows(),iSize),
487  vecPicks,
488  filterKernel,
489  bUseThreads);
490 
491  if(i == 0) {
492  matDataOut.block(0,0,mataData.rows(),sliceFiltered.cols()) += sliceFiltered;
493  } else {
494  matDataOut.block(0,from,mataData.rows(),sliceFiltered.cols()) += sliceFiltered;
495  }
496 
497  if(bFilterEnd && (i == 0)) {
498  matDataOut.block(0,0,matDataOut.rows(),iOrder) += m_matOverlapBack;
499  } else if (!bFilterEnd && (i == numSlices-1)) {
500  matDataOut.block(0,matDataOut.cols()-iOrder,matDataOut.rows(),iOrder) += m_matOverlapFront;
501  }
502 
503  from += iSize;
504  }
505  } else {
506  matDataOut = filterDataBlock(mataData,
507  vecPicks,
508  filterKernel,
509  bUseThreads);
510 
511  if(bFilterEnd) {
512  matDataOut.block(0,0,matDataOut.rows(),iOrder) += m_matOverlapBack;
513  } else {
514  matDataOut.block(0,matDataOut.cols()-iOrder,matDataOut.rows(),iOrder) += m_matOverlapFront;
515  }
516  }
517 
518  // Refresh the overlap matrix with the new calculated filtered data
519  m_matOverlapBack = matDataOut.block(0,matDataOut.cols()-iOrder,matDataOut.rows(),iOrder);
520  m_matOverlapFront = matDataOut.block(0,0,matDataOut.rows(),iOrder);
521 
522  if(bKeepOverhead) {
523  return matDataOut;
524  } else {
525  return matDataOut.block(0,0,matDataOut.rows(),mataData.cols());
526  }
527 }
528 
529 //=============================================================================================================
530 
532 {
533  m_matOverlapBack.resize(0,0);
534  m_matOverlapFront.resize(0,0);
535 }
Eigen::MatrixXd calculate(const Eigen::MatrixXd &matData, int type, double dCenterfreq, double dBandwidth, double dTransition, double dSFreq, int iOrder=1024, int designMethod=FilterKernel::m_designMethods.indexOf(FilterParameter("Cosine")), const Eigen::RowVectorXi &vecPicks=Eigen::RowVectorXi(), bool bFilterEnd=true, bool bUseThreads=true, bool bKeepOverhead=false)
#define FIFF_FIRST_SAMPLE
Definition: fiff_file.h:461
void applyFftFilter(Eigen::RowVectorXd &vecData, bool bKeepOverhead=false)
RTPROCESINGSHARED_EXPORT bool filterFile(QIODevice &pIODevice, QSharedPointer< FIFFLIB::FiffRawData > pFiffRawData, int type, double dCenterfreq, double dBandwidth, double dTransition, double dSFreq, int iOrder=4096, int designMethod=FilterKernel::m_designMethods.indexOf(FilterParameter("Cosine")), const Eigen::RowVectorXi &vecPicks=Eigen::RowVectorXi(), bool bUseThreads=true)
The FilterKernel class provides methods to create/design a FIR filter kernel.
Definition: filterkernel.h:132
RTPROCESINGSHARED_EXPORT Eigen::MatrixXd filterData(const Eigen::MatrixXd &mataData, const RTPROCESSINGLIB::FilterKernel &filterKernel, const Eigen::RowVectorXi &vecPicks=Eigen::RowVectorXi(), bool bUseThreads=true, bool bKeepOverhead=false)
RTPROCESINGSHARED_EXPORT Eigen::MatrixXd filterDataBlock(const Eigen::MatrixXd &mataData, const Eigen::RowVectorXi &vecPicks, const RTPROCESSINGLIB::FilterKernel &filterKernel, bool bUseThreads=true)
RTPROCESINGSHARED_EXPORT Eigen::MatrixXd filterData(const Eigen::MatrixXd &matData, int type, double dCenterfreq, double dBandwidth, double dTransition, double dSFreq, int iOrder=1024, int designMethod=FilterKernel::m_designMethods.indexOf(FilterParameter("Cosine")), const Eigen::RowVectorXi &vecPicks=Eigen::RowVectorXi(), bool bUseThreads=true, bool bKeepOverhead=false)
RTPROCESINGSHARED_EXPORT bool filterFile(QIODevice &pIODevice, QSharedPointer< FIFFLIB::FiffRawData > pFiffRawData, const RTPROCESSINGLIB::FilterKernel &filterKernel, const Eigen::RowVectorXi &vecPicks=Eigen::RowVectorXi(), bool bUseThreads=false)
FiffRawData class declaration.
MNEMath class declaration.
void prepareFilter(int iDataSize)
RTPROCESINGSHARED_EXPORT void filterChannel(FilterObject &channelDataTime)
Definition: filter.cpp:375
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
Filter declarations.