ML Estimation and MAP Estimation 差異

剛才看了 BAYESIAN ESTIMATION,   這篇文章把 MAP, ML, MMSE, MAVE 解釋和比較的很清楚 under the Bayesian estimation framework.  可以直接參考跳過本文。(Chap4 of “Advanced Digital Signal Processing and Noise Reduction”, Second Edition. by Saeed V. Vaseghi)

ML (Maximum Likelihood) and MAP (Maximum A posteriori Probability) estimation/detection 非常類似，也容易混淆。其實只有一點差異。

(2.1) 是 MAP detection;  (2.3) 是 ML detection.   MAP 和 ML 之間的關係是 (2.2) (Bayes theory). 

簡單來說，ML (eq. 2.3) 和 s 出現什麼的機率無關 P(s was sent).  意即完全不考慮 s 的 pdf; 或者假設所有的 s 有相同的機率。在這情況下，eq. 2.2 and eq. 2.3 are equivalent; ML estimation 等同 MAP estimation.  (P(v is observed) is a constant in eq. 2.2).

In summary,  ML (or MLE) 和 MAP estimator 的關係可以更簡化如下。唯一的差別就是 MAP 多乘了 theta 的  distribution.

 

舉例而言，s=+1 or –1 is a simple Bernoulli distribution. 

P(s=+1)=p  P(s=-1)=1-p=q           n ~ N(0, sig^2)

ML estimator:  s* = arg max P(v is observed | s was sent)

MAP estimator: s* = arg max P(v is observed | s was sent) * P(s was sent)

兩者的差異如下圖:

ML estimator:  v>0  s*=+1;  v<0  s*=-1   regardless p

MAP estimator: v>a  s*=+1;  v<a s*=-1   a depends on p (can be computed numerically)

如果 p=0.5,  a=0 and equivalent to ML estimator.  如果 p>0.5, a 就是負值，  s*=+1 的區間變大，反之亦然。直覺來說，p>0.5 代表 P(s=+1) 的機率比較大，因此在 estimate 時，比較 favor s=+1, 也就是 threshold a<0.

MMSE estimator 和 ML estimator 差異

ML estimator 是 Maximum Likelihood estimator (最有可能)

MMSE estimator 是 Minimum Mean Square Error estimator (最小平方誤差)

嚴格來說，ML estimation 是  classic estimation.  MMSE 是 Bayesian estimation.  請另文參考。本文主要是以應用來比較。

這兩種 estimator 在實際應用上常遇到。 到底差異在那? 以及 when to apply what?
ML example: match filter (maximize decision SNR), digital communication detection (minimum BER- Bit Error Rate), etc.    ML 一般是在 no prior information 時的最佳選擇。

MMSE example: Wiener filter (minimize error, maximize total SNR), Kalman filter (minimize error), etc.   MMSE = min E(|x-x*|^2) 直譯為 minimum mean square.  通常 reserve 給 Bayesian estimation quardratic cost function with prior information.

簡單來說，ML estimator 主要應用是 estimate parameters based on what is most likely to create the oberved samples (但 parameter 本身 no prior information).

MMSE estimator 主要差別是 estimate waveform (Wiener filter) or sequences (Kalman filter) based on samples to minimize a cost function (mean square error function).

當然以上只是粗略的區分。也有 sequence ML estimator (Viterbi decoding algorithm 就是有名的例子)。反之 MMSE 也常用來 estimate parameters, 請 reference 前文 “Kalman filter: MMSE estimator”.

簡單例子 (類比通信和數位通信)

類比通信如 AM 或 FM, 接收時會被 noise 污染，理想的類比接收器就是 MMSE filter (如 Wiener filter), 儘可能還原原始的類比信號。實際上為了成本的考量，一般都用簡單的 LPF.  ML estimator 在此就無用武之地。

數位通信如 PSK 或 FSK, 接收時也會被 noise 污染，但理想的接收器不是儘可能還原原始信號。如下圖所示，重點在 maximize decision moment 時的 SNR, 才能準確判斷收到是 0 or 1. 即使 waveform 有 distortion 也無妨，因為數位通信並不 care waveform. 理想的數位接收器就是 ML estimator (match filer), 能準確判斷收到(或送出)是 0 or 1.   MMSE estimator and ML estimator 在這個例子顯然不同。

一個反例 (MMSE and ML estimator 相同結果): Linear Regression

==> 以下的例子並非是 MMSE, least square function or least square estimation (LSE) .  常常容易混淆 MMSE and least square optimzation, 以下誤導大家。

Least square estimation (LSE) and ML estimation 的差異可見本文。

在 linear regression estimation, the goal is to estimate y = b_0 + b_1 x 中的 b_0 and b_1. 

注意這是 estimate parameters b_0 and b_1, 不是如 Kalman filter estimate waveform.

先 define square error function 如下:

 

 

ML 的方法是先 define likelihood function 如下:

Then find what theta is most likely given Xi, i.e. find maximum of likelihood function.

更常用的是 log-likelihood function, take log of the likelihood function 再微分。

 

Assuming epi_i is Gaussian with i.i.d:

How to maximize the likelihood function?  It tunrs out equivalent to minimize

 

以下為錯誤結論 :

In this case, MMSE and ML estimation are the same!!

到底  MMSE 和 ML 在 parameter estimation 有什麼不同，還是每次都會有相同的結果?  從上面的例子，MMSE estimator 和 error (epson_i) 的 pdf 沒有關係，但 ML estimator 明顯 depends on pdf. 

舉例而言，如果 eps_i 不是 Gaussian, 而是 exponential distribution, ML 要 maximize likelihood function 就不會 equivalent to minimize “mean absolute error” (沒有 square). 

In general, MMSE estimator 不管 error 的 pdf, 基本上只 care 1st order and 2nd order statistics (mean and variance).   換言之，就是假設 error pdf 是 Gaussian, MMSE estimator 就是 optimal estimator and achieve the Cramer-Rao lower bound (CRLB).

反之，ML estimator 會針對 error pdf (prior information) 做 optimization.  理論上會比 MMSE 更優?  如果 error pdf 是 Gaussian,  MMSE and ML 就會得到相同的結果。

以下為正確的結論:

1. 上述 MMSE 應改為 least square estimation (LSE), 並沒有 mean (prior pdf) 的觀念。

2. ML estimator depends on pdf (稱為 likelihood function). 但 parameter 本身 no prior information.

3. Depending on the likelihood function, ML 需要 optimize 的 function 不同。For Gaussian likelihood function, 即為 least square function.   For other likelihood function,  就會有不同的 minimal function.

 

另一個例子 (Sine Wave Estimation)

Sine wave estimation 有兩種 flavors: (i) Parameter estimation (假設 amplitude, frequency, phase are fixed); (ii) Waveform estimation (amplitude, frequency, phase can be time varying).

For (ii) waveform estimation, 很顯然只適合 MMSE estimator.  For (i) parameter estimation, ML and MMSE 皆可用。

Sine Wave Parameter Estimation

這個例子在 digial signal processing 常常用到。Estimate 一個 sine wave 的頻率和 phase，不過這個 sine wave 被 additive Gaussian noise 所污染。這個問題有很多變形:sine wave 的 IQ imbalance; sine wave carrier synchronization; sine wave 包含了 harmonics; 或是 harmonic level 的 estimation, etc.

基本上有兩種做法：(i) Time domain or (ii) Frequency domain. 

Time Domain Linear Method

Time domain method is mainly based on linear regression estimation.  It shows accurate result at moderate to high SNRs.

Tretter suggests the following method:

以上 eq. (3), estimate wo and theta using x_t,  t=0, 1, . ,N-1.  可以証明 ML estimator and MMSE estimator are equivalent.   wo and theta are obtained by the equations in linear regression example.  

There are two problems: (i) eq. (2) is only approximate model and accurate only at moderate to high SNRs.; (ii) (3) needs to do phase unwrapped and it can be error prone.

One variation of the sine wave estimation is coherent sampling.  That is, the sampling frequency is integer multiple of the sine wave.  In this case, the sine wave frequency is known, the only parameter to be estimated is the phase.  Therefore, the phase is inversely proportional to the N^2.  

The following is the math of problem statement, CRLB, and Tretter’s model

 

 

 

The following is the Tretter’s model.

 

 

 

Tretter’s algorithm involves unwrap(arctan) operation and achieve CRLB at a high SNR.  It can be used for IQ imbalance estimation withour prior knowledge of the frequency.  (IQ is arctan function).

In summary, Tretter’s algorithm

Pros: (i) No FFT, lower complexity;  (ii) Sequence estimation, low latency

Cons: (ii) phase unwrap add additional complexity and may fail at low SNR

 

Kay improves the algorithm by using the differenced phase of two adjacent samples. 

It is shown that Tretter and Kay’s algorithms are equivalent via Itoh’s phase unwrapping algorithm.

 

 

Time Domain Non-Linear Method

Zero crossing method, not surprising since arctan is nonlinear function.  We can avoid the arctan by using zero crossing

Can it generalize to non-sine waveform?  square wave or pulse like waveform?  It may be useful for asynchronous eye diagram analysis???

Frequency Domain

D.C. Rife

 

Waveform example:

MMSE is used due to close resamble the original waveform!

It is NOT required to be a fixed sine wave.  Some parameter could be time varying.

Example is like FM demodulation.  How about zero crossing MMSE for FM demodulation?

 

A sin(wt + theta)

For MMSE estimator

For ML estimator, estimate A, w, or theta

Kalman Filter : MMSE Estimator

常常說 Kalman filter 是 optimal estimate (in linear dynamic model), 如何解釋? 或者更基本，Kalman filter 的理論基礎為何? Kalman filter 和一般常用的 LPF (Butterworth, Cheby filter, etc.) 有什麼本質上的差異?

Kalman filter 和常用 LPF 差異

一般信號處理常用的 filter (FIR, IIR, Butterworth, Cheby) based on Laplace or Fourier tranform , (i) 都是 linear time invariant (LTI) filter;  (ii) filter 的特性 (e.g.頻寛)和輸入信號基本上 independent.

要計算輸出信號，我們通常先 derive (or define) filter 的 transfer function and impulse response, 再用 frequency domain (s/z-domain) multiplication 或 time domain convolution 來獲得。

Kalman filter, on the other hand, 是由 stochastic signal (random sequence 隨機信號) estimation 出發 (i) 輸入信號和相關雜訊的 statistics 決定了 Kalman filter 的特性;  (ii) Kalman filter 基本上是 time varying filter.  (在很多情況下 Kalman filter 會收歛到 steady state filter, 近似 time invariant filter).  以上兩點和常用 LPF 完全不同。

Kalman filter 和常用 LPF 可以類比嗎?

前文 ( Kalman filter for Engineer) 提到 Kalman filter 可以視為一個 feedback controller with time varying Kalman gain (Kg).  因此可以視為一個 time varying IIR filter, 頻寛會隨時間改變。

在 steady state (如果 steady state 存在)時，Kg 基本上近似常數，趨近一個 time invariant filter).  不過 Kg 由輸入信號和相關雜訊決定，和一般 LPF 仍有不同。

Kalman filter 和 Wiener Filter 的差異

可以 reference Binghamton 的 slides.

Wiener filter assumes WSS (wide sense stationary) signal + noise

Kalman filter assumes Dynamic model + observation noise

Wiener filter 是 time invariant filter, 但 filter 特性由輸入信號和雜訊決定，這點和 Kalman filter 相似。

If we let consider only after much time has elapsed: given IIR Wiener case and steady state Kalman filter, then Kalman = Wiener.

Step1: Random Variable MMSE Estimator

先從最簡單的 one variable 開始，X 的 MMSE (Minimum Mean-Square Error) estimator 是 find x* to minimize E(|X-x*|^2).  很容易可証明: x* = E(X),  MSE=E(|X-E(X)|^2)=σ^2.        Any other estimate E(X)+a 的 MSE=σ^2 + a^2 > σ^2.

當有 two variables X, Y, the MMSE (Minimum Mean-Squared Error) estimator:
Let y be observed.   Fnd x* an estimate of the unobserved x to minimize  E(|X-x*|^2).  Hint: E(|X-x*|^2) = Ex Ey (|X-x*|^2|y)

很明顯如果 X, Y are independent or uncorreleated. y is observed or not 對 X estimation 並沒有任何幫助:

x* = E(X) and MSE = σ_x^2.

如果 X and Y are correlated, y is observed 對 X estimation 多少有幫助。可以証明 x* = E(X|Y=y) 是 MMSE estimator. 並且 MSE error 會比沒有 y is observed 小。直覺上 make sense, 多了 y 的資訊對減小 X esimation error 多少有幫助。

注意 x* = E(X|Y=y) 並沒有假設任何 X and Y 的統計特性 (如 Gaussian)!  這是 Baysian estimator 的基礎。

X, Y Gaussian

再來假設 X, Y Gaussian 來進一步說明。

[x, y] ~ N( [ E(X), E(Y)], [σ_x^2, σ_xy; σ_xy, σ_y^2])

x* = E(X|Y=y) = E(X) + σ_xy (y-E(y)) / σ_y^2

MSE = σ_x^2 – |σ_xy|^2/σ_y^2  < σ_x^2

(i) x* = E(X|Y=y) = a + b y  是一個 linear function of y

(ii) MSE is reduced with additional observable y

(iii) If X, Y are uncorrelated, σ_xy = 0.  x* = E(X); MSE = σ_x^2 as expected

X, Y Non-Gaussian

在 X and Y 非 Gaussian 情況下要得到 x* = E(X|Y=y) 就比較困難。上述 Gaussian 結論需要修正:

(i) x* = E(X|Y=y) = g(y) in general in NOT a linear function of y.  不過可以証明在所有的 linear function (a + b y) 中，之前 Gaussian 所得的公式:

x* = E(X) + σ_xy (y-E(y)) / σ_y^2

MSE 是最小的，稱為 LMSE and LMMSE (Linear Minimum Mean Square Error).   同時可以証明 LMMSE 的確是 minimum (at least local minimum) and unbias.

但是可能存在其他非線性的  x* = g(y) 的 MSE 更小 (global minimum), some example are here.

(ii) and (iii) 仍然是正確。

Step2: MMSE with linear measurement

Consider the following special case (Y is linear transformation x plus another Gaussian noise)

y = H x + v,  x~N(E(X), σ_x^2),  v~N(E(V), σ_v^2)

x* = EX+ K (y-EY)   where EY = H EX + EV, 

K = σ_x^2 H’ (H σ_x^2 H’ + σ_v^2) ^-1  from equation of last section.  K 事實上就是 Kalman gain.

MSE = σ_x^2 – K * H σ_x^2 = σ_x^2 (1 – K H)   Therefore, MSE is reduced due to y is observed.

The above equation forms the base of Kalman filter. 

 

Step3: MMSE with sequences of linear measurement

y = Hx + v    given y1, y2, y3, … yk

x*(k) = E(X|y1, y2, y3, ..yk) =  E(X|y1|y2| .. |yk)

-> x*(1) = E(X|y1) = EX + K(1) (y1 – EY)

EY=H EX+EV  MSE(1)=P(1)=sig_x^2 H’(H sig_x^2 H’+ sig_v^2)^-1

K(1) = sig_x^2 H’ (H sig_x^2 H’ + sig_v^2) ^-1

-> x*(2) = E(X|y1|y2) = E(X|y1) + K(2)(y2 – EY(2))

                = x*(1) + K(2)(y2 – EY(2))

E(X|y1)=x*(1)     EY(2)=Hx*(1)+EV   MSE(2)=P(2)=P(1)^2 H’(HP(1)^2H’+sig_v^2)^-1  K(2) = sig_x^2 H’ (H sig_x^2 H’ + sig_v^2) ^-1

-> x*(3) = x*(2) + K(2)(y2 – EY(2))

其實上述公式就是 Kalman filter 的雛型。只要再加上 state  dynamic model  x(k+1) = A x(k) + W  就變成 Kalman filter, 這部份先帶過。  反之 Kalman filter 的 state equation 簡化為 x(k+1) = x(k), 就變回 sequencial MMSE estimator. 

Sequencial MMSE estimation 和 Kalman filter 一個很大的不同是沒有 W (state process noise).  因此愈多的 measurements,  MSE error 會一直減小，最後會接近 0, 也不會有 steady state. 

在 Kalman filter 有 W (state process noise).  在 prediction stage 時，MSE noise 會增加，update stage 時，MSE noise 才會減少。因此 Kalman filter MSE 最後會收歛到一個非零的 steady state. 圖示如下。

 

MMSE with nonlinear measurement

y = f(x) + v

直覺的 approach, 就是用 Taylor expansion 來 approximate f(x).   However, there are two problems:

(i) In general, y is no longer Gaussian 上述的公式不適合。

(ii) For some common nonlinear function f(x) may introduce bias to disturb the estimation.

除了 Taylor expansion 外，另一方法是用 Gaussian 近似 f(x) 的 pdf.  一個方法是 unscented transform (UT).  之後再介紹。

Kalman Filter Example: Water Tank

以下的 water tank example is from Cornell

Case 1: Constant water level with some measurement noise.

If we use the following Kalman filter with q=0.0001 and r=0.1, we get the following result that match with the true water level.

Case 2: 如果實際的 water level 是隨時間增加，但仍然用以上的 Kalman filter。意即 model error, 會發生什麼狀況?

可以看到 Kalman filter output 的誤差愈來愈大而發散，原因是 Kalman filter model 和實際有落差。如何解決這個問題: (1) 放大 q, state prediction 有更多的自由度; (2) 用 2D model 包含速度。

先看 (1), q=0.01 and r=0.1

有改善但仍然不夠。Try q=r=0.1, 可以追上但 (a) noise 增加; (b) q 需要 try and error, 不實用。

q=1 and r=0.1, 基本上 follow measurement (Kg~1).

正解是用 2D model 包含速度:

qf = 0.00001 and r=0.1 注意 Kalman filter 會自動收斂到對的速度, regardless initial condition.

qf = 0.00001 and r=0.1 for constant tank level. (velocity=0)

Case 3: 如果 water level 是 oscillation, 2D with position and velocity 會如何？

1. Estimated state is smoother but lags behind true value. 這和前面用 static model constant velocity 產生的 lag 類似。代表 model 有問題。

2. The amplitude of the estimated state is getting smaller and smaller.

This is because the model is slowly converging to what it thinks is the truth... a constant level, which is accurate over time.

需要 extended Kalman filter 才能解以上的問題。

Kalman Filter for Engineer

Kalman filter 最初的發展主要在導航應用，之後在很多領域都有新的應用。本文想從不同的角度以及例子了解 Kalman filter,  更直覺 catch Kalman filter physical insight.  在日益重要的 real time 數位信號處理，可以找到廣泛的應用。

一般對於所測量的 data (包含量測的誤差或雜訊)，傳統上就是取 mean (1st order) and variance (2nd order statistics) 如下圖。可以進一步做六個標準差或信心區間的 分析。

更進一步分析就是用 LMS (least mean square) curve fit 如下圖。一般用直線 (linear) 或拋物線 (quadratic) 做 least mean square error curve fit.  這是 Gauss 在研究行星軌道所發展的方法。為什麼用 linear or quadratic curve:  一是簡單 (視覺和 math)，二是在局部近似的效果很好 (like Taylor expansion).

 

 

Kalman (或是某位俄國數學家)不凡之處就是在於更進一步 (i) 探索 data 的內在結構 (state equation and noise covariance matrix), 以及 (ii) 發展出 recursive equation 來 estimate 內在的結構 (state or sometimes noise).

其實 (i) and (ii) 是一體兩面:  Kalman 假設 Markovian state model, 意即目前的 state 只 depends on 前一刻 state, 和更早的 state uncorrelated.  因此在解這個包含所有 time step 的 LMS matrix equation，這個 matrix 是 blocked tri-diagonal.  使用Ricatti recursion 求解時就會是 recursive form.  以上的說法非常 hand waving, 可以 refer 到 MIT 的 video.

(ii) Recursive equation 的好處顯而易見: estimation 不需要儲存所有的 data。事實上，只需要前一刻的必要 data, 加上目前的 measurement 就足夠。這非常適合 real time 的應用 (如導航)。

至於 (i) 探索 data 的內在結構，理想上可以更精確而且有效率的 estimate or even predict data.  實際上遇到的困難是 (a) 這樣的 state equation 是否存在? (b) 如果存在如何找到正確的 model and 參數? and (c) 如果 model or 參數不準確，Kalman filter estimate 的結果是否 robust?

Kalman 考慮的 State equation:

1. X(k)=F X(k-1)+G U(k)+W(k)   (covariance Q)

Measurement:

2. Z(k)=H X(k)+V(k)   (covariance R)

基本上是 linear, Markovian, H may be partially observable --> Hidden Markovian.  圖示如 Fig. 1.

 

 

Kalman Filter Implementation:

X(k|k-1)=F X(k-1|k-1)+G U(k) ……….. (1)
P(k|k-1)=F P(k-1|k-1) F'+Q ……… (2) 
X(k|k)= X(k|k-1)+Kg(k) (Z(k)-H X(k|k-1)) ……… (3) 
Kg(k)= P(k|k-1) H' / (H P(k|k-1) H' + R) ……… (4) 
P(k|k)=（I-Kg(k) H）P(k|k-1) ……… (5)

上述公式看來比起一般 LPF 複雜很多。其實重點是 Kalman filter 是一個 Gaussian pdf estimator.  我們要 estimate 是 Prob(Xk | Y1, Y2,  … Yk) 的 mean and variance.  mean 就是 X(k|k),  variance 就是 P(k|k).  分為兩步: 第一步是 estimate Prob(Xk | Y1, Y2, .. Yk-1), mean and variance are (1) and (2).  再來計算 Kalman gain (4), 以及 Prob(Xk | Y1, Y2, … Yk), mean and variance are (3) and (5).

Note that Kg is NOT a constant, changing based on equation (4).  F, H, and G are time invariant matrix most of time.  Usually, Kg and P can be divided into two regions: (i) transient region; and (ii) steady-state region.

 

 

 

 

 

 

 

 

From Wiki

Because the certainty of the measurements is often difficult to measure precisely, it is common to discuss the filter's behavior in terms of gain. The Kalman gain is a function of the relative certainty of the measurements and current state estimate, and can be "tuned" to achieve particular performance. With a high gain, the filter places more weight on the measurements, and thus follows them more closely. With a low gain, the filter follows the model predictions more closely, smoothing out noise but decreasing the responsiveness. At the extremes, a gain of one causes the filter to ignore the state estimate entirely, while a gain of zero causes the measurements to be ignored.

• Clearly from the figure: K=0 ignores measurement z completely.
• How to explain K=1 ignore state prediction from the figure?

Hand waving explanation of Kalman filter:

1. State prediction based on equation.  This prediction contains noise W.

2. Measurement contains noise V.

3. Trust state or observation?  based on the covariance of state prediction; and covariance of observation.   If lots of noise in state prediction, lean more to observation; or vice versa.   The ratio of covariances determines Kalman gain.

Feedback Explanation

Assuming G=0,  F=A, H=C,  Lt = Kg, Kalman filter 可以簡化為一個 feedback loop 如下圖。先考慮最簡單的 case: 

A is a scaler (1x1).  The plant 本身包含一個 integrator with feedback (Kalman low pass filter). The feedback controller Lt (Kg) behaves as a P (Proportion) controller.  
The close loop response = Lt Z^-1/[1-(1-Lt)Z^-1] ~ Lt/(Lt+s).
Lt (Kg) 大則 loop gain and loop bandwidth 大。yt 沒有太多的 filtering, trust more on yt (i.e. measurement).  (Lt=Kg=1 Kalman filter = Z^-1 no filtering at all)
Lt (Kg) 小則 loop gain and loop bandwidth 小。yt 會被 filtering, trust more on state prediction.  (Lt=Kg --> 0 Kalman filter becomes an integrator)

更複雜的 case: A = [1, 1; 1, 0] (as case 3 below), the feedback controller behaves as a PI (Proportion and Integrator) controller.   The detailed analysis 可參考前文: Kalman filter and PLL.

 

 

 

 

 

 

 

 

 

Kalman filter 和一般 feedback controller (P/PI/PID) 一個非常重要的差別是 Kg 是由 Eqs (1)-(5) 所控制，Kg 基本是 time varying and depends on the 內在結構 (state and noise covariance matrix), 後面有一些例子 Kg 是如何 time varying.   傳統的 P/PI/PID controller 基本上是 (i) 固定值或是 (ii) 分兩段 (acquisition and tracking), roughly 對應到 Kalman filter transient and steady-state 兩階段，當然只是簡化版。

Handwaving supplement

 

1.  Given measurement z(0), z(1), z(2), .... z(k-1), estimate x(0), x(1), ..., x(k-1),      x(k) and update after z(k) is received.

2a.  KF is Markovian: X(k) depends on X(k-1), not X(k-2) ... and W(k), V(k).  Therefore X(k|k-1, ..1, 0) = X(k|k-1)
2b.  KF is recursive: statistics of X(k) (covariance) depends on k-1.
2c.   Observable Z = H X usually has less degree of freedom than X; therefore hidden Markov.

3.  KF is linear:  X(k) = A X(k-1) + B U(k) + W(k)

Matrix Explanation

 

http://videolectures.net/mit18085f07_strang_lec16/

Some examples

 

Trivial case 0: (perfect measurement)

X(k)=A X(k-1)+B U(k)+W(k)

Z(k) = X(k)       Covariance matrix R = 0

 

Given Z(0), Z(1), ..., Z(k) without measurement errors

The best estimation of X(k) is simply Z(k)

H Kg(k) = I;      if H=I  then Kg(k) = 1;     P(k|k) = 0

Therefore, Kg=1 means perfect measurement; trust measurements and ignore state prediction completely.

 

 

Case 1: (constant state/no state noise, with measurement noise)

X(k) = X(k-1)

Z(k) = X(k) + V(k)

 

Given Z(0), Z(1), ..., Z(k) and V(k) is i.i.d. with covariance matrix R

The best estimation is:  Z(0) + ... + Z(k) / (k+1)

More average reduces the variance of the estimation

P(0|0) = R

P(k|k-1) = P(k-1|k-1) = R / k

P(k|k) = P(k|k-1) * k / (k+1)

Kg(k) = P(k|k-1) / (P(k|k-1) + R) = 1 / (k+1)

P(k|k) = P(0|0) / (k+1) = R / (k+1)

 

When k becomes large, the state prediction becomes more trustable and the new measurement contains very little new information, Kg --> 0 and P --> 0.  從 feedback loop bandwidth 角度來看，就是愈來愈小。

Let's estimate a constant signal of 1. The noise is normal distributed with a variance of 25. Therefore the SNR is only 1/25=0.04 (-28dB), which is very hard to give a good estimate of the value. The Kalman filter doesn't seem to care much, after a few ten's samples he gives a very good estimate (Fig.1).  The Kalman gain, Kg(k) (Fig.2) and estimation error, P(k|k) (Fig.3) are inversely proportional to k as expected. Fig.4 shows the measurement error ,Z(k)-X(k), in green and estimation error, X^(k)-X(k), in red after Kalman filter.  Again, Kalman filter is very effective in this case.

 

Fig. 1

 

 

 

 

Fig.2

 

 

Fig.3

Fig.4

 

 

 

 

 

 

 

 

 

 

 

Case 2: (State with noise accumulation like random walk or phase noise, with measurement noise)

X(k) = X(k-1) + W(k)

Z(k) = X(k) + V(k)

Given Z(0), Z(1), ..., Z(k) and W(k), V(k) is i.i.d. with covariance matrix Q, R

X(1) = X(0) + W(1)

X(2) = X(0) + W(1) + W(2)

.. X(k) = X(0) + sum(W(k))

==>  When k is large, the state behaves like random walk.   Moreover, there are measurement noise on the top of the random walk.

The best estimation is: ?

P(k|k-1)=P(k-1|k-1) +Q

H = I

Kg(k)= P(k|k-1) / (P(k|k-1) + R) = (P(k-1|k-1) + Q ) / (P(k-1|k-1) + R + Q)

P(k|k) = (I-Kg(k)) P(k|k-1) = RP(k|k-1)/(P(k|k-1)+R)

=> 1/P(k|k) = 1/P(k|k-1) + 1/R  (new sample is useful to reduce variance!)

      1/P(k|k) = 1 / (P(k-1|k-1) + Q) + 1/R

     P(k|k) = (P(k-1|k-1)+Q) R / ( P(k-1|k-1) + Q + R)

P(0|0)

P(1|0) = P(0|0) + Q

P(1|1) =  (P(0|0)+Q)R / (P(0|0) + Q + R)

In the steady state, P(k|k) --> P = 1/2(sqrt(Q^2+4QR)-Q) ;  Kg converges to (P+Q)/(P+R+Q)

If R is small,  P --> R

Kg and P have no close forms.  Kg and P are monotonic descending and converging to a constant (steady-state).  從 loop bandwidth 角度來看，就是變小最後收歛到一個固定值。這個固定值 depends on Q and R. 

Fig. 5 shows the random walk of state with variance 0.1 (in blue), with measurement noise of variance 25 as case 1 (in green), and after Kalman filter (in red).   Both Kg and P converge to steady state quickly shown in Fig. 6 and Fig. 7.   Fig. 8 shows the measurement error ,Z(k)-X(k), in green and estimation error, X^(k)-X(k), in red after Kalman filter.  Apparently, the estimate error is no longer zero as in the Case 1 even with a small perturbation state noise of 0.1.  P converges to 1.53; Kg converges to 0.061.

Fig.5

Fig.6

 

Fig.7

Fig.8

Even though there is no close forms of Kg and P (close form only for steady-state), a simple resistor ladder illustration is useful.   Po is the end resistor, serial with Q, and parallel with R, and so on.

If Q = 0; P(k) = R/k // P(0)
If R -->0;  P(k) --> 0















Case 3: (position and velocity estimation, no force with noise)

 

 

 

 

 

 

 

p''(t) = 0  (no acceleration) => [p, p']' = [0 1; 0 0][p, p'] = A [p,p']

Transition matrix = exp(A dt) for discrete time
exp(A dt) = I + A dt + A^2/2! + A^3/3! + ...    A^2=A^3=..=0
exp(A dt) = [1 dt; 0 1]

To verify the equation, assuming no noise:

p(k+1) = p(k) + p'(k) dt

p'(k+1) = p'(k)    ---> constant velocity

H = [1 0] only observe position

z(k) = p(k)  --> observable

這一點非常重要!  不需要事先知道 velocity (no prior knowledge of velocity), 就可以從 position measurement estimate 出來!

 

Case 4: (phase and frequency estimation, no cycle slip)

Case 4a: (phase and frequency estimation, with cycle slip, mod)

Case 5: Oscillation etimation
Case 6: Nonlinear estimation (sign estimation)

 

2.  Kalman filter is a special case of Baysian recursive filter! 

 

Some applications

DCOC (dc offset cancellation)

Baseline wander

sine wave estimation (extended Kalman filter)
multi-bit SAR (successive approximation ADC)