MNE-CPP  0.1.9
A Framework for Electrophysiology
rtnoise.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //=============================================================================================================
37 // INCLUDES
38 //=============================================================================================================
39 
40 #include "rtnoise.h"
41 
42 #include <iostream>
43 #include <fiff/fiff_cov.h>
44 
45 //=============================================================================================================
46 // QT INCLUDES
47 //=============================================================================================================
48 
49 #include <QDebug>
50 
51 //=============================================================================================================
52 // USED NAMESPACES
53 //=============================================================================================================
54 
55 using namespace RTPROCESSINGLIB;
56 using namespace FIFFLIB;
57 using namespace Eigen;
58 using namespace UTILSLIB;
59 
60 //=============================================================================================================
61 // DEFINE MEMBER METHODS
62 //=============================================================================================================
63 
64 RtNoise::RtNoise(qint32 p_iMaxSamples,
65  FiffInfo::SPtr p_pFiffInfo,
66  qint32 p_dataLen,
67  QObject *parent)
68 : QThread(parent)
69 , m_iFftLength(p_iMaxSamples)
70 , m_pFiffInfo(p_pFiffInfo)
71 , m_dataLength(p_dataLen)
72 , m_bIsRunning(false)
73 , m_iNumOfBlocks(0)
74 , m_iBlockSize(0)
75 , m_iSensors(0)
76 , m_iBlockIndex(0)
77 {
78  qRegisterMetaType<Eigen::MatrixXd>("Eigen::MatrixXd");
79  //qRegisterMetaType<QVector<double> >("QVector<double>");
80 
81  m_Fs = m_pFiffInfo->sfreq;
82 
83  m_bSendDataToBuffer = true;
84 
85  m_fWin.clear();
86 
87  //create a hanning window
88  m_fWin = hanning(m_iFftLength,0);
89 
90  qDebug()<<"Hanning window is created.";
91 }
92 
93 //=============================================================================================================
94 
96 {
97  if(this->isRunning()){
98  stop();
99  }
100 }
101 
102 //=============================================================================================================
103 
104 QVector <float> RtNoise::hanning(int N, short itype)
105 {
106  int half, i, idx, n;
107  QVector <float> w(N,0.0);
108 
109  if(itype==1) //periodic function
110  n = N-1;
111  else
112  n = N;
113 
114  if(n%2==0)
115  {
116  half = n/2;
117  for(i=0; i<half; i++) //CALC_HANNING Calculates Hanning window samples.
118  w[i] = 0.5 * (1 - cos(2*3.14159265*(i+1) / (n+1)));
119 
120  idx = half-1;
121  for(i=half; i<n; i++) {
122  w[i] = w[idx];
123  idx--;
124  }
125  }
126  else
127  {
128  half = (n+1)/2;
129  for(i=0; i<half; i++) //CALC_HANNING Calculates Hanning window samples.
130  w[i] = 0.5 * (1 - cos(23.14159265*(i+1) / (n+1)));
131 
132  idx = half-2;
133  for(i=half; i<n; i++) {
134  w[i] = w[idx];
135  idx--;
136  }
137  }
138 
139  if(itype==1) //periodic function
140  {
141  for(i=N-1; i>=1; i--)
142  w[i] = w[i-1];
143  w[0] = 0.0;
144  }
145  return(w);
146 }
147 
148 //=============================================================================================================
149 
150 void RtNoise::append(const MatrixXd &p_DataSegment)
151 {
152  if(!m_pCircularBuffer)
154 
155  if (m_bSendDataToBuffer)
156  m_pCircularBuffer->push(p_DataSegment);
157 }
158 
159 //=============================================================================================================
160 
162 {
163  //Check if the thread is already or still running. This can happen if the start button is pressed immediately after the stop button was pressed. In this case the stopping process is not finished yet but the start process is initiated.
164  if(this->isRunning())
165  QThread::wait();
166 
167  m_bIsRunning = true;
168  QThread::start();
169 
170  return true;
171 }
172 
173 //=============================================================================================================
174 
176 {
177  m_bIsRunning = false;
178 
179  m_pCircularBuffer->clear();
180 
181  qDebug()<<" RtNoise Thread is stopped.";
182 
183  return true;
184 }
185 
186 //=============================================================================================================
187 
189 {
190  #ifdef EIGEN_FFTW_DEFAULT
191  fftw_make_planner_thread_safe();
192  #endif
193 
194  bool FirstStart = true;
195  MatrixXd block;
196 
197  while(m_bIsRunning) {
198  if(m_pCircularBuffer) {
199  if(m_pCircularBuffer->pop(block)) {
200  if(FirstStart){
201  //init the circ buffer and parameters
202  if(m_dataLength < 0) m_dataLength = 10;
203  m_iNumOfBlocks = m_dataLength;//60;
204  m_iBlockSize = block.cols();
205  m_iSensors = block.rows();
206 
207  m_matCircBuf.resize(m_iSensors,m_iNumOfBlocks*m_iBlockSize);
208 
209  m_iBlockIndex = 0;
210  FirstStart = false;
211  }
212  //concate blocks
213  for (int i=0; i< m_iSensors; i++)
214  for (int j=0; j< m_iBlockSize; j++)
215  m_matCircBuf(i,j+m_iBlockIndex*m_iBlockSize) = block(i,j);
216 
217  m_iBlockIndex ++;
218  if (m_iBlockIndex >= m_iNumOfBlocks){
219 
220  //m_pCircularBuffer.clear(); //empty the buffer
221 
222  m_bSendDataToBuffer = false;
223  //stop collect block and start to calculate the spectrum
224  m_iBlockIndex = 0;
225 
226  MatrixXd sum_psdx = MatrixXd::Zero(m_iSensors,m_iFftLength/2+1);
227 
228  int nb = floor(m_iNumOfBlocks*m_iBlockSize/m_iFftLength)+1;
229  qDebug()<<"nb"<<nb<<"NumOfBlocks"<<m_iNumOfBlocks<<"BlockSize"<<m_iBlockSize;
230  MatrixXd t_mat(m_iSensors,m_iFftLength);
231  MatrixXd t_psdx(m_iSensors,m_iFftLength/2+1);
232  for (int n = 0; n<nb; n++){
233  //collect a data block with data length of m_iFftLength;
234  if(n==nb-1)
235  {
236  for(qint32 ii=0; ii<m_iSensors; ii++)
237  for(qint32 jj=0; jj<m_iFftLength; jj++)
238  if(jj+n*m_iFftLength<m_iNumOfBlocks*m_iBlockSize)
239  t_mat(ii,jj) = m_matCircBuf(ii,jj+n*m_iFftLength);
240  else
241  t_mat(ii,jj) = 0.0;
242 
243  }
244  else
245  {
246  for(qint32 ii=0; ii<m_iSensors; ii++)
247  for(qint32 jj=0; jj<m_iFftLength; jj++)
248  t_mat(ii,jj) = m_matCircBuf(ii,jj+n*m_iFftLength);
249  }
250 
251  //FFT calculation by row
252  for(qint32 i = 0; i < t_mat.rows(); i++){
253  RowVectorXd data;
254 
255  data = t_mat.row(i);
256 
257  //zero-pad data to m_iFftLength
258  RowVectorXd vecDataZeroPad = RowVectorXd::Zero(m_iFftLength);
259  vecDataZeroPad.head(data.cols()) = data;
260 
261  for (qint32 lk = 0; lk<m_iFftLength; lk++)
262  vecDataZeroPad[lk] = vecDataZeroPad[lk]*m_fWin[lk];
263 
264  //generate fft object
265  Eigen::FFT<double> fft;
266  fft.SetFlag(fft.HalfSpectrum);
267 
268  //fft-transform data sequence
269  RowVectorXcd vecFreqData(m_iFftLength/2+1);
270  fft.fwd(vecFreqData,vecDataZeroPad);
271 
272  // calculate spectrum from FFT
273  for(qint32 j=0; j<m_iFftLength/2+1;j++)
274  {
275  double mag_abs = sqrt(vecFreqData(j).real()* vecFreqData(j).real() + vecFreqData(j).imag()*vecFreqData(j).imag());
276  double spower = (1.0/(m_Fs*m_iFftLength))* mag_abs;
277  if (j>0&&j<m_iFftLength/2) spower = 2.0*spower;
278  sum_psdx(i,j) = sum_psdx(i,j) + spower;
279  }
280  }//row computing is done
281  }//nb
282 
283  //DB-calculation
284  for(qint32 ii=0; ii<m_iSensors; ii++)
285  for(qint32 jj=0; jj<m_iFftLength/2+1; jj++)
286  t_psdx(ii,jj) = 10.0*log10(sum_psdx(ii,jj)/nb);
287 
288  qDebug()<<"Send spectrum to Noise Estimator";
289  emit SpecCalculated(t_psdx); //send back the spectrum result
290 
291  m_bSendDataToBuffer = true;
292  }
293  }
294  }
295 
296  }
297 }
298 
QSharedPointer< FiffInfo > SPtr
Definition: fiff_info.h:87
RtNoise class declaration.
virtual bool stop()
Definition: rtnoise.cpp:175
void append(const Eigen::MatrixXd &p_DataSegment)
Definition: rtnoise.cpp:150
virtual bool start()
Definition: rtnoise.cpp:161
virtual void run()
Definition: rtnoise.cpp:188
RtNoise(qint32 p_iMaxSamples, FIFFLIB::FiffInfo::SPtr p_pFiffInfo, qint32 p_dataLen, QObject *parent=0)
Definition: rtnoise.cpp:64
FiffCov class declaration.
The TEMPLATE CIRCULAR BUFFER provides a template for thread safe circular buffers.
void SpecCalculated(Eigen::MatrixXd)
QSharedPointer< CircularBuffer > SPtr