MNE-CPP  0.1.9
A Framework for Electrophysiology
spectral.cpp
Go to the documentation of this file.
1 //=============================================================================================================
40 //=============================================================================================================
41 // INCLUDES
42 //=============================================================================================================
43 
44 #include "spectral.h"
45 #include "math.h"
46 
47 //=============================================================================================================
48 // EIGEN INCLUDES
49 //=============================================================================================================
50 
51 #include <unsupported/Eigen/FFT>
52 
53 //=============================================================================================================
54 // QT INCLUDES
55 //=============================================================================================================
56 
57 #include <QtMath>
58 #include <QtConcurrent>
59 #include <QVector>
60 
61 //=============================================================================================================
62 // USED NAMESPACES
63 //=============================================================================================================
64 
65 using namespace UTILSLIB;
66 using namespace Eigen;
67 
68 //=============================================================================================================
69 // DEFINE MEMBER METHODS
70 //=============================================================================================================
71 
72 MatrixXcd Spectral::computeTaperedSpectraRow(const RowVectorXd &vecData,
73  const MatrixXd &matTaper,
74  int iNfft)
75 {
76  //qDebug() << "Spectral::computeTaperedSpectra Matrixwise";
77  FFT<double> fft;
78  fft.SetFlag(fft.HalfSpectrum);
79 
80  //Check inputs
81  if (vecData.cols() != matTaper.cols() || iNfft < vecData.cols()) {
82  return MatrixXcd();
83  }
84 
85  //FFT for freq domain returning the half spectrum
86  RowVectorXd vecInputFFT;
87  RowVectorXcd vecTmpFreq;
88  MatrixXcd matTapSpectrum(matTaper.rows(), int(floor(iNfft / 2.0)) + 1);
89  for (int i=0; i < matTaper.rows(); i++) {
90  vecInputFFT = vecData.cwiseProduct(matTaper.row(i));
91  fft.fwd(vecTmpFreq, vecInputFFT, iNfft);
92  matTapSpectrum.row(i) = vecTmpFreq;
93  }
94 
95  return matTapSpectrum;
96 }
97 
98 //=============================================================================================================
99 
100 QVector<MatrixXcd> Spectral::computeTaperedSpectraMatrix(const MatrixXd &matData,
101  const MatrixXd &matTaper,
102  int iNfft,
103  bool bUseThreads)
104 {
105  #ifdef EIGEN_FFTW_DEFAULT
106  fftw_make_planner_thread_safe();
107  #endif
108 
109  QVector<MatrixXcd> finalResult;
110 
111  if(!bUseThreads) {
112  // Sequential
113 // QElapsedTimer timer;
114 // int iTime = 0;
115 // int iTimeAll = 0;
116 // timer.start();
117 
118  FFT<double> fft;
119  fft.SetFlag(fft.HalfSpectrum);
120 
121  RowVectorXd vecInputFFT, rowData;
122  RowVectorXcd vecTmpFreq;
123 
124  MatrixXcd matTapSpectrum(matTaper.rows(), int(floor(iNfft / 2.0)) + 1);
125  int j;
126  for (int i = 0; i < matData.rows(); ++i) {
127  rowData = matData.row(i);
128 
129  //FFT for freq domain returning the half spectrum
130  for (j = 0; j < matTaper.rows(); j++) {
131  vecInputFFT = rowData.cwiseProduct(matTaper.row(j));
132  fft.fwd(vecTmpFreq, vecInputFFT, iNfft);
133  matTapSpectrum.row(j) = vecTmpFreq;
134  }
135 
136  finalResult.append(matTapSpectrum);
137 
138 // iTime = timer.elapsed();
139 // qDebug() << QThread::currentThreadId() << "Spectral::computeTaperedSpectraMatrix - Row-wise computation:" << iTime;
140 // iTimeAll += iTime;
141 // timer.restart();
142  }
143 
144 // qDebug() << QThread::currentThreadId() << "Spectral::computeTaperedSpectraMatrix - Complete computation:" << iTimeAll;
145  } else {
146  // Parallel
147  QList<TaperedSpectraInputData> lData;
148  TaperedSpectraInputData dataTemp;
149  dataTemp.matTaper = matTaper;
150  dataTemp.iNfft = iNfft;
151 
152  for (int i = 0; i < matData.rows(); ++i) {
153  dataTemp.vecData = matData.row(i);
154 
155  lData.append(dataTemp);
156  }
157 
158  QFuture<QVector<MatrixXcd> > result = QtConcurrent::mappedReduced(lData,
159  compute,
160  reduce,
161  QtConcurrent::OrderedReduce);
162  result.waitForFinished();
163  finalResult = result.result();
164  }
165 
166  return finalResult;
167 }
168 
169 //=============================================================================================================
170 
171 MatrixXcd Spectral::compute(const TaperedSpectraInputData& inputData)
172 {
173  //qDebug() << "Spectral::compute";
174  return computeTaperedSpectraRow(inputData.vecData,
175  inputData.matTaper,
176  inputData.iNfft);
177 }
178 
179 //=============================================================================================================
180 
181 void Spectral::reduce(QVector<MatrixXcd>& finalData,
182  const MatrixXcd& resultData)
183 {
184  //qDebug() << "Spectral::reduce";
185  finalData.append(resultData);
186 }
187 
188 //=============================================================================================================
189 
190 //void Spectral::reduce(std::vector<MatrixXcd>& finalData,
191 // const MatrixXcd& resultData)
192 //{
193 // finalData.push_back(resultData);
194 //}
195 
196 //=============================================================================================================
197 
198 Eigen::RowVectorXd Spectral::psdFromTaperedSpectra(const Eigen::MatrixXcd &matTapSpectrum,
199  const Eigen::VectorXd &vecTapWeights,
200  int iNfft,
201  double dSampFreq)
202 {
203  //Check inputs
204  if (matTapSpectrum.rows() != vecTapWeights.rows()) {
205  return Eigen::RowVectorXd();
206  }
207 
208  //Compute PSD (average over tapers if necessary)
209  //Normalization via sFreq
210  //multiply by 2 due to half spectrum
211  double denom = vecTapWeights.cwiseAbs2().sum() * dSampFreq;
212  Eigen::RowVectorXd vecPsd = 2.0 * (vecTapWeights.asDiagonal() * matTapSpectrum).cwiseAbs2().colwise().sum() / denom;
213 
214  vecPsd(0) /= 2.0;
215  if (iNfft % 2 == 0){
216  vecPsd.tail(1) /= 2.0;
217  }
218 
219  return vecPsd;
220 }
221 
222 //=============================================================================================================
223 
224 Eigen::RowVectorXcd Spectral::csdFromTaperedSpectra(const Eigen::MatrixXcd &vecTapSpectrumSeed,
225  const Eigen::MatrixXcd &vecTapSpectrumTarget,
226  const Eigen::VectorXd &vecTapWeightsSeed,
227  const Eigen::VectorXd &vecTapWeightsTarget,
228  int iNfft,
229  double dSampFreq)
230 {
231 // QElapsedTimer timer;
232 // int iTime = 0;
233 // timer.start();
234 
235  //Check inputs
236  if (vecTapSpectrumSeed.rows() != vecTapSpectrumTarget.rows()) {
237  return Eigen::MatrixXcd();
238  }
239  if (vecTapSpectrumSeed.cols() != vecTapSpectrumTarget.cols()) {
240  return Eigen::MatrixXcd();
241  }
242  if (vecTapSpectrumSeed.rows() != vecTapWeightsSeed.rows()) {
243  return Eigen::MatrixXcd();
244  }
245  if (vecTapSpectrumTarget.rows() != vecTapWeightsTarget.rows()) {
246  return Eigen::MatrixXcd();
247  }
248 
249 // iTime = timer.elapsed();
250 // qDebug() << QThread::currentThreadId() << "Spectral::csdFromTaperedSpectra timer - Prepare:" << iTime;
251 // timer.restart();
252 
253  // Compute PSD (average over tapers if necessary)
254  // Multiply by 2 due to half spectrum
255  // Normalize via sFreq
256  double denom = sqrt(vecTapWeightsSeed.cwiseAbs2().sum()) * sqrt(vecTapWeightsTarget.cwiseAbs2().sum()) * dSampFreq;
257  Eigen::RowVectorXcd vecCsd = 2.0 * (vecTapWeightsSeed.asDiagonal() * vecTapSpectrumSeed).cwiseProduct((vecTapWeightsTarget.asDiagonal() * vecTapSpectrumTarget).conjugate()).colwise().sum() / denom;
258 
259 // iTime = timer.elapsed();
260 // qDebug() << QThread::currentThreadId() << "Spectral::csdFromTaperedSpectra timer - compute PSD:" << iTime;
261 // timer.restart();
262 
263  //multiply first and last element by 2 due to half spectrum
264  vecCsd(0) /= 2.0;
265  if (iNfft % 2 == 0){
266  vecCsd.tail(1) /= 2.0;
267  }
268 
269 // iTime = timer.elapsed();
270 // qDebug() << QThread::currentThreadId() << "Spectral::csdFromTaperedSpectra timer - half spectrum:" << iTime;
271 // timer.restart();
272 
273  return vecCsd;
274 }
275 
276 //=============================================================================================================
277 
278 VectorXd Spectral::calculateFFTFreqs(int iNfft, double dSampFreq)
279 {
280  //Compute FFT frequencies
281  RowVectorXd vecFFTFreqs;
282  if (iNfft % 2 == 0){
283  vecFFTFreqs = (dSampFreq / iNfft) * RowVectorXd::LinSpaced(iNfft / 2.0 + 1, 0.0, iNfft / 2.0);
284  } else {
285  vecFFTFreqs = (dSampFreq / iNfft) * RowVectorXd::LinSpaced((iNfft - 1) / 2.0 + 1, 0.0, (iNfft - 1) / 2.0);
286  }
287  return vecFFTFreqs;
288 }
289 
290 //=============================================================================================================
291 
292 QPair<MatrixXd, VectorXd> Spectral::generateTapers(int iSignalLength, const QString &sWindowType)
293 {
294  QPair<MatrixXd, VectorXd> pairOut;
295  if (sWindowType == "hanning") {
296  pairOut.first = hanningWindow(iSignalLength);
297  pairOut.second = VectorXd::Ones(1);
298  } else if (sWindowType == "ones") {
299  pairOut.first = MatrixXd::Ones(1, iSignalLength) / double(iSignalLength);
300  pairOut.second = VectorXd::Ones(1);
301  } else {
302  pairOut.first = hanningWindow(iSignalLength);
303  pairOut.second = VectorXd::Ones(1);
304  }
305  return pairOut;
306 }
307 
308 //=============================================================================================================
309 
310 std::pair<MatrixXd, VectorXd> Spectral::generateTapers(int iSignalLength, const std::string &sWindowType)
311 {
312  std::pair<MatrixXd, VectorXd> pairOut;
313  if (sWindowType == "hanning") {
314  pairOut.first = hanningWindow(iSignalLength);
315  pairOut.second = VectorXd::Ones(1);
316  } else if (sWindowType == "ones") {
317  pairOut.first = MatrixXd::Ones(1, iSignalLength) / double(iSignalLength);
318  pairOut.second = VectorXd::Ones(1);
319  } else {
320  pairOut.first = hanningWindow(iSignalLength);
321  pairOut.second = VectorXd::Ones(1);
322  }
323  return pairOut;
324 }
325 
326 //=============================================================================================================
327 
328 MatrixXd Spectral::hanningWindow(int iSignalLength)
329 {
330  MatrixXd matHann = MatrixXd::Zero(1, iSignalLength);
331 
332  //Main step of building the hanning window
333  for (int n = 0; n < iSignalLength; n++) {
334  matHann(0, n) = 0.5 - 0.5 * cos(2.0 * M_PI * n / (iSignalLength - 1.0));
335  }
336  matHann.array() /= matHann.row(0).norm();
337 
338  return matHann;
339 }
static Eigen::RowVectorXd psdFromTaperedSpectra(const Eigen::MatrixXcd &matTapSpectrum, const Eigen::VectorXd &vecTapWeights, int iNfft, double dSampFreq=1.0)
Definition: spectral.cpp:198
static QVector< Eigen::MatrixXcd > computeTaperedSpectraMatrix(const Eigen::MatrixXd &matData, const Eigen::MatrixXd &matTaper, int iNfft, bool bUseThreads=true)
Definition: spectral.cpp:100
static Eigen::RowVectorXcd csdFromTaperedSpectra(const Eigen::MatrixXcd &vecTapSpectrumSeed, const Eigen::MatrixXcd &vecTapSpectrumTarget, const Eigen::VectorXd &vecTapWeightsSeed, const Eigen::VectorXd &vecTapWeightsTarget, int iNfft, double dSampFreq=1.0)
Definition: spectral.cpp:224
static Eigen::MatrixXcd computeTaperedSpectraRow(const Eigen::RowVectorXd &vecData, const Eigen::MatrixXd &matTaper, int iNfft)
Definition: spectral.cpp:72
Declaration of Spectral class.
static Eigen::VectorXd calculateFFTFreqs(int iNfft, double dSampFreq)
Definition: spectral.cpp:278
static void reduce(QVector< Eigen::MatrixXcd > &finalData, const Eigen::MatrixXcd &resultData)
Definition: spectral.cpp:181
static Eigen::MatrixXcd compute(const TaperedSpectraInputData &inputData)
Definition: spectral.cpp:171
static QPair< Eigen::MatrixXd, Eigen::VectorXd > generateTapers(int iSignalLength, const QString &sWindowType="hanning")
Definition: spectral.cpp:292