Wiener Filter

假设分别有两个WSS process:$x[n]$,$y[n]$,这两个process之间存在某种关系,并且我们也了解这种关系。现在我们手头上有process $x[n]$,目的是要设计一个LTI系统,使得系统输出$y[n]$,不过$y[n]$是一个WSS process,我们不可能准确得到随机过程上的值,因此实际输出并不是$y[n]$,而是$hat{y}[n]$,此时我们就能通过引入MSE来判断实际输出$hat{y}[n]$与$y[n]$之间的差距。当我们所设计的系统使得$hat{y}[n]$与$y[n]$之间有Minimun MSE时,所得到的输出就是最优的,即

$e[n] riangleq hat{y}[n] – y[n]$

当满足Minimum MSE时,有

$epsilon = MMSE = E{e^2 [n]}$

那么这个LTI系统就被称为Wiener filter

image

本文的主要目的就是求Wiener filter的脉冲响应,所求得的脉冲响应需要使得系统有Minimun MSE。

※在开始之前,先说明本文一些计算所需的前提条件:

  1. 由于输入以及输出的随机过程都是实数,因此$R_{xx}[m],R_{yy}[m],R_{xy}[m],R_{yx}[m]$也会是实数
  2. 由于$R_{xx}[m] = R_{xx}[-m]$,因此$S_{xx}(e^{jOmega})$是实数,同理有$S_{yy}(e^{jOmega})$是实数 
  3. 本文假设$S_{xx}(e^{jOmega})$以及$S_{yy}(e^{jOmega})$在所有频率上都大于0
  4. 为了方便,本文假设随机过程$x[n],y[n]$的mean都为0
  5. 连续时间随机过程有上述同样条件

NONCAUSAL DT WIENER FILTER

DT表示的是离散时间系统,noncausl即非因果,因果性的系统的脉冲响应的$n<0$部分为$0$,非因果系统则无此要求。

对于DT系统,有两种求解Wiener filter脉冲响应的方式。这是因为DT系统的脉冲响应可以分为FIR以及IIR两种,分别对应不同的求解方法。

FIR

FIR即有限脉冲响应,脉冲响应$h[n]$的长度是有限的,这表明我们可以准确得到脉冲响应序列中的各个元素的值。

DT LTI系统以卷积来表征,即

$displaystyle{hat{y}[n] = sum_{k=-infty}^{infty}h[k]x[n-k] }$

那么MSE就可以表示为

$displaystyle{epsilon = Eleft{ left(sum_{k=-infty}^{infty}h[k]x[n-k]-y[n]  ight)^2 ight}}$

由于是一个FIR filter,因此假设$h[n]$是一个长度为N的脉冲响应,也就是说上面的式子有N个未知数,分别对应$h[0],h[1],cdotcdotcdot,h[N-1]$。为了得到该N元二次函数的最小值,需要分别对这N个未知数求偏导,当偏导数为0时可以得到最小的MSE。以序列其中的一个未知数$h[m]$为例

$egin{align*}
frac{partial epsilon}{partial h[m]} &= frac{displaystyle{partial Eleft{left(sum_{k=-infty}^{infty}h[k]x[n-k]-y[n] ight)^2 ight} }}{partial h[m]}\
&= frac{displaystyle{partial Eleft{left(sum_{k=0}^{N-1}h[k]x[n-k]-y[n] ight)^2 ight} }}{partial h[m]} qquad h[k] is FIR\
&= Eleft{ 2underbrace{left(sum_{k=0}^{N-1}h[k]x[n-k]-y[n] ight )}_{e[n]}x[n-m] ight}qquad chain rule\
&= 2EBig{e[n]x[n-m]Big}\
&= 2EBig{ig(hat{y}[n]-y[n]ig)x[n-m]Big}\
&= 2EBig{hat{y}[n]x[n-m]-y[n]x[n-m]Big}\
&= 2Big(R_{hat{y}x}[m]-R_{yx}[m]Big)
end{align*}$

当偏导数为0时有最小MSE,此时

$color{red}{R_{hat{y}x}[m] = R_{yx}[m]}$

根据correlation相关定理,有

$color{red}{R_{yx}[m] = h[m]*R_{xx}[m]}$

因此,只要我们能根据随机过程$x[n],y[n]$之间的关系求得$R_{yx}[m],R_{xx}[m]$,就可以得到最佳的脉冲响应$h[m]$。卷积的计算方式如下

$egin{align*}
R_{yx}[m] &= h[m]*R_{xx}[m] = R_{xx}[m]*h[m]\
&=egin{bmatrix}
R_{xx}[0] &R_{xx}[-1]  &cdots  &R_{xx}[1-N] \
R_{xx}[1] &R_{xx}[0]  &cdots  &R_{xx}[2-N] \
vdots  &vdots  &ddots   &vdots \
R_{xx}[N-1] &R_{xx}[N-2]  &cdots  &R_{xx}[0]
end{bmatrix}
egin{bmatrix}
h[0]\
h[1]\
vdots\
h[N-1]
end{bmatrix}
=egin{bmatrix}
R_{yx}[0]\
R_{yx}[1]\
vdots\
R_{yx}[N-1]
end{bmatrix} end{align*}$

IIR

IIR即无限脉冲响应,脉冲响应$h[n]$的长度是无限的,我们无法得到脉冲响应序列中的所有元素的值,因此主要求解方式需要先从频域获得系统的频率响应$H(e^{jOmega})$,然后根据需要转换为脉冲响应$h[n]$。

如前一小节所述,当MSE式子的偏导数为0时可以得到MMSE,此时系统能输出最佳的预测值。按照此思路,同理可得到

$egin{align*}
time-domain: &qquad R_{hat{y}x}[n] = R_{yx}[n]\
z-transform: &qquad S_{hat{y}x}(z) = S_{yx}(z)\
fourier-transform: &qquad S_{hat{y}x}(e^{jOmega}) = S_{yx}(e^{jOmega})
end{align*}$

根据correlation相关定理,又有

$egin{align*}
time-domain: &qquad h[n]*R_{xx}[n] = R_{hat{y}x}[n] = R_{yx}[n]\
z-transform: &qquad H(z)S_{xx}(z)=S_{hat{y}x}(z) = S_{yx}(z)\
fourier-transform: &qquad H(e^{jOmega})S_{xx}(e^{jOmega})=S_{hat{y}x}(e^{jOmega}) = S_{yx}(e^{jOmega})
end{align*}$

那么最佳的系统函数以及最佳的频率响应分别为

$egin{align*}
z-transform: &qquad H(z)= S_{yx}(z)/S_{xx}(z)\
fourier-transform: &qquad H(e^{jOmega})=S_{yx}(e^{jOmega})/S_{xx}(e^{jOmega})
end{align*}$

因此,只要我们能根据随机过程$x[n],y[n]$之间的关系求得它们在频域上的表示$S_{yx}(e^{jOmega}),S_{xx}(e^{jOmega})$,就可以得到最佳的频率响应$H(e^{jOmega})$

MMSE

当系统为最佳滤波器的时候,MSE为最小值,既有

$MMSE = R_{ee}[0]$

其中$R_{ee}[m]$可以通过以下方法得到

$egin{align*}
R_{ee}[m] &= E{e[n+m]e[n]}\
&= EBig{ig(y[n+m]-hat{y}[n+m]ig)ig(y[n]-hat{y}[n]ig)Big}\
&= EBig{y[n+m]y[n]-hat{y}[n+m]y[n]-y[n+m]hat{y}[n]+hat{y}[n+m]hat{y}[n]Big}\
&= R_{yy}[m]-R_{hat{y}y}[m]-R_{yhat{y}}[m]+R_{hat{y}hat{y}}[m]\
&= R_{yy}[m]-R_{hat{y}y}[m]-h[-m]*R_{yx}[m]+h[-m]*R_{hat{y}x}[m]\
&= R_{yy}[m]-R_{hat{y}y}[m]-h[-m]*R_{yx}[m]+h[-m]*R_{yx}[m]\
&= R_{yy}[m]-R_{hat{y}y}[m]\ &= R_{yy}[m]-h[m]*R_{xy}[m] qquad href{https://www.cnblogs.com/TaigaCon/p/9127534.html#DT}{correlation equation} end{align*}$

MMSE在频域上则能有如下表达

$egin{align*}
MMSE = R_{ee}[0]&=frac{1}{2pi}int_{-pi}^{pi}S_{ee}(e^{jOmega})e^{jOmega 0}dOmega\
&=frac{1}{2pi}int_{-pi}^{pi}S_{ee}(e^{jOmega})dOmega\
&=frac{1}{2pi}int_{-pi}^{pi}Big(S_{yy}-HS_{xy}Big)dOmega qquad drop (e^{jOmega})\
&=frac{1}{2pi}int_{-pi}^{pi}Big(S_{yy}-frac{S_{yx}}{S_{xx}}S_{xy}Big)dOmega\
&=frac{1}{2pi}int_{-pi}^{pi}left(S_{yy}left{1-frac{S_{yx}}{S_{xx}}frac{S_{xy}}{S_{yy}} ight} ight)dOmega\
&=frac{1}{2pi}int_{-pi}^{pi}Big(S_{yy}(1- ho_{yx} ho_{yx}^*)Big)dOmega
end{align*}$

其中$ ho_{yx}(e^{jOmega})$可以看作是频域的correlation coefficient(相关系数)。定义如下

$displaystyle{ ho_{yx}(e^{jOmega}) = frac{S_{yx}(e^{jOmega})}{sqrt{S_{xx}(e^{jOmega})S_{yy}(e^{jOmega})}}}$

由于我们前面以及讲过$R_{yx}[m]=R_{xy}[m]$是实数,因此有$href{https://www.cnblogs.com/TaigaCon/p/9125855.html}{S_{xy}(e^{jOmega}) = S_{yx}^*(e^{jOmega})}$,又有$S_{xx}(e^{jOmega}),S_{yy}(e^{jOmega})$在任意频域上都大于0,有了这两个条件容易得到

$displaystyle{ ho_{yx}^*(e^{jOmega}) = frac{S_{yx}^*(e^{jOmega})}{sqrt{S_{xx}(e^{jOmega})S_{yy}(e^{jOmega})}}}$

NONCAUSAL CT WIENER FILTER

不同于DT中离散的随机过程,CT中的随机过程是连续的,因此不能像DT一样对某一个采样点求导来得到最佳系统,需要采用其他计算方式。

image

与DT时一样基于MSE的预测,要得到最佳系统,那么该系统应该使得输出的随机过程$hat{y}(t)$与目标随机过程$y(t)$的MSE最小。

$MSE = EBig{ e^2(t) Big} = EBig{ ig(hat{y}(t)-y(t)ig)^2 Big}$

我们注意到

$MSE = EBig{e^2{t}Big} = R_{ee}(0)$

因此可以通过$R_{ee}( au)$来推导如何得到最佳系统。

$egin{align*}
R_{ee}( au) &= EBig{ig(hat{y}(t)-y(t)ig)ig(hat{y}(t+ au)-y(t+ au)ig)Big}\
&= EBig{y(t)y(t+ au)+hat{y}(t)hat{y}(t+ au)-y(t)hat{y}(t+ au)-hat{y}(t)y(t+ au)Big}\
&= R_{yy}( au)+R_{hat{y}hat{y}}( au)-R_{yhat{y}}( au)-R_{hat{y}y}( au)
end{align*}$

根据这条式子可以得到

$S_{ee}(jomega) = S_{yy}(jomega)+S_{hat{y}hat{y}}(jomega)-S_{yhat{y}}(jomega)-S_{hat{y}y}(jomega)$

下面来推导如何得到Minimun MSE

$egin{align*}
MSE = R_{ee}(0) &= frac{1}{2pi}int_{-infty}^{infty}S_{ee}(jomega)e^{jomega 0}domega\
&=frac{1}{2pi}int_{-infty}^{infty}S_{ee}(jomega)domega\
&=frac{1}{2pi}int_{-infty}^{infty}Big(S_{yy}+S_{hat{y}hat{y}}-S_{yhat{y}}-S_{hat{y}y}Big)domega qquad drop (jomega )\
&=frac{1}{2pi}int_{-infty}^{infty}Big(S_{yy}+HH^*S_{xx}-H^*S_{yx}-HS_{xy} Big)domega\
&=frac{1}{2pi}int_{-infty}^{infty}Big(S_{yy}+HH^*S_{xx}-H^*S_{yx}-HS_{yx}^* Big)domegaqquad R_{xy}(t)=R_{yx}(-t) is real\
&=frac{1}{2pi}int_{-infty}^{infty}left( S_{yy}+left( Hsqrt{S_{xx}}-frac{S_{yx}}{sqrt{S_{xx}}} ight )left( H^*sqrt{S_{xx}}-frac{S_{yx}^*}{sqrt{S_{xx}}} ight )-frac{S_{yx}S_{yx}^*}{S_{xx}} ight)domega\
&=frac{1}{2pi}int_{-infty}^{infty}left| Hsqrt{S_{xx}}-frac{S_{yx}}{sqrt{S_{xx}}} ight|^2domega+frac{1}{2pi}int_{-infty}^{infty}left(S_{yy}-frac{S_{yx}S_{yx}^*}{S_{xx}} ight)domega
end{align*}$

上面的式子采用了提取式子中的平方的方法来使式子得到最小值,当系统的频率响应满足以下条件时,即可得到Minimun MSE

$color{red}{displaystyle{ H(jomega) = frac{S_{yx}(jomega)}{S_{xx}(jomega)} }}$

这也表明只要我们能根据随机过程$x[n],y[n]$之间的关系求得它们在频域上的表示$S_{yx}(jomega),S_{xx}(jomega)$,就可以得到最佳的频率响应$H(jomega)$

此时式子中第一项为0,那么MMSE为

$egin{align*}
MMSE = R_{ee}(0)
&=frac{1}{2pi}int_{-infty}^{infty}left( S_{yy}-frac{S_{yx}S_{yx}^*}{S_{xx}} ight)domega\
&=frac{1}{2pi}int_{-infty}^{infty}S_{yy}left(1-frac{S_{yx}S_{yx}^*}{S_{xx}S_{yy}} ight )domega\
&=frac{1}{2pi}int_{-infty}^{infty}S_{yy}(1- ho ho^*)domega
end{align*}$

其中频域的correlation coefficient $ ho(jomega)$为

$displaystyle{ ho(jomega) = frac{S_{yx}(jomega)}{sqrt{S_{xx}(jomega)S_{yy}(jomega)}}}$

CAUSAL CT WIENER FILTER

观察上一小节推导得到的式子

$displaystyle{MSE =frac{1}{2pi}int_{-infty}^{infty}left| Hsqrt{S_{xx}}-frac{S_{yx}}{sqrt{S_{xx}}} ight|^2domega+frac{1}{2pi}int_{-infty}^{infty}left(S_{yy}-frac{S_{yx}S_{yx}^*}{S_{xx}} ight)domega}$

当所设计的系统的频率响应使得第一项为0的时候就能得到MMSE,不过如果我们对该系统有因果性的要求,也就意味着对频率响应$H(jomega)$有所限制,此时的$H(jomega)$可能就无法使得上述等式的第一项为0,因此对于因果系统,我们需要寻求其它的求解方法。

WSS Process On Causal LTI System中,我们得到了一个结论:一个Casual LTI系统可以把white noise建模成WSS process。我们这里假设一个频率响应为$M_{xx}(jomega)$的因果系统把white noise建模成了WSS process $x(t)$,那么就可以把$S_{xx}(jomega)$分解为

$S_{xx}(jomega) = M_{xx}(jomega)M^*_{xx}(jomega)$

那么前面的式子就可以变形为

$displaystyle{MSE =frac{1}{2pi}int_{-infty}^{infty}left| HM_{xx}-frac{S_{yx}}{M_{xx}^*} ight|^2domega+frac{1}{2pi}int_{-infty}^{infty}left(S_{yy}-frac{S_{yx}S_{yx}^*}{S_{xx}} ight)domega}$

其中$H(jomega)$为因果系统,$M_{xx}(jomega)$也是因果系统,$HM_{xx}$是这两个因果系统的级联,因此也是一个因果的。此时,我们可以把$displaystyle{frac{S_{yx}(jomega)}{M_{xx}^*(jomega)}}$分成两部分之和

$displaystyle{frac{S_{yx}(jomega)}{M_{xx}^*(jomega)}=left[frac{S_{yx}(jomega)}{M_{xx}^*(jomega)} ight]_+ + left[frac{S_{yx}(jomega)}{M_{xx}^*(jomega)} ight]_-}$

其中$displaystyle{left[frac{S_{yx}(jomega)}{M_{xx}^*(jomega)} ight]_+}$代表的是因果部分,$displaystyle{left[frac{S_{yx}(jomega)}{M_{xx}^*(jomega)} ight]_-}$代表的是非因果部分。当$H(jomega)$使得积分中的因果部分为0时,就能得到MMSE。此时有

$color{red}{displaystyle{H(jomega) = frac{1}{M_{xx}(jomega)}left[frac{S_{yx}(jomega)}{M_{xx}^*(jomega)} ight]_+ }}$

MMSE为

$displaystyle{MMSE=frac{1}{2pi}int_{-infty}^{infty}left|left[frac{S_{yx}}{M_{xx}^*} ight ]_- ight|^2domega+frac{1}{2pi}int_{-infty}^{infty}left( S_{yy}-frac{S_{yx}S_{yx}^*}{S_{xx}} ight)domega}$

CAUSAL DT WIENER FILTER

Causal DT的求解方式与Causal CT的求解方式基本一样

$color{red}{displaystyle{H(e^{jOmega}) = frac{1}{M_{xx}(e^{jOmega})}left[frac{S_{yx}(e^{jOmega})}{M_{xx}^*(e^{jOmega})} ight]_+ }}$

Dealing with Nonzero Means

前文所讨论都是基于LTI系统,而LTI系统的数学定义如下(CT):

$displaystyle{y(t) = int_{-infty}^{infty}h(t-k)x(k)dk}$

不过,有时随机过程$y(t)$除了跟$x(t)$有关,还有可能有常数偏差,比如我们之前在讨论LMMSE预测的时候会假设$y = ax+b$。对比上面LTI系统的定义,可以发现这种假设并不符合LTI系统,这也是我们在前面的讨论中把随机过程$x(t),y(t)$的mean都设为0的原因。

对于mean不为0的情况,我们可以假设LTI系统的输入是$x(t)-mu_x$,输出的是$hat{y}(t)-mu_y$。

image

为了输出最佳的预测process,此时的LTI系统有如下频率响应

$displaystyle{H(jomega) = frac{D_{yx}(jomega)}{D_{xx}(jomega)}}$

其中$D_{yx}(jomega),D_{xx}(jomega)$分别为covariance $C_{yx}( au),C_{xx}( au)$的傅里叶变换。

Reference:

Alan V. Oppenheim: Signals, Systems and Inference, Chapter 11:Wiener Filtering

原文地址:https://www.cnblogs.com/TaigaCon/p/9290484.html