MNE-CPP  0.1.9
A Framework for Electrophysiology
cosinefilter.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //=============================================================================================================
37 // INCLUDES
38 //=============================================================================================================
39 
40 #include "cosinefilter.h"
41 
42 #define _USE_MATH_DEFINES
43 #include <math.h>
44 
45 //=============================================================================================================
46 // EIGEN INCLUDES
47 //=============================================================================================================
48 
49 //#ifndef EIGEN_FFTW_DEFAULT
50 //#define EIGEN_FFTW_DEFAULT
51 //#endif
52 
53 #include <unsupported/Eigen/FFT>
54 
55 //=============================================================================================================
56 // USED NAMESPACES
57 //=============================================================================================================
58 
59 using namespace RTPROCESSINGLIB;
60 using namespace Eigen;
61 
62 //=============================================================================================================
63 // DEFINE MEMBER METHODS
64 //=============================================================================================================
65 
67 : m_iFilterOrder(0)
68 {
69 }
70 
71 //=============================================================================================================
72 
74  float lowpass,
75  float lowpass_width,
76  float highpass,
77  float highpass_width,
78  double sFreq,
79  TPassType type)
80 {
81  #ifdef EIGEN_FFTW_DEFAULT
82  fftw_make_planner_thread_safe();
83  #endif
84 
85  m_iFilterOrder = fftLength;
86 
87  int highpasss,lowpasss;
88  int highpass_widths,lowpass_widths;
89  int k,s,w;
90  int resp_size = fftLength/2+1; //Take half because we are not interested in the conjugate complex part of the spectrum
91 
92  double pi4 = M_PI/4.0;
93  float mult,add,c;
94 
95  RowVectorXcd filterFreqResp = RowVectorXcd::Ones(resp_size);
96 
97  //Transform frequencies into samples
98  highpasss = ((resp_size-1)*highpass)/(0.5*sFreq);
99  lowpasss = ((resp_size-1)*lowpass)/(0.5*sFreq);
100 
101  lowpass_widths = ((resp_size-1)*lowpass_width)/(0.5*sFreq);
102  lowpass_widths = (lowpass_widths+1)/2;
103 
104  if (highpass_width > 0.0) {
105  highpass_widths = ((resp_size-1)*highpass_width)/(0.5*sFreq);
106  highpass_widths = (highpass_widths+1)/2;
107  }
108  else
109  highpass_widths = 3;
110 
111  //Calculate filter freq response - use cosine
112  //Build high pass filter
113  if(type != LPF) {
114  if (highpasss > highpass_widths + 1) {
115  w = highpass_widths;
116  mult = 1.0/w;
117  add = 3.0;
118 
119  for (k = 0; k < resp_size; k++)
120  filterFreqResp(k) = 0.0;
121 
122  for (k = -w+1, s = highpasss-w+1; k < w; k++, s++) {
123  if (s >= 0 && s < resp_size) {
124  c = cos(pi4*(k*mult+add));
125  filterFreqResp(s) = filterFreqResp(s).real()*c*c;
126  }
127  }
128  }
129  }
130 
131  //Build low pass filter
132  if(type != HPF) {
133  if (lowpass_widths > 0) {
134  w = lowpass_widths;
135  mult = 1.0/w;
136  add = 1.0;
137 
138  for (k = -w+1, s = lowpasss-w+1; k < w; k++, s++) {
139  if (s >= 0 && s < resp_size) {
140  c = cos(pi4*(k*mult+add));
141  filterFreqResp(s) = filterFreqResp(s).real()*c*c;
142  }
143  }
144 
145  for (k = s; k < resp_size; k++)
146  filterFreqResp(k) = 0.0;
147  }
148  else {
149  for (k = lowpasss; k < resp_size; k++)
150  filterFreqResp(k) = 0.0;
151  }
152  }
153 
154  m_vecFftCoeff = filterFreqResp;
155 
156  //Generate windowed impulse response - invert fft coeeficients to time domain
157  Eigen::FFT<double> fft;
158  fft.SetFlag(fft.HalfSpectrum);
159 
160  //invert to time domain and
161  fft.inv(m_vecCoeff, filterFreqResp);/*
162  m_vecCoeff = m_vecCoeff.segment(0,1024).eval();
163 
164  //window/zero-pad m_vecCoeff to m_iFftLength
165  RowVectorXd vecCoeffZeroPad = RowVectorXd::Zero(fftLength);
166  vecCoeffZeroPad.head(m_vecCoeff.cols()) = m_vecCoeff;
167 
168  //fft-transform filter coeffs
169  m_vecFftCoeff = RowVectorXcd::Zero(fftLength);
170  fft.fwd(m_vecFftCoeff,vecCoeffZeroPad);*/
171 }
172 
173 //=============================================================================================================
Declaration of the CosineFilter class.
Eigen::RowVectorXd m_vecCoeff
Definition: cosinefilter.h:101
Eigen::RowVectorXcd m_vecFftCoeff
Definition: cosinefilter.h:100