WelchPsd

Namespace: RTPROCESSINGLIB · Library: DSP Library

Python equivalent

mne.time_frequency.psd_welch in MNE-Python.

#include <dsp/welch_psd.h>

class UTILSLIB::WelchPsd

Welch's averaged-periodogram power spectral density estimator.

Divides each channel into overlapping windowed segments, computes the FFT of each segment, and averages the squared magnitudes. The result is a one-sided PSD normalised so that integrating over frequency recovers the mean signal power.

// 600 Hz data, 512-sample FFT, 50 % overlap, Hann window
WelchPsdResult r = WelchPsd::compute(matData, 600.0, 512);
// r.matPsd   → n_channels × 257
// r.vecFreqs → [0, 1.17, 2.34, …, 300] Hz

Static Methods

compute(matData, dSFreq, iNfft, dOverlap, window, vecPicks)

Compute Welch PSD for every (selected) channel of a data matrix.

Parameters:

matData : const Eigen::MatrixXd & Data matrix (n_channels × n_samples).
dSFreq : double Sampling frequency in Hz.
iNfft : int FFT length / segment length in samples (default 256).
dOverlap : double Fractional overlap between adjacent segments, in [0, 1) (default 0.5).
window : WindowType Window function (default Hann).
vecPicks : const Eigen::RowVectorXi & Channel row indices to process; empty = all channels.

Returns:

WelchPsdResult — WelchPsdResult with matPsd (n_picks × iNfft/2+1) and vecFreqs.

computeVector(vecData, dSFreq, iNfft, dOverlap, window)

Compute Welch PSD for a single channel row vector.

Parameters:

vecData : const Eigen::RowVectorXd & Single-channel data (1 × n_samples row vector).
dSFreq : double Sampling frequency in Hz.
iNfft : int FFT / segment length in samples.
dOverlap : double Fractional segment overlap in [0, 1).
window : WindowType Window function.

Returns:

Eigen::RowVectorXd — One-sided PSD row vector of length iNfft/2+1.

freqAxis(iNfft, dSFreq)

Build the frequency axis for a given FFT length and sampling frequency.

Parameters:

iNfft : int FFT length.
dSFreq : double Sampling frequency in Hz.

Returns:

Eigen::RowVectorXd — Row vector of length iNfft/2+1 with frequencies in Hz.

Example

Source: src/examples/ex_spectral/main.cpp

#include <math.h>
#include <disp/plots/plot.h>
#include <math/spectral.h>
#include <utils/generics/mne_logger.h>
#include <mne/mne.h>

//=============================================================================================================
// Eigen
//=============================================================================================================

#include <Eigen/Core>

//=============================================================================================================
// QT INCLUDES
//=============================================================================================================

#include <QApplication>
#include <QtMath>

//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================

using namespace Eigen;
using namespace DISPLIB;
using namespace UTILSLIB;

//=============================================================================================================
// MAIN
//=============================================================================================================

//=============================================================================================================
/**
 * The function main marks the entry point of the program.
 * By default, main has the storage class extern.
 *
 * @param[in] argc (argument count) is an integer that indicates how many arguments were entered on the command line when the program was started.
 * @param[in] argv (argument vector) is an array of pointers to arrays of character objects. The array objects are null-terminated strings, representing the arguments that were entered on the command line when the program was started.
 * @return the value that was set to exit() (which is 0 if exit() is called via quit()).
 */
int main(int argc, char *argv[])
{
    qInstallMessageHandler(MNELogger::customLogWriter);
    QApplication a(argc, argv);

    // Generate input data
    int iNSamples = 500;
    double dSampFreq = 23.0;
    MatrixXd inputData = MatrixXd::Random(2, iNSamples);
    for (int n = 0; n < iNSamples; n++) {
        inputData(0, n) += 10.0 * sin(2.0 * M_PI *  10. * n / iNSamples );
        inputData(0, n) += 10.0 * sin(2.0 * M_PI * 100. * n / iNSamples);
        inputData(0, n) += 10.0 * sin(2.0 * M_PI * 200. * n / iNSamples);

        inputData(1, n) += 10.0 * sin(2.0 * M_PI *  50. * n / iNSamples);
        inputData(1, n) += 10.0 * sin(2.0 * M_PI * 100. * n / iNSamples);
        inputData(1, n) += 10.0 * sin(2.0 * M_PI * 150. * n / iNSamples);
    }

    //Generate hanning window
    QPair<MatrixXd, VectorXd> tapers = Spectral::generateTapers(iNSamples, QString("hanning"));
    MatrixXd matTaps = tapers.first;
    VectorXd vecTapWeights = tapers.second;

    //Plot hanning window
    VectorXd hann = matTaps.row(0).transpose();
    Plot plotWindow(hann);
    plotWindow.setTitle("Hanning window");
    plotWindow.setXLabel("X Axes");
    plotWindow.setYLabel("Y Axes");
    plotWindow.setWindowTitle("Corresponding function to MATLABs plot");
    plotWindow.show();

    //Compute Spectrum of first row of input data
    int iNfft = iNSamples;
    MatrixXcd matTapSpectrumSeed = Spectral::computeTaperedSpectraRow(inputData.row(0), matTaps, iNfft);

    //Compute PSD
    RowVectorXd vecPsd = Spectral::psdFromTaperedSpectra(matTapSpectrumSeed, vecTapWeights, iNfft, dSampFreq);

    //Plot PSD
    VectorXd psd = vecPsd.transpose();
    Plot plotPsd(psd);
    plotPsd.setTitle("PSD of Row 0");
    plotPsd.setXLabel("X Axes");
    plotPsd.setYLabel("Y Axes");
    plotPsd.setWindowTitle("Corresponding function to MATLABs plot");
    plotPsd.show();

    //Check PSD
    VectorXd psdTest = matTapSpectrumSeed.row(0).cwiseAbs2().transpose();
    psdTest *= 2.0;
    psdTest(0) /= 2.0;
    if (iNfft % 2 == 0){
        psdTest.tail(1) /= 2.0;
    }
    //Normalization
    psdTest /= dSampFreq;

    //Plot PSDTest
    Plot plotPsdTest(psdTest);
    plotPsdTest.setTitle("PSD of Row 0");
    plotPsdTest.setXLabel("X Axes");
    plotPsdTest.setYLabel("Y Axes");
    plotPsdTest.setWindowTitle("Corresponding function to MATLABs plot");
    plotPsdTest.show();

    //Compute CSD of matTapSpectrumSeed and matTapSpectrumSeed
    //The real part should be equivalent to the PSD of matTapSpectrumSeed)
    RowVectorXcd vecCsdSeed = Spectral::csdFromTaperedSpectra(matTapSpectrumSeed, matTapSpectrumSeed, vecTapWeights, vecTapWeights, iNfft, dSampFreq);
    VectorXd psdTest2 = vecCsdSeed.real();

    //Plot PSDTest2
    Plot plotPsdTest2(psdTest2);
    plotPsdTest2.setTitle("PSD of Row 0");
    plotPsdTest2.setXLabel("X Axes");
    plotPsdTest2.setYLabel("Y Axes");
    plotPsdTest2.setWindowTitle("Corresponding function to MATLABs plot");
    plotPsdTest2.show();

    //Check that sums of different psds of the same signal are equal
    qDebug()<<psd.sum();
    qDebug()<<psdTest.sum();
    qDebug()<<psdTest2.sum();
    RowVectorXd data_hann = inputData.row(0).cwiseProduct(matTaps.row(0));
    qDebug()<<data_hann.row(0).cwiseAbs2().sum() * double(iNSamples) / dSampFreq;

    //Compute Spectrum of second row of input data
    MatrixXcd matTapSpectrumTarget = Spectral::computeTaperedSpectraRow(inputData.row(1), matTaps, iNfft);

    //Compute CSD between seed and target
    RowVectorXcd vecCsd = Spectral::csdFromTaperedSpectra(matTapSpectrumSeed, matTapSpectrumTarget, vecTapWeights, vecTapWeights, iNfft, dSampFreq);

    //Plot real part of CSD
    VectorXd csd = vecCsd.transpose().real();
    Plot plotCsd(csd);
    plotCsd.setTitle("Real part of the CSD of Row 0 and 1");
    plotCsd.setXLabel("X Axes");
    plotCsd.setYLabel("Y Axes");
    plotCsd.setWindowTitle("Corresponding function to MATLABs plot");
    plotCsd.show();

    return a.exec();
}

Authors of this file

Christoph Dinh <christoph.dinh@mne-cpp.org>

Static Methods​

compute(matData, dSFreq, iNfft, dOverlap, window, vecPicks)​

computeVector(vecData, dSFreq, iNfft, dOverlap, window)​

freqAxis(iNfft, dSFreq)​

Example​

Authors of this file​

Static Methods

compute(matData, dSFreq, iNfft, dOverlap, window, vecPicks)

computeVector(vecData, dSFreq, iNfft, dOverlap, window)

freqAxis(iNfft, dSFreq)

Example

Authors of this file