传感器读数偶尔跳一下，监控曲线锯齿很多，或者一段趋势被单点尖峰打断，都会让人想到“加个滤波”。但平滑后的曲线更顺眼，不代表它离真实状态更近：突变可能被拖慢，真正的异常也可能被抹掉。

本文放在数据分析分类，承接 NumPy 常用统计与数据可视化。先比较四种常见方法，再手算一维卡尔曼的预测与更新。所有示例都只使用当前及过去数据，不把离线看过未来的结果冒充实时滤波。

先问任务

去尖峰、看趋势，还是估计一个隐藏状态？

再写约定

窗口方向、初始值、采样间隔与缺失规则。

最后评估

误差、突变响应和延迟必须一起观察。

同一合成信号含一个8的尖峰和从0到4的阶跃，比较三点移动平均、指数平滑、中值及一维卡尔曼的计算输出
去掉尖峰与及时跟踪突变，不是同一目标。 曲线由下方代码的九个合成观测计算，虚线是人为指定的参考状态；不是传感器实测，也不是算法排名。查看原图

进阶可继续阅读位置—速度卡尔曼与协方差，把标量预测扩展为联合状态，检查位置观测怎样通过交叉协方差修正速度，以及变化采样间隔下的过程噪声。

先区分四类方法

设观测序列为 y_t，希望输出更适合当前任务的估计 s_t。滤波器不能凭空知道哪些变化来自噪声，哪些变化是真实信号。

方法 利用的信息 常见用途 主要代价
移动平均 最近 W 个观测等权汇总 看局部均值、抑制独立零均值波动 滞后，尖峰会影响多个输出
指数平滑 当前观测与上一次估计 低状态成本的在线平滑 参数影响记忆长度与响应速度
中值滤波 窗口内排序后取中间值 对付孤立的极端读数 非线性，也可能消除真实短脉冲
卡尔曼滤波 动态模型、观测与误差方差 递归估计隐藏状态 依赖模型和噪声假设，普通版本不抗任意异常值

线性差分滤波可以用 FIR、IIR 的语言描述；SciPy lfilter 官方说明给出了输入输出的差分形式。中值不是线性滤波，不能把同一个线性频率响应直接套给它。本文只实现一维标量例子，不覆盖图像二维滤波和完整频域设计。

因果移动平均：窗口越大，响应越慢

满窗口时，尾随移动平均定义为：

1
s_t = (y_t + y_(t−1) + ... + y_(t−W+1)) / W

若窗口内真实状态恒为 μ，观测噪声独立、均值为零且方差同为 σ²，则平均仍无偏，噪声方差降为 σ²/W。这里的条件很重要：窗口内真实状态改变时，把不同状态平均在一起会引入偏差；噪声相关时，均值方差还要加入协方差项。

对于无噪声线性斜坡 y_t=a t+b，满窗口输出恰为 a(t−(W−1)/2)+b。这说明它相对斜坡落后了 (W−1)/2 个采样间隔，并不是所有信号都可以用一个固定时间移位解释。

本文在序列刚开始、不足 W 项时只平均已有数据，不做零填充。因此启动阶段的分母会变化，不能把整段实现都直接当作固定系数、时不变的 W 点滤波器。

指数平滑：把历史压成一个状态

取 0<α≤1：

1
2
s_0 = y_0
s_t = α y_t + (1−α) s_(t−1)

展开递推，较旧观测的权重按几何级数衰减；初始值的残留权重为 (1−α)^t。α 大时更跟随当前观测，α 小时历史记忆更长。

若真实状态恒定、观测噪声独立同方差，且忽略已衰减的启动影响，稳态噪声方差为：

1
σ² α² Σ_(j≥0) (1−α)^(2j) = σ² α/(2−α)

这同样是特定噪声模型下的结论，不是所有时间序列上的效果保证。α=1 时原样输出；本文不接受 α=0，以免把“不再响应新数据”混入默认接口。

中值滤波：去尖峰也可能去掉真实短事件

对最近 W 项排序后取中位数。本文使用奇数 W；序列启动时仍使用已有项，若已有项为偶数，取中间两值的平均。

三点窗口 [0,8,0] 的均值为 8/3，中值为 0；[0,4,4] 的中值为 4。这解释了它为何能抑制一个孤立尖峰，但“尖峰是噪声”必须来自任务判断。如果真实事件只持续一个采样点，同样可能被去掉。

窗口方向也是定义的一部分。本文是尾随窗口；SciPy medfilt还明确规定了零填充边界，不能不对齐方向和边界就逐项比较输出。离线居中窗口可能用到未来样本。

同一组观测怎样手算

使用观测 [0,0,8,0,0,4,4,4,4]，人为参考状态为 [0,0,0,0,0,4,4,4,4]。8 是注入的单点干扰，随后是持续阶跃；这只是用于解释行为的合成构造。

t 观测 参考状态 三点均值 指数平滑 α=0.5 尾随三点中值
0—1 0 0 0 0 0
2 8 0 8/3 4 0
3 0 0 8/3 2 0
4 0 0 8/3 1 0
5 4 4 4/3 2.5 0
6 4 4 8/3 3.25 4
7 4 4 4 3.625 4
8 4 4 4 3.8125 4

中值去掉了 t=2 的尖峰，却直到 t=6 才跟上阶跃；均值在 t=5、t=6 分别为 4/3、8/3，也展示了响应滞后。评估应查看完整轨迹，不能只挑一个看起来正确的点。

一维卡尔曼：先写模型，再讨论增益

采用等间隔随机游走和直接观测：

1
2
x_t = x_(t−1) + w_t       Var(w_t)=Q
y_t = x_t + v_t           Var(v_t)=R

假定过程噪声、观测噪声与初始误差零均值，且在需要的时间和相互关系上不相关；Q≥0，R>0，初始误差方差 P≥0。若初始状态及噪声进一步满足线性高斯条件，递推给出相应的后验均值和方差；仅有二阶矩条件时，不应自动宣称得到了任意分布下的完整后验。

状态空间递归估计的原始依据见 Kalman 1960 论文。本文只推导其标量特例，Q、R 都是方差，单位为观测单位的平方；Q 对应一个采样间隔，不能不加调整就移给另一种采样频率。

1
2
3
4
5
预测均值：m⁻ = m
预测方差：P⁻ = P + Q
增益：K = P⁻ / (P⁻ + R)
更新均值：m = m⁻ + K(y−m⁻)
更新方差：P = (1−K)P⁻

为什么是这个 K？用线性修正 m=m⁻+K(y−m⁻)，在预测误差与观测噪声不相关的条件下，修正误差方差是 (1−K)²P⁻+K²R。对 K 求导并令其为零，得到上式；二阶导数 2(P⁻+R)>0，因而是唯一最小值。

P⁻ 大意味着预测更不确定，K 较大；R 大意味着观测更不可靠，K 较小。这不是“测量一定比模型可信”或相反。代码用等价的标量 Joseph 形式 (1−K)²P⁻+K²R 更新方差，避免直接相减带来的部分数值问题。

三步卡尔曼怎样核对

在第一条观测之前设置 m=0、P=1，Q=1、R=1，每个时间步都先预测，第一步也先加 Q。依次收到 [1,2,3]：

步骤 观测 预测均值 P⁻ K 更新均值 更新方差
1 1 0 2 2/3 2/3 2/3
2 2 2/3 5/3 5/8 3/2 5/8
3 3 3/2 13/8 13/21 17/7 13/21

若某时刻没有观测，仍做 P←P+Q，但不做观测修正。把缺失值填成 0 再更新，相当于凭空告诉滤波器“确实测到了零”，含义完全不同。

完整标准库实现

均值和中值接口要求完整有限观测；卡尔曼额外接受 None 作为缺失标记。NaN、无穷和布尔值均不作为普通数值吞入。所有方法返回的列表用于展示；在线使用时可只保留相应的窗口或状态。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
import json
import math
import statistics
from collections import deque


def real(value, name):
    if isinstance(value, bool) or not isinstance(value, (int, float)):
        raise ValueError(name + " must be real")
    try:
        result = float(value)
    except OverflowError as error:
        raise ValueError(name + " is too large") from error
    if not math.isfinite(result):
        raise ValueError(name + " must be finite")
    return result


def series(values, missing=False):
    if not isinstance(values, (list, tuple)) or not values:
        raise ValueError("a nonempty series is required")
    return [None if missing and v is None else real(v, "observation") for v in values]


def window_size(w, odd=False):
    if type(w) is not int or w < 1 or (odd and w % 2 == 0):
        raise ValueError("window must be positive, and odd for median")
    return w


def checked_sum(values):
    try:
        return real(math.fsum(values), "sum")
    except OverflowError as error:
        raise ValueError("intermediate sum overflow") from error


def moving_mean(values, window=3):
    values, window = series(values), window_size(window)
    history, out, total = deque(), [], 0.0
    for value in values:
        removed = history.popleft() if len(history) == window else 0.0
        history.append(value)
        total = checked_sum((total, value, -removed))
        out.append(real(total/len(history), "mean"))
    return out


def exponential(values, alpha=0.5):
    values, alpha = series(values), real(alpha, "alpha")
    if not 0 < alpha <= 1:
        raise ValueError("require 0 < alpha <= 1")
    state, out = values[0], [values[0]]
    for value in values[1:]:
        state = checked_sum((alpha*value, (1-alpha)*state))
        out.append(state)
    return out


def moving_median(values, window=3):
    values, window = series(values), window_size(window, odd=True)
    history, out = deque(), []
    for value in values:
        if len(history) == window:
            history.popleft()
        history.append(value)
        out.append(real(statistics.median(history), "median"))
    return out


def kalman(values, q=0.25, r=1.0, mean=0.0, variance=1.0):
    values = series(values, missing=True)
    q, r = real(q, "Q"), real(r, "R")
    mean, variance = real(mean, "mean"), real(variance, "variance")
    if q < 0 or r <= 0 or variance < 0:
        raise ValueError("require Q>=0, R>0 and variance>=0")
    out = []
    for value in values:
        prediction = mean
        prior = checked_sum((variance, q))
        gain = 0.0
        mean, variance = prediction, prior
        if value is not None:
            gain = prior/checked_sum((prior, r))
            mean = checked_sum(((1-gain)*prediction, gain*value))
            variance = checked_sum(((1-gain)**2*prior, gain*gain*r))
        out.append({"observation": value, "prediction": prediction,
                    "prior_variance": prior, "gain": gain,
                    "estimate": mean, "variance": variance})
    return out


if __name__ == "__main__":
    observed = [0, 0, 8, 0, 0, 4, 4, 4, 4]
    print(json.dumps({"observed": observed,
                      "reference": [0, 0, 0, 0, 0, 4, 4, 4, 4],
                      "mean": moving_mean(observed),
                      "ema": exponential(observed),
                      "median": moving_median(observed),
                      "kalman": kalman(observed),
                      "hand": kalman([1, 2, 3], q=1, r=1),
                      "missing": kalman([1, None, 3], q=1, r=1)}))

滑动和实现的均值为 O(N) 时间、O(W) 工作空间；浮点滑动和仍可能随长序列累积舍入误差，必要时定期重算或使用更稳定的流式统计。指数平滑与标量卡尔曼均为 O(N) 时间、O(1) 状态。这里的中值每步排序，时间 O(N W log W)、工作空间 O(W)，不是高性能在线中位数数据结构。输入副本和返回列表另占 O(N) 空间。

验证程序 tests/verify-filtering-article.cjs 以独立的有理数窗口与显式历史权重核对 600 组序列的 23,489 个位置；对 220 组带缺失的随机模型，用批量高斯二次型消元核对 985 次递归后验均值和方差，并检查因果前缀、平移、Q=0、确定初态及 23 类非法输入。图中折线与本文实际计算输出一致；这些检查不证明某个真实设备的噪声满足模型假设。

验证程序 tests/verify-filtering-article.cjs 以独立的有理数窗口与显式历史权重核对 600 组序列的 23,489 个位置；对 220 组带缺失的随机模型，用批量高斯二次型消元核对 985 次递归后验均值和方差，并检查因果前缀、平移、Q=0、确定初态及 23 类非法输入。图中折线与本文实际计算输出一致；这些检查不证明某个真实设备的噪声满足模型假设。

选用与失败边界

情形 应检查什么
需要真实在线报警 当前输出有没有读到未来；平滑后报警是否延后
看到孤立尖峰 是传感器毛刺还是要检测的真实短事件
卡尔曼跟踪迟缓 Q、R 的大小与单位是否合理，状态是否真能用随机游走描述
卡尔曼被异常值拉走 普通高斯/二阶矩模型不自动具有鲁棒性
输入有缺失或不等间隔 缺失要跳过观测更新，不同时间间隔需要对应的状态和 Q
与某库的结果不同 先对齐窗口方向、填充方式、初始化和缺失规则
只看平滑后的图 同时看原始数据、真实事件和独立验证指标，不能只评“好看”

SciPy filtfilt采用前后双向处理并处理边界，不能当作不看未来的实时实现。本文的四个实现也不是全部滤波算法：频率选择任务还要分析采样率、目标频带和滤波器设计，状态模型复杂时则要扩展状态和观测方程。

自检：缺失观测为什么不能简单填成 0？

缺失表示没有新测量信息，随机游走模型只推进预测及不确定性；零则是一条有明确数值、会影响增益修正的观测。两者在概率模型中不是同一事件。