Add Hampel outlier filter for noisy powermeter sources

CHANGELOG.md

@@ -1,6 +1,7 @@
 # Changelog
 
 ## Next
+- Added optional Hampel outlier filter (`HAMPEL_WINDOW`, `HAMPEL_N_SIGMA`, `HAMPEL_MIN_THRESHOLD`) for rejecting wild samples from noisy powermeter sources (MQTT/HTTP/WiFi); configurable globally or per powermeter section, disabled by default
 - Added `DEDUPE_TIME_WINDOW` support to the Shelly emulator to drop burst-repeat requests from the same battery IP; the value can be set under `[GENERAL]` to apply regardless of which device type is emulated
 - Added opt-in web-based configuration editor (`WEB_CONFIG_ENABLED = True` in `[GENERAL]`) accessible at `http://<host>:52500/config`; supports editing all config sections and keys with type-appropriate inputs, comment preservation, and a Save & Restart button
 - **Breaking:** Rebrand project from "B2500 Meter" to "AstraMeter" (formerly b2500-meter). Package renamed to `astrameter`, CLI commands are now `astrameter` and `astra-sim`. Docker image moved from `ghcr.io/tomquist/b2500-meter` to `ghcr.io/tomquist/astrameter` (the legacy `ghcr.io/tomquist/b2500-meter` image is still published in parallel for backward compatibility). Home Assistant users must update their app repository URL to `https://github.com/tomquist/astrameter#main`.

README.md

@@ -202,6 +202,15 @@ Per-powermeter options (apply in any powermeter section, e.g. `[TASMOTA]` or `[H
 - **DEADBAND** (default 0 = disabled, W) — When the absolute reading is below this
   value, the wrapper emits zeros instead of chasing noise. Keeps batteries from
   hunting around the zero-crossing; 10–30 W is a sensible range.
+- **HAMPEL_WINDOW** (default 0 = disabled) — Rolling window size for
+  median-based outlier rejection. Typical values 5–7. Useful for MQTT/HTTP
+  sources that occasionally emit wild samples; applied after throttling and
+  before EMA smoothing.
+- **HAMPEL_N_SIGMA** (default 3.0) — Rejection threshold in MAD-derived sigmas.
+  Lower values reject more aggressively.
+- **HAMPEL_MIN_THRESHOLD** (default 0, W) — Minimum rejection threshold in
+  watts. Prevents spikes from passing through during long periods of constant
+  readings (the MAD=0 degenerate case); 50 W is a reasonable starting value.
 
 CT002/CT003 active-steering options (all under `[CT002]` or `[CT003]`):
 - **ACTIVE_CONTROL** — When true (default), the emulator smooths the grid reading, splits

config.ini.example

@@ -357,6 +357,26 @@ THROTTLE_INTERVAL = 0
 #PID_OUTPUT_MAX = 800
 #PID_MODE = bias
 
+## --- Hampel outlier filter (optional) ---
+## Stateful rolling-median rejection of wild samples (e.g. MQTT/WiFi glitches).
+## Outliers are replaced with the window median; normal samples pass through
+## unchanged.  Disabled by default.
+##
+## When HAMPEL_WINDOW > 0, the wrapper is inserted after throttling and before
+## EMA smoothing.  This is orthogonal to the EMA: Hampel rejects impulse noise,
+## the EMA smooths the cleaned signal.
+##
+## Recommended starting point for flaky HTTP/MQTT sources:
+##   HAMPEL_WINDOW        = 5    # odd sizes recommended; 5–7 works well
+##   HAMPEL_N_SIGMA       = 3.0  # reject beyond ~3σ (MAD-derived)
+##   HAMPEL_MIN_THRESHOLD = 50   # watts; floor for constant-signal (MAD=0) case
+##
+## Parameters can be set globally in [GENERAL] or per powermeter section.
+## Per-section values override the global ones.
+#HAMPEL_WINDOW = 5
+#HAMPEL_N_SIGMA = 3.0
+#HAMPEL_MIN_THRESHOLD = 50
+
 ## --- MQTT Insights (optional) ---
 ## Publishes internal state (grid power, targets, saturation, consumer topology)
 ## to MQTT with Home Assistant Device Discovery support.

src/astrameter/config/config_loader.py

@@ -37,6 +37,7 @@
     VZLogger,
     parse_sml_obis_config,
 )
+from astrameter.powermeter.wrappers.hampel import HampelPowermeter
 from astrameter.powermeter.wrappers.smoothing import (
     DeadbandPowermeter,
     SmoothedPowermeter,
@@ -162,6 +163,11 @@ def read_all_powermeter_configs(
     )
     global_max_smooth_step = config.getfloat("GENERAL", "MAX_SMOOTH_STEP", fallback=0.0)
     global_deadband = config.getfloat("GENERAL", "DEADBAND", fallback=0.0)
+    global_hampel_window = config.getint("GENERAL", "HAMPEL_WINDOW", fallback=0)
+    global_hampel_n_sigma = config.getfloat("GENERAL", "HAMPEL_N_SIGMA", fallback=3.0)
+    global_hampel_min_threshold = config.getfloat(
+        "GENERAL", "HAMPEL_MIN_THRESHOLD", fallback=0.0
+    )
     global_pid_kp = config.getfloat("GENERAL", "PID_KP", fallback=0.0)
     global_pid_ki = config.getfloat("GENERAL", "PID_KI", fallback=0.0)
     global_pid_kd = config.getfloat("GENERAL", "PID_KD", fallback=0.0)
@@ -208,6 +214,38 @@ def read_all_powermeter_configs(
                 )
                 powermeter = ThrottledPowermeter(powermeter, section_throttle_interval)
 
+            section_hampel_window = config.getint(
+                section, "HAMPEL_WINDOW", fallback=global_hampel_window
+            )
+            if section_hampel_window > 0:
+                section_hampel_n_sigma = config.getfloat(
+                    section, "HAMPEL_N_SIGMA", fallback=global_hampel_n_sigma
+                )
+                section_hampel_min_threshold = config.getfloat(
+                    section,
+                    "HAMPEL_MIN_THRESHOLD",
+                    fallback=global_hampel_min_threshold,
+                )
+                hampel_source = (
+                    "section-specific"
+                    if config.has_option(section, "HAMPEL_WINDOW")
+                    else "global"
+                )
+                logger.info(
+                    "Applying %s Hampel outlier filter (window=%d, n_sigma=%.2f, min_threshold=%.0fW) to %s",
+                    hampel_source,
+                    section_hampel_window,
+                    section_hampel_n_sigma,
+                    section_hampel_min_threshold,
+                    section,
+                )
+                powermeter = HampelPowermeter(
+                    powermeter,
+                    window=section_hampel_window,
+                    n_sigma=section_hampel_n_sigma,
+                    min_threshold=section_hampel_min_threshold,
+                )
+
             section_smooth_alpha = config.getfloat(
                 section, "SMOOTH_TARGET_ALPHA", fallback=global_smooth_alpha
             )

src/astrameter/powermeter/__init__.py

@@ -18,6 +18,7 @@
 from .vzlogger import VZLogger
 from .wrappers import (
     DeadbandPowermeter,
+    HampelPowermeter,
     PidPowermeter,
     PowermeterWrapper,
     SmoothedPowermeter,
@@ -30,6 +31,7 @@
     "DeadbandPowermeter",
     "ESPHome",
     "Emlog",
+    "HampelPowermeter",
     "HomeAssistant",
     "HomeWizardPowermeter",
     "IoBroker",

src/astrameter/powermeter/wrappers/__init__.py

@@ -1,11 +1,13 @@
 from .base import PowermeterWrapper
+from .hampel import HampelPowermeter
 from .pid import PidPowermeter
 from .smoothing import DeadbandPowermeter, SmoothedPowermeter
 from .throttling import ThrottledPowermeter
 from .transform import TransformedPowermeter
 
 __all__ = [
     "DeadbandPowermeter",
+    "HampelPowermeter",
     "PidPowermeter",
     "PowermeterWrapper",
     "SmoothedPowermeter",

src/astrameter/powermeter/wrappers/hampel.py

@@ -0,0 +1,88 @@
+"""Hampel outlier-rejection powermeter wrapper."""
+
+from __future__ import annotations
+
+import statistics
+from collections import deque
+
+from astrameter.config.logger import logger
+from astrameter.powermeter.base import Powermeter
+
+from .base import PowermeterWrapper
+
+
+class HampelPowermeter(PowermeterWrapper):
+    """Rolling-median outlier filter for sum-of-phases power readings.
+
+    Maintains a rolling window of the most recent ``window`` totals. When the
+    next total lies more than ``n_sigma * 1.4826 * MAD`` away from the window
+    median (with a floor of ``min_threshold`` watts to handle the constant-
+    signal MAD=0 degenerate case), the sample is treated as an outlier: the
+    reported total is replaced by the median and per-phase values are
+    redistributed proportionally (equal split when ``|raw_total|`` is near
+    zero). The window entry itself is mutated to the median so a single spike
+    does not poison future detections — this is the canonical Hampel
+    identifier formulation used in control literature.
+
+    Operates on the sum of phases, mirroring :class:`SmoothedPowermeter`.
+    A phase-cancelling outlier (e.g. +1000 W on L1 and -1000 W on L2) is
+    therefore invisible to this filter; that is acceptable because every
+    downstream wrapper (EMA, deadband, PID) also operates on sum-of-phases.
+    """
+
+    MAD_SCALE = 1.4826
+
+    def __init__(
+        self,
+        wrapped_powermeter: Powermeter,
+        window: int,
+        n_sigma: float = 3.0,
+        min_threshold: float = 0.0,
+    ) -> None:
+        if window < 1:
+            raise ValueError(f"Hampel window must be >= 1, got {window}")
+        if n_sigma < 0:
+            raise ValueError(f"Hampel n_sigma must be >= 0, got {n_sigma}")
+        if min_threshold < 0:
+            raise ValueError(f"Hampel min_threshold must be >= 0, got {min_threshold}")
+        super().__init__(wrapped_powermeter)
+        self._window: deque[float] = deque(maxlen=window)
+        self._window_size = window
+        self._n_sigma = n_sigma
+        self._min_threshold = min_threshold
+
+    def reset(self) -> None:
+        super().reset()
+        logger.debug("HampelPowermeter: reset (window size=%d)", len(self._window))
+        self._window.clear()
+
+    async def get_powermeter_watts(self) -> list[float]:
+        raw_values = await self.wrapped_powermeter.get_powermeter_watts()
+        if not raw_values:
+            return []
+
+        raw_total = sum(raw_values)
+        self._window.append(raw_total)
+
+        if len(self._window) < self._window_size:
+            return list(raw_values)
+
+        median = statistics.median(self._window)
+        mad = statistics.median(abs(x - median) for x in self._window)
+        threshold = max(self._n_sigma * self.MAD_SCALE * mad, self._min_threshold)
+
+        if threshold <= 0 or abs(raw_total - median) <= threshold:
+            return list(raw_values)
+
+        self._window[-1] = median
+        logger.debug(
+            "HampelPowermeter: outlier rejected raw=%.2f median=%.2f threshold=%.2f",
+            raw_total,
+            median,
+            threshold,
+        )
+
+        if abs(raw_total) < 1e-9:
+            return [median / len(raw_values)] * len(raw_values)
+        ratio = median / raw_total
+        return [v * ratio for v in raw_values]