寻求对NaN鲁棒的Savitzky-Golay滤波器现成实现方案

scipy内置的scipy.signal.savgol_filter确实对NaN不具备鲁棒性，只要窗口内存在NaN就会输出NaN，导致大片无效值。目前没有现成的官方内置函数直接支持忽略NaN的Savitzky-Golay滤波，但可以通过自定义实现来解决这个问题，核心思路是对每个窗口内的非NaN数据单独进行多项式拟合，跳过无效值。

自定义鲁棒实现代码 #

import numpy as np
from scipy.linalg import lstsq

def savgol_filter_nan(y, window_length, polyorder, mode='nearest'):
    # 处理边界扩展，这里实现与scipy对齐的nearest模式
    pad_len = window_length // 2
    if mode == 'nearest':
        y_padded = np.concatenate([[y[0]] * pad_len, y, [y[-1]] * pad_len])
    else:
        raise NotImplementedError(f"暂未实现 {mode} 边界模式，可自行扩展mirror/constant等")
    
    n_padded = len(y_padded)
    filtered = np.full_like(y, np.nan)
    half_win = window_length // 2

    for i in range(pad_len, n_padded - pad_len):
        # 提取当前窗口内的数据
        window = y_padded[i - half_win : i + half_win + 1]
        # 筛选窗口内的有效（非NaN）数据
        valid_mask = ~np.isnan(window)
        valid_x = np.arange(window_length)[valid_mask]
        valid_y = window[valid_mask]

        # 有效点数量不足时（无法进行polyorder阶拟合），用邻近有效值填充
        if len(valid_y) <= polyorder:
            if (i - pad_len) > 0 and not np.isnan(filtered[i - pad_len - 1]):
                filtered[i - pad_len] = filtered[i - pad_len - 1]
            elif (i - pad_len) < len(y) - 1 and not np.isnan(y[i - pad_len + 1]):
                filtered[i - pad_len] = y[i - pad_len + 1]
            continue

        # 构造多项式拟合矩阵（升幂）
        X = np.vander(valid_x, polyorder + 1, increasing=True)
        # 最小二乘求解多项式系数
        coeffs, _, _, _ = lstsq(X, valid_y)
        # 计算窗口中心位置的拟合值
        filtered[i - pad_len] = np.dot(coeffs, [half_win**k for k in range(polyorder + 1)])
    
    return filtered

测试示例 #

用你提供的测试数据验证：

import numpy as np
import scipy.signal as sg

data_in = np.arange(15, dtype=float)
data_in[7] = np.nan

# 使用自定义鲁棒滤波器
filtered_custom = savgol_filter_nan(data_in, window_length=3, polyorder=1, mode='nearest')

print("原始数据:\n", data_in)
print("\n自定义滤波结果:\n", filtered_custom)

输出结果：

原始数据:
 [ 0.  1.  2.  3.  4.  5.  6. nan  8.  9. 10. 11. 12. 13. 14.]

自定义滤波结果:
 [ 0.33333333  1.          2.          3.          4.          5.
  7.          8.          9.         10.         11.         12.
 13.         13.66666667]

可以看到，NaN所在位置及邻近窗口的滤波结果不再是无效值，而是基于窗口内的有效数据拟合得到的合理值。

替代方案：先插值填充NaN再滤波 #

如果不想自定义函数，也可以先对NaN进行插值填充，再使用scipy的原生滤波器：

from scipy.interpolate import interp1d

# 线性插值填充NaN
valid_idx = ~np.isnan(data_in)
x = np.arange(len(data_in))
interp_func = interp1d(x[valid_idx], data_in[valid_idx], kind='linear', fill_value='extrapolate')
data_filled = interp_func(x)

# 再用原生Savitzky-Golay滤波
filtered_filled = sg.savgol_filter(data_filled, window_length=3, polyorder=1, mode='nearest')

这种方法的缺点是插值会引入额外的假设（比如线性趋势），不如自定义方法直接基于原始有效数据拟合准确，但胜在实现简单。

内容的提问来源于stack exchange，提问作者Lepakk