Hampel Filter 算法原理与 C 语言实现

Hampel Filter 是一种基于中值绝对偏差（Median Absolute Deviation, MAD）的鲁棒离群点检测算法。与均值滤波、高斯滤波不同，Hampel Filter 对脉冲噪声和尖刺异常值具有天然免疫力，适用于传感器信号预处理场景。

算法原理

滑动窗口与中值

对于输入信号 $x[n]$，取以当前样本为中心、半径为 $k$ 的滑动窗口：

窗口长度为 $2k+1$。计算窗口内所有样本的中值：

窗口 $\{102, 98, 101, 500, 99\}$ 排序后 $\{98, 99, 101, 102, 500\}$，中值 $\tilde{x} = 101$。

中值绝对偏差（MAD）

MAD 定义为窗口内各点与中值偏差的绝对值的中值：

MAD 比标准差更鲁棒——单个异常值不会像对方差那样拉偏整个统计量。

异常判定与替换

将 MAD 乘以缩放系数 $\sigma$（通常取 1.4826，使 MAD 在高斯分布下与标准差一致），设定阈值 $t$（通常取 3）。$\sigma$ 在此语境下是 MAD 的缩放因子，作用是把 MAD 校正到与标准差等价的尺度，使得阈值 $t$ 可以直接沿用统计学中 "$3\sigma$ 准则"的直觉——$t=3$ 对应约 99.7% 的高斯置信区间：

判定规则：

• 若 $|x[n] - \tilde{x}| > t \cdot \sigma_{\text{MAD}}$，则将 $x[n]$ 视为离群点，替换为中值 $\tilde{x}$
• 否则保留原始值 $x[n]$

C 语言实现

以下实现适用于嵌入式环境，无动态内存分配，使用固定大小缓冲区。

#include <stdint.h>
#include <string.h>
#include <stdlib.h>

/* 比较函数，用于 qsort */
static int cmp_int16(const void *a, const void *b)
{
    int16_t va = *(const int16_t *)a;
    int16_t vb = *(const int16_t *)b;
    return (va > vb) - (va < vb);
}

/*
 * Hampel Filter - 基于中值绝对偏差的离群点检测与替换
 *
 * @param input:   输入信号数组
 * @param output:  输出信号数组（可与 input 相同，原地处理）
 * @param len:     信号长度
 * @param k:       窗口半径（窗口长度 = 2k+1）
 * @param t:       阈值系数（推荐 3.0）
 * @param sigma:   MAD 缩放系数（推荐 1.4826）
 */
void hampel_filter(const int16_t *input, int16_t *output,
                   uint32_t len, uint16_t k,
                   float t, float sigma)
{
    const uint16_t win_size = 2 * k + 1;

    /* 临时缓冲区，存放窗口数据的副本 */
    int16_t window[2 * 16 + 1]; /* 最大支持 k=16，按需调整 */
    int16_t abs_dev[2 * 16 + 1];

    for (uint32_t i = 0; i < len; i++) {
        /* 计算窗口边界，处理边界情况 */
        uint32_t start = (i >= k) ? i - k : 0;
        uint32_t end = (i + k < len) ? i + k : len - 1;
        uint16_t actual_size = (uint16_t)(end - start + 1);

        /* 复制窗口数据 */
        for (uint16_t j = 0; j < actual_size; j++) {
            window[j] = input[start + j];
        }

        /* 计算中值 */
        qsort(window, actual_size, sizeof(int16_t), cmp_int16);
        int16_t median = window[actual_size / 2];

        /* 计算各点与中值的绝对偏差 */
        for (uint16_t j = 0; j < actual_size; j++) {
            abs_dev[j] = (window[j] > median) ?
                         (window[j] - median) :
                         (median - window[j]);
        }

        /* 计算 MAD */
        qsort(abs_dev, actual_size, sizeof(int16_t), cmp_int16);
        int16_t mad = abs_dev[actual_size / 2];

        /* 异常判定 */
        float threshold = t * sigma * (float)mad;
        int16_t deviation = (input[i] > median) ?
                            (input[i] - median) :
                            (median - input[i]);

        if ((float)deviation > threshold) {
            output[i] = median;  /* 离群点，替换为中值 */
        } else {
            output[i] = input[i]; /* 正常点，保留原值 */
        }
    }
}

参数选择

使用示例

#include <stdio.h>

int main(void)
{
    int16_t signal[] = {100, 102, 98, 101, 500, 99, 103, 100, -200, 98};
    uint32_t len = sizeof(signal) / sizeof(signal[0]);
    int16_t filtered[len];

    hampel_filter(signal, filtered, len,
                  /*k=*/3, /*t=*/3.0f, /*sigma=*/1.4826f);

    printf("原始信号\t滤波后\n");
    for (uint32_t i = 0; i < len; i++) {
        printf("%d\t\t%d\n", signal[i], filtered[i]);
    }
    return 0;
}

输出：

原始信号        滤波后
100             100
102             102
98              98
101             101
500             101      ← 尖刺被替换为中值
99              99
103             103
100             100
-200            99       ← 尖刺被替换为中值
98              98

嵌入式应用注意事项

1. qsort 开销：每次窗口移动执行两次排序。若窗口较小（$k \leq 5$，窗口 $\leq 11$），排序开销可接受。对于实时性要求高的场景，可用插入排序替代。

2. 整数溢出：MAD 计算中的绝对偏差使用 int16_t，若信号动态范围接近 INT16_MAX，需改用 int32_t。

3. 边界处理：代码中边界处窗口自动缩小。若需要严格等长窗口，可在信号两端做镜像填充或重复边界值。

4. 浮点依赖：阈值比较使用 float。若 MCU 无 FPU，可将 $t \cdot \sigma$ 合并为常数，乘法改用定点运算。

适用场景

• 传感器数据预处理（温度、压力、心率等突变尖刺）
• 工业控制系统中的信号清洗
• 嵌入式系统中替代中值滤波器，增加自适应阈值能力
• 任何需要保留信号边缘同时剔除脉冲干扰的场景

本文由 AI 辅助生成，可能存在错误或遗漏，请以实际资料和官方文档为准。