From 40f934901da3455d57dc20885fb67271eafd57f2 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Wed, 19 Aug 2026 19:14:46 +0100 Subject: [PATCH] Retune the Kalman process noise, and fuse the two thermometers Two changes to the same signal path, one large and one small. The large one: all three filters were tuned to track one to three decades faster than their signals move. In a still room the temperature filter reported a median rate of 12.4 C/h while the air moved 0.4 C/h, and it overshot a real -36 C/h event by 77%. Sweeping q against the RMSE of the reported rate versus the true rate, using noise measured on the board (temperature 0.088 C, pressure 0.022 hPa, humidity 0.40 %): temperature 6.45 -> 0.37 C/h RMSE 2e-6 -> 1e-9 pressure 2.15 -> 0.24 hPa/h RMSE 1e-5 -> 1e-8 humidity 27.94 -> 3.55 %/h RMSE 5e-5 -> 2e-8 Tracking does not suffer. Lag against a genuine 2 C/h ramp is 0.003 C at both the old and new values, and the peak response to a five-minute event moves closer to the truth rather than further from it, because the overshoot goes away. What is given up is response to sub-minute transients, which for a station forecasting fifteen minutes to a day ahead is noise to reject. This matters most for pressure, whose tendency drives the precipitation forecast, and which was the worst tuned of the three. config.yaml shadowed kalman_q_temp, so editing the dataclass alone changed nothing. All six values are now listed there with that hazard spelled out, because a silent shadow cost real time here. The small one: temp_raw was the plain average of two thermometers whose white-noise sds differ by 7x (LPS25HB 0.007 C, HTS221 0.049 C), which throws the quiet one away. Inverse-variance weighting cuts the raw noise 3.5x. The trap is that the chips do not agree. They sit at different distances from the SoC and stand about 1.3 C apart, so weighting by variance alone drags temp_raw 0.48 C onto the LPS25HB, which after the 1.55x gain of the inverse compensator is 0.75 C of silent bias on every reading, since k was fitted against the mean of the two. The gradient is therefore tracked and removed before weighting and only the deviations are fused: measured mean shift 0.0001 C, noise still 3.5x lower. The tracked gradient is retained because it is a second observation of self-heating. Also corrected: the earlier claim that the HTS221 was the quieter channel was wrong, taken from twelve samples at a cadence slow enough that real drift dominated. At 0.5 s over 120 samples the LPS25HB is quieter by 7x and takes 98% of the weight. --- ashvale/config.py | 24 +++++++++-- ashvale/sensors.py | 89 +++++++++++++++++++++++++++++++++++++++- config.yaml | 13 +++++- tests/test_estimation.py | 85 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 205 insertions(+), 6 deletions(-) diff --git a/ashvale/config.py b/ashvale/config.py index 05c955c..d579482 100644 --- a/ashvale/config.py +++ b/ashvale/config.py @@ -110,11 +110,29 @@ class SensorConfig: hum_offset_min: float = -35.0 hum_offset_max: float = 35.0 # Kalman process/measurement noise (per-signal) - kalman_q_temp: float = 2.0e-6 + kalman_q_temp: float = 1.0e-9 kalman_r_temp: float = 0.02 - kalman_q_press: float = 1.0e-5 + # Process noise, retuned against measured sensor noise rather than guessed. + # + # The originals tracked far faster than any of these signals move. In a + # still room the temperature filter reported a median rate of 12.4 C/h + # while the air moved 0.37 C/h, and it overshot a real -36 C/h event by + # 77%. Sweeping q against the RMSE of the reported rate versus the true + # rate, using noise measured on the board (temp 0.088 C, press 0.022 hPa, + # hum 0.40 %), puts the minimum about two to three decades lower: + # + # temperature 6.45 -> 0.37 C/h RMSE at 2e-6 -> 1e-9 + # pressure 2.15 -> 0.24 hPa/h RMSE at 1e-5 -> 1e-8 + # humidity 27.94 -> 3.55 %/h RMSE at 5e-5 -> 2e-8 + # + # Tracking does not suffer: lag against a genuine 2 C/h ramp is 0.003 C at + # both the old and new values, and peak response to a 5-minute event is + # closer to the truth, not further from it. What is lost is response to + # sub-minute transients, which for a station forecasting 15 minutes to a + # day ahead is noise to reject rather than signal to chase. + kalman_q_press: float = 1.0e-8 kalman_r_press: float = 0.05 - kalman_q_hum: float = 5.0e-5 + kalman_q_hum: float = 2.0e-8 kalman_r_hum: float = 0.60 diff --git a/ashvale/sensors.py b/ashvale/sensors.py index e26e2d6..3500892 100644 --- a/ashvale/sensors.py +++ b/ashvale/sensors.py @@ -65,7 +65,41 @@ def read_cpu_temperature() -> float: # therefore unchanged and only the per-channel detail is new, which matters # because that gradient is a second observation of self-heating. K_HTS221, K_LPS25HB = 0.6164, 0.4889 -SD_HTS221, SD_LPS25HB = 0.060, 0.443 +SD_HTS221, SD_LPS25HB = 0.049, 0.007 + + +class _ChannelNoise: + """Running white-noise variance of one thermometer. + + Taken from the first difference rather than a windowed variance. Over one + 2 s sample the air moves far less than either chip's own jitter, so + var(diff)/2 is the noise and is blind to the weather underneath it. A + windowed variance would measure the weather instead and would rise, not + fall, on a calm day. + """ + + def __init__(self, prior_sd: float, lam: float = 0.995, warmup: int = 200): + self.var = float(prior_sd) ** 2 + self.prior = self.var + self.lam = float(lam) + self.warmup = int(warmup) + self.last: Optional[float] = None + self.n = 0 + + def update(self, value: float) -> float: + if not math.isfinite(value): + return max(self.var, 1e-8) + if self.last is not None: + d = value - self.last + self.var = self.lam * self.var + (1.0 - self.lam) * (d * d / 2.0) + self.n += 1 + self.last = value + if self.n < self.warmup: + # Blend toward the prior while the estimate is young, so one quiet + # minute cannot hand a channel 100% of the weight on no evidence. + w = self.n / float(self.warmup) + return max(w * self.var + (1.0 - w) * self.prior, 1e-8) + return max(self.var, 1e-8) class SimulatedBoard: @@ -228,6 +262,13 @@ class SenseBoard: self.bus = None self.tcs_addr = tcs_addr self._sim = SimulatedBoard(latitude, longitude) + self._noise_h = _ChannelNoise(SD_HTS221) + self._noise_p = _ChannelNoise(SD_LPS25HB) + # Slow EWMA of the standing gradient between the two chips. About a + # 10-minute time constant at the 2 s cadence: long enough to ignore + # per-sample noise, short enough to follow a real change in SoC load. + self._gradient: Optional[float] = None + self._gradient_lam = 0.9967 try: from sense_hat import SenseHat # type: ignore @@ -264,6 +305,48 @@ class SenseBoard: except Exception: return {"clear": 0, "red": 0, "green": 0, "blue": 0, "hex": "#334155", "cct": None} + def _fuse(self, t_h: float, t_p: float) -> tuple[float, float]: + """Combine the two thermometers by inverse variance. + + A plain average of a quiet sensor and a noisy one throws the quiet one + away. Measured on the board at 0.5 s: the LPS25HB carries a white-noise + sd of 0.007 C against the HTS221's 0.049 C, so optimal weighting is + about 98/2 and cuts the raw noise by roughly 3.7x. + + The trap is that the two chips do not agree. They sit at different + distances from the SoC and stand about 1.3 C apart, so weighting them + by variance would drag temp_raw most of the way onto the LPS25HB and + shift it by more than half a degree. The compensator's k was fitted + against the mean of the two, and after the 1.55x gain of the inverse + model that is a full degree of silent bias on every reading and every + forecast built from it. + + So the gradient is tracked and removed before weighting, and only the + deviations are fused. The mean is left exactly where the average put + it, k stays valid, and the noise still falls. The gradient itself is + kept because it is a second observation of self-heating and is what + would let k be identified without a reference thermometer. + """ + if not (math.isfinite(t_h) and math.isfinite(t_p)): + good = [v for v in (t_h, t_p) if math.isfinite(v)] + return (good[0] if good else float("nan")), float("nan") + + var_h = self._noise_h.update(t_h) + var_p = self._noise_p.update(t_p) + + gap = t_h - t_p + if self._gradient is None: + self._gradient = gap + else: + lam = self._gradient_lam + self._gradient = lam * self._gradient + (1.0 - lam) * gap + + # Centre both channels on what the plain average would have reported. + half = self._gradient / 2.0 + w_h, w_p = 1.0 / var_h, 1.0 / var_p + fused = (w_h * (t_h - half) + w_p * (t_p + half)) / (w_h + w_p) + return float(fused), float(1.0 / (w_h + w_p)) + def read(self) -> Dict[str, Any]: """One full multi-sensor sample. Raw, uncompensated, untouched.""" if not self.available: @@ -276,6 +359,7 @@ class SenseBoard: s = self.sense t_h = s.get_temperature_from_humidity() t_p = s.get_temperature_from_pressure() + temp_raw, temp_var = self._fuse(t_h, t_p) orientation = s.get_orientation_degrees() accel = s.get_accelerometer_raw() gyro = s.get_gyroscope_raw() @@ -285,7 +369,8 @@ class SenseBoard: return v - 360.0 if v > 180.0 else v return { - "temp_raw": (t_h + t_p) / 2.0, + "temp_raw": temp_raw, + "temp_var": temp_var, "temp_h": t_h, "temp_p": t_p, "hum": s.get_humidity(), diff --git a/config.yaml b/config.yaml index 2052b5f..a8774b8 100644 --- a/config.yaml +++ b/config.yaml @@ -14,8 +14,19 @@ sensor: persist_period_s: 30.0 rotation_deg: 90 cpu_heat_k: 0.55 # starting point only, calibrate from the dashboard - kalman_q_temp: 2.0e-6 # raise to track faster, lower to smooth harder + # Process noise. Raise to track faster, lower to smooth harder. These were + # retuned against noise measured on a real board by sweeping each q against + # the RMSE of the reported rate versus the true rate. The originals tracked + # two to three decades faster than any of these signals move: in a still room + # the temperature filter reported a median rate of 12.4 C/h while the air + # moved 0.4 C/h. All three are listed because this file shadows the defaults + # in ashvale/config.py, and a value present here silently wins. + kalman_q_temp: 1.0e-9 kalman_r_temp: 0.02 + kalman_q_press: 1.0e-8 + kalman_r_press: 0.05 + kalman_q_hum: 2.0e-8 + kalman_r_hum: 0.60 model: grid_s: 300 diff --git a/tests/test_estimation.py b/tests/test_estimation.py index 18bd68b..b1079bc 100644 --- a/tests/test_estimation.py +++ b/tests/test_estimation.py @@ -248,3 +248,88 @@ def test_forecast_head_migrates_state_from_before_the_setpoint_member(): # discover a migration bug. assert back.member_mae.size == len(MEMBERS) back.learn(np.zeros(4), 20.0, 20.5, 0.1, 0.2) # must not raise + + +# ------------------------------------------------- dual-thermometer fusion + +def _bare_board(): + from ashvale.sensors import SD_HTS221, SD_LPS25HB, SenseBoard, _ChannelNoise + b = SenseBoard.__new__(SenseBoard) + b._noise_h = _ChannelNoise(SD_HTS221) + b._noise_p = _ChannelNoise(SD_LPS25HB) + b._gradient = None + b._gradient_lam = 0.9967 + return b + + +def _two_channels(n=4000, seed=5): + from ashvale.sensors import K_HTS221, K_LPS25HB + rng = np.random.default_rng(seed) + cpu = 43.0 + 0.5 * np.sin(np.arange(n) / 500.0) + th = (24.0 + K_HTS221 * cpu) / (1 + K_HTS221) + 0.049 * rng.normal(size=n) + tp = (24.0 + K_LPS25HB * cpu) / (1 + K_LPS25HB) + 0.007 * rng.normal(size=n) + return th, tp + + +def test_fusion_does_not_move_the_mean(): + """The whole point of removing the gradient first. + + The two chips stand about 1.3 C apart, so weighting them by variance drags + temp_raw onto the quieter one. k was fitted against the mean of the two, and + after the 1.55x gain of the inverse model that shift becomes about a degree + of silent bias on every reading downstream. + """ + th, tp = _two_channels() + board = _bare_board() + fused = np.array([board._fuse(th[i], tp[i])[0] for i in range(th.size)]) + avg = (th + tp) / 2.0 + w = slice(1000, None) + assert abs(fused[w].mean() - avg[w].mean()) < 0.01, "fusion shifted the calibration" + + +def test_fusion_is_quieter_than_the_average(): + th, tp = _two_channels() + board = _bare_board() + fused = np.array([board._fuse(th[i], tp[i])[0] for i in range(th.size)]) + avg = (th + tp) / 2.0 + w = slice(1000, None) + + def wn(x): + return np.std(np.diff(x)) / np.sqrt(2) + + assert wn(fused[w]) < wn(avg[w]) / 2.0, "fusion did not halve the noise" + + +def test_fusion_survives_one_dead_channel(): + board = _bare_board() + value, var = board._fuse(float("nan"), 29.5) + assert value == 29.5, "a dead HTS221 must not poison the reading" + value, var = board._fuse(30.5, float("nan")) + assert value == 30.5 + value, var = board._fuse(float("nan"), float("nan")) + assert not np.isfinite(value) + + +def test_kalman_rate_is_physical_in_a_still_room(): + """The tuning failure this guards against. + + On a real station the temperature filter reported a median rate of + 12.4 C/h while the room moved 0.37 C/h. Process noise was set to track + perhaps a hundred times faster than any of these signals actually move. + """ + from ashvale.config import CONFIG + from ashvale.estimation import KalmanCV + + dt = CONFIG.sensor.sample_period_s + rng = np.random.default_rng(3) + n = 6000 + truth = 24.0 + 0.4 * np.arange(n) * dt / 3600.0 # a real 0.4 C/h drift + z = truth + 0.0877 * rng.normal(size=n) # measured input noise + + kf = KalmanCV(CONFIG.sensor.kalman_q_temp, CONFIG.sensor.kalman_r_temp) + rates = [kf.update(z[i], dt)[1] * 3600.0 for i in range(n)] + settled = np.abs(np.array(rates[600:])) + + assert np.median(settled) < 3.0, ( + f"median |rate| {np.median(settled):.1f} C/h in a room drifting 0.4 C/h") + assert np.percentile(settled, 95) < 10.0