MNE-CPP  0.1.9
A Framework for Electrophysiology
spectrogram.cpp
Go to the documentation of this file.
1 //=============================================================================================================
34 //=============================================================================================================
35 // INCLUDES
36 //=============================================================================================================
37 
38 #include "spectrogram.h"
39 
40 //=============================================================================================================
41 // EIGEN INCLUDES
42 //=============================================================================================================
43 
44 #include <Eigen/SparseCore>
45 #include <unsupported/Eigen/FFT>
46 
47 //=============================================================================================================
48 // QT INCLUDES
49 //=============================================================================================================
50 
51 #include <QDebug>
52 #include <QElapsedTimer>
53 #include <QThread>
54 #include <QtConcurrent>
55 
56 //=============================================================================================================
57 // USED NAMESPACES
58 //=============================================================================================================
59 
60 using namespace UTILSLIB;
61 using namespace Eigen;
62 
63 //=============================================================================================================
64 // DEFINE MEMBER METHODS
65 //=============================================================================================================
66 
67 MatrixXd Spectrogram::makeSpectrogram(VectorXd signal, qint32 windowSize = 0)
68 {
69  //QElapsedTimer timer;
70  //timer.start();
71 
72  signal.array() -= signal.mean();
73  QList<SpectogramInputData> lData;
74  int iThreadSize = QThread::idealThreadCount()*2;
75  int iStepsSize = signal.rows()/iThreadSize;
76  int iResidual = signal.rows()%iThreadSize;
77 
78  SpectogramInputData dataTemp;
79  dataTemp.vecInputData = signal;
80  dataTemp.window_size = windowSize;
81  if(dataTemp.window_size == 0) {
82  dataTemp.window_size = signal.rows()/15;
83  }
84 
85  for (int i = 0; i < iThreadSize; ++i) {
86  dataTemp.iRangeLow = i*iStepsSize;
87  dataTemp.iRangeHigh = i*iStepsSize+iStepsSize;
88  lData.append(dataTemp);
89  }
90 
91  dataTemp.iRangeLow = iThreadSize*iStepsSize;
92  dataTemp.iRangeHigh = iThreadSize*iStepsSize+iResidual;
93  lData.append(dataTemp);
94 
95  QFuture<MatrixXd> resultMat = QtConcurrent::mappedReduced(lData,
96  compute,
97  reduce);
98  resultMat.waitForFinished();
99 
100  //qDebug() << "Spectrogram::make_spectrogram - timer.elapsed()" << timer.elapsed();
101  return resultMat.result();
102 }
103 
104 //=============================================================================================================
105 
106 VectorXd Spectrogram::gaussWindow(qint32 sample_count, qreal scale, quint32 translation)
107 {
108  VectorXd gauss = VectorXd::Zero(sample_count);
109 
110  for(qint32 n = 0; n < sample_count; n++)
111  {
112  qreal t = (qreal(n) - translation) / scale;
113  gauss[n] = exp(-3.14 * pow(t, 2))*pow(sqrt(scale),(-1))*pow(qreal(2),(0.25));
114  }
115 
116  return gauss;
117 }
118 
119 //=============================================================================================================
120 
121 MatrixXd Spectrogram::compute(const SpectogramInputData& inputData)
122 {
123  #ifdef EIGEN_FFTW_DEFAULT
124  fftw_make_planner_thread_safe();
125  #endif
126 
127  Eigen::FFT<double> fft;
128  MatrixXd tf_matrix = MatrixXd::Zero(inputData.vecInputData.rows()/2, inputData.vecInputData.rows());
129  VectorXd envelope, windowed_sig, real_coeffs;
130  VectorXcd fft_win_sig;
131  qint32 window_size = inputData.window_size;
132 
133  for(quint32 translate = inputData.iRangeLow; translate < inputData.iRangeHigh; translate++) {
134  envelope = gaussWindow(inputData.vecInputData.rows(), window_size, translate);
135 
136  windowed_sig = VectorXd::Zero(inputData.vecInputData.rows());
137  fft_win_sig = VectorXcd::Zero(inputData.vecInputData.rows());
138 
139  windowed_sig = inputData.vecInputData.array() * envelope.array();\
140 
141  fft.fwd(fft_win_sig, windowed_sig);
142 
143  real_coeffs = fft_win_sig.segment(0,inputData.vecInputData.rows()/2).array().abs2();
144 
145  tf_matrix.col(translate) = real_coeffs;
146  }
147 
148  return tf_matrix;
149 }
150 
151 //=============================================================================================================
152 
153 void Spectrogram::reduce(MatrixXd &resultData,
154  const MatrixXd &data)
155 {
156  if(resultData.size() == 0) {
157  resultData = data;
158  } else {
159  resultData += data;
160  }
161 }
Declaration of spectrogram class.
static Eigen::MatrixXd makeSpectrogram(Eigen::VectorXd signal, qint32 windowSize)
Definition: spectrogram.cpp:67