1 基于能量约束最优化的维纳滤波推导
维纳滤波是一种对于平稳随机信号的均方误差最优的线性滤波器，即添加了平稳随机噪声（如高斯白噪声）的确定性信号经过维纳滤波系统后的输出与原始信号的均方误差最小。具体的维纳滤波推导方式有很多，因为维纳滤波本身并没有给出具体的滤波器形式，本节主要翻译自以下链接：
并从线性问题的最优化角度来对维纳滤波进行推导，在理解上会比较直观和简单，可以作为理解维纳滤波的一种参考，如果希望更加深入理解可查看第二节的内容，或者在网上查阅相关的文献或博客。
维纳滤波一般会在频域中进行分析，因为时域或空域中的卷积在频域中可表示为乘积，为了方便推导一般将信号的频谱展开为一维的列向量，而频响函数则表示为对角矩阵，那么系数的输出就可表示为频响矩阵与信号向量的相乘。对于某个原始信号的频谱 s \bf{s} s，假设添加的噪声的频谱为 n \bf{n} n，因为线性变换满足线性叠加性质，所以添加了噪声的信号的频谱可表示为
f = s + n . (1.1) {\mathbf{f}} = {\mathbf{s}} + {\mathbf{n}}.\tag{1.1} f=s+n.(1.1)
实际上，在添加噪声之前，一般认为信号 s \bf{s} s 首先会经过一个退化系统，例如一个低通滤波器，或者图像中由于运动产生的动态模糊等等，这些可以表示为一个线性移不变（Linear Shift-Invariant, LSI）系统，记其频响函数为 H d \bf{H_d} Hd，那么信号 f \bf{f} f 可表示为
f = H d s + n . (1.2) {\mathbf{f}} = {{\mathbf{H}}_{\mathbf{d}}}{\mathbf{s}} + {\mathbf{n}}.\tag{1.2} f=Hds+n.(1.2)
这时，我们希望能找到一个 LSI 系统响应 H \bf{H} H，作为对信号 f \bf{f} f 的维纳滤波，并得到相应的最优估计
g = H f . (1.3) {\mathbf{g}} = {\mathbf{Hf}}.\tag{1.3} g=Hf.(1.3)
维纳滤波希望输出的信号 g \bf{g} g 与原始信号 s \bf{s} s 的均方误差最小，但是直接求解 H \bf{H} H 通常是比较困难的，因为随机噪声 n \bf{n} n 的存在使得原始信号 s \bf{s} s 的值难以得知，而且本身要推导 H \bf{H} H 的解析解也不是一件易事。所以，维纳滤波的基本分析方法是将系统输出 g \bf{g} g 作为需要优化的量，然后间接地求解出相应的系统响应 H \bf{H} H。
根据式 (1.2)，我们可以得到噪声 n \bf{n} n 的能量为
∣ H d s − f ∣ 2 = ∣ n ∣ 2 . (1.4) {\left| {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{s}} - {\mathbf{f}}} \right|^2} = {\left| {\mathbf{n}} \right|^2}.\tag{1.4} ∣Hds−f∣2=∣n∣2.(1.4)
因为我们希望系统的输出 g \bf{g} g 尽可能地与原始信号 s \bf{s} s 相同，所以可以认为其满足以下的噪声能量约束，即
∣ H d g − f ∣ 2 = ∣ n ∣ 2 . (1.5) {\left| {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{g}} - {\mathbf{f}}} \right|^2} = {\left| {\mathbf{n}} \right|^2}.\tag{1.5} ∣Hdg−f∣2=∣n∣2.(1.5)
如果噪声 n \bf{n} n 的能量很小，这时我们可以将问题简化为最小化式 (1.5) 的左边部分，也就是使得一阶导数为零，即
∂ ( H d g − f ) ∗ T ( H d g − f ) ∂ g = 2 H d ∗ T ( H d g − f ) = 0. (1.6) \frac{{\partial {{\left( {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{g}} - {\mathbf{f}}} \right)}^{*T}}\left( {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{g}} - {\mathbf{f}}} \right)}}{{\partial {\mathbf{g}}}} = 2{\mathbf{H}}_{\mathbf{d}}^{*T}\left( {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{g}} - {\mathbf{f}}} \right) = {\mathbf{0}}.\tag{1.6} ∂g∂(Hdg−f)∗T(Hdg−f)=2Hd∗T(Hdg−f)=0.(1.6)
注意式中的转置应该为共轭转置，因为这里的值都为复数。那么可得
g = ( H d ∗ T H d ) − 1 H d ∗ T f = H d − 1 f . (1.7) {\mathbf{g}} = {\left( {{\mathbf{H}}_{\mathbf{d}}^{*T}{{\mathbf{H}}_{\mathbf{d}}}} \right)^{ - 1}}{\mathbf{H}}_{\mathbf{d}}^{*T}{\mathbf{f}} = {\mathbf{H}}_{\mathbf{d}}^{ - 1}{\mathbf{f}}.\tag{1.7} g=(Hd∗THd)−1Hd∗Tf=Hd−1f.(1.7)
注意，正如我们前面提到的，频率响应函数 H d \bf{H_d} Hd 为对角复数矩阵。所以，在不考虑噪声 n \bf{n} n 的情况下，为了使维纳滤波的输出 g \bf{g} g 与原始信号 f \bf{f} f 的均方误差最小，其频响函数刚好为退化系统响应 H d \bf{H_d} Hd 的逆。其实这也很容易理解，当不考虑噪声的影响时，信号的改变都来自于退化系统，而退化系统的响应函数是确定的，所以我们能够很容易地无失真地恢复到原始的信号。然而，很多时候噪声的能量都是不可忽略的，这时前面的假设就不再适合了。例如假设退化系统为低通滤波器，那么其高频响应幅值一般很小，信号 s \bf{s} s 相应的高频分量经过退化系统后也会变得很小，那么对于信号 g \bf{g} g 的功率谱，有
∣ g ( u ) ∣ 2 = ∣ H d − 1 ( u ) ⋅ f ( u ) ∣ 2 ≈ ∣ H d ( u ) ∣ − 2 ∣ n ( u ) ∣ 2 . (1.8) {\left| {g(u)} \right|^2} = {\left| {H_d^{ - 1}(u) \cdot f(u)} \right|^2} \approx {\left| {{H_d}(u)} \right|^{ - 2}}{\left| {n(u)} \right|^2}.\tag{1.8} ∣g(u)∣2= Hd−1(u)⋅f(u) 2≈∣Hd(u)∣−2∣n(u)∣2.(1.8)
所以，这时高频分量中噪声的能量将会被极大地放大。不过，如果不考虑退化系统，即 H d = I \bf{H_d=I} Hd=I，那么噪声的能量不会被放大，当然，也不会被缩小。
通常来说，噪声的引入会导致信号的能量增加。所以，在滤波系统输出 g \bf{g} g 满足式 (1.5) 的噪声能量约束的同时，我们希望其本身的能量也是尽可能低的，或者说至少某些分量上的能量尽可能低。为了解决以上的问题，我们引入一个评价矩阵 Q \bf{Q} Q，用于赋予信号 g \bf{g} g 的不同分量以不同的权值。注意，评价矩阵 Q \bf{Q} Q 也是一个对角矩阵，虽然在变换域中每个对角元素只作用于信号 g \bf{g} g 的单个分量，但在时域或空域中由于卷积的缘故可能涉及多个时刻或位置的输入。那么，我们就可得到以下的最优化问题，即
min J ( g ) = ∣ Q g ∣ 2 , s . t . ∣ H d g − f ∣ 2 = ∣ n ∣ 2 . (1.9) \begin{aligned} & \min \quad J\left( {\mathbf{g}} \right) = {\left| {{\mathbf{Qg}}} \right|^2}, \\ & s.t.\quad {\left| {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{g}} - {\mathbf{f}}} \right|^2} = {\left| {\mathbf{n}} \right|^2}. \\ \end{aligned} \tag{1.9} minJ(g)=∣Qg∣2,s.t.∣Hdg−f∣2=∣n∣2.(1.9)
为了解决这个最优化问题，我们可以引入拉格朗日乘子 λ \lambda λ，根据有约束极值问题的拉格朗日一阶必要条件，有
∂ ∂ g [ J ( g ) + λ ( ∣ H d g − f ∣ 2 − ∣ n ∣ 2 ) ] = ∂ ∂ g { g ∗ T Q ∗ T Q g + λ [ ( H d g − f ) ∗ T ( H d g − f ) − n ∗ T n ] } = 2 Q ∗ T Q g + 2 λ H d ∗ T ( H d g − f ) = 0. (1.10) \begin{aligned} & \frac{\partial }{{\partial {\mathbf{g}}}}\left[ {J\left( g \right) + \lambda \left( {{{\left| {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{g}} - {\mathbf{f}}} \right|}^2} - {{\left| {\mathbf{n}} \right|}^2}} \right)} \right] \\ & = \frac{\partial }{{\partial {\mathbf{g}}}}\left\{ {{{\mathbf{g}}^{*T}}{{\mathbf{Q}}^{*T}}{\mathbf{Qg}} + \lambda \left[ {{{\left( {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{g}} - {\mathbf{f}}} \right)}^{*T}}\left( {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{g}} - {\mathbf{f}}} \right) - {{\mathbf{n}}^{*T}}{\mathbf{n}}} \right]} \right\} \\ & = 2{{\mathbf{Q}}^{*T}}{\mathbf{Qg}} + 2\lambda {\mathbf{H}}_{\mathbf{d}}^{*T}\left( {{{\mathbf{H}}_{\mathbf{d}}}{\mathbf{g}} - {\mathbf{f}}} \right) = {\mathbf{0}}. \\ \end{aligned} \tag{1.10} ∂g∂[J(g)+λ(∣Hdg−f∣2−∣n∣2)]=∂g∂{g∗TQ∗TQg+λ[(Hdg−f)∗T(Hdg−f)−n∗Tn]}=2Q∗TQg+2λHd∗T(Hdg−f)=0.(1.10)
那么可得
g = ( H d ∗ T H d + γ Q ∗ T Q ) − 1 H d ∗ T f , γ = 1 / λ . (1.11) {\mathbf{g}} = {\left( {{\mathbf{H}}_{\mathbf{d}}^{*T}{{\mathbf{H}}_{\mathbf{d}}} + \gamma {{\mathbf{Q}}^{*T}}{\mathbf{Q}}} \right)^{ - 1}}{\mathbf{H}}_{\mathbf{d}}^{*T}{\mathbf{f}},\quad \gamma = 1/\lambda .\tag{1.11} g=(Hd∗THd+γQ∗TQ)−1Hd∗Tf,γ=1/λ.(1.11)
可以看到，通过这种方法，我们也可以间接得到从有噪声信号f到滤波输出 g \bf{g} g 的系统频率响应 H \bf{H} H。然而，我们很难得到拉格朗日乘子 λ \lambda λ 的解析解，虽然也可以通过迭代的方法来逐步接近期望的值，但在实际应用中一般都是直接按照经验来设置 γ \gamma γ 的值。因为式中的矩阵都为对角矩阵，我们可以按标量写为
H ( u ) = H d ∗ ( u ) ∣ H d ( u ) ∣ 2 + γ ∣ Q ( u ) ∣ 2 . (1.12)) H(u) = \frac{{H_d^ * (u)}}{{{{\left| {{H_d}(u)} \right|}^2} + \gamma {{\left| {Q(u)} \right|}^2}}}.\tag{1.12)} H(u)=∣Hd(u)∣2+γ∣Q(u)∣2Hd∗(u).(1.12))
对于不同的 γ \gamma γ 和 Q \bf{Q} Q，我们将得到具有不同频率响应特性的维纳滤波器。当两者至少其一接近零的时候，所得的频率响应函数与式 (1.7) 一致，也就是退化系统的逆，即我们只考虑退化系统对原始信号的影响，而不考虑噪声所带来的失真，所以这时的维纳滤波并不能去除信号中的噪声。除此以外常用的选择有以下三种：
（1） Q = I \bf{Q=I} Q=I。这时，我们相当于在最小化输出信号 g \bf{g} g 的能量。前面也有提到，一般来说噪声会增加信号的总能量，所以在满足式 (1.5) 的噪声能量约束的前提下，尽可能降低输出信号的能量有助于减少噪声的影响，使得输出更加平滑。并且这种情况不需要任何关于原始信号 s \bf{s} s 和噪声 n \bf{n} n 的信息，适用于很多两者频谱难以获得的情况。类似于式 (1.8)，同样考虑退化系统为一个低通滤波器，即高频响应幅值很小，那么对于高频分量来说，有
H ( u ) ≈ γ − 1 H d ∗ ( u ) . (1.13) H(u) \approx {\gamma ^{ - 1}}H_d^ * (u).\tag{1.13} H(u)≈γ−1Hd∗(u).(1.13)
∣ g ( u ) ∣ 2 ≈ ∣ γ − 1 H d ∗ ( u ) ⋅ f ( u ) ∣ 2 ≈ γ − 2 ∣ H d ( u ) ∣ 2 ∣ n ( u ) ∣ 2 . (1.14) {\left| {g(u)} \right|^2} \approx {\left| {{\gamma ^{ - 1}}H_d^ * (u) \cdot f(u)} \right|^2} \approx {\gamma ^{ - 2}}{\left| {{H_d}(u)} \right|^2}{\left| {n(u)} \right|^2}.\tag{1.14} ∣g(u)∣2≈ γ−1Hd∗(u)⋅f(u) 2≈γ−2∣Hd(u)∣2∣n(u)∣2.(1.14)
所以，通过增加约束以及选择合适的 γ \gamma γ 值，我们并不会将信号f的高频分量中的噪声的能量放大，甚至还可以进行缩小。注意，这里的分析只适用于退化系统为低通滤波器等高频响应幅值很小的情况。
（2）不同于类型 (1)，类型 (2) 需要用到原始信号 s \bf{s} s 和噪声 n \bf{n} n 的频谱信息，并且评价矩阵 Q \bf{Q} Q 不再是单位矩阵，即有
γ ∣ Q ( u ) ∣ 2 = P n ( u ) P s ( u ) . (1.15) \gamma {\left| {Q(u)} \right|^2} = \frac{{{P_n}(u)}}{{{P_s}(u)}}.\tag{1.15} γ∣Q(u)∣2=Ps(u)Pn(u).(1.15)
其中 P n ( u ) P_n(u) Pn(u) 和 P s ( u ) P_s(u) Ps(u) 分别代表噪声 n \bf{n} n 和原始信号 s \bf{s} s 的功率谱，这时的滤波称之为维纳反卷积滤波（Wiener Deconvolution Filter），相应地式 (1.12) 可改写为
H ( u ) = H d ∗ ( u ) ∣ H d ( u ) ∣ 2 + P n ( u ) P s ( u ) . (1.16) H(u) = \frac{{H_d^ * (u)}}{{{{\left| {{H_d}(u)} \right|}^2} + \frac{{{P_n}(u)}}{{{P_s}(u)}}}}.\tag{1.16} H(u)=∣Hd(u)∣2+Ps(u)Pn(u)Hd∗(u).(1.16)
回顾式 (1.9)，我们通过评价矩阵 Q \bf{Q} Q 来控制不同分量的权值，当噪声的强度相比于原始信号越高时，应该赋予该分量更高的权重，因为该分量的滤波结果更有可能包含更多的噪声。关于式 (1.15) 的详细推导可查看第二节的内容。同样考虑退化系统为低通滤波器的情况，不同于式 (1.14)，对于高频分量有
∣ g ( u ) ∣ 2 ≈ [ P s ( u ) P n ( u ) ] 2 ⋅ ∣ H d ( u ) ∣ 2 ⋅ ∣ n ( u ) ∣ 2 . (1.17) {\left| {g(u)} \right|^2} \approx {\left[ {\frac{{{P_s}(u)}}{{{P_n}(u)}}} \right]^2} \cdot {\left| {{H_d}(u)} \right|^2} \cdot {\left| {n(u)} \right|^2}.\tag{1.17} ∣g(u)∣2≈[Pn(u)Ps(u)]2⋅∣Hd(u)∣2⋅∣n(u)∣2.(1.17)
当噪声能量越高，而原始信号能量越低时，输出信号 g \bf{g} g 的高频分量的噪声能量将会被衰减得越多，所以相比于类型 (1)，类型 (2) 能够获得更加平滑的输出。但是，类型 (2) 的缺点也很明显，即需要知道噪声 n \bf{n} n 和原始信号 s \bf{s} s 的功率谱，但这两者通常是不可获得的，极大地限制了其应用范围。
（3）不同于前两者直接最小化输出信号 g \bf{g} g 的能量，类型 (3) 实际上是最小化其时域或空域二阶微分的能量，即令
q ( n ) = [ − 1 2 − 1 ] , q(n) = \left[ {\begin{array}{c} { - 1}&2&{ - 1} \end{array}} \right], q(n)=[−12−1],
而矩阵 Q \bf{Q} Q 为 q ( n ) q(n) q(n) 的频谱对角化的结果。最小化二阶微分也就是说希望输出信号 g \bf{g} g 的一阶微分基本为常数，使得其变化相对平稳，例如图像中的渐变色。而且 q ( n ) q(n) q(n) 在频域中的幅值类似于二次曲线，随着频率的提升而接近二次方地增长。所以，相比于类型 (2)，类型 (3) 甚至能够获得更加平滑的输出，并且不需要知道关于噪声和原始信号的任何信息。
图 1 展示了一个例子，其中原始图像首先经过高斯低通滤波器模糊（即退化系统）后，再添加均值为零的加性高斯白噪声，从而获得相应的退化图像。右上角为类型 (1) 的维纳滤波结果，其中 γ = 0.02 \gamma = 0.02 γ=0.02，可以看到其平滑效果效果是最差的，但是能恢复较多的细节。左下角为类型 (3) 的滤波，其平滑效果最好，但是也抹去了很多细节。右下角为类型 (2) 的滤波，因为我们知道原始图像和噪声，所以在获得适中的平滑效果的同时能够恢复较多的细节。
图1 维纳滤波效果示例
2 维纳滤波与最小二乘法
通常我们只能获取信号受噪声影响后的结果，而无法得知其真实的取值，毕竟如果能够直接获取真实信号，就没必要对受噪声影响的信号进行处理了，因此维纳滤波通常基于噪声的数学期望进行分析，包括第一节中最小化能量的先验假设亦是如此。然而在某些场景中，我们确实可以同时得到原始信号与噪声信号，并且非常有必要对噪声信号进行维纳滤波处理，这些场景中的典型就是视频图像的编解码。图像有损压缩的过程不可避免会引入失真噪声信号，但我们也不可能把原始图片传输到解码端。因为维纳滤波基于卷积实现，而实际使用的卷积核相对于图像尺寸通常差了几个数量级，所以我们可以在编码端基于维纳滤波的最小化均方误差准则求得相应的卷积核，并编码传输至解码端，就可以以尽可能小的码率代价最大化解码图像的重建质量，这也是 H.266 中自适应环路滤波的基本思路。因为我们可以同时得到原始信号与噪声信号，此时的维纳滤波实际上是一种确定性的优化问题，我们也能发现确定性的维纳滤波本质上就是最小二乘法，因为它们两者都是基于最小化均方误差准则出发的。
以离散信号的时域或空域滤波为例，记真实信号为 s i {s_i} si，受噪声影响后的信号为 f i {f_i} fi，其中信号的长度并不要求为有限值，为了方便暂且记为 I I I。同理，记长度为 N N N 的滤波核为 h n {h_n} hn，用 p n {p_n} pn 来表示滤波核各个点相对于滤波中心的坐标偏移值。那么可得卷积滤波后的均方误差为
M S E = E [ ( ∑ n = 0 N − 1 h n ⋅ f i + p n − s i ) 2 ] . (2.1) MSE = E\left[ {{{\left( {\sum\limits_{n = 0}^{N - 1} {{h_n} \cdot {f_{i + {p_n}}}} - {s_i}} \right)}^2}} \right]. \tag{2.1} MSE=E (n=0∑N−1hn⋅fi+pn−si)2 .(2.1)
为了实现均方误差的最优化，对 h n {h_n} hn 分别求偏导可得以下方程组
∂ M S E ∂ h m = E [ 2 ( ∑ n = 0 N − 1 h n ⋅ f i + p n − s i ) ⋅ f i + p m ] = 0. (2.2) \frac{{\partial MSE}}{{\partial {h_m}}} = E\left[ {2\left( {\sum\limits_{n = 0}^{N - 1} {{h_n} \cdot {f_{i + {p_n}}}} - {s_i}} \right) \cdot {f_{i + {p_m}}}} \right] = 0. \tag{2.2} ∂hm∂MSE=E[2(n=0∑N−1hn⋅fi+pn−si)⋅fi+pm]=0.(2.2)
通过变形可得
E ( ∑ n = 0 N − 1 h n ⋅ f i + p n ⋅ f i + p m ) = E ( s i ⋅ f i + p m ) ⇒ ∑ n = 0 N − 1 h n ⋅ E ( f i + p n ⋅ f i + p m ) = E ( s i ⋅ f i + p m ) . (2.3) \begin{array}{c} E\left( {\sum\limits_{n = 0}^{N - 1} {{h_n} \cdot {f_{i + {p_n}}}} \cdot {f_{i + {p_m}}}} \right) = E\left( {{s_i} \cdot {f_{i + {p_m}}}} \right)\\ \Rightarrow \sum\limits_{n = 0}^{N - 1} {{h_n} \cdot E\left( {{f_{i + {p_n}}} \cdot {f_{i + {p_m}}}} \right)} = E\left( {{s_i} \cdot {f_{i + {p_m}}}} \right). \end{array}\tag{2.3} E(n=0∑N−1hn⋅fi+pn⋅fi+pm)=E(si⋅fi+pm)⇒n=0∑N−1hn⋅E(fi+pn⋅fi+pm)=E(si⋅fi+pm).(2.3)
以上方程组表示为矩阵形式可得
[ E ( f i + p 0 ⋅ f i + p 0 ) E ( f i + p 0 ⋅ f i + p 1 ) ⋯ E ( f i + p 0 ⋅ f i + p N − 1 ) E ( f i + p 1 ⋅ f i + p 0 ) E ( f i + p 1 ⋅ f i + p 1 ) ⋯ E ( f i + p 1 ⋅ f i + p N − 1 ) ⋮ ⋮ ⋱ ⋮ E ( f i + p N − 1 ⋅ f i + p 0 ) E ( f i + p N − 1 ⋅ f i + p 1 ) ⋯ E ( f i + p N − 1 ⋅ f i + p N − 1 ) ] ⋅ [ h 0 h 1 ⋮ h N − 1 ] = [ E ( s i ⋅ f i + p 0 ) E ( s i ⋅ f i + p 1 ) ⋮ E ( s i ⋅ f i + p N − 1 ) ] . (2.4) \left[ {\begin{array}{c} {E\left( {{f_{i + {p_0}}} \cdot {f_{i + {p_0}}}} \right)}&{E\left( {{f_{i + {p_0}}} \cdot {f_{i + {p_1}}}} \right)}& \cdots &{E\left( {{f_{i + {p_0}}} \cdot {f_{i + {p_{N - 1}}}}} \right)}\\ {E\left( {{f_{i + {p_1}}} \cdot {f_{i + {p_0}}}} \right)}&{E\left( {{f_{i + {p_1}}} \cdot {f_{i + {p_1}}}} \right)}& \cdots &{E\left( {{f_{i + {p_1}}} \cdot {f_{i + {p_{N - 1}}}}} \right)}\\ \vdots & \vdots & \ddots & \vdots \\ {E\left( {{f_{i + {p_{N - 1}}}} \cdot {f_{i + {p_0}}}} \right)}&{E\left( {{f_{i + {p_{N - 1}}}} \cdot {f_{i + {p_1}}}} \right)}& \cdots &{E\left( {{f_{i + {p_{N - 1}}}} \cdot {f_{i + {p_{N - 1}}}}} \right)} \end{array}} \right] \cdot \left[ {\begin{array}{c} {{h_0}}\\ {{h_1}}\\ \vdots \\ {{h_{N - 1}}} \end{array}} \right] = \left[ {\begin{array}{c} {E\left( {{s_i} \cdot {f_{i + {p_0}}}} \right)}\\ {E\left( {{s_i} \cdot {f_{i + {p_1}}}} \right)}\\ \vdots \\ {E\left( {{s_i} \cdot {f_{i + {p_{N - 1}}}}} \right)} \end{array}} \right].\tag{2.4} E(fi+p0⋅fi+p0)E(fi+p1⋅fi+p0)⋮E(fi+pN−1⋅fi+p0)E(fi+p0⋅fi+p1)E(fi+p1⋅fi+p1)⋮E(fi+pN−1⋅fi+p1)⋯⋯⋱⋯E(fi+p0⋅fi+pN−1)E(fi+p1⋅fi+pN−1)⋮E(fi+pN−1⋅fi+pN−1) ⋅ h0h1⋮hN−1 = E(si⋅fi+p0)E(si⋅fi+p1)⋮E(si⋅fi+pN−1) .(2.4)
当信号的长度 I I I 足够大时，令
R = [ f i 0 + p 0 f i 0 + p 1 ⋯ f i 0 + p N − 1 f i 1 + p 0 f i 1 + p 1 ⋯ f i 1 + p N − 1 ⋮ ⋮ ⋱ ⋮ f i I − 1 + p 0 f i I − 1 + p 1 ⋯ f i I − 1 + p N − 1 ] ∈ R I × N . (2.5) {\mathbf{R}} = \left[ {\begin{array}{c} {{f_{{i_0} + {p_0}}}}&{{f_{{i_0} + {p_1}}}}& \cdots &{{f_{{i_0} + {p_{N - 1}}}}} \\ {{f_{{i_1} + {p_0}}}}&{{f_{{i_1} + {p_1}}}}& \cdots &{{f_{{i_1} + {p_{N - 1}}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{f_{{i_{I - 1}} + {p_0}}}}&{{f_{{i_{I - 1}} + {p_1}}}}& \cdots &{{f_{{i_{I - 1}} + {p_{N - 1}}}}} \end{array}} \right] \in {\mathbb{R}^{I \times N}}. \tag{2.5} R= fi0+p0fi1+p0⋮fiI−1+p0fi0+p1fi1+p1⋮fiI−1+p1⋯⋯⋱⋯fi0+pN−1fi1+pN−1⋮fiI−1+pN−1 ∈RI×N.(2.5)
则式 (2.4) 近似可表示为
R T R h = R T s . (2.6) {{\bf{R}}^T}{\bf{Rh}} = {{\bf{R}}^T}{\bf{s}}.\tag{2.6} RTRh=RTs.(2.6)
于是可得方程组的解为
h ∗ = ( R T R ) − 1 R T s . (2.7) {{\bf{h}}^*} = {\left( {{{\bf{R}}^T}{\bf{R}}} \right)^{ - 1}}{{\bf{R}}^T}{\bf{s}}. \tag{2.7} h∗=(RTR)−1RTs.(2.7)
可以发现，式 (2.7) 所表示的解与求解超定方程组
R h = s (2.8) {\bf{Rh}} = {\bf{s}} \tag{2.8} Rh=s(2.8)
所使用最小二乘法求得的解是一致的，因此也证明了基于最小化均方误差准则的确定性维纳滤波本质上就是最小二乘法。
然而，我们需要注意一个比较容易忽略的问题。文章的开头着重强调了维纳滤波主要用于去除平稳随机噪声，其一种简单的理解是噪声的分布如方差等不随时间或空间而改变，但这通常是不能满足的。例如在图像中，通常认为像素的亮度因为光量子的性质而符合泊松分布，因此在跨越纹理边缘时，噪声的分布特性会因为像素亮度的急剧变化而产生较大的波动。这就意味着使用一个简单的卷积核并不能很好地去除图像中的噪声。为此，我们需要尝试增大卷积核的长度，使其甚至能够与信号的长度相当，以此来增强对不同噪声分布的适应性。但这时我们又面临着更多的问题。当卷积核的长度 N N N 变得很大时，式 (2.6) 所表示的线性方程组将是一个病态问题，求解结果极易受到噪声的影响，而式 (2.6) 又是式 (2.4) 的近似结果，其必然包含一定的噪声扰动。因此，如果直接按照前面的最小二乘法对一个很大的卷积核进行求解，所得的结果通常是比较不稳定的，更不要说在这过程中所需的极高计算复杂度。为了避免单个较小的卷积核所带来的局限，H.266 实际使用了多达 25 组的小卷积核进行自适应选择，以此来平衡滤波性能与计算复杂度的矛盾。
既然这样，我们似乎没有办法求解一个稳定的超大卷积核了，但事实并非如此。要注意，前面的结论其实是在时域或空域推导所得的，但我们应该可以发现，式 (2.6) 其实可以在频域中有等价的形式。首先我们需要简单地了解相关函数，即对于序列 { x i } \{ {x_i}\} {xi} 与 { y i } \{ {y_i}\} {yi}，其相关性可表示为
R x , y ( n ) = ∑ i x i ⋅ y i + n . (2.9) {R_{x,y}}\left( n \right) = \sum\limits_i {{x_i} \cdot {y_{i + n}}} .\tag{2.9} Rx,y(n)=i∑xi⋅yi+n.(2.9)
当两个序列相同时，称为自相关；否则称为互相关。相关函数具有一些良好的性质，包括
R x y ( n ) = ∑ i x i ⋅ y i + n = ∑ j x j − n ⋅ y j = R y x ( − n ) , R x x ( − n ) = ∑ i x i ⋅ x i − n = ∑ j x j + n ⋅ x j = R x x ( n ) , R x x ( n − m ) = R x x ( m − n ) = ∑ i x i ⋅ x i − n + m = ∑ j x j + n ⋅ x j + m . (2.10) \begin{array}{c} {R_{xy}}\left( n \right) = \sum\limits_i {{x_i} \cdot {y_{i + n}}} = \sum\limits_j {{x_{j - n}} \cdot {y_j}} = {R_{yx}}\left( { - n} \right),\\ {R_{xx}}\left( { - n} \right) = \sum\limits_i {{x_i} \cdot {x_{i - n}}} = \sum\limits_j {{x_{j + n}} \cdot {x_j}} = {R_{xx}}\left( n \right),\\ {R_{xx}}\left( {n - m} \right) = {R_{xx}}\left( {m - n} \right) = \sum\limits_i {{x_i} \cdot {x_{i - n + m}}} = \sum\limits_j {{x_{j + n}} \cdot {x_{j + m}}} . \end{array}\tag{2.10} Rxy(n)=i∑xi⋅yi+n=j∑xj−n⋅yj=Ryx(−n),Rxx(−n)=i∑xi⋅xi−n=j∑xj+n⋅xj=Rxx(n),Rxx(n−m)=Rxx(m−n)=i∑xi⋅xi−n+m=j∑xj+n⋅xj+m.(2.10)
那么，式 (2.6) 可作以下变形，即
R T R h = R T s ⇔ ∑ n = 0 N − 1 h ( n ) ⋅ R f f ( m − n ) = h ( m ) ∗ R f f ( m ) = R s f ( m ) . (2.11) \begin{array}{c} {{\bf{R}}^T}{\bf{Rh}} = {{\bf{R}}^T}{\bf{s}} \Leftrightarrow \\ \sum\limits_{n = 0}^{N - 1} {h\left( n \right) \cdot {R_{ff}}\left( {m - n} \right)} = h\left( m \right) * {R_{ff}}\left( m \right) = {R_{sf}}\left( m \right). \end{array}\tag{2.11} RTRh=RTs⇔n=0∑N−1h(n)⋅Rff(m−n)=h(m)∗Rff(m)=Rsf(m).(2.11)
我们知道，时域卷积在频域中就是乘积，实际上相关也具有类似的性质，即
∫ − ∞ ∞ [ ∫ − ∞ ∞ f ( τ ) g ( t + τ ) d τ ] e − j ω t d t = ∫ − ∞ ∞ f ( τ ) [ ∫ − ∞ ∞ g ( t + τ ) e − j ω t d t ] d τ = ∫ − ∞ ∞ f ( τ ) e j ω τ d τ ∫ − ∞ ∞ g ( t + τ ) e − j ω ( t + τ ) d ( t + τ ) = F ∗ ( ω ) G ( ω ) . (2.12) \begin{array}{l} \int\limits_{ - \infty }^\infty {\left[ {\int\limits_{ - \infty }^\infty {f\left( \tau \right)g\left( {t + \tau } \right)d\tau } } \right]{e^{ - j\omega t}}dt} = \int\limits_{ - \infty }^\infty {f\left( \tau \right)\left[ {\int\limits_{ - \infty }^\infty {g\left( {t + \tau } \right){e^{ - j\omega t}}dt} } \right]d\tau } \\ = \int\limits_{ - \infty }^\infty {f\left( \tau \right){e^{j\omega \tau }}d\tau \int\limits_{ - \infty }^\infty {g\left( {t + \tau } \right){e^{ - j\omega \left( {t + \tau } \right)}}d\left( {t + \tau } \right)} } = {F^*}\left( \omega \right)G\left( \omega \right). \end{array} \tag{2.12} −∞∫∞[−∞∫∞f(τ)g(t+τ)dτ]e−jωtdt=−∞∫∞f(τ)[−∞∫∞g(t+τ)e−jωtdt]dτ=−∞∫∞f(τ)ejωτdτ−∞∫∞g(t+τ)e−jω(t+τ)d(t+τ)=F∗(ω)G(ω).(2.12)
也就是说，由于相关与卷积相差了一个负号，它们在频域就相差了一个共轭。于是，式 (2.11) 在频域中可表示为
H ( ω ) F ∗ ( ω ) F ( ω ) = S ∗ ( ω ) F ( ω ) . (2.13) H\left( \omega \right){F^*}\left( \omega \right)F\left( \omega \right) = {S^*}\left( \omega \right)F\left( \omega \right).\tag{2.13} H(ω)F∗(ω)F(ω)=S∗(ω)F(ω).(2.13)
于是，我们就得到了形式非常简洁的解
H ( ω ) = S ∗ ( ω ) F ∗ ( ω ) . (2.14) H\left( \omega \right) = \frac{{{S^*}\left( \omega \right)}}{{{F^*}\left( \omega \right)}}.\tag{2.14} H(ω)=F∗(ω)S∗(ω).(2.14)
由于离散信号的傅里叶变换可通过快速傅里叶变换算法求解，同时也不需要求解任何矩阵的逆，所以以上的计算复杂度并不高，相比于时域的求解无疑是十分可行的。而且此时的卷积核的大小并不需要作为一个预设的参数，其可以是任意且无穷的，当然由于离散傅里叶变换的采样限制，一般可以认为与信号的长度一致。不过必须要承认的是，由于式 (2.6) 只是式 (2.3) 的近似估计，其并没有利用任何噪声的先验分布参数，所以在实际的应用中，式 (2.14) 由于噪声的偏差往往具有较差的效果。为此，虽然我们已经把维纳滤波转换到频域，我们也必须从噪声分布的角度来进行进一步的分析。
正如第一节的内容，我们假设噪声与原始信号无关，且可表示线性叠加的关系，即如
f i = s i + n i . (2.15) {f_i} = {s_i} + {n_i}.\tag{2.15} fi=si+ni.(2.15)
一般假设噪声的均值为零，那么有
E ( f i + n ⋅ f i + m ) = E ( s i + n ⋅ s i + m ) + E ( n i + n ⋅ n i + m ) , E ( s i ⋅ f i + n ) = E ( s i + n ⋅ s i + m ) . (2.16) \begin{array}{c} E\left( {{f_{i + n}} \cdot {f_{i + m}}} \right) = E\left( {{s_{i + n}} \cdot {s_{i + m}}} \right) + E\left( {{n_{i + n}} \cdot {n_{i + m}}} \right),\\ E\left( {{s_i} \cdot {f_{i + n}}} \right) = E\left( {{s_{i + n}} \cdot {s_{i + m}}} \right). \end{array}\tag{2.16} E(fi+n⋅fi+m)=E(si+n⋅si+m)+E(ni+n⋅ni+m),E(si⋅fi+n)=E(si+n⋅si+m).(2.16)
于是从数学期望的角度来看，式 (2.13) 可表示为
H ( ω ) [ S ∗ ( ω ) S ( ω ) + N ∗ ( ω ) N ( ω ) ] = S ∗ ( ω ) S ( ω ) . (2.17) H\left( \omega \right)\left[ {{S^*}\left( \omega \right)S\left( \omega \right) + {N^*}\left( \omega \right)N\left( \omega \right)} \right] = {S^*}\left( \omega \right)S\left( \omega \right).\tag{2.17} H(ω)[S∗(ω)S(ω)+N∗(ω)N(ω)]=S∗(ω)S(ω).(2.17)
注意复数乘以其共轭，得到的是其模的平方，是一个实数，在频域中通常称为功率谱。因此，记噪声的功率谱为 P n ( ω ) {P_n}\left( \omega \right) Pn(ω)，原始信号的功率谱为 P s ( ω ) {P_s}\left( \omega \right) Ps(ω)，则基于以上噪声先验所得的维纳滤波响应可表示为
H ( ω ) = P s ( ω ) P s ( ω ) + P n ( ω ) = 1 1 + 1 S N R , S N R = P s ( ω ) P n ( ω ) . (2.18) H\left( \omega \right) = \frac{{{P_s}\left( \omega \right)}}{{{P_s}\left( \omega \right) + {P_n}\left( \omega \right)}} = \frac{1}{{1 + {\textstyle{1 \over {SNR}}}}},\;SNR = \frac{{{P_s}\left( \omega \right)}}{{{P_n}\left( \omega \right)}}.\tag{2.18} H(ω)=Ps(ω)+Pn(ω)Ps(ω)=1+SNR11,SNR=Pn(ω)Ps(ω).(2.18)
这样，我们就得到第一节中前二类维纳滤波器形式的由来。要注意的是，虽然当我们知道真实信号时，噪声的具体采样取值也可得知，但此时的噪声的功率谱 P n ( ω ) {P_n}\left( \omega \right) Pn(ω) 通常不会直接通过该采样取值的傅里叶变换求解，这是因为单次采样的结果并不能反映噪声的真实功率谱。例如当噪声为理想的高斯白噪声时，其功率谱在整个频域中都应该为常数，当傅里叶变换未正交化时即为高斯分布的方差乘以信号的长度，这时我们只需直接将其代入其功率谱即可。对于大部分情况，信号与噪声的功率谱都难以获取，因此更常用的方式是直接配置信噪比。应该看到，式 (2.18) 总是小于 1 的。一般信号的主要能量都集中在低频区域，而高频区往往具有很小的值。因此，基于最小化均方误差准则的维纳滤波对于信号高频的衰减要比低频严重，例如在图像中会导致一些细小纹理丢失，以及边缘的锐度会有所下降，总体上会变得相对模糊，这也是低通滤波器的基本特性。
扫一扫
711
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
