diff --git a/plugins/processing/signal-processing/src/algorithms/connectivity/connectivityMeasure.cpp b/plugins/processing/signal-processing/src/algorithms/connectivity/connectivityMeasure.cpp
index 3026f7f0324ed13b2e970d857556b1754e838772..7536ec1e4ce330864355eae3c52351faf7d1b04a 100644
--- a/plugins/processing/signal-processing/src/algorithms/connectivity/connectivityMeasure.cpp
+++ b/plugins/processing/signal-processing/src/algorithms/connectivity/connectivityMeasure.cpp
@@ -113,18 +113,19 @@ bool ConnectivityMeasure::initialize(const EConnectMetric metric, uint64_t sampR
m_nbChannels = nbChannels;
m_fftSize = fftSize;
m_dcRemoval = dcRemoval;
+
return true;
}
-bool ConnectivityMeasure::process(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix, const EPsdMode mode)
+bool ConnectivityMeasure::process(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix, const EPsdMode mode, const std::vector<size_t>& vLookup)
{
switch(m_psdMode) {
- case EPsdMode::Welch: return processWelch(samples, connectivityMatrix);
- case EPsdMode::Burg: return processBurg(samples, connectivityMatrix);
+ case EPsdMode::Welch: return processWelch(samples, connectivityMatrix, vLookup);
+ case EPsdMode::Burg: return processBurg(samples, connectivityMatrix, vLookup);
}
}
-bool ConnectivityMeasure::processWelch(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix)
+bool ConnectivityMeasure::processWelch(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup)
{
// Windowing
const size_t windowShift = m_window.size() - m_windowOverlapSamples;
@@ -152,24 +153,26 @@ bool ConnectivityMeasure::processWelch(const std::vector<Eigen::VectorXd>& sampl
}
// Cross-SPECTRA
- for (size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
- for (size_t chan2 = chan1; chan2 < m_nbChannels; chan2++) {
- crossSpectralDensity(m_dft[chan1], m_dft[chan2], m_cpsd[chan1][chan2], nbWindows);
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+ for (size_t chan2 = 0; chan2 < m_nbChannels; chan2++) {
+ if (vLookup[chan1] != chan2) {
+ crossSpectralDensity(m_dft[vLookup[chan1]], m_dft[chan2], m_cpsd[vLookup[chan1]][chan2], nbWindows);
+ }
}
}
// Use PSDs & x-spectra to compute the connectivity matrix
std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>> placeHolder; // empty vector, to use same fcts...
switch (m_metric) {
- case EConnectMetric::Coherence: return coherence(placeHolder, connectivityMatrix);
- case EConnectMetric::MagnitudeSquaredCoherence: return magnitudeSquaredCoherence(placeHolder, connectivityMatrix);
- case EConnectMetric::ImaginaryCoherence: return imaginaryCoherence(placeHolder, connectivityMatrix);
- case EConnectMetric::AbsImaginaryCoherence: return absImaginaryCoherence(placeHolder, connectivityMatrix);
- default: return coherence(placeHolder, connectivityMatrix);
+ case EConnectMetric::Coherence: return coherence(placeHolder, connectivityMatrix, vLookup);
+ case EConnectMetric::MagnitudeSquaredCoherence: return magnitudeSquaredCoherence(placeHolder, connectivityMatrix, vLookup);
+ case EConnectMetric::ImaginaryCoherence: return imaginaryCoherence(placeHolder, connectivityMatrix, vLookup);
+ case EConnectMetric::AbsImaginaryCoherence: return absImaginaryCoherence(placeHolder, connectivityMatrix, vLookup);
+ default: return coherence(placeHolder, connectivityMatrix, vLookup);
}
}
-bool ConnectivityMeasure::processBurg(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix)
+bool ConnectivityMeasure::processBurg(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup)
{
// NOT MULTITHREADED VERSION - for reference / debugging purpose
/*for (size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
@@ -181,25 +184,38 @@ bool ConnectivityMeasure::processBurg(const std::vector<Eigen::VectorXd>& sample
}*/
// THREAD POOL (MAX THREADS / STANDALONE without init)
- for (size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
- for (size_t chan2 = chan1 + 1; chan2 < m_nbChannels; chan2++) {
- m_threadPoolMax.QueueJob([&, chan1, chan2] {
- // No need to use mutexes : samples data can be read from multiple threads without issue.
- // and writing occurs to indices of m_mvarCoeffsMatrix that are independent for every threads
- mvarSpectralEstimation(samples[chan1], samples[chan2], m_mvarCoeffsMatrix[chan1][chan2], m_autoRegOrder, m_fftSize, m_sampRate, m_dcRemoval, m_expMat);
- });
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+ size_t chan1idx = vLookup[chan1];
+ size_t startingIdxChan2 = 0;
+
+ // if we consider all channels, override this channel idx nonsense
+ if (vLookup.size() == m_nbChannels) {
+ startingIdxChan2 = chan1 + 1;
+ chan1idx = chan1;
+ }
+
+ for (size_t chan2 = startingIdxChan2; chan2 < m_nbChannels; chan2++) {
+
+ if (chan1idx != chan2) {
+ m_threadPoolMax.QueueJob([&, chan1idx, chan2] {
+ // No need to use mutexes : samples data can be read from multiple threads without issue.
+ // and writing occurs to indices of m_mvarCoeffsMatrix that are independent for every threads
+ mvarSpectralEstimation(samples[chan1idx], samples[chan2], m_mvarCoeffsMatrix[chan1idx][chan2], m_autoRegOrder, m_fftSize, m_sampRate, m_dcRemoval, m_expMat);
+ });
+ }
}
+
}
// Wait until all jobs have been processed
m_threadPoolMax.waitUntilCompleted();
// Use MVAR results to compute the connectivity matrix
switch (m_metric) {
- case EConnectMetric::Coherence: return coherence(m_mvarCoeffsMatrix, connectivityMatrix);
- case EConnectMetric::MagnitudeSquaredCoherence: return magnitudeSquaredCoherence(m_mvarCoeffsMatrix, connectivityMatrix);
- case EConnectMetric::ImaginaryCoherence: return imaginaryCoherence(m_mvarCoeffsMatrix, connectivityMatrix);
- case EConnectMetric::AbsImaginaryCoherence: return absImaginaryCoherence(m_mvarCoeffsMatrix, connectivityMatrix);
- default: return coherence(m_mvarCoeffsMatrix, connectivityMatrix);
+ case EConnectMetric::Coherence: return coherence(m_mvarCoeffsMatrix, connectivityMatrix, vLookup);
+ case EConnectMetric::MagnitudeSquaredCoherence: return magnitudeSquaredCoherence(m_mvarCoeffsMatrix, connectivityMatrix, vLookup);
+ case EConnectMetric::ImaginaryCoherence: return imaginaryCoherence(m_mvarCoeffsMatrix, connectivityMatrix, vLookup);
+ case EConnectMetric::AbsImaginaryCoherence: return absImaginaryCoherence(m_mvarCoeffsMatrix, connectivityMatrix, vLookup);
+ default: return coherence(m_mvarCoeffsMatrix, connectivityMatrix, vLookup);
}
}
@@ -251,55 +267,90 @@ bool ConnectivityMeasure::crossSpectralDensity(const Eigen::MatrixXcd& dft1, con
return true;
}
-bool ConnectivityMeasure::coherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix)
+bool ConnectivityMeasure::coherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup)
{
// TODO : coherence is a complex matrix. Here we also used MSC... but with sqrt in the denominator.
if(m_psdMode == EPsdMode::Welch) {
// PSDs and CPSD are directly available as class members m_psd and m_cpsd
- for (size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
- connectivityMatrix[chan1].row(chan1) = Eigen::VectorXd::Zero(m_psd[0].size());
- for (size_t chan2 = chan1 + 1; chan2 < m_nbChannels; chan2++) {
- connectivityMatrix[chan1].row(chan2) = m_cpsd[chan1][chan2].cwiseAbs2().cwiseQuotient(m_psd[chan1].cwiseProduct(m_psd[chan2])).cwiseSqrt();
- connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2);
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+ size_t chan1idx = vLookup[chan1];
+ size_t startingIdxChan2 = 0;
+ if (vLookup.size() == m_nbChannels) {
+ chan1idx = chan1;
+ startingIdxChan2 = chan1 + 1;
}
- }
- } else { // m_psdMode == EPsdMode::Burg
+ connectivityMatrix[chan1idx].row(chan1) = Eigen::VectorXd::Zero(m_psd[0].size());
- // For each channel pair : Coh = S(0,1) / sqrt(S(0,0).*S(1,1))
- // Input matrix dimensions : chan1 x chan2 x nfft x Matrix2cd (S11,S12),(S21,S22))
- // mvarMatrix[chan1][chan2] is matrix S
- // mvarMatrix[chan1][chan2][0][0] is psd of chan1
- // mvarMatrix[chan1][chan2][1][1] is psd of chan2
- // mvarMatrix[chan1][chan2][0][1] is cpsd of pair of channels
+ for (size_t chan2 = startingIdxChan2; chan2 < m_nbChannels; chan2++) {
+ if (chan1idx != chan2) {
+ connectivityMatrix[chan1].row(chan2) = m_cpsd[chan1idx][chan2].cwiseAbs2().cwiseQuotient(m_psd[chan1idx].cwiseProduct(m_psd[chan2])).cwiseSqrt();
+ if (vLookup.size() == m_nbChannels) {
+ connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2); // diagonal matrix
+ }
+ }
+ }
+ }
+
+ }
+ else { // m_psdMode == EPsdMode::Burg
+
+ // For each channel pair : Coh = S(0,1) / sqrt(S(0,0).*S(1,1))
+ // Input matrix dimensions : chan1 x chan2 x nfft x Matrix2cd (S11,S12),(S21,S22))
+ // mvarMatrix[chan1][chan2] is matrix S
+ // mvarMatrix[chan1][chan2][0][0] is psd of chan1
+ // mvarMatrix[chan1][chan2][1][1] is psd of chan2
+ // mvarMatrix[chan1][chan2][0][1] is cpsd of pair of channels
+
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+
+ size_t chan1idx = vLookup[chan1];
+ size_t startingIdxChan2 = 0;
+ if (vLookup.size() == m_nbChannels) {
+ chan1idx = chan1;
+ startingIdxChan2 = chan1 + 1;
+ }
- for(size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
connectivityMatrix[chan1].row(chan1) = Eigen::VectorXd::Zero(m_fftSize);
- for (size_t chan2 = chan1 + 1; chan2 < m_nbChannels; chan2++)
- {
- connectivityMatrix[chan1].row(chan2) = mvarMatrix[chan1][chan2][0][1].cwiseAbs2().cwiseQuotient(mvarMatrix[chan1][chan2][0][0].real().cwiseProduct(mvarMatrix[chan1][chan2][1][1].real())).cwiseSqrt();
- connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2);
+
+ for (size_t chan2 = startingIdxChan2; chan2 < m_nbChannels; chan2++) {
+ {
+ if (chan1idx != chan2) {
+ connectivityMatrix[chan1].row(chan2) = mvarMatrix[chan1idx][chan2][0][1].cwiseAbs2().cwiseQuotient(mvarMatrix[chan1idx][chan2][0][0].real().cwiseProduct(mvarMatrix[chan1idx][chan2][1][1].real())).cwiseSqrt();
+ if (vLookup.size() == m_nbChannels) {
+ connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2); // diagonal matrix
+ }
+ }
+ }
}
}
-
}
return true;
}
-
-bool ConnectivityMeasure::magnitudeSquaredCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix)
+bool ConnectivityMeasure::magnitudeSquaredCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup)
{
if(m_psdMode == EPsdMode::Welch) {
// PSDs and CPSD are directly available as class members m_psd and m_cpsd
- for (size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+ size_t chan1idx = vLookup[chan1];
+ size_t startingIdxChan2 = 0;
+ if (vLookup.size() == m_nbChannels) {
+ chan1idx = chan1;
+ startingIdxChan2 = chan1 + 1;
+ }
connectivityMatrix[chan1].row(chan1) = Eigen::VectorXd::Zero(m_psd[0].size());
- for (size_t chan2 = chan1 + 1; chan2 < m_nbChannels; chan2++) {
- connectivityMatrix[chan1].row(chan2) = m_cpsd[chan1][chan2].cwiseAbs2().cwiseQuotient(m_psd[chan1].cwiseProduct(m_psd[chan2]) );
- connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2);
+ for (size_t chan2 = startingIdxChan2; chan2 < m_nbChannels; chan2++) {
+ if (chan1 != chan2) {
+ connectivityMatrix[chan1].row(chan2) = m_cpsd[chan1idx][chan2].cwiseAbs2().cwiseQuotient(m_psd[chan1idx].cwiseProduct(m_psd[chan2]));
+ if (vLookup.size() == m_nbChannels) {
+ connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2); // diagonal matrix
+ }
+ }
}
}
@@ -311,11 +362,22 @@ bool ConnectivityMeasure::magnitudeSquaredCoherence(const std::vector<std::vecto
// mvarMatrix[chan1][chan2][0][0] is psd of chan1
// mvarMatrix[chan1][chan2][1][1] is psd of chan2
// mvarMatrix[chan1][chan2][0][1] is cpsd of pair of channels
- for (size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+ size_t chan1idx = vLookup[chan1];
+ size_t startingIdxChan2 = 0;
+ if (vLookup.size() == m_nbChannels) {
+ chan1idx = chan1;
+ startingIdxChan2 = chan1 + 1;
+ }
+
connectivityMatrix[chan1].row(chan1) = Eigen::VectorXd::Zero(m_fftSize);
- for (size_t chan2 = chan1 + 1; chan2 < m_nbChannels; chan2++) {
- connectivityMatrix[chan1].row(chan2) = mvarMatrix[chan1][chan2][0][1].cwiseAbs2().cwiseQuotient(mvarMatrix[chan1][chan2][0][0].real().cwiseProduct(mvarMatrix[chan1][chan2][1][1].real()));
- connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2);
+ for (size_t chan2 = startingIdxChan2; chan2 < m_nbChannels; chan2++) {
+ if (chan1idx != chan2) {
+ connectivityMatrix[chan1].row(chan2) = mvarMatrix[chan1idx][chan2][0][1].cwiseAbs2().cwiseQuotient(mvarMatrix[chan1idx][chan2][0][0].real().cwiseProduct(mvarMatrix[chan1idx][chan2][1][1].real()));
+ if (vLookup.size() == m_nbChannels) {
+ connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2); // diagonal matrix
+ }
+ }
}
}
@@ -325,16 +387,28 @@ bool ConnectivityMeasure::magnitudeSquaredCoherence(const std::vector<std::vecto
return true;
}
-bool ConnectivityMeasure::imaginaryCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix)
+bool ConnectivityMeasure::imaginaryCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup)
{
if(m_psdMode == EPsdMode::Welch) {
// PSDs and CPSD are directly available as class members m_psd and m_cpsd
- for (size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
- connectivityMatrix[chan1].row(chan1) = Eigen::VectorXd::Zero(m_psd[0].size());
- for (size_t chan2 = chan1 + 1; chan2 < m_nbChannels; chan2++) {
- connectivityMatrix[chan1].row(chan2) = m_cpsd[chan1][chan2].imag().cwiseQuotient(m_psd[chan1].cwiseProduct(m_psd[chan2]).cwiseSqrt());
- connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2);
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+ size_t chan1idx = vLookup[chan1];
+ size_t startingIdxChan2 = 0;
+ if (vLookup.size() == m_nbChannels) {
+ chan1idx = chan1;
+ startingIdxChan2 = chan1 + 1;
+ }
+
+ connectivityMatrix[chan1idx].row(chan1) = Eigen::VectorXd::Zero(m_psd[0].size());
+
+ for (size_t chan2 = startingIdxChan2; chan2 < m_nbChannels; chan2++) {
+ if (chan1idx != chan2) {
+ connectivityMatrix[chan1].row(chan2) = m_cpsd[chan1idx][chan2].imag().cwiseQuotient(m_psd[chan1idx].cwiseProduct(m_psd[chan2]).cwiseSqrt());
+ if (vLookup.size() == m_nbChannels) {
+ connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2); // diagonal matrix
+ }
+ }
}
}
@@ -347,11 +421,22 @@ bool ConnectivityMeasure::imaginaryCoherence(const std::vector<std::vector<std::
// mvarMatrix[chan1][chan2][1][1] is psd of chan2
// mvarMatrix[chan1][chan2][0][1] is cpsd of pair of channels
- for(size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+ size_t chan1idx = vLookup[chan1];
+ size_t startingIdxChan2 = 0;
+ if (vLookup.size() == m_nbChannels) {
+ chan1idx = chan1;
+ startingIdxChan2 = chan1 + 1;
+ }
+
connectivityMatrix[chan1].row(chan1) = Eigen::VectorXd::Zero(m_fftSize);
- for (size_t chan2 = chan1 + 1; chan2 < m_nbChannels; chan2++) {
- connectivityMatrix[chan1].row(chan2) = mvarMatrix[chan1][chan2][0][1].imag().cwiseQuotient(mvarMatrix[chan1][chan2][0][0].real().cwiseProduct(mvarMatrix[chan1][chan2][1][1].real()).cwiseSqrt());
- connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2);
+ for (size_t chan2 = startingIdxChan2; chan2 < m_nbChannels; chan2++) {
+ if (chan1idx != chan2) {
+ connectivityMatrix[chan1].row(chan2) = mvarMatrix[chan1idx][chan2][0][1].imag().cwiseQuotient(mvarMatrix[chan1idx][chan2][0][0].real().cwiseProduct(mvarMatrix[chan1idx][chan2][1][1].real()).cwiseSqrt());
+ if (vLookup.size() == m_nbChannels) {
+ connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2); // diagonal matrix
+ }
+ }
}
}
}
@@ -359,16 +444,28 @@ bool ConnectivityMeasure::imaginaryCoherence(const std::vector<std::vector<std::
return true;
}
-bool ConnectivityMeasure::absImaginaryCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix)
+bool ConnectivityMeasure::absImaginaryCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup)
{
if(m_psdMode == EPsdMode::Welch) {
// PSDs and CPSD are directly available as class members m_psd and m_cpsd
- for (size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
- connectivityMatrix[chan1].row(chan1) = Eigen::VectorXd::Zero(m_psd[0].size());
- for (size_t chan2 = chan1 + 1; chan2 < m_nbChannels; chan2++) {
- connectivityMatrix[chan1].row(chan2) = m_cpsd[chan1][chan2].imag().cwiseAbs().cwiseQuotient(m_psd[chan1].cwiseProduct(m_psd[chan2]).cwiseSqrt());
- connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2);
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+ size_t chan1idx = vLookup[chan1];
+ size_t startingIdxChan2 = 0;
+ if (vLookup.size() == m_nbChannels) {
+ chan1idx = chan1;
+ startingIdxChan2 = chan1 + 1;
+ }
+
+ connectivityMatrix[chan1idx].row(chan1) = Eigen::VectorXd::Zero(m_psd[0].size());
+
+ for (size_t chan2 = startingIdxChan2; chan2 < m_nbChannels; chan2++) {
+ if (chan1idx != chan2) {
+ connectivityMatrix[chan1].row(chan2) = m_cpsd[chan1idx][chan2].imag().cwiseAbs().cwiseQuotient(m_psd[chan1idx].cwiseProduct(m_psd[chan2]).cwiseSqrt());
+ if (vLookup.size() == m_nbChannels) {
+ connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2); // diagonal matrix
+ }
+ }
}
}
@@ -381,11 +478,22 @@ bool ConnectivityMeasure::absImaginaryCoherence(const std::vector<std::vector<st
// mvarMatrix[chan1][chan2][1][1] is psd of chan2
// mvarMatrix[chan1][chan2][0][1] is cpsd of pair of channels
- for(size_t chan1 = 0; chan1 < m_nbChannels; chan1++) {
+ for (size_t chan1 = 0; chan1 < vLookup.size(); chan1++) {
+ size_t chan1idx = vLookup[chan1];
+ size_t startingIdxChan2 = 0;
+ if (vLookup.size() == m_nbChannels) {
+ chan1idx = chan1;
+ startingIdxChan2 = chan1 + 1;
+ }
+
connectivityMatrix[chan1].row(chan1) = Eigen::VectorXd::Zero(m_fftSize);
- for (size_t chan2 = chan1 + 1; chan2 < m_nbChannels; chan2++) {
- connectivityMatrix[chan1].row(chan2) = mvarMatrix[chan1][chan2][0][1].imag().cwiseAbs().cwiseQuotient(mvarMatrix[chan1][chan2][0][0].real().cwiseProduct(mvarMatrix[chan1][chan2][1][1].real()).cwiseSqrt());
- connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2);
+ for (size_t chan2 = startingIdxChan2; chan2 < m_nbChannels; chan2++) {
+ if (chan1idx != chan2) {
+ connectivityMatrix[chan1].row(chan2) = mvarMatrix[chan1idx][chan2][0][1].imag().cwiseAbs().cwiseQuotient(mvarMatrix[chan1idx][chan2][0][0].real().cwiseProduct(mvarMatrix[chan1idx][chan2][1][1].real()).cwiseSqrt());
+ if (vLookup.size() == m_nbChannels) {
+ connectivityMatrix[chan2].row(chan1) = connectivityMatrix[chan1].row(chan2); // diagonal matrix
+ }
+ }
}
}
}
diff --git a/plugins/processing/signal-processing/src/algorithms/connectivity/connectivityMeasure.hpp b/plugins/processing/signal-processing/src/algorithms/connectivity/connectivityMeasure.hpp
index 14f9a23ab18600367aba8326441fa8484dbc81a2..ce709b628cfdad6bff66f5b55522c955204bbd50 100644
--- a/plugins/processing/signal-processing/src/algorithms/connectivity/connectivityMeasure.hpp
+++ b/plugins/processing/signal-processing/src/algorithms/connectivity/connectivityMeasure.hpp
@@ -108,7 +108,9 @@ inline std::string toString(const EConnectWindowMethod method)
class ConnectivityMeasure
{
public:
- ConnectivityMeasure(){m_threadPoolMax.Start( std::thread::hardware_concurrency());}
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
+ ConnectivityMeasure() {m_threadPoolMax.Start(std::thread::hardware_concurrency()); }
~ConnectivityMeasure(){m_threadPoolMax.Stop();}
/// \brief Initialize the connectivity measurement class, Welch specific parameters
@@ -136,11 +138,11 @@ public:
/// \brief Select the function to call for the connectivity measurement.
/// \param samples The input data set \f$x\f$. With \f$ nChan \f$ channels and \f$ nSamp \f$ samples per channel
/// \param connectivityMatrix The connectivity matrix
- bool process(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix, const EPsdMode mode);
+ bool process(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix, const EPsdMode mode, const std::vector<size_t>& vLookup);
protected:
- bool processWelch(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix);
- bool processBurg(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix);
+ bool processWelch(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup);
+ bool processBurg(const std::vector<Eigen::VectorXd>& samples, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup);
private:
/// \brief Initialize the connectivity measurement class, with common parameters
@@ -166,7 +168,7 @@ private:
/// \param connectivityMatrix The connectivity matrix as a vector of 2D Matrices of size nbChannels x (nbChannels x frequency_taps)
/// \return True if it succeeds, false if it fails.
/// \todo don't recompute the whole table, because connectivity measures overlap !
- bool coherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix);
+ bool coherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup);
/// \brief Calculation of the Magnitude Squared Coherence
/// \f$ Pxx = PSD of x \f$ // \f$ Pyy = PSD of y \f$ // \f$ Sxy = cross spectral density of x and y \f$
@@ -183,7 +185,7 @@ private:
/// \param connectivityMatrix The connectivity matrix as a vector of 2D Matrices of size nbChannels x (nbChannels x frequency_taps)
/// \return True if it succeeds, false if it fails.
/// \todo don't recompute the whole table, because connectivity measures overlap !
- bool magnitudeSquaredCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix);
+ bool magnitudeSquaredCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup);
/// \brief Calculation of the Imaginary part of the coherence
/// \f$ Pxx = PSD of x \f$ // \f$ Pyy = PSD of y \f$ // \f$ Sxy = cross spectral density of x and y \f$
@@ -199,7 +201,7 @@ private:
/// \param connectivityMatrix The connectivity matrix as a vector of 2D Matrices of size nbChannels x (nbChannels x frequency_taps)
/// \return True if it succeeds, false if it fails.
/// \todo don't recompute the whole table, because connectivity measures overlap !
- bool imaginaryCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix);
+ bool imaginaryCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup);
/// \brief Calculation of the absolute value of the Imaginary part of the coherence
/// \f$ Pxx = PSD of x \f$ // \f$ Pyy = PSD of y \f$ // \f$ Sxy = cross spectral density of x and y \f$
@@ -215,7 +217,7 @@ private:
/// \param connectivityMatrix The connectivity matrix as a vector of 2D Matrices of size nbChannels x (nbChannels x frequency_taps)
/// \return True if it succeeds, false if it fails.
/// \todo don't recompute the whole table, because connectivity measures overlap !
- bool absImaginaryCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix);
+ bool absImaginaryCoherence(const std::vector<std::vector<std::vector<std::vector<Eigen::VectorXcd>>>>& mvarMatrix, std::vector<Eigen::MatrixXd>& connectivityMatrix, const std::vector<size_t>& vLookup);
/// \brief Generates periodigram of signal
/// \param input The signal
diff --git a/plugins/processing/signal-processing/src/box-algorithms/connectivity/CBoxAlgorithmConnectivityMeasure.cpp b/plugins/processing/signal-processing/src/box-algorithms/connectivity/CBoxAlgorithmConnectivityMeasure.cpp
index f1e863187b5ed27be3f53ba9ea63a7364131ed5b..1b2ce73df2a17fbadabc70d148f0da6bfe5f1c1d 100644
--- a/plugins/processing/signal-processing/src/box-algorithms/connectivity/CBoxAlgorithmConnectivityMeasure.cpp
+++ b/plugins/processing/signal-processing/src/box-algorithms/connectivity/CBoxAlgorithmConnectivityMeasure.cpp
@@ -25,6 +25,7 @@
#include "CBoxAlgorithmConnectivityMeasure.hpp"
#include "eigen/convert.hpp"
#include <cmath>
+#include <chrono>
namespace OpenViBE {
namespace Plugins {
@@ -39,26 +40,28 @@ bool CBoxAlgorithmConnectivityMeasure::initialize()
m_oMatrix = m_matrixEncoder.getInputMatrix();
// Settings
- m_metric = EConnectMetric(uint64_t(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 0)));
- m_connectLengthSeconds = double(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 1));
- m_connectOverlap = int(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 2));
- m_dcRemoval = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 3);
- m_fftSize = size_t(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 4));
- m_psdMode = EPsdMode(uint64_t(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 5)));
-
+ m_selectedChannels = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 0);
+ m_metric = EConnectMetric(uint64_t(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 1)));
+ m_connectLengthSeconds = double(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 2));
+ m_connectOverlap = int(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 3));
+ m_dcRemoval = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 4);
+ m_fftSize = size_t(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 5));
+ m_psdMode = EPsdMode(uint64_t(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 6)));
+
OV_ERROR_UNLESS_KRF(m_psdMode == EPsdMode::Welch || m_psdMode == EPsdMode::Burg, "Invalid PSD mode", Kernel::ErrorType::BadSetting);
switch(m_psdMode) {
case EPsdMode::Welch:
- m_windowMethod = EConnectWindowMethod(uint64_t(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 6)));
- m_windowLengthSeconds = double(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 7));
- m_windowOverlap = int(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 8));
+ m_windowMethod = EConnectWindowMethod(uint64_t(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 7)));
+ m_windowLengthSeconds = double(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 8));
+ m_windowOverlap = int(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 9));
break;
case EPsdMode::Burg:
- m_autoRegOrder = int(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 6));
+ m_autoRegOrder = int(FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 7));
break;
}
+
return true;
}
@@ -101,13 +104,54 @@ bool CBoxAlgorithmConnectivityMeasure::process()
m_windowOverlapSamples = int(std::floor(double(m_windowLength * m_windowOverlap / 100.0)));
m_connectShiftSamples = m_connectLength - m_connectOverlapSamples;
+ if (m_selectedChannels == CString("")) {
+ // Keep all channels
+ for (size_t j = 0; j < m_iMatrix->getDimensionSize(0); ++j) {
+ m_vLookup.push_back(j);
+ }
+ } else {
+ std::vector<CString> token;
+ const size_t nToken = split(m_selectedChannels, Toolkit::String::TSplitCallback<std::vector<CString>>(token), OV_Value_EnumeratedStringSeparator);
+
+ for (size_t j = 0; j < m_iMatrix->getDimensionSize(0); ++j) {
+ for (size_t k = 0; k < nToken; ++k) {
+ if (Toolkit::String::isAlmostEqual(m_iMatrix->getDimensionLabel(0, j), token[k], false)) {m_vLookup.push_back(j); }
+ }
+ }
+
+ if (nToken != m_vLookup.size()) {
+ this->getLogManager() << Kernel::LogLevel_Warning << "Some channels in the subselection were not found in the signal stream.\n";
+ for (size_t k = 0; k < nToken; ++k) {
+ bool found = false;
+ for (size_t j = 0; j < m_vLookup.size(); j++) {
+ if (Toolkit::String::isAlmostEqual(m_iMatrix->getDimensionLabel(0, m_vLookup[j]), token[k], false)) {
+ found = true; break;
+ }
+ }
+
+ if (!found) {
+ this->getLogManager() << Kernel::LogLevel_Warning << token[k] << " was not found in the signal stream.\n";
+ }
+ }
+ }
+
+ }
+
+ m_nbChanOutput = m_vLookup.size();
+ OV_ERROR_UNLESS_KRF(m_nbChanOutput > 0, "Connectivity channel subselection error. Mismatch between the subselection and the channel names? ",
+ Kernel::ErrorType::BadArgument);
+
+
+
+ this->getLogManager() << Kernel::LogLevel_Info << "NB CHANNELS FOR CONNECTIVITY OUTPUT: " << m_nbChanOutput << ")\n";
+
m_timeIncrement = uint64_t(m_connectLengthSeconds * std::pow(2, 32));
this->getLogManager() << Kernel::LogLevel_Debug << "TIME INCREMENT: " << m_timeIncrement << " (" << (double)(m_timeIncrement)*std::pow(2,-32) << ")\n";
m_nbChannels = matrix->getDimensionSize(0);
const auto sampPerChan = matrix->getDimensionSize(1);
- this->getLogManager() << Kernel::LogLevel_Debug << "HEADER : nChannels " << m_nbChannels << ", " << sampPerChan << " samples, sampling rate " << sampling << "\n";
+ this->getLogManager() << Kernel::LogLevel_Info << "HEADER : nChannels " << m_nbChannels << ", " << sampPerChan << " samples, sampling rate " << sampling << "\n";
this->getLogManager() << Kernel::LogLevel_Debug << "Connectivity segments overlap percent/samples : " << m_connectOverlap << " " << m_connectOverlapSamples << "\n";
this->getLogManager() << Kernel::LogLevel_Debug << "Connectivity segment shift (samples): " << m_connectShiftSamples << "\n";
this->getLogManager() << Kernel::LogLevel_Debug << "Welch Windows overlap percent/samples : " << m_windowOverlap << " " << m_windowOverlapSamples << "\n";
@@ -127,15 +171,18 @@ bool CBoxAlgorithmConnectivityMeasure::process()
}
- m_oMatrix->resize({ m_fftSize, m_nbChannels, m_nbChannels }); // nbFreqs x nbChan x nbChan
+ m_oMatrix->resize({ m_fftSize, m_nbChannels, m_nbChanOutput }); // nbFreqs x nbChan x nbChanOutput
this->getLogManager() << Kernel::LogLevel_Debug << "output Matrix dimensions: " << m_oMatrix->getDimensionSize(0) << " x " << m_oMatrix->getDimensionSize(1) << "\n";
for (size_t aa = 0; aa < m_fftSize; aa++) {
m_oMatrix->setDimensionLabel(0, aa, std::to_string(aa));
}
for (size_t aa = 0; aa < m_nbChannels; aa++) {
m_oMatrix->setDimensionLabel(1, aa, matrix->getDimensionLabel(0, aa));
- m_oMatrix->setDimensionLabel(2, aa, matrix->getDimensionLabel(0, aa));
- this->getLogManager() << Kernel::LogLevel_Debug << "Channel " << aa << " : " << m_oMatrix->getDimensionLabel(1, aa) << "/" << m_oMatrix->getDimensionLabel(2, aa) << "\n";
+ //this->getLogManager() << Kernel::LogLevel_Debug << "Channel " << aa << " : " << m_oMatrix->getDimensionLabel(1, aa) << "/" << m_oMatrix->getDimensionLabel(2, aa) << "\n";
+ }
+ for (size_t aa = 0; aa < m_nbChanOutput; aa++) {
+ m_oMatrix->setDimensionLabel(2, aa, matrix->getDimensionLabel(0, m_vLookup[aa]));
+ this->getLogManager() << Kernel::LogLevel_Info << "Output matrix: Channel " << aa << " : " << m_oMatrix->getDimensionLabel(2, aa) << "\n";
}
m_matrixEncoder.encodeHeader(); // Pass the header to the next boxes
@@ -185,16 +232,16 @@ bool CBoxAlgorithmConnectivityMeasure::process()
m_signalChannelBuffers[aa].erase(m_signalChannelBuffers[aa].begin(), m_signalChannelBuffers[aa].begin() + m_connectShiftSamples);
}
- // 3D Matrix init, vector of size (chan) of Matrices (chan x fftsize)
- std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> connectivityMatrix(m_nbChannels, Eigen::MatrixXd::Zero(m_nbChannels, m_fftSize));
+ // 3D Matrix init, vector of size (m_nbChanOutput) of Matrices (chan x fftsize)
+ std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> connectivityMatrix(m_nbChanOutput, Eigen::MatrixXd::Zero(m_nbChannels, m_fftSize));
- OV_ERROR_UNLESS_KRF(connectivityMeasure.process(m_vectorXdBuffer, connectivityMatrix, m_psdMode), "Connectivity measurement error",
+ OV_ERROR_UNLESS_KRF(connectivityMeasure.process(m_vectorXdBuffer, connectivityMatrix, m_psdMode, m_vLookup), "Connectivity measurement error",
Kernel::ErrorType::BadProcessing);
-
+
this->getLogManager() << Kernel::LogLevel_Debug << "Exited connectivityMeasure.process() : connect size "
<< (size_t)connectivityMatrix.size() << " x " << (size_t)connectivityMatrix[0].rows() << " x " << (size_t)connectivityMatrix[0].cols() << "\n";
- // Convert output matrix (nchan x nchan x fftsize) to matrix (fftsize x nchan x nchan)
+ // Convert output matrix (nchanoutput x nchan x fftsize) to matrix (fftsize x nchan x nchanoutput)
MatrixConvert(connectivityMatrix, *m_oMatrix, 5);
m_matrixEncoder.encodeBuffer();
diff --git a/plugins/processing/signal-processing/src/box-algorithms/connectivity/CBoxAlgorithmConnectivityMeasure.hpp b/plugins/processing/signal-processing/src/box-algorithms/connectivity/CBoxAlgorithmConnectivityMeasure.hpp
index 2896ec199db5f67cde559fe718d8961539ef4b67..d33f094aada99dcfb597e53b12f2afd4296fa18d 100644
--- a/plugins/processing/signal-processing/src/box-algorithms/connectivity/CBoxAlgorithmConnectivityMeasure.hpp
+++ b/plugins/processing/signal-processing/src/box-algorithms/connectivity/CBoxAlgorithmConnectivityMeasure.hpp
@@ -31,7 +31,6 @@
#include <toolkit/ovtk_all.h>
#include "connectivityMeasure.hpp"
-
#define OV_AttributeId_Box_FlagIsUnstable OpenViBE::CIdentifier(0x666FFFFF, 0x666FFFFF)
namespace OpenViBE {
@@ -44,6 +43,9 @@ namespace SignalProcessing {
class CBoxAlgorithmConnectivityMeasure final : virtual public Toolkit::TBoxAlgorithm<IBoxAlgorithm>
{
public:
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
void release() override { delete this; }
bool initialize() override;
@@ -87,6 +89,11 @@ protected:
int m_autoRegOrder = 11; // auto-regressive model order, in samples
// Common
+ std::vector<size_t> m_vLookup;
+ size_t m_nbChanOutput = 0;
+ //uint64_t totalDuration = 0; //for measuring performance
+
+ CString m_selectedChannels; // subselection of channels, separated with ";". Blank = select all
int m_connectLength = 256; // size of one full connectivity estimation occurrence (samples)
double m_connectLengthSeconds = 0.5; // size of one full connectivity estimation occurrence (sec)
int m_connectOverlap = 50; // overlap btw connectivity measurements (%)
@@ -99,7 +106,7 @@ protected:
uint64_t m_timeIncrement = 0;
std::vector<uint64_t> m_startTimes;
-
+
};
///
@@ -187,6 +194,8 @@ public:
prototype.addInput("Input signal", OV_TypeId_Signal);
prototype.addOutput("Connectivity Matrix", OV_TypeId_StreamedMatrix);
+ prototype.addSetting("Channel subselection", OV_TypeId_String, "");
+
prototype.addSetting("Metric", OVP_TypeId_Connectivity_Metric, toString(EConnectMetric::Coherence).c_str());
prototype.addSetting("Connectivity Measure Length (in sec)", OV_TypeId_Float, "0.5"); //s