MNE-CPP  0.1.9
A Framework for Electrophysiology
filterkernel.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "filterkernel.h"
42 
43 #include <utils/mnemath.h>
44 
45 #include "parksmcclellan.h"
46 #include "cosinefilter.h"
47 
48 #include <iostream>
49 
50 //=============================================================================================================
51 // QT INCLUDES
52 //=============================================================================================================
53 
54 #include <QDebug>
55 
56 //=============================================================================================================
57 // EIGEN INCLUDES
58 //=============================================================================================================
59 
60 #include <Eigen/SparseCore>
61 //#ifndef EIGEN_FFTW_DEFAULT
62 //#define EIGEN_FFTW_DEFAULT
63 //#endif
64 #include <unsupported/Eigen/FFT>
65 
66 //=============================================================================================================
67 // USED NAMESPACES
68 //=============================================================================================================
69 
70 using namespace RTPROCESSINGLIB;
71 using namespace Eigen;
72 using namespace UTILSLIB;
73 
74 
75 
76 //=============================================================================================================
77 // INIT STATIC MEMBERS
78 //=============================================================================================================
79 
80 QVector<RTPROCESSINGLIB::FilterParameter> FilterKernel::m_designMethods ({
81  FilterParameter(QString("Cosine"), QString("A cosine filter")),
82  FilterParameter(QString("Tschebyscheff"), QString("A tschebyscheff filter"))
83 // FilterParameter(QString("External"), QString("An external filter"))
84 });
85 QVector<RTPROCESSINGLIB::FilterParameter> FilterKernel::m_filterTypes ({
86  FilterParameter(QString("LPF"), QString("An LPF filter")),
87  FilterParameter(QString("HPF"), QString("An HPF filter")),
88  FilterParameter(QString("BPF"), QString("A BPF filter")),
89  FilterParameter(QString("NOTCH"), QString("A NOTCH filter")),
90  FilterParameter(QString("UNKNOWN"), QString("An UNKNOWN filter"))
91 });
92 
93 //=============================================================================================================
94 // DEFINE MEMBER METHODS
95 //=============================================================================================================
96 
98 : m_sFreq(1000)
99 , m_dCenterFreq(0.5)
100 , m_dBandwidth(0.1)
101 , m_dParksWidth(0.1)
102 , m_dLowpassFreq(40)
103 , m_dHighpassFreq(4)
104 , m_iFilterOrder(80)
105 , m_iDesignMethod(m_designMethods.indexOf(FilterParameter("Cosine")))
106 , m_iFilterType(m_filterTypes.indexOf(FilterParameter("BPF")))
107 , m_sFilterName("Unknown")
108 , m_sFilterShortDescription("")
109 {
110  designFilter();
111 }
112 
113 //=============================================================================================================
114 
115 FilterKernel::FilterKernel(const QString& sFilterName,
116  int iFilterType,
117  int iOrder,
118  double dCenterfreq,
119  double dBandwidth,
120  double dParkswidth,
121  double dSFreq,
122  int iDesignMethod)
123 : m_sFreq(dSFreq)
124 , m_dCenterFreq(dCenterfreq)
125 , m_dBandwidth(dBandwidth)
126 , m_dParksWidth(dParkswidth)
127 , m_iFilterOrder(iOrder)
128 , m_iDesignMethod(iDesignMethod)
129 , m_iFilterType(iFilterType)
130 , m_sFilterName(sFilterName)
131 , m_sFilterShortDescription()
132 {
133  if(iOrder < 9) {
134  qWarning() << "[FilterKernel::FilterKernel] Less than 9 taps were provided. Setting number of taps to 9.";
135  }
136 
137  designFilter();
138 }
139 
140 //=============================================================================================================
141 
142 void FilterKernel::prepareFilter(int iDataSize)
143 {
144  int iFftLength, exp;
145 
146  iFftLength = iDataSize + m_vecCoeff.cols();
147  exp = ceil(MNEMath::log2(iFftLength));
148  iFftLength = pow(2, exp);
149 
150  // Transform coefficients anew if needed
151  if(m_vecCoeff.cols() != (iFftLength/2+1)) {
152  fftTransformCoeffs(iFftLength);
153  }
154 }
155 
156 //=============================================================================================================
157 
158 RowVectorXd FilterKernel::applyConvFilter(const RowVectorXd& vecData,
159  bool bKeepOverhead) const
160 {
161  //Do zero padding or mirroring depending on user input
162  RowVectorXd vecDataZeroPad = RowVectorXd::Zero(2*m_vecCoeff.cols() + vecData.cols());
163  RowVectorXd vecFilteredTime = RowVectorXd::Zero(2*m_vecCoeff.cols() + vecData.cols());
164 
165  vecDataZeroPad.segment(m_vecCoeff.cols(), vecData.cols()) = vecData;
166 
167  //Do the convolution
168  for(int i = m_vecCoeff.cols(); i < vecFilteredTime.cols(); i++) {
169  vecFilteredTime(i-m_vecCoeff.cols()) = vecDataZeroPad.segment(i-m_vecCoeff.cols(),m_vecCoeff.cols()) * m_vecCoeff.transpose();
170  }
171 
172  //Return filtered data
173  if(!bKeepOverhead) {
174  return vecFilteredTime.segment(m_vecCoeff.cols()/2, vecData.cols());
175  }
176 
177  return vecFilteredTime.head(vecData.cols()+m_vecCoeff.cols());
178 }
179 
180 //=============================================================================================================
181 
182 void FilterKernel::applyFftFilter(RowVectorXd& vecData,
183  bool bKeepOverhead)
184 {
185  #ifdef EIGEN_FFTW_DEFAULT
186  fftw_make_planner_thread_safe();
187  #endif
188 
189  // Make sure we always have the correct FFT length for the given input data and filter overlap
190  int iFftLength = vecData.cols() + m_vecCoeff.cols();
191  int exp = ceil(MNEMath::log2(iFftLength));
192  iFftLength = pow(2, exp);
193 
194  // Transform coefficients anew if needed
195  if(m_vecFftCoeff.cols() != (iFftLength/2+1)) {
196  fftTransformCoeffs(iFftLength);
197  }
198 
199  //generate fft object
200  Eigen::FFT<double> fft;
201  fft.SetFlag(fft.HalfSpectrum);
202 
203  // Zero padd if necessary. Please note: The zero padding in Eigen's FFT is only working for column vectors -> We have to zero pad manually here
204  int iOriginalSize = vecData.cols();
205  if (vecData.cols() < iFftLength) {
206  int iResidual = iFftLength - vecData.cols();
207  vecData.conservativeResize(iFftLength);
208  vecData.tail(iResidual).setZero();
209  }
210 
211  //fft-transform data sequence
212  RowVectorXcd vecFreqData;
213  fft.fwd(vecFreqData, vecData, iFftLength);
214 
215  //perform frequency-domain filtering
216  vecFreqData = m_vecFftCoeff.array() * vecFreqData.array();
217 
218  //inverse-FFT
219  fft.inv(vecData, vecFreqData);
220 
221  //Return filtered data
222  if(!bKeepOverhead) {
223  vecData = vecData.segment(m_vecCoeff.cols()/2, iOriginalSize).eval();
224  } else {
225  vecData = vecData.head(iOriginalSize + m_vecCoeff.cols()).eval();
226  }
227 }
228 
229 //=============================================================================================================
230 
231 QString FilterKernel::getName() const
232 {
233  return m_sFilterName;
234 }
235 
236 //=============================================================================================================
237 
238 void FilterKernel::setName(const QString& sFilterName)
239 {
240  m_sFilterName = sFilterName;
241 }
242 
243 //=============================================================================================================
244 
245 double FilterKernel::getSamplingFrequency() const
246 {
247  return m_sFreq;
248 }
249 
250 //=============================================================================================================
251 
252 void FilterKernel::setSamplingFrequency(double dSFreq)
253 {
254  m_sFreq = dSFreq;
255 }
256 
257 //=============================================================================================================
258 
259 int FilterKernel::getFilterOrder() const
260 {
261  return m_iFilterOrder;
262 }
263 
264 //=============================================================================================================
265 
266 void FilterKernel::setFilterOrder(int iOrder)
267 {
268  m_iFilterOrder = iOrder;
269 }
270 
271 //=============================================================================================================
272 
273 double FilterKernel::getCenterFrequency() const
274 {
275  return m_dCenterFreq;
276 }
277 
278 //=============================================================================================================
279 
280 void FilterKernel::setCenterFrequency(double dCenterFreq)
281 {
282  m_dCenterFreq = dCenterFreq;
283 }
284 
285 //=============================================================================================================
286 
287 double FilterKernel::getBandwidth() const
288 {
289  return m_dBandwidth;
290 }
291 
292 //=============================================================================================================
293 
294 void FilterKernel::setBandwidth(double dBandwidth)
295 {
296  m_dBandwidth = dBandwidth;
297 }
298 
299 //=============================================================================================================
300 
301 double FilterKernel::getParksWidth() const
302 {
303  return m_dParksWidth;
304 }
305 
306 //=============================================================================================================
307 
308 void FilterKernel::setParksWidth(double dParksWidth)
309 {
310  m_dParksWidth = dParksWidth;
311 }
312 
313 //=============================================================================================================
314 
315 double FilterKernel::getHighpassFreq() const
316 {
317  return m_dHighpassFreq;
318 }
319 
320 //=============================================================================================================
321 
322 void FilterKernel::setHighpassFreq(double dHighpassFreq)
323 {
324  m_dHighpassFreq = dHighpassFreq;
325 }
326 
327 //=============================================================================================================
328 
329 double FilterKernel::getLowpassFreq() const
330 {
331  return m_dLowpassFreq;
332 }
333 
334 //=============================================================================================================
335 
336 void FilterKernel::setLowpassFreq(double dLowpassFreq)
337 {
338  m_dLowpassFreq = dLowpassFreq;
339 }
340 
341 //=============================================================================================================
342 
343 Eigen::RowVectorXd FilterKernel::getCoefficients() const
344 {
345  return m_vecCoeff;
346 }
347 
348 //=============================================================================================================
349 
350 void FilterKernel::setCoefficients(const Eigen::RowVectorXd& vecCoeff)
351 {
352  m_vecCoeff = vecCoeff;
353 }
354 
355 //=============================================================================================================
356 
357 Eigen::RowVectorXcd FilterKernel::getFftCoefficients() const
358 {
359  return m_vecFftCoeff;
360 }
361 
362 //=============================================================================================================
363 
364 void FilterKernel::setFftCoefficients(const Eigen::RowVectorXcd& vecFftCoeff)
365 {
366  m_vecFftCoeff = vecFftCoeff;
367 }
368 
369 //=============================================================================================================
370 
371 bool FilterKernel::fftTransformCoeffs(int iFftLength)
372 {
373  #ifdef EIGEN_FFTW_DEFAULT
374  fftw_make_planner_thread_safe();
375  #endif
376 
377  if(m_vecCoeff.cols() > iFftLength) {
378  std::cout <<"[FilterKernel::fftTransformCoeffs] The number of filter taps is bigger than the FFT length."<< std::endl;
379  return false;
380  }
381 
382  //generate fft object
383  Eigen::FFT<double> fft;
384  fft.SetFlag(fft.HalfSpectrum);
385 
386  // Zero padd if necessary. Please note: The zero padding in Eigen's FFT is only working for column vectors -> We have to zero pad manually here
387  RowVectorXd vecInputFft;
388  if (m_vecCoeff.cols() < iFftLength) {
389  vecInputFft.setZero(iFftLength);
390  vecInputFft.block(0,0,1,m_vecCoeff.cols()) = m_vecCoeff;
391  } else {
392  vecInputFft = m_vecCoeff;
393  }
394 
395  //fft-transform filter coeffs
396  RowVectorXcd vecFreqData;
397  fft.fwd(vecFreqData, vecInputFft, iFftLength);
398  m_vecFftCoeff = vecFreqData;;
399 
400  return true;
401 }
402 
403 //=============================================================================================================
404 
405 void FilterKernel::designFilter()
406 {
407  // Make sure we only use a minimum needed FFT size
408  int iFftLength = m_iFilterOrder;
409  int exp = ceil(MNEMath::log2(iFftLength));
410  iFftLength = pow(2, exp);
411 
412  switch(m_iDesignMethod) {
413  case 1: {
414  ParksMcClellan filter(m_iFilterOrder,
415  m_dCenterFreq,
416  m_dBandwidth,
417  m_dParksWidth,
418  static_cast<ParksMcClellan::TPassType>(m_iFilterType));
419  m_vecCoeff = filter.FirCoeff;
420 
421  //fft-transform m_vecCoeff in order to be able to perform frequency-domain filtering
422  fftTransformCoeffs(iFftLength);
423 
424  break;
425  }
426 
427  case 0: {
428  CosineFilter filtercos;
429 
430  switch(m_iFilterType) {
431  case 0:
432  filtercos = CosineFilter(iFftLength,
433  (m_dCenterFreq)*(m_sFreq/2.),
434  m_dParksWidth*(m_sFreq/2),
435  (m_dCenterFreq)*(m_sFreq/2),
436  m_dParksWidth*(m_sFreq/2),
437  m_sFreq,
438  static_cast<CosineFilter::TPassType>(m_iFilterType));
439 
440  break;
441 
442  case 1:
443  filtercos = CosineFilter(iFftLength,
444  (m_dCenterFreq)*(m_sFreq/2),
445  m_dParksWidth*(m_sFreq/2),
446  (m_dCenterFreq)*(m_sFreq/2),
447  m_dParksWidth*(m_sFreq/2),
448  m_sFreq,
449  static_cast<CosineFilter::TPassType>(m_iFilterType));
450 
451  break;
452 
453  case 2:
454  filtercos = CosineFilter(iFftLength,
455  (m_dCenterFreq + m_dBandwidth/2)*(m_sFreq/2),
456  m_dParksWidth*(m_sFreq/2),
457  (m_dCenterFreq - m_dBandwidth/2)*(m_sFreq/2),
458  m_dParksWidth*(m_sFreq/2),
459  m_sFreq,
460  static_cast<CosineFilter::TPassType>(m_iFilterType));
461 
462  break;
463  }
464 
465  //This filter is designed in the frequency domain, hence the time domain impulse response need to be shortend by the users dependent number of taps
466  m_vecCoeff.resize(m_iFilterOrder);
467 
468  m_vecCoeff.head(m_iFilterOrder/2) = filtercos.m_vecCoeff.tail(m_iFilterOrder/2);
469  m_vecCoeff.tail(m_iFilterOrder/2) = filtercos.m_vecCoeff.head(m_iFilterOrder/2);
470 
471  //Now generate the fft version of the shortened impulse response
472  fftTransformCoeffs(iFftLength);
473 
474  break;
475  }
476  }
477 
478  switch(m_iFilterType) {
479  case 0:
480  m_dLowpassFreq = 0;
481  m_dHighpassFreq = m_dCenterFreq*(m_sFreq/2);
482  break;
483 
484  case 1:
485  m_dLowpassFreq = m_dCenterFreq*(m_sFreq/2);
486  m_dHighpassFreq = 0;
487  break;
488 
489  case 2:
490  m_dLowpassFreq = (m_dCenterFreq + m_dBandwidth/2)*(m_sFreq/2);
491  m_dHighpassFreq = (m_dCenterFreq - m_dBandwidth/2)*(m_sFreq/2);
492  break;
493  }
494  getShortDescription();
495 }
496 
497 //=============================================================================================================
498 
499 QString FilterKernel::getShortDescription() const
500 {
501  QString description(m_designMethods.at(m_iDesignMethod).getName() + " - " + \
502  QString::number(m_dHighpassFreq,'g',4) + "Hz to " + QString::number(m_dLowpassFreq,'g',4) + "Hz - " \
503  "Ord: " + QString::number(m_iFilterOrder));
504  return description;
505 }
506 
507 //=============================================================================================================
508 
509 RTPROCESSINGLIB::FilterParameter FilterKernel::getDesignMethod() const
510 {
511  if(m_iDesignMethod < 0){
512  return m_designMethods.at(0);
513  }
514  return m_designMethods.at(m_iDesignMethod);
515 }
516 
517 //=============================================================================================================
518 
519 RTPROCESSINGLIB::FilterParameter FilterKernel::getFilterType() const
520 {
521  if(m_iFilterType < 0){
522  return m_filterTypes.at(0);
523  }
524  return m_filterTypes.at(m_iFilterType);
525 }
526 
527 //=============================================================================================================
528 
529 void FilterKernel::setDesignMethod(int iDesignMethod)
530 {
531  if(iDesignMethod < 0){
532  m_iDesignMethod = 0;
533  } else {
534  m_iDesignMethod = iDesignMethod;
535  }
536 }
537 
538 //=============================================================================================================
539 
540 void FilterKernel::setFilterType(int iFilterType)
541 {
542  if(iFilterType < 0){
543  m_iFilterType = 0;
544  } else {
545  m_iFilterType = iFilterType;
546  }
547 }
548 
549 //=============================================================================================================
550 
552 :FilterParameter("Unknown", "")
553 {
554 }
555 
556 //=============================================================================================================
557 
559 :FilterParameter(sName,"")
560 {
561 }
562 
563 //=============================================================================================================
564 
566  QString sDescription)
567 : m_sName(sName)
568 , m_sDescription(sDescription)
569 {
570 
571 }
572 
573 //=============================================================================================================
574 
576 {
577  return m_sName;
578 }
579 
580 //=============================================================================================================
FilterKernel()
FilterKernel creates a default FilterKernel object.
The FilterKernel class represents a filter object that generates the FIR filter coefficients using Pa...
Declaration of the CosineFilter class.
The ParksMcClellan class provides the ParksMcClellan filter desing algorithm.
void applyFftFilter(Eigen::RowVectorXd &vecData, bool bKeepOverhead=false)
Eigen::RowVectorXd m_vecCoeff
Definition: cosinefilter.h:101
Creates a cosine filter response in the frequency domain.
Definition: cosinefilter.h:68
static QVector< FilterParameter > m_filterTypes
Definition: filterkernel.h:241
Eigen::RowVectorXd applyConvFilter(const Eigen::RowVectorXd &vecData, bool bKeepOverhead=false) const
MNEMath class declaration.
void prepareFilter(int iDataSize)
static QVector< FilterParameter > m_designMethods
Definition: filterkernel.h:240
The FilterParameter class.
Definition: filterkernel.h:83