Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 1 addition & 4 deletions numerical/analysis/ConvolutionCorrelation.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -127,11 +127,8 @@ namespace analysis
for (std::size_t i = 0; i < M + K - 1; ++i)
y[i] = yFull[i];
}
}

#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
namespace analysis
{
extern template void LinearConvolution<float, 3, 3>(
const infra::BoundedVector<float>::WithMaxSize<3>&,
const infra::BoundedVector<float>::WithMaxSize<3>&,
Expand Down Expand Up @@ -159,5 +156,5 @@ namespace analysis
const infra::BoundedVector<float>::WithMaxSize<3>&,
infra::BoundedVector<float>::WithMaxSize<5>&,
FastFourierTransform<float>&);
}
#endif
}
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,13 @@ namespace
};
}

TYPED_TEST(TestFastFourierTransform, log2_runtime_both_branches)
{
auto& fftInst = *this->fft;
EXPECT_EQ(fftInst.Log2(1), 0u);
EXPECT_EQ(fftInst.Log2(8), 3u);
}

TYPED_TEST(TestFastFourierTransform, zero_input_produces_zero_output)
{
this->timeDomain.clear();
Expand Down
24 changes: 24 additions & 0 deletions numerical/analysis/windowing/test/TestWindowing.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -123,3 +123,27 @@ TYPED_TEST(WindowingTest, WindowSymmetry)
}
}
}

TYPED_TEST(WindowingTest, HammingWindowPower)
{
windowing::HammingWindow<TypeParam> w;
EXPECT_NEAR(math::ToFloat(w.Power(8)), 0.397f, this->kEpsilon);
}

TYPED_TEST(WindowingTest, HanningWindowPower)
{
windowing::HanningWindow<TypeParam> w;
EXPECT_NEAR(math::ToFloat(w.Power(8)), 0.375f, this->kEpsilon);
}

TYPED_TEST(WindowingTest, BlackmanWindowPower)
{
windowing::BlackmanWindow<TypeParam> w;
EXPECT_NEAR(math::ToFloat(w.Power(8)), 0.305f, this->kEpsilon);
}

TYPED_TEST(WindowingTest, RectangularWindowPower)
{
windowing::RectangularWindow<TypeParam> w;
EXPECT_NEAR(math::ToFloat(w.Power(8)), 0.999f, this->kEpsilon);
}
5 changes: 4 additions & 1 deletion numerical/control_analysis/ControllabilityObservability.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -131,8 +131,11 @@ namespace control_analysis
std::size_t rank{ 0 };
std::size_t pivotRow{ 0 };

for (std::size_t col = 0; col < Cols && pivotRow < Rows; ++col)
for (std::size_t col = 0; col < Cols; ++col)
{
if (pivotRow >= Rows)
break;

std::size_t maxRow{ pivotRow };
T maxVal{ T(0) };
for (std::size_t r = pivotRow; r < Rows; ++r)
Expand Down
15 changes: 15 additions & 0 deletions numerical/control_analysis/test/TestFrequencyResponse.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -317,3 +317,18 @@ TEST_F(TestFrequencyResponseHighpass, magnitude_at_quarter_nyquist_is_minus_3db)

EXPECT_NEAR(bestMag, -3.0103f, 0.5f);
}

TEST_F(TestFrequencyResponseUnity, zero_denominator_coefficients_produce_finite_output)
{
std::array<float, 1> bz{ 1.0f };
std::array<float, 1> az{ 0.0f };
control_analysis::FrequencyResponse<float, 64> frZeroDenom{ bz, az, kSampleFrequency };

auto [frequencies, magnitudes, phases] = frZeroDenom.Calculate();

for (const auto& m : magnitudes)
{
EXPECT_FALSE(std::isnan(m));
EXPECT_FALSE(std::isinf(m));
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -35,13 +35,12 @@ namespace controllers

private:
std::array<SchedulePoint<T, GainSize>, N> table;
std::array<T, GainSize> active;
std::array<T, GainSize> active{};
};

template<typename T, std::size_t N, std::size_t GainSize>
GainScheduledController<T, N, GainSize>::GainScheduledController(std::array<SchedulePoint<T, GainSize>, N> scheduleTable)
: table{ scheduleTable }
, active{}
{
for (std::size_t i{ 0 }; i < N - 1; ++i)
assert(table[i].breakpoint < table[i + 1].breakpoint);
Expand Down
3 changes: 1 addition & 2 deletions numerical/controllers/implementations/LuenbergerObserver.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -45,15 +45,14 @@ namespace controllers

Plant plant;
GainMatrix L;
StateVector xhat;
StateVector xhat{};
};

template<typename T, std::size_t StateSize, std::size_t InputSize, std::size_t OutputSize>
LuenbergerObserver<T, StateSize, InputSize, OutputSize>::LuenbergerObserver(
const Plant& plantModel, const GainMatrix& observerGain)
: plant{ plantModel }
, L{ observerGain }
, xhat{}
{}

template<typename T, std::size_t StateSize, std::size_t InputSize, std::size_t OutputSize>
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -195,3 +195,32 @@ TEST_F(TestIntegralStateFeedbackLqi, control_uses_negative_feedback)
EXPECT_NEAR(u.at(0, 0), expected, math::Tolerance<float>());
EXPECT_LT(u.at(0, 0), 0.0f);
}

TEST_F(TestIntegralStateFeedbackLqi, direct_gain_constructor_matches_lqr_constructor)
{
const auto kx = controller.GetGainState();
const auto ki = controller.GetGainIntegral();

Controller directController{ kx, ki, 0.01f };

math::Vector<float, 2> x{ { 0.5f }, { 0.1f } };
math::Vector<float, 1> reference{ { 1.0f } };
math::Vector<float, 1> measured{ { 0.3f } };

auto u1 = controller.ComputeControl(x, reference, measured);
auto u2 = directController.ComputeControl(x, reference, measured);

EXPECT_NEAR(u1.at(0, 0), u2.at(0, 0), math::Tolerance<float>());
}

TEST_F(TestIntegralStateFeedbackLqi, direct_gain_constructor_get_gains)
{
math::Matrix<float, 1, 2> kx{ { 1.5f, 0.8f } };
math::Matrix<float, 1, 1> ki{ { 2.0f } };

Controller c{ kx, ki, 0.01f };

EXPECT_FLOAT_EQ(c.GetGainState().at(0, 0), 1.5f);
EXPECT_FLOAT_EQ(c.GetGainState().at(0, 1), 0.8f);
EXPECT_FLOAT_EQ(c.GetGainIntegral().at(0, 0), 2.0f);
}
3 changes: 1 addition & 2 deletions numerical/estimators/offline/ExpectationMaximization.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -104,8 +104,7 @@ namespace estimators
while (iter < maxIterations)
{
const auto smootherOutput = smoother_.Smooth(
currentParams.F, currentParams.H,
currentParams.Q, currentParams.R,
{ currentParams.F, currentParams.H, currentParams.Q, currentParams.R },
observations, numSteps,
currentParams.initialState, currentParams.initialCovariance);

Expand Down
18 changes: 18 additions & 0 deletions numerical/estimators/offline/test/TestPolynomialFitting.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -180,6 +180,24 @@ TEST_F(TestPolynomialFitting, determinism_identical_fit_produces_identical_coeff
EXPECT_FLOAT_EQ(c1.at(2, 0), c2.at(2, 0));
}

TEST_F(TestPolynomialFitting, predict_evaluates_fitted_polynomial)
{
math::Matrix<float, 8, 1> x;
math::Matrix<float, 8, 1> y;

for (std::size_t i = 0; i < 8; ++i)
{
float xi = static_cast<float>(i) * 0.25f;
x.at(i, 0) = xi;
y.at(i, 0) = 1.0f - 0.5f * xi + 0.25f * xi * xi;
}

fitter.Fit(x, y);

EXPECT_NEAR(fitter.Predict(0.5f), 1.0f - 0.5f * 0.5f + 0.25f * 0.25f, math::Tolerance<float>());
EXPECT_NEAR(fitter.Predict(1.0f), 1.0f - 0.5f * 1.0f + 0.25f * 1.0f, math::Tolerance<float>());
Comment thread
gabrielfrasantos marked this conversation as resolved.
}

TEST_F(TestPolynomialFitting, predict_all_coefficients_finite_on_far_from_origin_data)
{
math::Matrix<float, 8, 1> x;
Expand Down
2 changes: 0 additions & 2 deletions numerical/filters/active/AlphaBetaFilter.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,7 @@ namespace filters
: samplePeriod{ Ts }
, gainAlpha{ alpha }
, gainBeta{ beta }
, gainGamma{ T{} }
, betaOverTs{ beta / Ts }
, twoGammaOverTs2{ T{} }
{}

template<typename T, std::size_t Order>
Expand Down
48 changes: 22 additions & 26 deletions numerical/filters/active/KalmanSmoother.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -40,11 +40,16 @@ namespace filters

KalmanSmoother() = default;

struct SmootherParams
{
const StateMatrix& F;
const MeasurementMatrix& H;
const StateMatrix& Q;
const MeasurementCovariance& R;
};

OPTIMIZE_FOR_SPEED SmootherOutput Smooth(
const StateMatrix& F,
const MeasurementMatrix& H,
const StateMatrix& Q,
const MeasurementCovariance& R,
const SmootherParams& params,
const std::array<MeasurementVector, MaxSteps>& observations,
std::size_t numSteps,
const StateVector& initialState,
Expand All @@ -68,10 +73,7 @@ namespace filters
std::array<KalmanGainMatrix, MaxSteps> kalmanGains_{};

OPTIMIZE_FOR_SPEED float RunForwardPass(
const StateMatrix& F,
const MeasurementMatrix& H,
const StateMatrix& Q,
const MeasurementCovariance& R,
const SmootherParams& params,
const std::array<MeasurementVector, MaxSteps>& observations,
std::size_t numSteps,
const StateVector& initialState,
Expand All @@ -97,10 +99,7 @@ namespace filters
OPTIMIZE_FOR_SPEED
typename KalmanSmoother<StateSize, MeasurementSize, MaxSteps>::SmootherOutput
KalmanSmoother<StateSize, MeasurementSize, MaxSteps>::Smooth(
const StateMatrix& F,
const MeasurementMatrix& H,
const StateMatrix& Q,
const MeasurementCovariance& R,
const SmootherParams& params,
const std::array<MeasurementVector, MaxSteps>& observations,
std::size_t numSteps,
const StateVector& initialState,
Expand All @@ -109,18 +108,15 @@ namespace filters
really_assert(numSteps >= 2 && numSteps <= MaxSteps);

SmootherOutput output;
output.logLikelihood = RunForwardPass(F, H, Q, R, observations, numSteps, initialState, initialCovariance);
RunBackwardPass(F, H, numSteps, output);
output.logLikelihood = RunForwardPass(params, observations, numSteps, initialState, initialCovariance);
RunBackwardPass(params.F, params.H, numSteps, output);
return output;
}

template<std::size_t StateSize, std::size_t MeasurementSize, std::size_t MaxSteps>
OPTIMIZE_FOR_SPEED float
KalmanSmoother<StateSize, MeasurementSize, MaxSteps>::RunForwardPass(
const StateMatrix& F,
const MeasurementMatrix& H,
const StateMatrix& Q,
const MeasurementCovariance& R,
const SmootherParams& params,
const std::array<MeasurementVector, MaxSteps>& observations,
std::size_t numSteps,
const StateVector& initialState,
Expand All @@ -133,27 +129,27 @@ namespace filters

for (std::size_t t = 0; t < numSteps; ++t)
{
const auto nu = observations[t] - H * predictedMeans_[t];
const auto S = math::CongruenceTransform(H, predictedCovariances_[t]) + R;
const auto nu = observations[t] - params.H * predictedMeans_[t];
const auto S = math::CongruenceTransform(params.H, predictedCovariances_[t]) + params.R;

const auto K = solvers::SolveSystem<float, MeasurementSize, StateSize>(
S, H * predictedCovariances_[t])
S, params.H * predictedCovariances_[t])
.Transpose();

filteredMeans_[t] = predictedMeans_[t] + K * nu;

const auto IminusKH = StateMatrix::Identity() - K * H;
const auto IminusKH = StateMatrix::Identity() - K * params.H;
filteredCovariances_[t] =
math::CongruenceTransform(IminusKH, predictedCovariances_[t]) + math::CongruenceTransform(K, R);
math::CongruenceTransform(IminusKH, predictedCovariances_[t]) + math::CongruenceTransform(K, params.R);

kalmanGains_[t] = K;

logLikelihood += ComputeLogLikelihoodContribution(nu, S);

if (t < numSteps - 1)
{
predictedMeans_[t + 1] = F * filteredMeans_[t];
predictedCovariances_[t + 1] = math::CongruenceTransform(F, filteredCovariances_[t]) + Q;
predictedMeans_[t + 1] = params.F * filteredMeans_[t];
predictedCovariances_[t + 1] = math::CongruenceTransform(params.F, filteredCovariances_[t]) + params.Q;
}
}

Expand Down Expand Up @@ -244,7 +240,7 @@ namespace filters
const StateVector& initialState,
const StateMatrix& initialCovariance)
{
return Smooth(plant.A, plant.C, Q, R, observations, numSteps, initialState, initialCovariance);
return Smooth({ plant.A, plant.C, Q, R }, observations, numSteps, initialState, initialCovariance);
}

#ifdef NUMERICAL_TOOLBOX_COVERAGE_BUILD
Expand Down
12 changes: 6 additions & 6 deletions numerical/filters/active/test/TestKalmanSmoother.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ namespace

Smoother::SmootherOutput RunSmoother()
{
return smoother.Smooth(F, H, Q, R, observations, T, x0, P0);
return smoother.Smooth({F, H, Q, R}, observations, T, x0, P0);
Comment thread
gabrielfrasantos marked this conversation as resolved.
}
};

Expand Down Expand Up @@ -80,7 +80,7 @@ namespace

Smoother::SmootherOutput RunSmoother()
{
return smoother.Smooth(F, H, Q, R, obs, T, x0, P0);
return smoother.Smooth({F, H, Q, R}, obs, T, x0, P0);
Comment thread
gabrielfrasantos marked this conversation as resolved.
}
};
}
Expand Down Expand Up @@ -210,7 +210,7 @@ TEST_F(TestKalmanSmoother, smoothed_rmse_is_less_than_filtered_rmse)
}

Smoother synthSmoother;
const auto output = synthSmoother.Smooth(F, H, Q, R, syntheticObs, T, x0, P0);
const auto output = synthSmoother.Smooth({F, H, Q, R}, syntheticObs, T, x0, P0);
Comment thread
gabrielfrasantos marked this conversation as resolved.

float filteredMse = 0.0f;
float smoothedMse = 0.0f;
Expand Down Expand Up @@ -270,7 +270,7 @@ TEST_F(TestKalmanSmoother, smoothed_covariance_at_final_step_equals_filtered_cov
TEST_F(TestKalmanSmoother, minimum_steps_boundary_produces_finite_output)
{
constexpr std::size_t minSteps = 2;
const auto output = smoother.Smooth(F, H, Q, R, observations, minSteps, x0, P0);
const auto output = smoother.Smooth({F, H, Q, R}, observations, minSteps, x0, P0);
Comment thread
gabrielfrasantos marked this conversation as resolved.

EXPECT_TRUE(std::isfinite(output.logLikelihood));
for (std::size_t t = 0; t < minSteps; ++t)
Expand All @@ -281,7 +281,7 @@ TEST_F(TestKalmanSmoother, minimum_steps_boundary_produces_finite_output)
TEST_F(TestKalmanSmoother, near_zero_process_noise_produces_no_nan)
{
StateMatrix Qsmall{ { 1e-6f, 0.0f }, { 0.0f, 1e-6f } };
const auto output = smoother.Smooth(F, H, Qsmall, R, observations, T, x0, P0);
const auto output = smoother.Smooth({F, H, Qsmall, R}, observations, T, x0, P0);
Comment thread
gabrielfrasantos marked this conversation as resolved.

EXPECT_TRUE(std::isfinite(output.logLikelihood));
for (std::size_t t = 0; t < T; ++t)
Expand Down Expand Up @@ -321,7 +321,7 @@ TEST_F(TestKalmanSmoother, nees_within_chi_squared_band_on_simulated_system)
trialObs[t] = MeasurementVector{ { trialTrue[t].at(0, 0) + measNoise(rng) } };

Smoother trialSmoother;
const auto out = trialSmoother.Smooth(F, H, Q, R, trialObs, T, x0, P0);
const auto out = trialSmoother.Smooth({F, H, Q, R}, trialObs, T, x0, P0);
Comment thread
gabrielfrasantos marked this conversation as resolved.

for (std::size_t t = 0; t < T; ++t)
{
Expand Down
3 changes: 1 addition & 2 deletions numerical/filters/passive/MedianFilter.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,11 @@ namespace filters::passive
private:
std::array<T, N> window;
std::array<T, N> scratch;
std::size_t head;
std::size_t head{ 0 };
};

template<typename T, std::size_t N>
MedianFilter<T, N>::MedianFilter(T initial) noexcept
: head{ 0 }
{
window.fill(initial);
scratch.fill(T{});
Expand Down
Loading
Loading