假设分别有两个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。
本文的主要目的就是求Wiener filter的脉冲响应，所求得的脉冲响应需要使得系统有Minimun MSE。
※在开始之前，先说明本文一些计算所需的前提条件：
- 由于输入以及输出的随机过程都是实数，因此$R_{xx}[m],R_{yy}[m],R_{xy}[m],R_{yx}[m]$也会是实数
- 由于$R_{xx}[m] = R_{xx}[-m]$，因此$S_{xx}(e^{jOmega})$是实数，同理有$S_{yy}(e^{jOmega})$是实数
- 本文假设$S_{xx}(e^{jOmega})$以及$S_{yy}(e^{jOmega})$在所有频率上都大于0
- 为了方便，本文假设随机过程$x[n],y[n]$的mean都为0
- 连续时间随机过程有上述同样条件
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一样对某一个采样点求导来得到最佳系统，需要采用其他计算方式。
与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$。
为了输出最佳的预测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