这篇是理解Kalman filter的最后一道关卡。如果你发现这篇文章读的略吃力，可能你需要配合前面两篇文章一起使用。前两篇文章传送门：
三脚猫Frank：21. 卡尔曼滤波器详解——从零开始(1) Kalman Filter from Zero
三脚猫Frank：22. 卡尔曼滤波器详解——从零开始(2) Kalman Filter from Zero
学完这篇我们就不会因理论基础不够再徘徊在Kalman Filter的的大门之外了。甚至学完前面3篇基础，你都可以比较轻松地入门统计学习或者机器学习的一些内容。说到底还是因为，随机信号处理和参数估计和状态估计等问题本来就是统计学的领域中重要的课题。熟悉统计学基础知识之后，再去学ML或者DL就比较顺利，甚至你会觉得很多人的文章里没有什么特别难的数学，只不过是搞的比较复杂。
本篇我们来讲讲估计理论中的重要部分，也是KF重要的两大理论基础：Theory of Bayesian Estimation（贝叶斯估计理论) 与 Wiener Filter(维纳滤波器)。很多人对Wiener filter知之甚少，殊不知它是Kalman filter的“前身”。没有Wiener等人的工作，是没有后来这么有用的KF的。深入学习之后，我们在下一篇正式研究KF时，你会发现一切推导都是那么自然又简单。
本篇目录
1.Bayesian Estimation 贝叶斯估计
1.1 Understand Bayes' theorem 理解贝叶斯公式
1.2 Why Bayesian estimation 为什么采用贝叶斯估计？
1.3 MMSE estimator 最小均方差估计
1.4 MMSE estimator with Gaussian prior 高斯先验分布下的MMSE估计
1.5 General Bayesian Estimators 一般形式的贝叶斯估计
1.6 Maximum a posteriori estimator (MAP) 最大后验估计
1.7 Linear MMSE Estimator 线性MMSE估计
1.8 Geometrical Interpretations - Orthogornal Projections 几何解释 —— 正交投影
2. Wiener Filter 维纳滤波器
2.1 Some motivations 维纳滤波器的动机
2.2 Model of Measurements 信号测量模型
2.3 FIR Wiener Filter 有限冲激响应维纳滤波器
2.4 IIR Wiener Filter 无限冲激响应维纳滤波器
3. Summary 全文核心内容总结
1. Bayesian Estimation 贝叶斯估计
1.1 Understand Bayes' theorem 理解贝叶斯公式
频率学派处理参数估计问题是建立在参数真值是确定的假设下。而在贝叶斯学派的眼中，待估参数本身也是一个随机变量。要理解贝叶斯估计理论的核心思想，我们必须要对一个核心公式很熟悉，那就是著名的贝叶斯公式（Bayes' theorem，or Bayes's rule）:
(1)这样形式的Baye's theorem是可以从乘法公式(Multiplication rule)中直接得出因为：
当我们采用全概率公式(Law of total probability)改写 ——
发生的概率是在样本空间的一个分割
上的条件概率之和，即：
我们就可以写出样本空间上的一个事件 在
发生时的条件概率：
这就是贝叶斯公式的第二种形式。(1) 非常容易理解，因为 是
共同发生的概率，那么根据(2)我们很自然就有了(1)，从而很自然就有了(4)。
仅从(2)中我们就已经对贝叶斯公式有了基本认识了——共同发生的概率，等同于
先发生，
在
上述公式是离散的，我们要把推广到连续随机变量的应用上。看看(1)就知道我们其实是在研究两个随机变量之间的概率关系。我们可以定义两个随机变量 的联合分布函数和联合概率密度函数(joint pdf)。无论是联合分布函数还是joint pdf，它们都可以用来描述
两个变量一起取某值时的概率。
我们以joint pdf为例，记为 ，这是一个关于
的二元函数。那么什么是条件联合概率密度？那就是当我们确定其中一个变量的值时，另外一个变量的概率密度函数。比如当
已经确定时，
就变成了
，这是关于
的一个一元函数。条件概率密度函数的定义是通过求取极限得到的，我们直接写出最后的定义：
其中 是
的边际概率密度函数(marginal pdf)。由(5)和全概率公式，我们可以推出连续随机变量的贝叶斯公式：
可以看出(6)与(4)的形式非常一致，只不过是把离散的概率由连续的条件概率密度与边际概率密度替代。贝叶斯学派对于这个公式的解读与频率学派不同，这也是贝叶斯推断的理论基础。
无论是采用分布函数还是概率密度函数，我们只要知道贝叶斯公式中每一项的含义就可以了。下面我们都采用密度函数作为标准，因此下面所说的各种概率都是指概率密度。
关于双随机变量之间概率密度函数和分布的关系 可参考这份slides：Results for Two Random Variables，或者参考概率论教材。
回归到贝叶斯估计问题。贝叶斯学派认为未知参数 也是一个随机变量，写出样本和参数的联合概率密度函数
。当参数
已知时，我们有
，意思是样本的概率密度函数是在
取某值的条件下的条件概率密度。这个条件概率函数是样本的函数，考虑
取确定值时样本的条件联合概率。因此根据乘法公式，我们可以重写joint pdf为：
其中 是
的边际概率密度，在贝叶斯估计中被称为prior probability (density function) (先验概率，或者验前概率)，它对应了随机变量参数
自身的一个概率分布，称之为prior distribution（先验分布）。prior probability往往并不容易得知，因为很多时候我们也不知道真实的参数符合何种具体分布和具体分布参数，所以当我们对
情况知晓不多时，我们必须根据问题给出一种假设的概率分布——这被贝叶斯学派认为是一种主观上的对事件发生的可能性判断，即主观概率(subjective probability)。
的取值是最容易获得的，因为它直接来源于样本。虽然我们假设这个条件概率是一个样本的函数，
那一刻是确定的。但是当我们有了具体样本之后，代入
，我们就得到了一个关于
的函数。代入样本后的
正是我们在最大似然估计中利用
求出的似然函数(Likelihood function)。
在经典估计中表示的仍然是样本之间的联合概率密度，而只是把
当做一个预先确定的值。样本的产生是因为
的含义实际是等价的，因此我们说，代入样本后，它正是Likelihood function。
那么根据乘法公式我们又可以写出：
其中 是样本的边际概率密度，意思是说计算当
取所有可能值时，样本的联合概率密度。而
是出现某一组样本时， 此时
的条件概率密度。令(7)与(8)相等，我们得到了：
被叫做posterior probability(后验概率，或者验后概率），对应的分布叫posterior distribution（后验分布），它的含义是在观察到样本后我们得到的
的分布，与prior probability是相对的。
贝叶斯定理或者贝叶斯公式是Bayesian Inference的核心依据。从下面这张图（采用了离散的形式）中我们可以看到每一部分的组成，按照图我们再次解读一下：公式的左边要计算的是，一个假设的概率，不过它是一个条件概率，是得知某事件的观察证据
的情况下假设的
是我们假设出来的
是在这个
，
1.2 Why Bayesian estimation 为什么采用贝叶斯估计？
经典估计中采用确定性参数的假设有时并不一定合适，贝叶斯估计能够让我们把对参数的prior knowledge给加入到估计之中，去提高准确率。并不是所有的参数取值的可能性都是一样的，就比如你觉得上海有lady M的可能性大，还是十八线小城市有的可能性大？我们据主观初判可以给我提供一个合理的可能性假设。当然前提是我们如何得到prior distribution，这是一个重要课题。
在Bayesian estimation中我们能够找到最小化MSE的estimator，即MMSE estimator。而在classic estimation中我们在上篇中已经讨论过，一般是不存在MMSE的，只有可能找到UMVUE。还有很多特定问题采用贝叶斯估计有它们各自的理由，就不继续展开了。
1.3 MMSE estimator 最小均方差估计
我们在上一篇中详细介绍过在频率学派中使用的指标MSE（Mean Square Error）。当时我们得出过一个结论：一般情况下（指没有限定估计类型），是不存在最小化MSE的estimator存在的。但是在Bayesian estimation中，这样的MMSE estimator却是存在的,而且总是存在！在Bayesian estimation中有许多的优化指标可以使用，而MMSE是最常用的。
设随机变量参数 估计量为
，其MSE为：
这里采用了乘法公式把联合概率密度给拆成条件概率和边际概率了。在经典估计理论中，我们认为 其实是
。为了最小化(10)，我们看到
是一个只关于
的条件概率, 令花括号中的部分对
求偏导,并使其为0：
注意其中不是
我们得到令MSE最小化的解：
所以(12) 表明在MMSE这个标准下，最优估计是参数 的后验分布的期望，即mean of posterior distrubution。这是一个非常重要的结论！
试想在观测任何样本之前，经典估计无法给出估计结论。但贝叶斯估计中会把最开始的先验分布的均值当做是MMSE估计。而后如果加入了样本，我们就可以根据贝叶斯公式计算出后验分布，则最优估计就变成了（12）所列的。
所以这就完事了吗？ 谁告诉我一下 怎么求？ 上面已经说了，你需要一个后验分布，或者说后验概率函数
。虽然Bayes' theorem能够帮我们计算，但返回公式去看，我们需要
，即likelihood function。因此这种最优估计必须要求我们对概率模型的形式有一个清楚的了解。问题是就算得到了
，我们计算：
最后积分的结果会是closed form 吗？ 这暗示最后解的形式可能非常复杂。
1.4 MMSE estimator with Gaussian prior 高斯先验分布下的MMSE估计
先给出结论：若噪音和参数都满足高斯分布，那么(13) 具有closed form。
这里给一篇我认为写的非常不错的对于多元高斯分布的解释。那我就默认你已经熟悉多元高斯分布了。现在考虑样本 ，满足联合高斯分布，其joint pdf为：
其中 是样本
协方差矩阵，
是均值向量。我们知道
是充分统计量，完全描述了这个高维分布。
还记得协方差矩阵长什么样吗？方阵，对角线是各自分量的方差，其他为两两之间的协方差
。
我们还没有加入参数 到联合分布之中。现在让我们假设参数
是也是符合高斯分布的。首先要解决一个问题，
是符合联合高斯分布的，而
也符合联合高斯分布吗？ 答案是：不一定！
当 两个变量各自都是高斯分布的时候，它们并不一定是联合高斯分布！可以参考Berkeley EE 223 lecture note中的一些重要结论和例子：
既然不是在任何时候高斯边际分布可以反推出联合高斯分布，我们自然要问什么时候可以这么做？下面这个定理阐述了一个事实：
Theorem 1: Real valued random variablesare jointly Gaussian iff for all
, the real R.V.
is Gaussian.
这意味着联合高斯分布自动imply了它们的任意线性组合也是高斯分布；反过来，只要任意线性组合是高斯分布，那么它们之间就是联合高斯分布。这是一个非常重要的结论。[证明请参考EE223中的叙述。]
下面给出不加证明的两则定理。它们告诉了我们，如果样本和参数的确是满足jointly Gaussian的，那么后验概率密度函数有closed form。下面两则公式的结果，都可以自己很容易地推导出来，因此在此省略了。
Theorem 2(二元高斯分布): Let X and Y be random variables distributed jointly Gaussian with mean vectorand covariance matrix：
posterior pdfis also Gaussian with mean and variance given by:
Theorem 3(多元高斯分布): Let(k×1) and
(l×1) be random vectors distributed jointly Gaussian with mean vector
and covariance matrix
posterior pdfis also Gaussian with mean and variance given by:
我们能看到 MMSE的解 是可以写成closed form的，因为
和它们的方差都是来自高斯分布，是常数。
根据这两个定理，我们可以为下面这个线性量测模型在Bayesian estimation的框架下寻找MMSE estimator：
这个量测模型我们在上一篇文章中就已经提及了。 是
样本，
是已知矩阵，
是未知参数
随机向量，
是高斯白噪声服从
。现在在贝叶斯估计的框架下，我们有
服从高斯分布
，
为均值向量。我们发现：
与
它们是jointly Guassian的（什么原因，见Theorem 1）。这时根据 Theorem 2，我们就可以直接写出该参数向量的MMSE估计:
MMSE估计为：
估计的后验协方差为：
这张图大概标注了每一项的意义，其实和我们的recursive least square算法(上一篇中提到的)与Kalman filter的方程结构已经非常接近了。这里都是将预测误差乘以某个增益后，添加到上一次的后验（即这一次的先验）上去更新下一次的后验。
所以这里的关键就是jointly Gaussian distribution可以得到一些closed form solution of MMSE。但是一般的联合分布就不好说了，通常难以得到一个封闭解。
1.5 General Bayesian Estimators 一般形式的贝叶斯估计
Square Error只是其中一种cost function 或者 loss function。Loss function是把参数和估计量一起映射到一个实数，因为一般由被估计参数 与采用的估计量
之间的函数关系组成。
Loss function本身并不一定是个确定的值，就像这里源于样本 具有随机性。因此常常会采用expected loss function，或者叫risk function，来作为评价
这个选择的好坏，或者在决策理论中，我们也把这个称之为决策
。
频率学派常用的risk function为：
此时的 是一个确定的值，只不过有一个取值范围。
下面的
就是表示在此时的
的值。而整个积分的被积变量是
，是每一个观察值，因此这个期望是求取了整个样本空间上的loss function的均值。这就碰到了一个问题，risk function在不同的
取值时会有不同的表现——就像我们上一篇中提到的此时不存在MMSE estimator一样，对有些参数，某些estimators会表现更好，而其他参数，它们的risk值可能反而比较差。有很多妥协的办法可以考虑，比如限制estimator的类型，选择“worst case is the best”等等。
Bayesian对loss的处理考虑到了参数本身也具有随机性，而且不是所有参数取值可能性都一样。他们仅考虑在已经观察到的样本上去对参数 的概率函数求期望，定义posterior risk：
其中 在这里是参数
的prior pdf，
就是在
时
的后验概率函数。比较(18)和(19)，它们的不同在于，Bayesian risk是一个样本条件期望，在观察到的样本下对分布中所有可能的参数求平均，前者确是对样本空间求平均。因此(19)看起来照顾到了不同参数对loss的影响，而仅考虑观察到的样本对推断的影响，不考虑其他“没关系”的样本。在(10)中应该已经能看到(19)的影子了。
那么贝叶斯估计中，我们设计一个posterior risk之后，最佳估计就应该使得posterior risk最小化。比如，选择square error loss function，我们得到：
最小化(20)，我们已经在(11)中得到了结果，就是 。这是后验概率函数
的mean。
很有意思的是，当两种不同观点相互结合的时候，我们还可以定义出别的risk function来，其中一种就是Bayes Risk。
所以Bayes risk正是上面的posterior risk对样本 求期望。我们的(10)中求MMSE就是采用了Bayes risk。
现在考虑另外两种loss function，一种是绝对误差，一种是Hit-or-Miss：
利用这两种loss function可以得到最优Bayesian estimator分别是：
具体推导可以参考[5]。
1.6 Maximum a posteriori estimator (MAP) 最大后验估计
我们在最大似然估计中考虑的是使得likelihood function最大化时的参数 ，在贝叶斯学派是想要让posterior pdf最大化，即最大化
。因为贝叶斯公式中分母是与参数无关的，所以最大化后验概率函数等同于最大化
。关于MAP我在这里不想过多的展开讨论，本身概念也十分清晰和简单。关键是要拥有概率密度函数以及prior pdf。我们在上面采用Hit-or-Miss的loss function时就已经得到了MAP的结果，它就是取的后验概率的最大值（mode of probability density function）作为估计。
1.7 Linear MMSE Estimator 线性MMSE估计
在1.3节中我们给出了MMSE estimator是后验分布的期望。但是期望的求取可能涉及到复杂的积分过程。1.4节中我们发现只要后验分布是高斯的，那么显然我们能得到closed form solution。那么如果我们不假设高斯分布，是否依然能够的导closed form solution呢？答案跟我们上一篇文章中的BLUE一样，要满足线性估计，并且需要知道随机变量的一阶矩和二阶矩。
我们考虑只有一个参数 时的标量情况.
假定我们有 个样本序列
，其中的
满足某个概率密度函数但是是未知的，因此它与样本的联合概率密度
也是未知的，但是
与
的一阶和二阶矩是知道的。我们想要构造一个线性估计量：
使得估计量的MSE是最小化的，试找到这样一个估计。
直接列出：
第一步先求出最小化时的 。就直接令(23)对
求偏导并为0。得到：
第二步把(24)代回(23)得到：
我们要寻找使得(23)最小化的 ，所以我们整理式子（展开上面的平方，整合）：
其中 。
是
与
之间的协方差向量（或行向量）。
是
协方差矩阵，
是参数的方差。令它对
求偏导，求导的时候注意
因为是实数，所以也可以把两个向量互相转置
：
把它代回后，我们得到Linear MMSE：
把解带回去就可以得到此时的MSE，这里就不计算了。
当有 个参数
时，我们可以得到Vector form LMMSE：
是一个
互协方差矩阵。差别就在于这个矩阵上。此时最小MSE为：
(28) 的结果并不依赖于 与
的关系，而是限定了
与
的线性关系，因此是一个比较general的结果。如果我们采用线性量测模型，即假定参数与观察样本之间也是线性关系，那么我们就有类似于上一篇中BLUE的Bayesian Gauss-Markov Theorem:
Bayesian Gauss-Markov Theorem: 对于线性量测模型，其中
, 协方差矩阵为
,噪声为零均值，协方差矩阵为
。（无需要求具体分布类型）则线性最小均方差估计和其协方差由下式给出：
这里要求不仅是线性估计，而且要求量测模型也是线性的。但我们在本节中的结论(28)是没有限定量测模型类型的。这是差别。
1.8 Geometrical Interpretations - Orthogonal Projections
LMMSE的几何解释 —— 正交投影
我们可以从几何的角度去考虑得到Linear MMSE estimator。实际上在上一篇中我们提到的ordinary least square(最小二乘法)的计算公式也可以通过几何的方法得到。让我们以R.E. Kalman自己在[1]中所写的一点内容来介绍Orthogornal Principle。我会修改一下它paper中的符号（因为我不是很喜欢他文中的符号）。
考虑我们观察到一个有限数量的，均值为0的随机过程 的样本，一个时间从
到
的信号序列
，其中每一个时刻的观察值都是随机变量
。现在我们令
所有的线性组合构成一个vector space，让我们把它记为
。里面的任意一个向量
都可以表示为:
由于样本观察的数量是不定的，我们认为 是所有可能观察到的样本（目前还没观察到的）的线性组合成的向量空间的子空间。
的维度是有限维，维度为观察到的样本的数量。而回顾我们所谓的估计量，在满足线性估计的假设后，就是所有已得到的观察
的线性组合，因此它属于
。这里的向量空间和以下的定义都是抽象的，但符合向量空间和相关定义的。
对于任意给定的两个向量 ，我们可以定义它们的内积为：
如果内积为0，我们认为它们是正交的(orthogonal)，即 。现在考虑任意一个均值为0的随机变量
（它不一定在
中）。根据正交分解定理（Orthogonal Decomposition Theorem）（无论是linear space还是Hilbert space），随机变量
可以表示为两个部分，一部分为封闭子空间
中的向量，另一部分是正交于
的向量。
其中 。由于
是正交于整个
子空间的，因此它也正交于
中任意一个向量
，即：
自然满足 。(31)就是我们的orthogornal principal，可以解读为：估计误差
总是与用于构成估计量的观察样本
正交。此时我们有
是
在
上的投影(projection)，并且如果我们计算
与任意
之间的平均距离：
因此在 中只有
的投影
是与
的距离最小的，这正是我们要寻找的Linear MMSE estimator。由于
均值为0，所以线性估计的形式没有bias项（以下
均为向量）：
我们根据(31)知道 满足Orthogonal principle：
因此：
这个结果与(26)是一致的。我们就得到了均值为0时的LMMSE：
这就是根据正交投影得到的结果，可以很方便地推广到多个参数的情况，即和上面(28)一样，只不过是均值为0，故比较特殊。
(34)中不光满足这个正交条件，而是所有
和它们的任意线性组合都正交于误差。
2. Wiener Filter 维纳滤波器
2.1 Some motivations 维纳滤波器的动机
Norbert Wiener这个名字，对于学control的人来说应该都不陌生。一个令他出名的原因是他的著作：cybernetics（控制论）。另外一个原因就是Wiener-Kolmogorov filter。Wiener在第二次世界大战期间为军方设计高射火炮控制器。系统必须要在雷达信息的配合下预测飞机的位置，雷达信息又充满了噪声。他利用统计学和概率论的知识将信号作为随机过程来研究，利用信号和噪音的自相关函数来获得了least mean square error, LMSE（即minimum mean square error，MMSE）意义下的最优预测。他的研究工作直到战后才被解密，被记录在了《Extrapolation, interpolation and smoothing of stationary time series》一书中。当时苏联数学家Andrey Kolmogorov在1941年也推导出了离散系统的最优线性predictor。
所以从当初的动机来看，Wiener filter最初目的是信号预测（prediciton），而非滤波（filtering）。不过我们马上就可以看到，其实我们可以用之前已经推导过的LMMSE的结果，直接解决filtering, smoothing以及prediciton三种信号处理问题，从而使得Wiener's problem得以解决。
我们之前3篇，包括这篇，说了这么多基础知识，最后就是要在Kalman filter之前引出Wiener filter（维纳滤波器）。我们讨论了如何最小化MSE得到最优估计，而Wiener filter就是这样一个optimal filter。很多人学Kalman filter之前可能听说过，但从来没有认真研究过Wiener filter。今天我们就来细细地来讲一讲。等下一篇我们剖析Kalman filter后，再回过头来看看它们之间密切的联系。
2.2 Model of Measurements 信号测量模型
让我们考虑一个通用的信号估计问题，也就是Wiener filtering problem。我们考虑一系列信号序列，首先定义一系列符号：
第n个受到噪声干扰信号(已测量)：
第n个噪声采样：
第n个真实的待估计或者还原的信号(未知)：
通过最优滤波器后的输出信号(第
个）：
真实的个待估计或者还原的信号(未知)：
第n步估计误差：
最优滤波器的脉冲函数：
问题阐述： 真实的信号 受到噪音(假设均值为0)污染，因此我们只能测量到
（因而
均值也为0）。Wiener filtering的目标是设计一个linear optimal filter，使得给定过去的测量序列
通过optimal filter后得到第
步的估计值
，使其MSE的均值最小，即求MMSE estimator。设计linear optimal filter等同于设计其impulse response function（脉冲函数，即传递函数的Laplace逆变换。）或者传递函数，即设计函数
或者
。
s[n]和v[n]均值为0，是下面运用orthogonal principle的重要假设。对于非零均值问题，我们可以令信号减去均值后，再通过Wiener filter，最后加回均值，因此均值为0不是问题。
假设： 序列是WSS过程(宽平稳)
满足jointly WSS(联合宽平稳过程，即联合概率函数有常均值，自相关函数仅与时间延迟有关，见我写的本系列第一篇文章。)
问题分类：
当 时：filtering problem。我们要通过当前的测量值
和过去的所有测量数据，估计当前的真实信号值
。其估计为
。
当 时：prediciton problem。我们要通过当前的测量值
和过去的所有测量数据，估计未来的真实信号值
。其估计为
。filter必须是causal的。
当 时：smoothing problem。我们要通过当前的测量值
和过去的所有测量数据，来估计过去的真实信号值
。其估计为
。filter可以是noncausal的。
如果不清楚这个问题的定义，请返回我的系列文章的第一篇中。我这里不多讲了。讲一下这里的causal与noncausal，如果一个filter的输入受到未来的输出的影响，那就是noncausal的。我们一般的系统都是causal的。但是noncausal的实现，只要我们改变时间起点就可以——比如smoothing problem，我们站在过去的时间点上看，那么未来的观测值和这个时间点过去的观测值都将可以使用，因此在smoothing problem中我们可以使用noncausal filter。
我们下面就按照FIR filter和IIR filter分别讲讲Wiener filter的推导过程和重要结论。
2.3 FIR Wiener Filter 有限冲激响应维纳滤波器
2.3.1 Causal FIR Wiener filter
FIR filter（Finite impulse response filter, 有限冲激响应滤波器) 是一种常用的数字滤波器。它将过去有限个信号作为输入，输出它们的线性组合，因此是一种线性滤波器。我们考虑filtering problem，即 。我们现在设计一个滤波器，它的输出为现在和过去有限个( 总共
+1个)测量信号的线性组合，且
时
，即causal
阶FIR filter：
时域中impulse function和输入的卷积为线性系统的输出，还记得吗？
假设 ，我们就有
阶FIR filter：
我们希望 和
的误差能够满足minimum MSE，于是我们定义：
好了，我们不用继续算下去了。因为无论是(28)的结论还是orthogonal principle (31)都已经可以直接解答这个问题了。我们试着用orthogonal principle来解答这个问题——估计误差与所有测量值是正交的。选择任意一个测量值，记为 ：
根据cross-correlation function(互相关函数）的定义，(39)为 与
的cross-correlation
。我们展开(39)：
或者说我们从orthogonal principle 得到：
因为
这里代表了离散卷积。
因为 是任意的，对于(42)代表了一系列的等式：
我们把这个结果整理成矩阵：
注意，autocorrelation是偶函数，所以
。
注意我这里是阶FIR，很多书和资料喜欢用
阶，因此(43)看起来会少一维。
我们把方程：
称之为Wiener-Hopf Equation，其矩阵形式如(43)。求解此方程就可以得到此时的optimal filter ，它是由
这几个系数，以及系统的阶数
所决定的。我们可以选择随时间不断增加FIR的阶数，如此FIR Wiener是一个time-varying的filter，滤波器系数序列
是时变的。(43)中记左边的方阵为
，它是一个symmetric 并且Toeplitz的矩阵，因此(43)或者(44)可以采用Levinson–Durbin recursion算法求解。当然我们也可以固定只用最近的
个观察来构建滤波器，这样假设我们准确地知道
那么方程的解就是固定的，我们只需要直接解一次(43)就行了。
采用(27)的LMMSE的结论也可以帮我们得到Wiener filter：
（注意 均值为0）。如果我们把(45)与方程(43)一对照就会发现，其中的：
(27)的解形式，正是对应了Wiener-Hopf Equation的解，其中 是
的autocorrelation matrix，
是cros-correlation vector。
注意(45)中根据计算出来的的系数
的顺序要颠倒过来，这样定义出来的就是
，使得输出是它们之间的卷积。
那么此时的MMSE是多少呢？ 留给你计算了。
2.3.2 Non causal FIR Wiener filter
我们让Wiener filter成为non causal的，即filter的输出不仅依赖于过去和现在，还取决于未来。令 阶 non causal FIR filter的impulse response
满足:
其中整数 。这样的设定意味当前的输出还取决于未来的输入
。我们回顾一下上面的步骤，其实这影响了方程组中等式的个数。
想象一下做smoothing的时候，我们站在过去的时间点上看，未来和过去都是可以利用。所以这里的“未来”并不是真的是未知的，而是相对于某个时间点而言的未来，因此虽然non causal的filter听起来不可实现，但实际上在特定场景下是可以实现的。
2.4 IIR Wiener Filter 无限冲激响应维纳滤波器
有了FIR之后，我们再看IIR filter就很容易了。IIR就是采用无限的观察值来构造下一个输出。同样的，我们可以分为causal和non causal IIR，同样的causal是主要是为了filtering和prediciton，non causal IIR是为了smoother。
我们这次先考虑non causal IIR，令滤波器的输出为：
我们直接写出Wiener Hopf Equation，
我们可以把(49)写成积分形式，而左边正是无穷时域离散卷积：
如果我们把(50)或者(49)进行两边 变换：
其中 是
域中的互谱和自谱。诶，怎么有点熟悉？我们在第一篇中就讲过传递函数可以通过PSD和CPSD来求解，这里的Wiener filter正是这种形式，并且告诉了我们这是一种最优的Linear filter，满足假设条件的情况下可以最小化MSE。换句话说，所谓的Wiener filter实际上构造了一个从
到
的传递函数！
注意，这样得出的(51)的 是non causal的，因为进行了双边
变换。如果有必要，我们也可以用two-sided Fourier transformation来替换z-transformation。同理Laplace变换也可以代替。实际上可以证明连续频域中也有：
详细请见[1]。
Causal IIR是做prediciton必须满足的。它的推导略微有点复杂，我们令滤波器满足：
即过去所有的 都被用于滤波器构建。我们先假设
输入本身是白噪声，均值为0，单位方差，其自相关函数是
。于是有Wiener Hopf Equation如下（无限个方程的组合）：
而对于一个最优的noncausal IIR，根据Wiener Hopf方程——如果输入是白噪声，同样有：
我们发现如果 本身是白噪声的话，causal和non causal IIR的结果是一样的！只要把non causal IIR的non causal部分去除掉就行了。
因为我们认为 的随机性是由于噪声
的随机性造成的，信号
不具有随机性，所以如果滤波器的输入是白噪声的话，我们就可以认为causal IIR和non causal IIR的结果是完全一样的。job is done！
当滤波器的输入不是白噪声时，我们就要将输入进行whitening（白化处理）。我们在第二篇文章中已经讲过了。我们这里在这里讲大致思想：就是构造一个白化滤波器使得一个已知频谱的信号通过后能输出白噪声。我们需要知道输入的频谱 ，然后进行谱分解得到：
取稳定部分，我们得到白化过滤器： 。但这样的滤波器只是假想出来的，因为本身
是stable和causal的，而其倒数就不是causal的。原来的信号
可以由
表示。
当 通过
变会成为白噪声，记为
。红框里的部分都是causal的，是我们要求的causal IIR
与
的乘积。由于红框的输入是白噪声
，因此它的non causal IIR和causal IIR的结果是一致的，只要把non causal IIR的non causal部分去除掉就可以了。我们假设对
设计最优的non causal IIR (50)，可以直接使用结论：
让我们把 做逆z变换到时域，然后把non causal部分去掉，再使用z变换到频域得到，记为：
（这个操作没法用表达式，只能描述）。
就代表了红框中
和
的乘积。因此我们要求的causal IIR Wiener filter就是：
3. Summary 全文核心内容总结
贝叶斯估计的核心思想还是利用了贝叶斯公式中，后验概率函数和先验概率函数、似然函数之间的关系。我们已经知道，在贝叶斯框架下，最小化MMSE的estimator就是 。选择不同的loss function，会让我们最优化的estimator变得不同。对于MMSE estimator，我们讨论了当满足Gaussian的时候，解有封闭形式；或者满足线性估计的时候，我们也能得到封闭形式的解。在分布满足Gaussian时，线性估计器打败一切非线性，是最优的，这个结论应该是不仅限于quadratic loss function的（这个在Kalman 1960年那篇文章中有讨论）。当分布类型不是Gaussian时，只要限定在线性估计器，我们依旧能找到封闭形式的解，并且这与分布为Gaussian时的结果一模一样！（参看(28)和Theorem 3的结果）。LMMSE estimator和正交原则奠定了Kalman filter和Wiener filter的理论基础。
Wiener filter是在Kalman filter发现之前的一种以最小化MSE为优化标准的最优滤波器。用orthogonal principle或者LMMSE的结论，我们可以很方便地推出Wiener Hopf方程。我们讨论了几种不同的实现，满足因果和非因果的FIR和IIR滤波器。不过Wiener filter有其自身的短板，我们下一篇将会正式开启讨论Kalman filter，届时会进行比较。粗略地讲，Kalman filter就是一种可以recursively执行的，结合了线性系统动态方程的Wiener filter。
Reference
[1]MIT Introduction to Communication, Control, and Signal Processing 6.011
[2]Uppsala University：Optimal filter(1)
[3]Uppsala University：Optimal filter(2)
[4]A Brief Introduction to Hilbert Space
[5]SUNY-Binghamton： Lecture Notes
[6]Stanford EE264: Wiener Filtering
评论（0）
您还未登录，请登录后发表或查看评论