卡爾曼濾波器詳解——從零開始(3) Kalman Filter from Zero

來源 · 知乎

作者 · 三腳貓Frank

編輯 · 睿小妹

原文 · https://zhuanlan.zhihu.com/p/245728093

這篇是理解Kalman filter的最後一道關卡。如果你發現這篇文章讀的略吃力，可能你需要配合前面兩篇文章一起使用。前兩篇文章傳送門：

卡爾曼濾波器詳解——從零開始(1) Kalman Filter from Zero

卡爾曼濾波器詳解——從零開始(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 全文核心內容總結

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)的形式非常一致，只不過是把離散的概率由連續的條件概率密度與邊際概率密度替代。貝葉斯學派對於這個公式的解讀與頻率學派不同，這也是貝葉斯推斷的理論基礎。

無論是採用分布函數還是概率密度函數，我們只要知道貝葉斯公式中每一項的含義就可以了。下面我們都採用密度函數作為標準，因此下面所說的各種概率都是指概率密度。

回歸到貝葉斯估計問題。貝葉斯學派認為未知參數  也是一個隨機變量，寫出樣本和參數的聯合概率密度函數  。當參數  已知時，我們有  ，意思是樣本的概率密度函數是在  取某值的條件下的條件概率密度。這個條件概率函數是樣本的函數，考慮  取確定值時樣本的條件聯合概率。因此根據乘法公式，我們可以重寫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的核心依據。從下面這張圖（採用了離散的形式）中我們可以看到每一部分的組成，按照圖我們再次解讀一下：公式的左邊要計算的是，一個假設  的概率，不過它是一個條件概率，是得知某事件的觀察證據  的情況下假設的  概率。那麼這個概率怎麼算？公式的右邊，分子中的概率  是我們假設出來的  的概率，稱之為Prior，這個概率是在觀察到更多證據  之前的假定概率。  是在這個  假設下，事件發生收集到  的概率，稱之為Likelihood。分母則是Marginal Probabiliy，是在所有假設下(包括  在內的)事件的證據  被觀察到的概率。其實仔細一看分子的乘積不就是乘法公式裡的分子嘛，代表了  ，  也發生， 且證據 也被觀察到的概率。那麼根據乘法公式，我們自然可以算得，在證據e發生時，我們的假設  成立的概率。這是個很有意思的事情，因為我們往往容易拿到證據發生的概率，比如用頻率法取近似一個概率，然後如果我們假設一個  ，再計算所有可能導致e發生的概率P(e)，就可以推出Posterior，所謂後驗概率，它實際上是站在觀察到的證據的基礎上，更新了原來Prior中假設的主觀概率。可以說，我們更關心的是Posterior中的  ，prior可能是在信息不足的情況下的一個估計概率，但是隨著證據的增多，我們會發現Posterior的可信度會越來越接近事實。因此總結來說，Bayes' Theorem是計算觀察到證據的情況下，假設  的概率分布，做法是先自己假設一個prior，然後通過計算假設成立下證據被觀察到概率，除以證據被觀察到的所有概率（所有可能來源，包括假設 和其他假設 )。

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為：

其中  是樣本  協方差矩陣，  是均值向量。我們知道  是充分統計量，完全描述了這個高維分布。

還記得協方差矩陣長什麼樣嗎？方陣，對角線是各自分量的方差  ，其他為兩兩之間的協方差  。

我們還沒有加入參數  到聯合分布之中。現在讓我們假設參數  是也是符合高斯分布的。首先要解決一個問題，  是符合聯合高斯分布的，而  也符合聯合高斯分布嗎？ 答案是：不一定！

當  兩個變量各自都是高斯分布的時候，它們並不一定是聯合高斯分布！

既然不是在任何時候高斯邊際分布可以反推出聯合高斯分布，我們自然要問什麼時候可以這麼做？下面這個定理闡述了一個事實：

Theorem 1: Real valued random variables  are jointly Gaussian iff for all  , the real R.V.  is Gaussian.

這意味著聯合高斯分布自動imply了它們的任意線性組合也是高斯分布；反過來，只要任意線性組合是高斯分布，那麼它們之間就是聯合高斯分布。這是一個非常重要的結論。

下面給出不加證明的兩則定理。它們告訴了我們，如果樣本和參數的確是滿足jointly Gaussian的，那麼後驗概率密度函數有closed form。下面兩則公式的結果，都可以自己很容易地推導出來，因此在此省略了。

Theorem 2(二元高斯分布): Let X and Y be random variables distributed jointly Gaussian with mean vector  and covariance matrix：

posterior pdf  is 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 pdf  is 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的結果。如果我們採用線性量測模型，即假定參數與觀察樣本之間也是線性關係。

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.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 有限衝激響應維納濾波器

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)中根據計算出來的  的係數  的順序要顛倒過來，這樣定義出來的就是  ，使得輸出是它們之間的卷積。  的值就是filter的impulse response，正是因為其值是有限個的，所以才叫FIR。
那麼此時的MMSE是多少呢？留給你計算了。

我們讓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就是：

貝葉斯估計的核心思想還是利用了貝葉斯公式中，後驗概率函數和先驗概率函數、似然函數之間的關係。我們已經知道，在貝葉斯框架下，最小化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濾波器。

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

[7]GaTech: Bayesian Statistics: Handouts

[8]UC Berkeley: EECS260 Lecture Notes

如何領取參考資料

1) 關注睿慕課公眾號

2) 回復「卡爾曼」獲取下載連結