降噪：维纳滤波

降噪 Débruitage

在许多情况下，去噪信号是必要的。实际上，信号（如声音、图像、振动等）的测量不可避免地会受到不确定性的影响。例如以下噪声来源：

热噪声（例如，原子波动等）
脉冲噪声（如图像中的坏点像素）
量化噪声
干扰信号
等等

建模

将所得的测量建模为：

\[ y = x+b \]

其中，\(x\)是目标信号，\(b\)是噪音，\(y\)是测量值。我们目的是根据现有信息获得\(x\)的估计值：

\[ \hat{x}=f(y) \]

在本文中，我们将假设观察信号 \(x\) 是随机变量 \(\mathbf{X}\) 的一个实现，同样地，\(b\) 和 \(y\) 是随机变量 \(\mathbf{B}\) 和 \(\mathbf{Y}\) 的实现。

假设\(\mathbf{X}\) 和 \(\mathbf{B}\) 是独立的。
文中出现的所有随机变量都属于 \(L^2\) 空间且均值为零。

我们通过概率分布获得关于 \(\mathbf{X}\) 的信息，这属于贝叶斯框架

我们的目标是最小化估计误差，这被称为估计器的“风险”（Risque），定义如下：

\[ r(f, \mathbf{X})=E\left(\|\mathbf{X}-\hat{\mathbf{X}}\|^2\right) \]

这是估计误差的期望值，因此被\(\mathbf{X}\) 的概率分布加权。我们将称能最小化此风险的估计器为“最优估计 estimateur optimal”。

\(L^2\)空间

\(L^2\) 空间（也称为平方可积函数空间）是一个数学分析中非常重要的函数空间，\(L^2\) 空间是定义在某个测度空间上的所有平方可积函数的集合。

函数 \(f(x)\)属于\(L^2\)空间，如果满足以下条件：

\[ \int_{\Omega}|f(x)|^2 \mathrm{~d} \mu<\infty \]

其中:

\(\Omega\) 是函数的定义域；

\(\mu\) 是定义在 \(\Omega\) 上的测度；

一般情况下的最优估计 Estimateur optimal dans le cas général

从观测值 \(\mathbf{Y}\) 中估计 \(\mathbf{X}\) 的最佳估计器 \(\hat{\mathbf{X}}\) 是条件期望：

\[ \hat{\mathbf{X}}=E(\mathbf{X} \mid \mathbf{Y}) \]

实际上， \(E(\mathbf{X} \mid \mathbf{Y})\) 是 \(\mathbf{X}\) 在 \(\mathbf{Y}\)-可测变量上的正交投影。因此，它最小化 \(E\left(\|\mathbf{X}-\mathbf{Y}\|^2\right)\) 。此外，由于它是一个 \(\mathbf{Y}\)-可测变量，因此可以表示为 \(f(\mathbf{Y})\) 的形式。实际上，条件期望 \(E(\mathbf{X} \mid \mathbf{Y})\) 很少能够以解析形式写出，通常是通过概率分布的采样技术（如蒙特卡罗（MonteCarlo）方法等）近似计算得到。

Y-可测变量

\(Z\) 是 \(Y\)-可测变量，当且仅当 \(Z\) 是关于 \(\sigma(Y)\) 的可测函数。

给定一个随机变量 \(Y\) ，可以定义一个由 \(Y\) 所生成的 \(\sigma-\text{代数}\)，记作 \(\sigma(Y)\) 。这个 \(\sigma-代数\)是所有可以通过 \(Y\) 表达的事件的集合。例如，事件" \(Y\) 落在区间 \([a, b]\) "属于 \(\sigma(Y)\) 。如果\(Z\) 是 \(Y\)-可测变量，则有\(Z = f(Y)\)

高斯情况下的最优估计 Estimateur optimal dans le cas gaussien

假设以下条件成立:

\[ \mathbf{X} \sim \mathcal{N}\left(0, \Sigma_X\right), \mathbf{B} \sim \mathcal{N}\left(0, \Sigma_B\right) . \]

因此:

\[ \mathbf{Y} \sim \mathcal{N}\left(0, \Sigma_X+\Sigma_B\right) . \]

回顾假设， \(X\) 和 \(B\) 是相互独立的。

在这种情况下，最优估计器由下式给出:

\[ \hat{\mathbf{X}}=E(\mathbf{X} \mid \mathbf{Y})=\Sigma_X\left(\Sigma_X+\Sigma_B\right)^{-1} \mathbf{Y} \]

证明是一个经典的CIP问题。我们定义 \(\mathbf{U}=\mathbf{X}-\Sigma_X\left(\Sigma_X+\Sigma_B\right)^{-1} \mathbf{Y}\) 。可以验证 \(\mathbf{U}\) 和 \(\mathbf{Y}\)是联合高斯分布的，并且 \(E\left(\mathbf{U}^T \mathbf{Y}\right)=0\)

证明\(E\left(\mathbf{U}^T \mathbf{Y}\right)=0\)

首先计算：

\[ E\left[X_i X_j\right]=\Sigma_{X,ij} \]

\[ E\left[Y_i Y_j\right]=E\left[(X_i+B_j)(X_jB_j)\right]=E\left[X_i X_j\right]+E\left[X_i B_j\right]+E\left[B_i X_j\right]+E\left[B_i B_j\right]=\Sigma_{X,ij}+\Sigma_{B,ij} \]

然后观察到：

\[ U^t Y=\left(X-\Sigma_X\left(\Sigma_X+\Sigma_B\right)^{-1} Y\right)^t Y=X^t Y-Y^t \Sigma_X\left(\Sigma_X+\Sigma_B\right)^{-1} Y \]

其中：

\[ E\left[X_i Y_j\right]=\Sigma_{X,ij}，E\left[Y_i A_{jk} Y_k\right]=A_{jk} E\left[Y_i Y_k\right]=\Sigma_{X,jk}\left(\Sigma_{X,jk}+\Sigma_{B,jk}\right)^{-1}\left(\Sigma_{X,jk}+\Sigma_{B,jk}\right) = \Sigma_{X,jk} \]

故有：

\[ E\left(\mathbf{U}^T \mathbf{Y}\right)=0 \]

这意味着 \(\mathbf{U}\) 和 \(\mathbf{Y}\) 是独立的，且 \(E(\mathbf{U} \mid \mathbf{Y})=0\) 。因此:

\[ E(\mathbf{U} \mid \mathbf{Y})=E(\mathbf{X} \mid \mathbf{Y})-\Sigma_X\left(\Sigma_X+\Sigma_B\right)^{-1} E(\mathbf{Y} \mid \mathbf{Y})=E(\mathbf{X} \mid \mathbf{Y})-\Sigma_X\left(\Sigma_X+\Sigma_B\right)^{-1} \mathbf{Y} =0 \]

故有：

\[ E(\mathbf{X} \mid \mathbf{Y})=\Sigma_X\left(\Sigma_X+\Sigma_B\right)^{-1} \mathbf{Y} \]

此估计是线性估计。

线性最佳估计 Estimateur linéaire optimal

假设\(X\)和\(B\)的分布未知，假设其输入\(L^2\)空间，互相独立，均值为零。假设估计器形式如下：

\[ \hat{\mathbf{X}}=H \mathbf{Y} \]

其中\(H\in \mathbb R^{n\times n}\)。首先证明估计误差与数据本身无关：

\[ E\left(\left(X_n-\hat{X}_n\right) Y_k\right)=0, \forall n, k \]

风险函数写作：

\[ r=E\left(\|\mathbf{X}-\hat{\mathbf{X}}\|^2\right) \]

展开：

\[ r=E\left(\sum_{n=1}^N\left(X_n-\sum_{k=1}^N h_{n k} Y_k\right)^2\right) \]

为讨论风险函数最小的情况，对\(h_{nk}\)求偏导：

\[ \frac{\partial r}{\partial h_{n k}}=-2 E\left(Y_k\left(X_n-\sum_{l=1}^N h_{n l} Y_l\right)\right)=-2 E\left(Y_k\left(X_n-\hat{X}_n\right)\right) = 0 \]

又因为\(X\)和\(B\)相互独立，均值为零，有：

\[ E\left(Y_k X_n\right)=E\left(X_n\left(X_k+B_k\right)\right)=\Sigma_{X, n k}\\E\left(Y_k Y_l\right)=E\left(\left(X_k+B_k\right)\left(X_l+B_l\right)\right)=\Sigma_{X, k l}+\Sigma_{B, k l} \]

由上式即可推的：

\[ \begin{gathered} \Sigma_{X, n k}-\sum_{k=1}^N h_{n l}\left(\Sigma_{X, l k}+\Sigma_{B, l k}\right)=0 \\ \Sigma_X-H\left(\Sigma_X+\Sigma_B\right)=0 \end{gathered} \]

得到最佳线性估计的表达式：

\[ \hat{\mathbf{X}}=H\mathbf{Y} = \Sigma_X\left(\Sigma_X+\Sigma_B\right)^{-1} \mathbf{Y} \]

与高斯情况一致。

特殊情况之一：协方差矩阵为对角矩阵

在 \(\Sigma_X\) 和 \(\Sigma_B\) 为对角矩阵的情况下（例如噪声是白噪声（即协方差矩阵是对角阵，无论基如何），并且信号表示在其协方差矩阵为对角形式的基（如 Karhunen-Loève 基）上）。记 \(\sigma_{X, n}^2\) 和 \(\sigma_{B, n}^2\) 为这些协方差矩阵的对角元素。因此可得:

\[ \hat{X}n=\frac{\sigma{X, n}^2}{\sigma_{X, n}^2+\sigma_{B, n}^2} Y_n \]

在这种情况下，维纳滤波（Filtrage de Wiener）对应于对 \(Y_n\) 的估计 \(\hat{X}n\) ，通过信噪比的比例来调整对 \(Y_n\) 的观测强度。记 \(a_n=\frac{\sigma{X, n}^2}{\sigma_{X, n}^2+\sigma_{B, n}^2}\) ，我们可以总结出以下特性:

当 \(\sigma_{X, n}^2 \gg \sigma_{B, n}^2\) 时, \(a_n \approx 1\);
当 \(\sigma_{X, n}^2=\sigma_{B, n}^2\) 时， \(a_n=1 / 2\);
当 \(\sigma_{X, n}^2 \ll \sigma_{B, n}^2\) 时， \(a_n \approx 0\) 。

特殊情况之二：广义平稳过程 stationnaires au sens large

广义平稳（Stationarity au sens large）是随机过程的一种统计性质。它是对严格平稳性的一个放宽版本，仅要求随机过程的均值和协方差与时间无关。一个随机过程 \(\{X(t)\}\) 是广义平稳的，当且仅当它满足以下两个条件:

均直不随时间变化:

\[ E[X(t)]=\mu, \quad \forall t, \]

自协方当函数只与时间间隔（滞后）有关：

\[ \Sigma(X(t), X(t+\tau))=E[(X(t)-\mu)(X(t+\tau)-\mu)]=\gamma_X(\tau), \]

我们现在关注随机过程 \(\{X[n]\}{n \in \mathbb{Z}}\) **和 \(\{B[n]\}{n \in \mathbb{Z}}\) 。假设它们是广义平稳的、均值为零、相互独立，并且具有已知的自协方差函数：

\[ \begin{gathered}\gamma_X[k]=E(X[n] X[n+k]), \\\gamma_B[k]=E(B[n] B[n+k]) .\end{gathered} \]

这些过程是实数值的，且满足 \(\gamma_X[-k]=\gamma_X[k]\) 和 \(\gamma_B[-k]=\gamma_B[k]\) 。观察到:

\[ \begin{gathered}Y[n]=X[n]+B[n], \\\gamma_Y[k]=\gamma_X[k]+\gamma_B[k] .\end{gathered} \]

我们对 \(X[n]\) 的估计为:

\[ \hat{X}=h * Y, \]

在有限情况下，可用类似的方法得到：

\[ \gamma_X=h \star\left(\gamma_X+\gamma_B\right) \]

另一种方法是将估计误差写成功率谱密度（DSP）和h的频率相应形式：

\[ \begin{aligned}& r=E\left(|X[n]-\hat{X}[n]|^2\right)=\frac{1}{2 \pi} \int_{-\pi}^\pi \Gamma_{X-\hat{X}}(\nu) d \nu .\end{aligned} \]

\(P_x = E(x[n]^2) = \frac{1}{2 \pi} \int_{-\pi}^\pi \Gamma_{X}(\nu) d \nu .\) \(h \star X\)的DSP是\(|H(\nu)|^2 \Gamma_X(\nu)\)

其中，根据干扰公式formule des interférences：

\[ \begin{aligned}X-\hat{X} & =X-h \star X-h \star B \\& =(\delta-h) \star X-h \star B \\\Gamma_{X-\hat{X}}(\nu) & =|1-H(\nu)|^2 \Gamma_X(\nu)+|H(\nu)|^2 \Gamma_B(\nu)\end{aligned} \]

同样求导疑惑的最小化误差的情况，可以得到：

\[ H(\nu)=\frac{\Gamma_X(\nu)}{\Gamma_X(\nu)+\Gamma_B(\nu)} \]

这与维纳滤波在有限维协方差矩阵为对角阵情况下得到的系数形式类似。对 \(\gamma_X=h \star\left(\gamma_X+\gamma_B\right)\)做傅里叶变换（TFTD）可以得到相同的结果。

此时的估计误差为：

\[ r=\int_{-\pi}^\pi \frac{\Gamma_X(\nu) \Gamma_B(\nu)}{\Gamma_X(\nu)+\Gamma_B(\nu)} d \nu . \]

然而，这种滤波器在实践中难以直接应用。滤波器 \(H(\nu)\) 的傅里叶逆变换通常难以计算。此外，由于 \(H(\nu)\) 是实对称函数（与 DSP 类似），滤波器通常是非因果的，因此需要知道 \(Y[n+k], k>0\) 。这一点会导致稳定性问题需要解决。

特殊情况之三：RIF滤波

我们现在对估计滤波器 \(h\) 添加限制，使其为因果滤波器，并具有有限的冲激响应。记 \(h_p\) 为滤波器系数，满足 \(0 \leq p<P<\infty\) ，其系数非零。可以表示为:

\[ \hat{X}[n]=\sum_{p=0}^{P-1} h_p Y[n-p] \]

误差：

\[ r=E\left((X[n]-\hat{X}[n])^2\right)=E\left(\left(X[n]-\sum_{p=0}^{P-1} h_p Y[n-p]\right)^2\right) . \]

同样求导：

\[ \frac{\partial r}{\partial h_p}=-2 E\left(Y[n-p]\left(X[n]-\sum_{l=0}^{P-1} h_l Y[n-l]\right)\right)=-2\left(E(Y[n-p] X[n])-\sum_{l=0}^{P-1} h_l E(Y[n-p] Y[n-l])\right) \]

验证可得，与有限维情况类似，满足：

\[ E(X[n] Y[n+k])=\gamma_X[k], \quad E(Y[n] Y[n+k])=\gamma_X[k]+\gamma_B[k] \]

遂得：

\[ \gamma_X[p]=\sum_{l=0}^{P-1} h_l\left(\gamma_X[l-p]+\gamma_B[l-p]\right) \]

构建自相关矩阵矩阵 \(R_X 、 R_B\) 和互相关向量 \(r_X\) ，其中 \(0 \leq n, m<P\) :

\[ \begin{gathered}R_{X, n m}=\gamma_X[n-m], \\R_{B, n m}=\gamma_B[n-m], \\r_{X, n}=\gamma_X[n],\end{gathered} \]

最优 RIF 滤波器的系数向量\(h\)为：

\[ \mathbf{h}=\left(R_X+R_B\right)^{-1} \mathbf{r}_X \]

而已通过对\(h\)进行离散时间傅里叶变换（DTFT）得到\(H\)。