SciPy lfilter滤波技巧：正确设置初始状态zi

1. 数字滤波器与内部状态

数字滤波器，特别是无限脉冲响应（IIR）滤波器，在处理信号时需要利用过去的输入和输出值来计算当前的输出。为了实现这一点，滤波器内部会维护一个“状态”（state），这个状态包含了计算当前输出所需的历史信息。在SciPy中，这个内部状态由zi参数表示。

当对一个信号进行滤波时，滤波器的初始状态会显著影响其在信号开始部分的输出。如果滤波器从“静止”状态（即所有内部寄存器都为零，表示之前没有任何输入或输出）开始，那么在处理最初的几个样本时，其行为会与已经运行了一段时间并达到某种稳定状态的滤波器不同。

2. scipy.signal.lfilter：一次性滤波与迭代滤波

scipy.signal.lfilter是SciPy中用于应用数字IIR或FIR滤波器的核心函数。它既可以一次性处理整个信号数组，也可以在迭代模式下逐个样本地处理数据，这在实时数据流处理中非常有用。

2.1 一次性滤波（One-shot Filtering）

在一次性滤波模式下，我们通常将整个输入信号传递给lfilter函数。在这种情况下，如果未显式提供zi参数，lfilter默认会假定滤波器处于“初始静止”状态。

import scipy.signal
import numpy as np

# 滤波器参数
fc_bessel = 0.14  # 截止频率 [Hz]
ordre_bessel = 3  # 滤波器阶数
fs = 300          # 采样频率 [Hz]

# 生成示例输入数据
# 为了更好地观察初始状态的影响，我们使用一个从0开始的信号
t = np.arange(0, 10, 1/fs)
input_data = np.sin(2 * np.pi * 0.05 * t) + np.random.randn(len(t)) * 0.1
# 确保第一个值为0，方便对比
input_data[0] = 0.0

# 设计Bessel低通滤波器
b, a = scipy.signal.bessel(ordre_bessel, fc_bessel, 'low', analog=False, output='ba', fs=fs)

# 一次性滤波
filter_once = scipy.signal.lfilter(b, a, input_data)

print(f"一次性滤波的第一个值: {filter_once[0]:.6f}")

2.2 迭代滤波（Iterative Filtering）及其遇到的问题

在实时应用中，数据通常是逐个样本到达的。为了在这种场景下应用滤波器，我们需要在每次处理一个样本后保存并更新滤波器的内部状态zi。用户尝试的迭代滤波代码如下：

# 尝试的迭代滤波版本
# 初始化zi，使用了lfilter_zi
z_wrong = scipy.signal.lfilter_zi(b, a)

filter_iter_wrong = []
z_current_wrong = z_wrong # 每次迭代前复制初始状态
for input_value in input_data:
    # lfilter处理单个值，并返回新的状态
    filtered_value, z_current_wrong = scipy.signal.lfilter(b, a, [input_value], zi=z_current_wrong)
    filter_iter_wrong.append(filtered_value[0])

print(f"错误迭代滤波的第一个值: {filter_iter_wrong[0]:.6f}")

运行上述代码，会发现filter_iter_wrong[0]与filter_once[0]的值存在显著差异（例如，filter_once[0]可能是0，而filter_iter_wrong[0]可能是0.999...）。这表明两种滤波方式的初始行为不一致。

3. 问题根源：zi的初始化差异

造成这种差异的关键在于zi的初始化方式。

• lfilter的默认行为： 当zi参数未提供时（如一次性滤波），lfilter默认假定滤波器处于“初始静止”（initial rest）状态。这意味着它认为在处理当前输入之前，滤波器所有的内部状态寄存器都是零。
• lfilter_zi的行为： scipy.signal.lfilter_zi(b, a)函数用于构造lfilter的初始条件，但它旨在模拟“步进响应稳态”（step response steady-state）。这意味着它计算的zi值，会使得滤波器在接收到一个大的阶跃输入后，立即处于一种稳定的输出状态。这显然与“初始静止”的概念完全不同。当滤波器被初始化为这种“稳态”时，即使输入为零，它也会立即产生一个非零的输出，因为它“认为”之前已经有了一个阶跃输入，现在正在从该输入导致的稳定状态开始。

因此，使用lfilter_zi作为迭代滤波的初始状态，实际上是让滤波器从一个已经“预热”或“稳定”的状态开始，而不是从一个“冷启动”或“静止”的状态开始，这与一次性滤波的默认行为相悖。

4. 解决方案：正确初始化zi

为了使迭代滤波的结果与一次性滤波（假定初始静止）的结果一致，我们需要确保迭代滤波的zi也表示“初始静止”状态。

SciPy提供了scipy.signal.lfiltic(b, a, y0, x0)函数来构造lfilter的初始条件。当我们将y0和x0（分别代表滤波器输出和输入的初始条件）都设置为零时，它就能模拟“初始静止”状态。

# 正确初始化zi
# lfiltic(b, a, y0, x0) - 当y0和x0都为0时，表示初始静止
z_correct = scipy.signal.lfiltic(b, a, 0) # 0表示所有初始输出和输入都为0
# 或者，更简单地，直接使用全零数组，因为对于初始静止，zi就是全零
# z_correct = np.zeros(ordre_bessel) # 对于N阶滤波器，zi的长度是N

有了正确的zi，我们就可以实现与一次性滤波结果一致的迭代滤波：

# 正确的迭代滤波版本
filter_iter_correct = []
z_current_correct = z_correct # 每次迭代前复制初始状态
for input_value in input_data:
    filtered_value, z_current_correct = scipy.signal.lfilter(b, a, [input_value], zi=z_current_correct)
    filter_iter_correct.append(filtered_value[0])

print(f"正确迭代滤波的第一个值: {filter_iter_correct[0]:.6f}")

现在，filter_iter_correct[0]将与filter_once[0]的值非常接近，表明两种滤波方式的初始行为已经一致。

5. 完整示例与验证

下面的代码展示了三种滤波方式的对比，并验证了正确初始化zi后的结果一致性。

import scipy.signal
import numpy as np
import matplotlib.pyplot as plt

# 滤波器参数
fc_bessel = 0.14  # 截止频率 [Hz]
ordre_bessel = 3  # 滤波器阶数
fs = 300          # 采样频率 [Hz]

# 生成示例输入数据 (从0开始的信号，以便观察初始影响)
t = np.arange(0, 5, 1/fs)
input_data = np.sin(2 * np.pi * 0.05 * t) + np.random.randn(len(t)) * 0.1
input_data[0] = 0.0 # 确保第一个值为0

# 设计Bessel低通滤波器
b, a = scipy.signal.bessel(ordre_bessel, fc_bessel, 'low', analog=False, output='ba', fs=fs)

print("--- 滤波结果对比 ---")

# 1. 一次性滤波 (默认初始静止)
filter_once = scipy.signal.lfilter(b, a, input_data)
print(f"一次性滤波的第一个值: {filter_once[0]:.6f}")

# 2. 错误迭代滤波 (使用lfilter_zi初始化)
z_wrong_init = scipy.signal.lfilter_zi(b, a)
filter_iter_wrong = []
z_current_wrong = z_wrong_init
for input_value in input_data:
    filtered_value, z_current_wrong = scipy.signal.lfilter(b, a, [input_value], zi=z_current_wrong)
    filter_iter_wrong.append(filtered_value[0])
print(f"错误迭代滤波的第一个值: {filter_iter_wrong[0]:.6f}")

# 3. 正确迭代滤波 (使用lfiltic或np.zeros初始化为初始静止)
# 方法一：使用lfiltic
z_correct_init_lfiltic = scipy.signal.lfiltic(b, a, 0) # y0=0, x0=0 模拟初始静止
# 方法二：直接使用全零数组 (对于初始静止，zi就是全零)
# z_correct_init_zeros = np.zeros(ordre_bessel) # 滤波器阶数即为zi的长度

filter_iter_correct = []
z_current_correct = z_correct_init_lfiltic # 使用lfiltic初始化的状态
for input_value in input_data:
    filtered_value, z_current_correct = scipy.signal.lfilter(b, a, [input_value], zi=z_current_correct)
    filter_iter_correct.append(filtered_value[0])
print(f"正确迭代滤波的第一个值: {filter_iter_correct[0]:.6f}")

# 验证一致性
print(f"\n一次性滤波与正确迭代滤波是否几乎一致 (np.allclose): {np.allclose(filter_once, filter_iter_correct)}")
print(f"一次性滤波与错误迭代滤波是否几乎一致 (np.allclose): {np.allclose(filter_once, filter_iter_wrong)}")


# 绘制结果进行可视化对比
plt.figure(figsize=(12, 6))
plt.plot(t, input_data, label='原始输入信号', alpha=0.7)
plt.plot(t, filter_once, label='一次性滤波 (参考)', linewidth=2)
plt.plot(t, filter_iter_wrong, label='错误迭代滤波 (lfilter_zi)', linestyle='--', alpha=0.8)
plt.plot(t, filter_iter_correct, label='正确迭代滤波 (lfiltic/zeros)', linestyle=':', linewidth=2)
plt.title('SciPy lfilter 初始状态对滤波结果的影响')
plt.xlabel('时间 [s]')
plt.ylabel('幅值')
plt.legend()
plt.grid(True)
plt.xlim(0, 0.5) # 放大初始部分以便观察差异
plt.show()

# 绘制差异
plt.figure(figsize=(12, 4))
plt.plot(t, filter_iter_wrong - filter_once, label='错误迭代 - 一次性', color='red')
plt.plot(t, filter_iter_correct - filter_once, label='正确迭代 - 一次性', color='green')
plt.title('滤波结果与一次性滤波的差异')
plt.xlabel('时间 [s]')
plt.ylabel('差异')
plt.legend()
plt.grid(True)
plt.xlim(0, 0.5) # 放大初始部分以便观察差异
plt.show()

从输出和图中可以看出，filter_iter_correct的结果与filter_once几乎完全一致，而filter_iter_wrong在信号开始部分存在明显偏差。

6. 注意事项与总结

• zi的选择取决于应用场景：
  • 如果你希望滤波器从一个“冷启动”或“初始静止”状态开始处理数据，例如处理一个全新的传感器数据流，或者每次处理一个独立的数据块，那么应该使用scipy.signal.lfiltic(b, a, 0)或直接使用np.zeros(ordre_bessel)来初始化zi。这是最常见的实时滤波需求。
  • 如果你需要滤波器在处理数据时立即达到某种“稳态”响应，例如在已知信号已经达到稳定状态的情况下继续滤波，或者希望模拟阶跃响应后的稳定输出，那么scipy.signal.lfilter_zi(b, a)可能适用。但这种情况相对较少。
• 实时处理的连续性： 在处理连续的实时数据流时，正确地维护和传递zi至关重要。每次调用lfilter处理一个或一批新样本后，都必须将返回的新zi状态保存起来，并在下一次调用时作为输入参数传入，以确保滤波过程的连续性和准确性。
• 滤波器阶数与zi长度： 对于一个N阶的滤波器，b和a系数数组通常有N+1个元素，而zi数组的长度是N。因此，np.zeros(ordre_bessel)是初始化zi为全零的简洁方法。

总之，在使用scipy.signal.lfilter进行迭代滤波时，理解并正确初始化zi是确保滤波结果准确性和一致性的关键。对于大多数从“初始静止”状态开始的实时滤波场景，使用scipy.signal.lfiltic(b, a, 0)或np.zeros(ordre_bessel)来初始化zi是正确的做法。