dpss.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "dpss.h"
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
 
//=============================================================================================================
// STD INCLUDES
//=============================================================================================================
 
#include <cmath>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace UTILSLIB;
using namespace Eigen;
 
namespace
{
constexpr double DPSS_PI = 3.14159265358979323846;
}
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
DpssResult Dpss::compute(int nSamples, double halfBandwidth, int nTapers, bool sym, bool lowBias)
{
    if (nTapers < 0)
        nTapers = static_cast<int>(std::floor(2.0 * halfBandwidth - 1.0));
 
    const int N = sym ? nSamples : nSamples + 1;
 
    const double W = halfBandwidth / static_cast<double>(N);
 
    // Build the symmetric tridiagonal matrix T
    // Diagonal:     d(i) = ((N-1-2*i)/2)^2 * cos(2*pi*W)
    // Off-diagonal: e(i) = (i+1)*(N-i-1)/2
    VectorXd diag(N);
    VectorXd offdiag(N > 1 ? N - 1 : 0);
 
    const double cosW = std::cos(2.0 * DPSS_PI * W);
 
    for (int i = 0; i < N; ++i) {
        const double val = (static_cast<double>(N - 1) - 2.0 * static_cast<double>(i)) / 2.0;
        diag[i] = val * val * cosW;
    }
 
    for (int i = 0; i < N - 1; ++i) {
        offdiag[i] = static_cast<double>(i + 1) * static_cast<double>(N - i - 1) / 2.0;
    }
 
    // Construct tridiagonal matrix and solve eigenvalue problem
    MatrixXd T = MatrixXd::Zero(N, N);
    for (int i = 0; i < N; ++i)
        T(i, i) = diag[i];
    for (int i = 0; i < N - 1; ++i) {
        T(i, i + 1) = offdiag[i];
        T(i + 1, i) = offdiag[i];
    }
 
    SelfAdjointEigenSolver<MatrixXd> solver(T);
 
    // Eigenvalues are sorted ascending; the tapers are the eigenvectors of the largest nTapers
    const MatrixXd& allEigenvectors = solver.eigenvectors();
 
    DpssResult result;
    result.matTapers.resize(nTapers, N);
    result.vecEigenvalues.resize(nTapers);
 
    // Concentration kernel r(m) = sin(2 pi W m) / (pi m), r(0) = 2W (Percival & Walden 1993, p. 390)
    VectorXd kernel(N);
    kernel[0] = 2.0 * W;
    for (int m = 1; m < N; ++m)
        kernel[m] = std::sin(2.0 * DPSS_PI * W * m) / (DPSS_PI * m);
 
    for (int k = 0; k < nTapers; ++k) {
        const int idx = N - 1 - k; // largest eigenvalue first
 
        // Extract eigenvector as a row, normalize to unit L2 norm
        VectorXd taper = allEigenvectors.col(idx);
        const double norm = taper.norm();
        if (norm > 0.0)
            taper /= norm;
 
        // Sign convention: first element should be positive
        if (taper[0] < 0.0)
            taper = -taper;
 
        result.matTapers.row(k) = taper.transpose();
 
        // Fraction of the taper's energy inside [-W, W]: sum over lags of autocorrelation x kernel
        double ratio = kernel[0] * taper.squaredNorm();
        for (int m = 1; m < N; ++m)
            ratio += 2.0 * kernel[m] * taper.head(N - m).dot(taper.tail(N - m));
        result.vecEigenvalues[k] = ratio;
    }
 
    if (lowBias) {
        int nKeep = 1;
        while (nKeep < nTapers && result.vecEigenvalues[nKeep] > 0.9)
            ++nKeep;
        result.matTapers.conservativeResize(nKeep, Eigen::NoChange);
        result.vecEigenvalues.conservativeResize(nKeep);
    }
    if (!sym)
        result.matTapers.conservativeResize(Eigen::NoChange, nSamples);
 
    return result;
}