[AUD-027][Medium] Preserve pole semantics in frequency response

[gabrielsantosphilips]

Audit finding

• ID: AUD-027
• Status: Verified defect
• Severity: Medium
• Confidence: High
• Audited revision: 2f479320d805a1f9f35ebe4afaaeeded48913a94

Problem

When the evaluated transfer denominator is exactly zero, FrequencyResponse::Calculate() replaces it with 1+0i. An exact pole is therefore reported as the numerator magnitude instead of infinity, a documented saturation, or failure. The current test encodes this false finite result.

Source:

numerical-toolbox-cpp/numerical/control_analysis/FrequencyResponse.hpp

Lines 63 to 80 in 2f47932

	std::complex<T> numerator(static_cast<T>(0), static_cast<T>(0));
	std::complex<T> denominator(static_cast<T>(0), static_cast<T>(0));
	for (std::size_t i = 0; i < b.size(); ++i)
	numerator += b[i] * std::polar(static_cast<T>(1), -omega * static_cast<T>(i));
	for (std::size_t i = 0; i < a.size(); ++i)
	denominator += a[i] * std::polar(static_cast<T>(1), -omega * static_cast<T>(i));
	if (denominator == static_cast<T>(0))
	denominator.real(static_cast<T>(1));
	auto h = numerator / denominator;
	auto magnitude = std::max(std::abs(h), std::numeric_limits<T>::min());
	frequencies.emplace_back(f);
	response.emplace_back(static_cast<T>(20) * math::Log10(magnitude));
	phase.emplace_back(std::arg(h) * static_cast<T>(180) / static_cast<T>(std::numbers::pi));

Acceptance criteria

• Define pole and near-pole behavior in the API/documentation.
• Preserve infinity or return explicit status rather than substituting an unrelated denominator.
• Add a known-pole test with an independent expected outcome.
• Handle finite-range dB display separately from transfer-function evaluation.

[gabrielfrasantos]

Fixed on claude/repo-issues-triage-pgbb6e, together with #303 since both live in Calculate(). Two notes.

Confirmed the reachability question first, because it matters for how much this was worth doing. The substitution only fires on an exactly zero denominator, which a log-spaced frequency grid essentially never hits for a normal transfer function — the grid starts at fs/N and never samples omega = 0. It is reachable for a degenerate denominator such as a = {0}, and there the old code reported 13.98 dB, precisely the numerator magnitude 20·log10(5), exactly as described.

Behaviour now:

• result buffers are cleared at the start of each call, so Calculate() is reusable ([AUD-010][High] Make FrequencyResponse::Calculate reusable #303);
• an exactly-zero denominator stores +infinity in the magnitude vector rather than substituting 1;
• ClampDb(value, floorDb, ceilingDb) is added for finite-range display, keeping presentation separate from transfer-function evaluation.

Note 1 — API behaviour change. The magnitude vector can now contain +infinity. Anything plotting these values directly needs ClampDb, or it will get an unbounded axis. Nothing in the repo consumes FrequencyResponse outside its tests today, so nothing is broken right now, but it is a trap for the first consumer.

Note 2 — worth knowing repo-wide. Writing the test for this surfaced something that is not specific to this issue: std::isfinite returns true for a value whose bits are 0x7f800000, in these translation units. Every one of these headers opens with #pragma GCC optimize("O3", "fast-math"), and -ffinite-math-only lets the compiler fold the classification functions to constants. The stored value really is infinity — verified via std::bit_cast — but std::isfinite will not tell you so.

That means AGENTS.md's "unit tests call std:: directly, never math::" is a hazard for exactly these functions. This test asserts on the bit pattern instead. I also added a bit-inspecting math::IsFinite to Math.hpp for production use, since the same trap applies to every finiteness guard added across this batch of fixes. Might be worth a line in AGENTS.md.

Tests added: RepeatedCalculateOnSameObjectIsIdentical, ExactPoleReportsInfinityNotNumeratorMagnitude, ClampDbGivesFiniteDisplayRange.

---

Generated by Claude Code