csd.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "csd.h"
#include "dpss.h"
#include "morlet_tfr.h"
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Core>
#include <Eigen/Dense>
//#ifndef EIGEN_FFTW_DEFAULT
//#define EIGEN_FFTW_DEFAULT
//#endif
#include <unsupported/Eigen/FFT>
 
//=============================================================================================================
// STD INCLUDES
//=============================================================================================================
 
#include <cmath>
#include <complex>
#include <vector>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace UTILSLIB;
using namespace Eigen;
 
namespace
{
constexpr double CSD_PI = 3.14159265358979323846;
}
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
CsdResult Csd::computeMultitaper(const MatrixXd& matData,
                                 double sfreq,
                                 double fmin,
                                 double fmax,
                                 double halfBandwidth,
                                 int nTapers)
{
    const int nChannels = static_cast<int>(matData.rows());
    const int nTimes = static_cast<int>(matData.cols());
    const int nFreqsFull = nTimes / 2 + 1;
 
    if (fmax < 0.0)
        fmax = sfreq / 2.0;
 
    // 1. Compute DPSS tapers; by default those mne.time_frequency.csd_array_multitaper uses
    const bool lowBias = nTapers < 0;
    DpssResult dpss = Dpss::compute(nTimes, halfBandwidth, lowBias ? static_cast<int>(2.0 * halfBandwidth) : nTapers, false, lowBias);
    const int nTap = static_cast<int>(dpss.matTapers.rows());
 
    // Compute eigenvalue weight sum
    double weightSum = 0.0;
    for (int t = 0; t < nTap; ++t)
        weightSum += dpss.vecEigenvalues[t];
 
    // Build full frequency axis
    const double freqRes = sfreq / static_cast<double>(nTimes);
 
    // Determine frequency bin range (DC excluded, as in MNE-Python)
    const int iBinMin = std::max(1, static_cast<int>(std::ceil(fmin / freqRes)));
    const int iBinMax = std::min(nFreqsFull - 1, static_cast<int>(std::floor(fmax / freqRes)));
    const int nFreqsSel = std::max(0, iBinMax - iBinMin + 1);
 
    // Initialize per-frequency CSD accumulators
    std::vector<MatrixXcd> csdAccum(static_cast<std::size_t>(nFreqsSel),
                                    MatrixXcd::Zero(nChannels, nChannels));
 
    Eigen::FFT<double> fft;
 
    // 2. For each taper, compute FFT for all channels, then accumulate CSD
    for (int t = 0; t < nTap; ++t) {
        const double w = dpss.vecEigenvalues[t];
 
        // Compute FFT spectra for all channels with this taper
        MatrixXcd spectra(nChannels, nFreqsFull);
 
        for (int ch = 0; ch < nChannels; ++ch) {
            VectorXd tapered = (matData.row(ch).array() - matData.row(ch).mean()).transpose() * dpss.matTapers.row(t).transpose().array();
 
            VectorXcd spec;
            fft.fwd(spec, tapered);
 
            spectra.row(ch) = spec.head(nFreqsFull).transpose();
        }
        if (nTimes % 2 == 0)
            spectra.col(nFreqsFull - 1) /= std::sqrt(2.0);
 
        // Accumulate weighted CSD: X(:,f) * X(:,f)^H for selected frequency bins
        for (int fi = 0; fi < nFreqsSel; ++fi) {
            const int fBin = iBinMin + fi;
            VectorXcd col = spectra.col(fBin);
            csdAccum[static_cast<std::size_t>(fi)] += w * (col * col.adjoint());
        }
    }
 
    // 3. Assemble result
    CsdResult result;
    result.vecFreqs.resize(nFreqsSel);
    result.csdByFreq.resize(nFreqsSel);
 
    MatrixXcd meanCsd = MatrixXcd::Zero(nChannels, nChannels);
 
    for (int fi = 0; fi < nFreqsSel; ++fi) {
        // One-sided and per Hz, as mne.time_frequency.csd_array_multitaper
        csdAccum[static_cast<std::size_t>(fi)] *= 2.0 / (weightSum * sfreq);
        result.csdByFreq[fi] = csdAccum[static_cast<std::size_t>(fi)];
        result.vecFreqs[fi] = (iBinMin + fi) * freqRes;
        meanCsd += result.csdByFreq[fi];
    }
 
    if (nFreqsSel > 0)
        meanCsd /= static_cast<double>(nFreqsSel);
 
    result.matCsd = meanCsd;
    return result;
}
 
//=============================================================================================================
 
CsdResult Csd::computeFourier(const MatrixXd& matData,
                              double sfreq,
                              double fmin,
                              double fmax,
                              int nFft,
                              double overlap)
{
    const int nChannels = static_cast<int>(matData.rows());
    const int nTimes = static_cast<int>(matData.cols());
 
    if (nFft > nTimes)
        nFft = nTimes;
 
    const int nFreqsFull = nFft / 2 + 1;
 
    if (fmax < 0.0)
        fmax = sfreq / 2.0;
 
    const double freqRes = sfreq / static_cast<double>(nFft);
    const int iBinMin = std::max(1, static_cast<int>(std::ceil(fmin / freqRes)));
    const int iBinMax = std::min(nFreqsFull - 1, static_cast<int>(std::floor(fmax / freqRes)));
    const int nFreqsSel = std::max(0, iBinMax - iBinMin + 1);
 
    // Build Hann window
    VectorXd hannWin(nFft);
    for (int n = 0; n < nFft; ++n) {
        const double t = static_cast<double>(n) / static_cast<double>(nFft - 1);
        hannWin[n] = 0.5 * (1.0 - std::cos(2.0 * CSD_PI * t));
    }
 
    const int iStep = std::max(1, static_cast<int>(std::round(static_cast<double>(nFft) * (1.0 - overlap))));
 
    // Initialize per-frequency CSD accumulators
    std::vector<MatrixXcd> csdAccum(static_cast<std::size_t>(nFreqsSel),
                                    MatrixXcd::Zero(nChannels, nChannels));
 
    Eigen::FFT<double> fft;
    int nSeg = 0;
 
    // Process each segment
    for (int start = 0; start + nFft <= nTimes; start += iStep) {
        // FFT all channels for this segment
        MatrixXcd spectra(nChannels, nFreqsFull);
 
        for (int ch = 0; ch < nChannels; ++ch) {
            VectorXd seg = matData.row(ch).segment(start, nFft).transpose();
            seg = (seg.array() - seg.mean()) * hannWin.array();
 
            VectorXcd spec;
            fft.fwd(spec, seg);
 
            spectra.row(ch) = spec.head(nFreqsFull).transpose();
        }
        // One-sided spectrum: DC and an even-length Nyquist bin have no mirror image
        spectra.col(0) /= std::sqrt(2.0);
        if (nFft % 2 == 0)
            spectra.col(nFreqsFull - 1) /= std::sqrt(2.0);
 
        // Accumulate CSD for selected frequency bins
        for (int fi = 0; fi < nFreqsSel; ++fi) {
            const int fBin = iBinMin + fi;
            VectorXcd col = spectra.col(fBin);
            csdAccum[static_cast<std::size_t>(fi)] += col * col.adjoint();
        }
 
        ++nSeg;
    }
 
    // 3. Assemble result
    CsdResult result;
    result.vecFreqs.resize(nFreqsSel);
    result.csdByFreq.resize(nFreqsSel);
 
    MatrixXcd meanCsd = MatrixXcd::Zero(nChannels, nChannels);
 
    // As mne.time_frequency.csd_array_fourier: 2 for one-sided, 8/3 for the Hann power loss, per Hz
    const double dNorm = (nSeg > 0) ? 2.0 * 8.0 / 3.0 / (static_cast<double>(nSeg) * nFft * sfreq) : 0.0;
 
    for (int fi = 0; fi < nFreqsSel; ++fi) {
        csdAccum[static_cast<std::size_t>(fi)] *= dNorm;
        result.csdByFreq[fi] = csdAccum[static_cast<std::size_t>(fi)];
        result.vecFreqs[fi] = (iBinMin + fi) * freqRes;
        meanCsd += result.csdByFreq[fi];
    }
 
    if (nFreqsSel > 0)
        meanCsd /= static_cast<double>(nFreqsSel);
 
    result.matCsd = meanCsd;
    return result;
}
 
//=============================================================================================================
 
CsdResult Csd::computeMorlet(const MatrixXd& matData,
                             double sfreq,
                             const RowVectorXd& frequencies,
                             int nCycles)
{
    const int nChannels = static_cast<int>(matData.rows());
    const int nTimes = static_cast<int>(matData.cols());
    const int nFreqs = static_cast<int>(frequencies.size());
 
    CsdResult result;
    result.vecFreqs = frequencies;
    result.csdByFreq.resize(nFreqs);
 
    MatrixXcd meanCsd = MatrixXcd::Zero(nChannels, nChannels);
 
    for (int fi = 0; fi < nFreqs; ++fi) {
        int halfLen = 0;
        const VectorXcd wavelet = MorletTfr::buildWavelet(frequencies[fi], sfreq, nCycles, halfLen);
        const int wavLen = static_cast<int>(wavelet.size());
 
        // Same-length convolution; conj(w[k]) = w[wavLen - 1 - k] for a Morlet wavelet without offset
        MatrixXcd analytic(nChannels, nTimes);
 
        for (int ch = 0; ch < nChannels; ++ch) {
            for (int t = 0; t < nTimes; ++t) {
                std::complex<double> sum(0.0, 0.0);
                for (int k = 0; k < wavLen; ++k) {
                    const int idx = t - halfLen + k;
                    if (idx >= 0 && idx < nTimes) {
                        sum += matData(ch, idx) * std::conj(wavelet[k]);
                    }
                }
                analytic(ch, t) = sum;
            }
        }
 
        // CSD at this frequency = mean over time of C(:,t) * C(:,t)^H
        MatrixXcd csdAtFreq = MatrixXcd::Zero(nChannels, nChannels);
        for (int t = 0; t < nTimes; ++t) {
            VectorXcd col = analytic.col(t);
            csdAtFreq += col * col.adjoint();
        }
        // Scaled by the sampling frequency like mne.time_frequency.csd_array_morlet
        csdAtFreq /= static_cast<double>(nTimes) * sfreq;
 
        result.csdByFreq[fi] = csdAtFreq;
        meanCsd += csdAtFreq;
    }
 
    if (nFreqs > 0)
        meanCsd /= static_cast<double>(nFreqs);
 
    result.matCsd = meanCsd;
    return result;
}