滤波方法(三):维纳滤波

上一篇指出,移动平均和指数平滑都只能做单变量滤波,无法利用信号和噪声的统计特性来设计全局最优滤波器。维纳滤波正是从这个方向突破:它在频域中根据信号和噪声的功率谱密度来设计传递函数,使输出均方误差最小。这是历史上第一个在统计意义上最优的线性滤波器。

与前两种方法的对比

移动平均和指数平滑的滤波器设计完全靠经验调参(窗口宽度 $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)$,然后在频域应用该传递函数,再逆变换回时域。

滤波03-图1 维纳滤波:频域去噪过程

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.

Category
Tagcloud
Junck Kalman Tool Science Moving Average Exponential Smoothing UKF Nonlinear Filtering FckZhiHu Wiener Filter Code GIS Probability Code Generation Bayesian Story Architecture Non-Gaussian Nonlinear Hackintosh Hadoop Matplotlib OSX-KVM Optimal Estimation History Lens Cellular Automata Translate Unscented Transform Discuss Sequential Monte Carlo Bayesian Estimation Photo OpenWebUI Time Series Windows AI PVE Computability State Space Algorithm QGIS Kalman Filter Photography CUDA Signal Processing Pyenv RX590 PHD C Radio Communicate Prerequisites Ubuntu Sigma Point Simulation LLM Tools Geology Ollama Visualization GlumPy Programming University Scholar Particle Filter ChromeBook Hack RTL-SDR Kivy Guide Qwen3 Camera Math EKF Turing LlamaFactory AMD Jacobian Book Data Game Ventoy Hardware VirtualMachine Mac Frequency Domain Memory Computer GPT-OSS Learning Linux Filtering ML Python Life NumPy Mathematical Modeling Poem Microscope