Scipy滤波处理含NaN数组时返回全NaN值问题排查

问题描述 #

尝试用Scipy的FIR滤波器处理包含NaN值的数组，代码如下，但运行后结果全为NaN：

import numpy as np
import scipy.signal as sp

def apply_filter(x,fs,fc):
    l_filt = 2001
    b = sp.firwin(l_filt, fc, window='blackmanharris', pass_zero='lowpass', fs=fs)

    # zero-phase filter:
    xmean = np.nanmean(x)
    y = sp.filtfilt(b, 1, x - xmean, padlen=9)
    y += xmean
    return y

my_array = [13.049393453879606, 11.710994125276567, 15.39159227893492, 14.053192950331884, np.nan, np.nan, np.nan, np.nan, np.nan, 18.57029068436713, np.nan, np.nan, np.nan, np.nan, 15.893492027161058, 16.228091859311817, 15.558892195010298, np.nan, 8.866895551995118, 14.053192950331882]

tt = apply_filter(my_array,64,30)

错误原因 #

• NaN值传播：Scipy的filtfilt及底层的lfilter不支持直接处理含NaN的输入，计算过程中NaN会逐步扩散，最终导致整个输出数组变为NaN。
• 滤波器长度不合理：设置的滤波器长度l_filt=2001远大于输入数组的长度（仅20个元素），这种极端长度的滤波器会引发计算异常，再加上NaN的存在，直接导致全NaN输出。

解决方案 #

步骤1：预处理填充NaN值 #

先通过插值或有效值填充的方式消除数组中的NaN，确保滤波输入无缺失值。常用方法包括线性插值、前后值填充等。

步骤2：调整滤波器长度 #

滤波器长度需与输入信号长度匹配，一般建议设置为信号长度的1/10~1/5左右，且FIR滤波器通常使用奇数长度。

步骤3：优化滤波参数 #

调整padlen参数为滤波器长度的一半左右，避免边界处理异常。

修改后的代码示例 #

import numpy as np
import scipy.signal as sp

def apply_filter(x, fs, fc):
    x = np.array(x)
    # 用线性插值填充NaN
    mask = ~np.isnan(x)
    if not np.any(mask):
        return x
    
    xp = np.where(mask)[0]
    fp = x[mask]
    x_filled = np.interp(np.arange(len(x)), xp, fp)
    
    # 设置合理的滤波器长度（奇数，不超过信号长度的1/2）
    l_filt = min(21, len(x_filled) // 2)
    if l_filt % 2 == 0:
        l_filt += 1
    
    b = sp.firwin(l_filt, fc, window='blackmanharris', pass_zero='lowpass', fs=fs)

    # 零相位滤波，padlen设为滤波器长度的一半
    xmean = np.mean(x_filled)
    y = sp.filtfilt(b, 1, x_filled - xmean, padlen=l_filt//2)
    y += xmean
    return y

my_array = [13.049393453879606, 11.710994125276567, 15.39159227893492, 14.053192950331884, np.nan, np.nan, np.nan, np.nan, np.nan, 18.57029068436713, np.nan, np.nan, np.nan, np.nan, 15.893492027161058, 16.228091859311817, 15.558892195010298, np.nan, 8.866895551995118, 14.053192950331882]

tt = apply_filter(my_array, 64, 30)
print(tt)

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