savitzky_golay()

MAYAFLUX_API std::vector< double > MayaFlux::Kinesis::Discrete::savitzky_golay	(	size_t	window,
		size_t	poly_order,
		size_t	derivative = 0
	)

Savitzky-Golay FIR coefficients for smoothing or differentiation.

Fits a polynomial of degree poly_order to a sliding window by least squares and evaluates the requested derivative of that fit at the window centre. With derivative 0 this smooths while preserving peak shape better than a boxcar of the same length; with derivative 1 or 2 it estimates rate of change with far less noise amplification than a raw finite difference, since the differentiation acts on the fitted polynomial rather than on the samples.

This is the block and node form of the kernel that Kinesis::Differential::backward_difference computes per sample on glm types. The streaming form divides by an observed dt and handles irregular timing; this form assumes uniform spacing and returns coefficients in sample units, so a caller needing physical units scales by 1 / dt^derivative.

Returned in FIR tap order: index 0 multiplies the newest sample. Group delay is (window - 1) / 2 samples.

Parameters

window	Kernel length, forced odd, must exceed poly_order
poly_order	Degree of the fitted polynomial
derivative	Which derivative of the fit to evaluate; 0 smooths

Returns

Coefficient vector of length window, or empty if window <= poly_order

Definition at line 216 of file Coefficients.cpp.

{
    const size_t n = force_odd(window);
    if (n <= poly_order || derivative > poly_order)
        return {};
 
    const auto rows = static_cast<Eigen::Index>(n);
    const auto cols = static_cast<Eigen::Index>(poly_order) + 1;
    const double half = static_cast<double>(n - 1) * 0.5;
 
    Eigen::MatrixXd vandermonde(rows, cols);
    for (Eigen::Index i = 0; i < rows; ++i) {
        const double z = static_cast<double>(i) - half;
        double p = 1.0;
        for (Eigen::Index j = 0; j < cols; ++j) {
            vandermonde(i, j) = p;
            p *= z;
        }
    }
 
    const Eigen::MatrixXd pseudo = vandermonde
                                       .completeOrthogonalDecomposition()
                                       .pseudoInverse();
 
    double factorial = 1.0;
    for (size_t k = 2; k <= derivative; ++k)
        factorial *= static_cast<double>(k);
 
    const Eigen::VectorXd row = pseudo.row(static_cast<Eigen::Index>(derivative)) * factorial;
 
    std::vector<double> h(n);
    for (size_t i = 0; i < n; ++i)
        h[i] = row(static_cast<Eigen::Index>(n - 1 - i));
 
    return h;
}

References h, and k.