上一篇指出,移动平均和指数平滑都只能做单变量滤波,无法利用信号和噪声的统计特性来设计全局最优滤波器。维纳滤波正是从这个方向突破:它在频域中根据信号和噪声的功率谱密度来设计传递函数,使输出均方误差最小。这是历史上第一个在统计意义上最优的线性滤波器。
与前两种方法的对比
移动平均和指数平滑的滤波器设计完全靠经验调参(窗口宽度 $W$ 或平滑因子 $\alpha$),没有任何"最优"保证。维纳滤波则不同——它需要知道信号和噪声的功率谱密度,然后直接算出使均方误差最小的传递函数。代价是:它要求信号和噪声都是平稳随机过程,且需要离线处理整段数据。
问题设定
给定含噪观测序列 $z(t) = s(t) + n(t)$,其中 $s(t)$ 是信号,$n(t)$ 是加性噪声,二者均为平稳随机过程且互相独立。维纳滤波的任务是设计一个线性时不变滤波器,其脉冲响应为 $h(t)$,使得滤波输出 $\hat{s}(t) = h(t) * z(t)$ 与真实信号 $s(t)$ 之间的均方误差最小:
$$\min_h \; E\left[ \left( s(t) - \int_0^\infty h(\tau) z(t - \tau) \, d\tau \right)^2 \right]$$
由正交性原理,最优滤波器的脉冲响应满足 Wiener–Hopf 积分方程:
$$\int_0^\infty h(\tau) \, R_{zz}(t - \tau) \, d\tau = R_{sz}(t), \quad t \geq 0$$
其中 $R_{zz}(\tau) = E[z(t) z(t+\tau)]$ 是观测的自相关函数,$R_{sz}(\tau) = E[s(t) z(t+\tau)]$ 是信号与观测的互相关函数。
频域解
Wiener–Hopf 方程在时域中不易直接求解。转换到频域后,最优滤波器的传递函数具有极其简洁优美的形式。设 $S_{ss}(\omega)$ 为信号的功率谱密度,$S_{nn}(\omega)$ 为噪声的功率谱密度,且信号与噪声不相关,那么非因果维纳滤波器的传递函数为:
$$H(\omega) = \frac{S_{ss}(\omega)}{S_{ss}(\omega) + S_{nn}(\omega)}$$
这个公式揭示了一种深刻的直觉。在信噪比高的频率分量上($S_{ss} \gg S_{nn}$),$H(\omega) \approx 1$,滤波器让信号几乎无损通过;在信噪比低的频率分量上($S_{ss} \ll S_{nn}$),$H(\omega) \approx 0$,滤波器几乎完全抑制该分量。维纳滤波器在每一个频率上独立地做信噪比加权。
因果维纳滤波器(只能使用过去和当前的观测)的表达式稍复杂,涉及频谱分解——将 $\frac{S_{ss}(\omega)}{S_{ss}(\omega) + S_{nn}(\omega)}$ 分解为因果部分和反因果部分,但核心思想完全相同。
可视化
下图的 GIF 展示了维纳滤波的工作过程:在频域中,根据信号功率谱和噪声功率谱计算出最优传递函数 $H(\omega)$,然后在频域应用该传递函数,再逆变换回时域。

Python 实现
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
os.makedirs('images', exist_ok=True)
np.random.seed(42)
N = 512
k = np.arange(N)
true = np.sin(2 * np.pi * k / 50) + 0.5 * np.sin(2 * np.pi * k / 8)
noise = np.random.randn(N) * 0.8
obs = true + noise
Z = np.fft.fft(obs)
S = np.fft.fft(true)
Nn = np.fft.fft(noise)
freqs = np.fft.fftfreq(N)
S_ss = np.abs(S)**2 / N
S_nn = np.abs(Nn)**2 / N
S_ss_smooth = np.convolve(S_ss, np.ones(10)/10, mode='same')
S_nn_smooth = np.convolve(S_nn, np.ones(10)/10, mode='same')
H = S_ss_smooth / (S_ss_smooth + S_nn_smooth + 1e-12)
Z_filtered = Z * H
filtered = np.real(np.fft.ifft(Z_filtered))
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
axes[0].plot(k, true, 'k-', label='True Signal', linewidth=1.2)
axes[0].plot(k, obs, '.', color='gray', alpha=0.3, label='Observations')
axes[0].plot(k, filtered, 'r-', label='Wiener Filtered', linewidth=1.2)
axes[0].set_xlabel('Time Step k')
axes[0].set_ylabel('Value')
axes[0].set_title('Time Domain')
axes[0].legend(loc='upper right')
axes[0].grid(True, alpha=0.3)
pos = freqs >= 0
axes[1].plot(freqs[pos]*N, S_ss_smooth[pos], label='Signal PSD', linewidth=1.2)
axes[1].plot(freqs[pos]*N, S_nn_smooth[pos], label='Noise PSD', linewidth=1.2)
axes[1].plot(freqs[pos]*N, H[pos], label='Transfer Function H(omega)', linewidth=1.5)
axes[1].set_xlabel('Frequency Bin')
axes[1].set_ylabel('Magnitude')
axes[1].set_title('Frequency Domain')
axes[1].legend(loc='upper right')
axes[1].grid(True, alpha=0.3)
fig.tight_layout()
fig.savefig('images/wiener_filter.svg')
fig.savefig('images/wiener_filter.png', dpi=150)
plt.close(fig)
fig, ax = plt.subplots(figsize=(10, 5))
line_true, = ax.plot(k, true, 'k-', label='True Signal', linewidth=1.2)
line_obs, = ax.plot(k, obs, '.', color='gray', alpha=0.2, label='Observations')
line_filt, = ax.plot(k, filtered, 'r-', label='Wiener Filtered', linewidth=1.2)
title = ax.set_title('Wiener Filter - Time Domain')
ax.set_xlabel('Time Step k')
ax.set_ylabel('Value')
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_ylim(-3.5, 3.5)
modes = ['time'] * 6 + ['freq'] * 6 + ['time'] * 6
def update(i):
mode = modes[i]
if mode == 'time':
ax.clear()
ax.plot(k, true, 'k-', label='True Signal', linewidth=1.2)
ax.plot(k, obs, '.', color='gray', alpha=0.2, label='Observations')
ax.plot(k, filtered, 'r-', label='Wiener Filtered', linewidth=1.2)
ax.set_xlabel('Time Step k')
ax.set_ylabel('Value')
ax.set_title('Wiener Filter - Time Domain Reconstruction')
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_ylim(-3.5, 3.5)
else:
ax.clear()
pos = freqs >= 0
ax.plot(freqs[pos]*N, S_ss_smooth[pos], label='Signal PSD', linewidth=1.2)
ax.plot(freqs[pos]*N, S_nn_smooth[pos], label='Noise PSD', linewidth=1.2)
ax.plot(freqs[pos]*N, H[pos], label='H(omega)', linewidth=1.5)
ax.set_xlabel('Frequency Bin')
ax.set_ylabel('Magnitude')
ax.set_title('Wiener Filter - Frequency Domain Analysis')
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_ylim(0, max(S_ss_smooth.max(), S_nn_smooth.max()) * 1.2)
return []
ani = animation.FuncAnimation(fig, update, frames=len(modes), interval=600, blit=False)
writer = animation.PillowWriter(fps=1.5)
ani.save('images/滤波03-图1.gif', writer=writer)
plt.close(fig)
import csv
with open('images/wiener_filter.csv', 'w', newline='') as f:
w = csv.writer(f)
w.writerow(['k', 'true', 'obs', 'filtered'])
for i in range(N):
w.writerow([i, true[i], obs[i], filtered[i]])
演化与局限
维纳滤波相比移动平均和指数平滑有了质的飞跃:它第一次引入了统计最优性——在已知信号和噪声功率谱的前提下,给出均方误差最小的线性滤波器。但它的三个致命局限阻碍了实际应用:要求信号和噪声平稳、需要离线批量处理、无法处理多维耦合状态。
1960 年,鲁道夫·卡尔曼提出了一个完全不同的框架:用状态空间模型描述系统内部动力学,在时域中以递推方式实现最优估计。这个框架不要求平稳性,可以实时处理,天然支持多维状态——这就是卡尔曼滤波,下一篇的主题。
参考文献
[1] arXiv:2112.03533 — A Time-domain Real-valued Generalized Wiener Filter for Multi-channel Neural Separation.
[2] arXiv:2606.07876 — Optimal Wiener-Filter Solutions for Denoising of Graph Signals on Directed Graphs.
[3] arXiv:2410.21184 — Shannon-like Interpolation with Spectral Priors and Weighted Hilbert Spaces.
CycleUser