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