MarkovProcess.hpp

Go to the documentation of this file.

// First-Order Gauss-Markov Process
// Exponentially correlated noise model
#pragma once
 
#include <cmath>
#include <janus/janus.hpp>
 
namespace vulcan::markov {
 
// =============================================================================
// Process State
// =============================================================================

template <typename Scalar> struct State {
    Scalar x; 
};

template <typename Scalar> State<Scalar> init_state() {
    return State<Scalar>{.x = Scalar(0)};
}

template <typename Scalar>
State<Scalar> init_state(const Scalar &initial_value) {
    return State<Scalar>{.x = initial_value};
}
 
// =============================================================================
// Discretized Process Parameters
// =============================================================================

struct Coeffs {
    double phi; 
    double q;   
};

inline Coeffs discretize(double tau, double sigma, double dt) {
    double phi = std::exp(-dt / tau);
    double q = sigma * std::sqrt(1.0 - phi * phi);
    return Coeffs{.phi = phi, .q = q};
}
 
// =============================================================================
// Continuous-Time Parameters
// =============================================================================

inline double process_noise_psd(double tau, double sigma) {
    return 2.0 * sigma * sigma / tau;
}

inline double process_noise_intensity(double tau, double sigma) {
    return std::sqrt(2.0 / tau) * sigma;
}
 
// =============================================================================
// Simulation Step
// =============================================================================

template <typename Scalar>
Scalar step(State<Scalar> &state, const Coeffs &coeffs,
            const Scalar &noise_input) {
    state.x = coeffs.phi * state.x + coeffs.q * noise_input;
    return state.x;
}

template <typename Scalar>
Scalar step(State<Scalar> &state, double tau, double sigma, double dt,
            const Scalar &noise_input) {
    Coeffs coeffs = discretize(tau, sigma, dt);
    return step(state, coeffs, noise_input);
}
 
// =============================================================================
// Analysis Utilities
// =============================================================================

inline double autocorrelation(double tau, double sigma, double tau_lag) {
    return sigma * sigma * std::exp(-std::abs(tau_lag) / tau);
}

inline double psd_at_frequency(double tau, double sigma, double f) {
    double omega_tau = 2.0 * M_PI * f * tau;
    return sigma * sigma * 2.0 * tau / (1.0 + omega_tau * omega_tau);
}
 
} // namespace vulcan::markov