> For the complete documentation index, see [llms.txt](https://aparande.gitbook.io/berkeley-notes/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://aparande.gitbook.io/berkeley-notes/eecs225a-0/eecs225a-4.md).

# Filtering

If we think of our signal as a discrete time random process, then like a normal deterministic signal, we can try filtering our random process.

![Figure 1: Filtering a Disrete Time Random Process with an LTI system with transfer function $H(z)$](/files/g2NEXA6QPAYNFO6ETPmD)

Filtering can either be accomplished with an LTI system or some other non-linear/non-time-invariant system just like with deterministic signals.

## LTI Filtering on WSS Processes

If we use an LTI filter on a WSS process, then we can easily compute how the filter impacts the spectrum of the signal.

{% hint style="info" %}

#### Theorem 13

When $$Y(n)$$ is formed by passing a WSS process $$X\_n$$ through a stable LTI system with impulse response $$h\[n]$$ and transfer function $$H(z)$$, then $$S\_Y(z) = H(z)S\_X(z)H^*(z^{-*})$$ and $$S\_{YX}(z) = H(z)S\_X(z)$$. If we have a third process $$Z\_n$$ that is jointly WSS with $$(Y\_n, X\_n)$$, then $$S\_{ZY}(z) = S\_{ZX}(z)H^*(z^{-*})$$.
{% endhint %}

This gives us an interesting interpretation of the spectral factorization (Definition 14) since it essentially passing a WSS process with auto-correlation $$R\_W(k) = r\_e\delta\[n]$$ through a minimum-phase filter with transfer function $$L(z)$$.

## Wiener Filter

Suppose we have a stochastic WSS process $$Y\_n$$ that is jointly WSS with $$X\_n$$ and that we want to find the best linear estimator of $$X\_n$$ using $$Y\_n$$. The best linear estimator of $$X\_n$$ given the observations $$Y\_n$$ can be written as

$$\hat{X}*n = \sum*{m\in\mathbb{Z}}h(m)Y\_{n-m} = h\[n] \* Y\_n.$$

This is identical to passing $$Y\_n$$ through an LTI filter. If we restrict ourselves to using $${Y\_i}\_{i=-\infty}^{n}$$ to estimate $$X\_n$$, then the best linear estimator can be written as

$$\hat{X}*n = \sum*{m=0}^\infty h(m)Y\_{n-m} = h\[n] \* Y\_n.$$

It is identical to passing $$Y\_n$$ through a causal LTI filter. Since we are trying to find a best linear estimator, it would be nice if each of the random variables we are using for estimating were uncorrelated with each other. In other words, instead of using $$Y$$ directly, we want to transform $$Y$$ into a new process $$W$$ where $$R\_W(k) = \delta\[k]$$. This transformation is known as whitening. From the spectral factorization of $$Y$$, we know if we use the filter $$G(z) =\frac{1}{S\_Y^+(z)}$$ then

$$S\_W(z) = \frac{S\_Y(z)}{S\_Y^+(z)S\_Y^{+*}(z^{-*})} = \frac{S\_Y(z)}{S\_Y^+(z)S\_Y^-(z)} = 1.$$

Now we want to find the best linear estimator of $$X$$ using our new process $$W$$ by designing an LTI filter $$Q(z)$$.

![Figure 2: Finding the best linear estimator of $X$ using $W$ with a two-stage filter that first whitens the input.](/files/MnLA2VskhNSzch1ioBq8)

### Non-Causal Case

Starting with noncausal case, we can apply the orthogonality principle,

$$\begin{aligned} \mathbb{E}\left\[(X\_n-\hat{X}*n)W*{n-k}^*\right] = 0 &\implies \mathbb{E}\left\[X\_nW^**{n-k}\right] = \sum*{m\in\mathbb{Z}}q(m)\mathbb{E}\left\[W\_{n-m}W^**{n-k}\right] \ \therefore R*{XW}(k) = \sum\_{m\in\mathbb{Z}}q(m)R\_W(k-m) &\implies S\_{XW}(z) = Q(z)S\_W(z)\ \therefore Q(z) = \frac{S\_{XW}(z)}{S\_W(z)} = S\_{XW}(z) &= S\_{XY}(z)(S\_Y^+(z^{-*}))^{-\*} = \frac{S\_{XY}(z)}{S\_Y^-(z)}\end{aligned}$$

When we cascade these filters,

$$H(z) = Q(z)G(z) = \frac{S\_{XY}(z)}{S\_Y^-(z)} \frac{1}{S\_Y^+(z)}= \frac{S\_{XY}(z)}{S\_Y(z)}.$$

{% hint style="info" %}

#### Definition 21

The best linear estimator of $$X\_n$$ using $$Y\_n$$ where $$(X\_n, Y\_n)$$ is jointly WSS is given by the non-causal Wiener filter.

$$H(z) = \frac{S\_{XY}(z)}{S\_Y(z)}$$
{% endhint %}

If we interpret Definition 21 in the frequency domain, for a specific $$\omega$$, we can understand $$H(e^{j\omega})$$ as an optimal linear estimator for $$F\_X(\omega)$$ where $$F\_X(\omega)$$ is the the stochastic process given by the Cramer-Khinchin decomposition (Theorem 7). More specifically, we can use the Cramer-Khinchin decomposition of $$Y\_n$$.

$$\begin{aligned} \hat{X}*n &= \sum*{i\in\mathbb{Z}}h\[i]\int\_{-\pi}^\pi e^{j\omega(n-i)}dF\_Y(\omega)\ &= \int\_{-\pi}^{\pi}\left(\sum\_{i\in\mathbb{Z}}h\[i]e^{-j\omega i}\right)e^{j\omega n}dF\_Y(\omega) \ &= \int\_{-\pi}^\pi H(e^{j\omega})e^{j\omega n}dF\_Y(\omega)\\\end{aligned}$$

Since $$F\_X$$ and $$F\_Y$$ have jointly orthogonal increments, this tells us that $$H(e^{j\omega})$$ is just the optimal linear estimator of $$dF\_X(\omega)$$ using $$dF\_Y(\omega)$$. $$dF\_X(\omega)$$ and $$dF\_Y(\omega)$$ exist on a Hilbert space, meaning we are essentially projecting each frequency component of $$X\_n$$ onto the corresponding frequency component of $$Y\_n$$.

### Causal Case

First, note that in the causal case, whitening doesn’t break causality because $$\frac{1}{S\_Y^+(z)}$$ is causal. When we apply the orthogonality principle,

$$\begin{aligned} \mathbb{E}\left\[(X\_n-\hat{X}*n)W*{n-k}^*\right] = 0 &\implies \mathbb{E}\left\[X\_nW^**{n-k}\right] = \sum*{m=0}^\infty q(m)\mathbb{E}\left\[W\_{n-m}W^\**{n-k}\right] \ \therefore R*{XW}(k) &= \sum\_{m = 0}^\infty q\[m]R\_W(k-m) \qquad k \geq 0\end{aligned}$$

We can’t take the Z-transform of both sides because the equation is not necessarily true for $$k < 0$$. Instead, we can look at the function

$$f(k) = R\_{XW}(k) - \sum\_{m=0}^\infty R\_W(k-m)q\[m] = \begin{cases} 0 & k\geq 0,\ ? & \text{ else.}\end{cases}$$

Taking the unilateral Z-transform of both sides,

$$\begin{aligned} \left\[F(z)\right]*+ &= \left\[S*{XW}(z) - S\_W(z)Q(z)\right]*+ = \left\[S*{XW}(z)\right]*+ - Q(z) = 0\ Q(z) &= \left\[S*{XW}(z)\right]*+ = \left\[\frac{S*{XY}(z)}{S\_Y^-(z)}\right]\_+ \end{aligned}$$

Thus the filter $$H$$ which gives the causal best linear estimator of $$X$$ using $$Y$$ is

$$H(z) = Q(z) G(z)= \left\[\frac{S\_{XY}(z)}{S\_Y^-(z)}\right]\_+ \frac{1}{S\_Y^+(z)}.$$

{% hint style="info" %}

#### Definition 22

The best linear estimator of $$X\_n$$ using $${Y\_i}\_{i=-\infty}^{n}$$ is given by the causal Wiener filter.

$$H(z) = Q(z)G(z) = \left\[\frac{S\_{XY}(z)}{S\_Y^-(z)}\right]\_+ \frac{1}{S\_Y^+(z)}.$$
{% endhint %}

Intuitively, this should make sense because we are using the same $$W$$ process as in the non-causal case, but only the ones which we are allowed to use, hence use the unilateral Z-transform of the non-causal Wiener filter, which amounts to truncated the noncausal filter to make it causal.

{% hint style="info" %}

#### Theorem 14

If $$\hat{X}*{NC}(n)$$ is the non-causal Wiener filter of $$X$$, then the causal wiener filter of $$X$$ given $$Y$$ is the same as the causal wiener filter of $$\hat{X}*{NC}$$ given $$Y$$, and if $$Y$$ is white noise, then

$$\hat{X}*C(n) = \sum*{i=0}^{\infty}h\[i]Y\_{n-i}$$
{% endhint %}

### Vector Case

Suppose that instead of a Wide-Sense Stationary process, we an $$N$$ length signal $$\boldsymbol{X}$$ which we want to estimate with another $$N$$ length signal $$\boldsymbol{Y}$$. We can represent both $$\boldsymbol{X}$$ and $$\boldsymbol{Y}$$ as vectors in $$\mathbb{C}^N$$. If we are allowed to use all entries of $$\boldsymbol{Y}$$ to estimate $$\boldsymbol{X}$$, this is identical to linear estimation.

{% hint style="info" %}

#### Definition 23

The non-causal Wiener filter of a finite length $$N$$ signal $$\boldsymbol{Y}$$ is given by

$$K\_s = R\_{\boldsymbol{X}\boldsymbol{Y}}R\_{\boldsymbol{Y}}^{-1}.$$
{% endhint %}

Note that this requires $$R\_{\boldsymbol{Y}} \succ 0$$. Suppose that we wanted to design a causal filter for the vector case, so $$\hat{X}*i$$ only depends on $${Y\_j}*{j=1}^i$$. By the orthogonality principle,

$$\forall 1 \leq l \leq i,\ \mathbb{E}\left\[X\_i - \sum\_{j=1}^iK\_{f, ij}Y\_jY\_l^\*\right] = 0 \implies R\_{\boldsymbol{XY}}(i, l) = \sum\_{j=1}^i K\_{f, ij}R\_{\boldsymbol{Y}}(j, l)$$

In matrix form, this means

$$R\_{\boldsymbol{XY}} - K\_fR\_{\boldsymbol{Y}} = U^+$$

where $$U^+$$ is strictly upper triangular.

{% hint style="info" %}

#### Theorem 15

If matrix $$H \succ 0$$, then there exists a unique lower-diagonal upper triangular factorization of $$H=LDL^\*$$ where $$L$$ is lower diagonal and invertible with unit diagonal entries and $$D$$is diagonal with positive entries.
{% endhint %}

Applying the LDL decomposition, we see that

$$\begin{aligned} R\_{\boldsymbol{XY}} - K\_fLDL^\* = U^+ \implies R\_{\boldsymbol{XY}}L^{-*}D^{-1} -K\_f L = U^+L^{-*}D^{-1}\ \therefore \[R\_{\boldsymbol{XY}}L^{-\*}D^{-1}]\_L -K\_f L = 0\end{aligned}$$

where $$\[\cdot]\_L$$ represent the lower triangular part of a matrix.

{% hint style="info" %}

#### Definition 24

The causal Wiener filter of a finite length $$N$$ signal $$\boldsymbol{Y}$$ is given by

$$K\_f = \[R\_{\boldsymbol{XY}}L^{-\*}D^{-1}]\_LL^{-1}$$
{% endhint %}

## Hidden Markov Model State Estimation

Suppose we have a Hidden Markov Process $${Y\_n}*{n\geq1}$$. We can think of determining the state $${X\_n}*{n\geq1}$$ as filtering $${Y\_n}\_{n\geq1}$$.

### Causal Distribution Estimation

Suppose we want to know the distribution of $$X\_t$$ after we have observered $$Y^t$$.

$$\begin{aligned} p(x\_t|y^t) &= \frac{p(x\_t,y^t)}{p(y^t)} = \frac{p(x\_t)p(y\_t, y^{t-1}|x\_t)}{\sum\_{x}p(y\_t,y^{t-1}|x\_t=x)p(x\_t=x)}\ &= \frac{p(x\_t)p(y\_t|x\_t)p(y^{t-1}|x\_t)}{\sum\_xp(y\_t|x\_t=x)p(y^{t-1}|x\_t=x)p(x\_t=x)} = \frac{p(y\_t|x\_t)p(y^{t-1})p(x\_t|y^{t-1})}{\sum\_{x}p(y\_t|x\_t=x)p(y^{t-1})p(x\_t=x|y^{t-1})}\ &=\frac{p(y\_t|x\_t)p(x\_t|y^{t-1})}{\sum\_{x}p(y\_t|x\_t=x)p(x\_t=x|y^{t-1})}\end{aligned}$$

Now if we know $$p(x\_t|y^{t-1})$$, then we are set.

$$\begin{aligned} p(x\_t|y^{t-1}) &= \sum\_xp(x\_t,x\_{t-1}=x|y^{t-1}) = \sum\_x p(x\_{t-1}=x|y^{t-1})p(x\_t|x\_{t-1}=x,y^{t-1}) \ &= \sum\_x p(x\_{t-1}=x|y^t)p(x\_t|x\_{t-1}=x)\end{aligned}$$

Now we have a recursive algorithm for computing the distribution of $$x\_t$$.

![](/files/UCmfFn9lXH5odR7mDoWx)

### Non-Causal Distribution Estimation

Suppose we are allowed to non-causally filter our signal and we care about the distribution of $$X\_t$$ after we have observed $$Y^n$$. In other words, for $$t \geq n$$, we want to find $$\gamma\_t(x\_t) = p(x\_t|y^n)$$. When $$t=n$$, $$\gamma\_n(x\_n) = \alpha\_n(x\_n)$$. If we continue expanding backwards, then

$$\begin{aligned} p(x\_t|y^n) &= \sum\_x p(x\_t,x\_{t+1}=x|y^n) = \sum\_x p(x\_{t+1}=x|y^n)p(x\_t|x\_{t+1}=x,y^t,y\_{t+1}^n)\ &= \sum\_x p(x\_{t+1}=x|y^n)p(x\_t|x\_{t+1},y^t) = \sum\_x p(x\_{t+1}=x|y^n)\frac{p(x\_t|y^t)p(x\_{t+1}=x|x\_t,y^t)}{p(x\_{t+1}=x|y^t)}\ &= \sum\_x \gamma\_{t+1}(x)\frac{\alpha\_t(x\_t)p(x\_{t+1}=x|x\_t)}{\beta\_{t+1}(x)}\end{aligned}$$

This gives us a clear algorithm for non-causally computing the distribution of $$x\_t$$.

![](/files/v6mKXNNvLHbbh2wnw8gO)

### State Sequence Estimation

Suppose we want to find the most likely sequence of states given our observations. This means we should compute

$$\hat{X}^n = \text{argmax}\_{X^n}p(x^n|y^n)$$

$$\begin{aligned} p(x^t, y^t) &= p(x^{t-1}, y^{t-1})p(x\_t, y\_t|x^{t-1},y^{t-1})\ &= p(x^{t-1}, y^{t-1})p(x\_t|x^{t-1},y^{t-1})p(y\_t|x\_t,x^{t-1},y^{t-1}) \ &= p(x^{t-1},y^{t-1})p(x\_t|x\_{t-1})p(y\_t|x\_t)\end{aligned}$$

We see that there is a recursion in the joint distribution, so if we let $$V\_t(x\_t) = \max\_{x^{t-1}}p(x^t,y^t)$$, then

$$\begin{aligned} V\_t(x\_t) &= \max\_{x^{t-1}} p(x^t, y^t) = p(y\_t|x\_t)\max\_{x^{t-1}}p(x^{t-1},y^{t-1})p(x\_t|x\_{t-1})\ &= p(y\_t|x\_t)\max\_{x^{t-1}}\left\[p(x\_t|x\_{t-1}) \max\_{x^{t-2}} p(x^{t-1},y^{t-1})\right]\ &= p(y\_t|x\_t)\max\_{x^{t-1}}p(x\_t|x\_{t-1}) V\_{t-1}(x\_{t-1})\end{aligned}$$

The base case is that $$V\_1(x\_1) = p(x\_1)p(y\_1|x\_1)$$. $$V\_t$$ is useful because $$\hat{x}*n = \text{argmax}*{x\_n}V\_n(x\_n)$$. This is because we can first maximize over $$\hat{X}^{n-1}$$ and $$Y^n$$, so the only thing left to maximize is $$\hat{x}\_n$$. Once we have $$\hat{x}*t$$, then we can comptue $$\hat{x}*{t-1}$$ by

$$\hat{x}*{t-1} = \text{argmax}*{x\_{t-1}}p(\hat{x}*t|x*{t-1})V\_{t-1}(x\_{t-1}).$$

Putting these equations gives us the Viterbi algorithm.

![](/files/5D603pHECL6gssn8ZA11)

## Kalman Filtering

In the Kalman Filter setup, we assume that the signal we would like to filter can be represented by a state-space model. We want to predict the state vectors $$\boldsymbol{\hat{X}}\_i$$ using some linear combination of the observations $$\boldsymbol{Y}\_i$$.

### Kalman Prediction Filter

Suppose that we want to compute the one-step prediction. In other words, given $$\boldsymbol{Y}^i$$, we want to predict $$\boldsymbol{\hat{X}}\_{i+1}$$. Our observations $$\boldsymbol{Y}$$ are the only thing which give us information about the state, so it would be nice if we could de-correlate all of the $$\boldsymbol{Y}$$. To do this, we can define the innovation process

$$\boldsymbol{e}\_i = \boldsymbol{Y}*i - \boldsymbol{\hat{Y}*{i|i-1}} = \boldsymbol{Y}*i - H\_i\boldsymbol{\hat{X}}*{i|i-1}$$

The last equality follows from the state-space modela and that past observation noises are uncorrelated with the current one. Now, to compute the one-step prediction, we just need to project $$\boldsymbol{\hat{X}}\_i$$ onto the innovations.

$$\begin{aligned} \boldsymbol{\hat{X}}*{i+1|i} &= \sum*{j=0}^i\langle \boldsymbol{X}*{i+1}, \boldsymbol{e}*j \rangle R*{\boldsymbol{e},j}^{-1}\boldsymbol{e}*j \ &= \langle \boldsymbol{X}*{i+1}, \boldsymbol{e}*{i} \rangle R\_{\boldsymbol{e},i}^{-1}\boldsymbol{e}*i + \sum*{j=0}^{i-1}\langle \boldsymbol{X}*{i+1}, \boldsymbol{e}*j \rangle R*{\boldsymbol{e},j}^{-1}\boldsymbol{e}*j\ &= \langle \boldsymbol{X}*{i+1}, \boldsymbol{e}*{i} \rangle R\_{\boldsymbol{e},i}^{-1}\boldsymbol{e}*i + \boldsymbol{\hat{X}}*{i+1|i-1} = \langle \boldsymbol{X}*{i+1}, \boldsymbol{e}*{i} \rangle R\_{\boldsymbol{e},i}^{-1}\boldsymbol{e}*i + \boldsymbol{\hat{X}}*{i+1|i}\ &= \langle \boldsymbol{X}*{i+1}, \boldsymbol{e}*{i} \rangle R\_{\boldsymbol{e},i}^{-1}\boldsymbol{e}*i + F\_i\boldsymbol{\hat{X}}*{i|i-1}\end{aligned}$$

The second to last equality follows from the Wide-Sense Markovity of state space models, and the last equality is due to the state evolution noises being uncorrelated. If we let $$K\_{p,i} = \langle \boldsymbol{X}\_{i+1}, \boldsymbol{e}*i \rangle R*{\boldsymbol{e},i}^{-1}$$ (called the prediction gain), then we have a recursive estimate of the optimal one-step predictor.

$$\boldsymbol{\hat{X}}*{i+1|i} = F\_i\boldsymbol{\hat{X}}*{i|i-1} + K\_{p,i}\boldsymbol{e}\_i.$$

Now, we just need to find a recursive formulation for $$K\_{p,i}$$ and $$R\_{\boldsymbol{e},i}$$. Starting with $$R\_{\boldsymbol{e},i}$$, notice that we can write $$\boldsymbol{e}\_i = \boldsymbol{Y}*i-H\_i\boldsymbol{\hat{X}}*{i|i-1} = H\_i(\boldsymbol{X}*i-\boldsymbol{\hat{X}}*{i|i-1})+\boldsymbol{V}\_i$$.

$$\begin{aligned} R\_{\boldsymbol{e},i} &= \langle H\_i(\boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1})+ \boldsymbol{V\_i}, H\_i(\boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1})+\boldsymbol{V\_i} \rangle \ &= H\_i\langle \boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1}, \boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1} \rangle H\_i^\* + R\_i\end{aligned}$$

To find $$K\_{p,i}$$, we should first find $$\langle \boldsymbol{X}\_{i+1}, \boldsymbol{e}\_i \rangle$$.

$$\begin{aligned} \langle \boldsymbol{X}\_{i+1}, \boldsymbol{e}\_i \rangle &= F\_i\langle \boldsymbol{X}\_i, \boldsymbol{e}\_i \rangle + G\_i\langle \boldsymbol{U}\_i, \boldsymbol{e}\_i \rangle \ &= F\_i\langle \boldsymbol{X}\_i, H\_i(\boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1})+\boldsymbol{V}\_i \rangle + \langle \boldsymbol{U}\_i, H\_i(\boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1})+\boldsymbol{V}\_i \rangle \ &= F\_i\langle \boldsymbol{X}*i, \boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1} \rangle H\_i^\* + G\_iS\_i\ &= F\_i\langle (\boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1}) + \boldsymbol{\hat{X}}*{i|i-1}, \boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1} \rangle H\_i^\* + G\_iS\_i\ &= F\_i\langle \boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1}, \boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1} \rangle H\_i^\* + G\_iS\_i\end{aligned}$$

Notice that the matrix $$P\_i = \langle \boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1}, \boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1} \rangle$$is the auto-correlation of the estimation error, and it shows up in both $$K\_{p,i}$$ and $$R\_{\boldsymbol{e}\_i}$$. It would be useful to have a recursive solution for this matrix as well.

$$\begin{aligned} P\_{i+1} &= \Pi\_{i+1}-\langle \boldsymbol{\hat{X}}*{i+1|i}, \boldsymbol{\hat{X}}*{i+1|i} \rangle \ &= F\_i\Pi\_iF\_i^\* + G\_iQ\_iG\_i^\* - \langle F\_i\boldsymbol{\hat{X}}*{i|i-1}+K*{p,i}\boldsymbol{e}*i, F\_i\boldsymbol{\hat{X}}*{i|i-1}+K\_{p,i}\boldsymbol{e}*i \rangle \ &= F\_i\Pi\_iF\_i^\* + G\_iQ\_iG\_i^\* - F\_i\langle \boldsymbol{\hat{X}}*{i|i-1}, \boldsymbol{\hat{X}}*{i|i-1} \rangle F^\*+K*{p,i}R\_{\boldsymbol{e},i}K\_{p,i}^*\ &= F\_iP\_iF\_i^* + G\_iQ\_iG\_i^\* - K\_{p,i}R\_{\boldsymbol{e},j}K\_{p,i}^\*\end{aligned}$$

Putting this into a concrete algorithm, we get the Kalman Prediction Filter.

### Schmidt’s Modification of the Kalman Filter

The predictive Kalman filter goes directly from $$\boldsymbol{\hat{X}}*{i|i-1}$$ to $$\boldsymbol{\hat{X}}*{i+1|i}$$ without ever determining $$\boldsymbol{\hat{X}}\_{i|i}$$. The Schmidt Modification of the Kalman filter separates the predictive kalman filter into two steps, allowing us to estimate the current state.

1. Measurement Update: Find $$\boldsymbol{\hat{X}}*{i|i}$$ given the latest observation $$\boldsymbol{Y}*{i}$$ and $$\boldsymbol{\hat{X}}\_{i|i-1}$$.
2. State Evolution (Time) Update: Find $$\boldsymbol{\hat{X}}\_{i+1|i}$$ using what we know about the state evolution.

This mimics the approach of the forward algorithm for Hidden Markov Models, which separated updates to the distribution using a time update and a measurement update. Using our innovation process,

$$\begin{aligned} \boldsymbol{\hat{X}}*{i|i} &= \sum*{j=0}^i \langle \boldsymbol{X}*i, \boldsymbol{e}*j \rangle R*{\boldsymbol{e},j}^{-1}\boldsymbol{e}*j \ &= \boldsymbol{\hat{X}}*{i|i-1} + \langle \boldsymbol{X}*i, \boldsymbol{e}*j \rangle R*{\boldsymbol{e},i}^{-1}\boldsymbol{e}^j\ &= \boldsymbol{\hat{X}}*{i|i-1} + \langle (\boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1}) + \boldsymbol{\hat{X}}*{i|i-1}, H\_i(\boldsymbol{X}*i - \boldsymbol{\hat{X}}*{i|i-1}) + \boldsymbol{V}*i \rangle R*{\boldsymbol{e},i}^{-1}\boldsymbol{e}^j\ &= \boldsymbol{\hat{X}}*{i|i-1} + P\_iH\_i^\*R*{\boldsymbol{e},i}^{-1}\boldsymbol{e}\_i\end{aligned}$$

The gain on the coefficient of the innovation $$K\_{f,i}=P\_iH\_i^\*R\_{\boldsymbol{e},i}$$ is called the Kalman Filter Gain. The error of our estimator $$P\_{i|i} = \langle \boldsymbol{\hat{X}}*{i|i}, \boldsymbol{\hat{X}}*{i|i} \rangle$$is given by

$$P\_{i|i} = P\_i - P\_iH\_i^\*R\_{\boldsymbol{e},i}^{-1}H\_iP\_i.$$

For the time update,

$$\begin{aligned} \boldsymbol{\hat{X}}*{i+1|i} &= F\_i\boldsymbol{\hat{X}}*{i|i} + G\_i\boldsymbol{\hat{U}}*{i|i} \ &= F\_i\boldsymbol{\hat{X}}*{i|i} + G\_i\langle \boldsymbol{U}\_i, \boldsymbol{e}*i \rangle R*{\boldsymbol{e},i}^{-1}\boldsymbol{e}*i = F\_i\boldsymbol{\hat{X}}*{i|i} + G\_i\langle \boldsymbol{U}\_i, \boldsymbol{e}*i \rangle R*{\boldsymbol{e},i}^{-1}\boldsymbol{e}*i\ &= F\_i\boldsymbol{\hat{X}}*{i|i} + G\_i\langle \boldsymbol{U}*i, H\boldsymbol{X}*i+\boldsymbol{V}*i - H\_i\boldsymbol{\hat{X}}*{i|i-1} \rangle R*{\boldsymbol{e},i}^{-1}\boldsymbol{e}*i\ &= F\_i\boldsymbol{\hat{X}}*{i|i} + G\_iS\_iR*{\boldsymbol{e},i}^{-1}\boldsymbol{e}\_i\end{aligned}$$

We can re-write the error if this estimator $$P\_{i+1}$$ as

$$P\_{i+1} = F\_iP\_{i|i}F\_i^\* + G\_i(Q\_i - S\_iR\_{\boldsymbol{e},i}^{-1}S\_i^*)G\_i^* - F\_iK\_{f,i}S\_i^*G\_i^* - G\_iS\_iK\_{f,i}^*F\_i^*$$

Writing this as an algorithm,

### Kalman Smoother

The Kalman Prediction Filter and Schmidt’s modification of the Kalman filter are both causal filters. The Kalman Smoother provides a non-causal estimate in the same way that the Backward Recursion algorithm does for Hidden Markov Processes. In other words, the Kalman Smoother predicts $$\hat{\boldsymbol{X}}\_{i|n}$$, the best linear estimator of $$\hat{\boldsymbol{X}}\_i$$ from $${\boldsymbol{Y}\_0, \cdots, \boldsymbol{Y}\_N}$$. As before, we can start with the innovation process $$\boldsymbol{e}\_j = \boldsymbol{Y}*j - H\_j\hat{\boldsymbol{X}}*{j|j-1}$$.

$$\hat{\boldsymbol{X}}*{i|N} = \sum*{j=0}^N \langle \boldsymbol{X}*i, \boldsymbol{e}*j \rangle R*{e,j}^{-1}\boldsymbol{e}*j = \hat{\boldsymbol{X}}*{i|i-1} + \sum*{j=i}^N \langle X\_i, \boldsymbol{e}*j \rangle R*{e,j}^{-1}\boldsymbol{e}\_j.$$

We can compute $$\hat{\boldsymbol{X}}\_{i|i-1}$$ using the Predictive Kalman filter, so we just need to compute the $$\langle \boldsymbol{X}\_i, \boldsymbol{e}*j \rangle R*{e,j}^{-1}\boldsymbol{e}\_j$$.

$$\begin{aligned} \langle \boldsymbol{X}\_i, \boldsymbol{e}\_j \rangle &= \langle \boldsymbol{X}\_i, H\_j(\boldsymbol{X}*i - \hat{\boldsymbol{X}}*{j|j-1}) + \boldsymbol{V}\_j \rangle = \langle \boldsymbol{X}*i, H\_j(\boldsymbol{X}*j - \hat{\boldsymbol{X}}*{j|j-1}) \rangle \ &= \langle \hat{\boldsymbol{X}}*{i|i-1} + (\boldsymbol{X}*i - \hat{\boldsymbol{X}}*{i|i-1}), \boldsymbol{X}*j - \hat{\boldsymbol{X}}*{j|j-1} \rangle H\_j^\*\end{aligned}$$

$$\boldsymbol{X}*j - \hat{\boldsymbol{X}}*{j|j-1}$$ is orthgonal to any linear function of $${\boldsymbol{Y}*0, \cdots, \boldsymbol{Y}*{j-1}}$$, so when $$j \geq i$$, it must be orthgonal to $$\hat{\boldsymbol{X}}\_{i|i-1}$$ since it is a function of $${\boldsymbol{Y}\_k}^{i-1}\_0$$. Thus, for $$j\geq i$$,

$$\langle \boldsymbol{X}\_i, \boldsymbol{e}\_j \rangle = \langle \boldsymbol{X}*i - \hat{\boldsymbol{X}}*{i|i-1}, \boldsymbol{X}*j - \hat{\boldsymbol{X}}*{j|j-1} \rangle H\_j^\*$$

If we denote $$P\_{ij} = \langle \boldsymbol{X}*i - \hat{\boldsymbol{X}}*{i|i-1}, \boldsymbol{X}*j - \hat{\boldsymbol{X}}*{j|j-1} \rangle$$, then

$$P\_{ij} = P\_{i}\Phi\_{p}^\*(j, i) = P\_{i} \begin{cases} I & \text{ if } j = i\ \prod\_{k=i}^{j-1} F\_{p,k} & \text{ if } j > i \end{cases}$$

where $$F\_{p,k} = (F\_k - K\_{p, k}H\_k)$$. This gives us the expression

$$\hat{\boldsymbol{X}}*{i|N} = \hat{\boldsymbol{X}}*{i|i-1} + P\_i \sum\_{j=i}^N \Phi\_p^\*(j, i)H\_i^\*R\_{e,j}^{-1}\boldsymbol{e}\_j.$$

If we let $$\lambda\_{i|N} = \sum\_{j=i}^N \Phi\_p^\*(j, i)H\_i^\*R\_{e,j}^{-1}\boldsymbol{e}\_j$$, then we get the recursion

$$\lambda{i|N} = F\_{p,i}^\*\lambda\_{i+1|N} + H\_i^\*R\_{e, i}^{-1}\boldsymbol{e}\_i.$$

If we want to look at the error of this estimator, we see that

$$P\_{i|N} = P\_i - P\_i \Lambda\_{i|N} P\_i \text{ where } \Lambda\_{i|N} = \langle \lambda\_{i|N}, \lambda\_{i|N} \rangle = F\_{p,i}^\*\Lambda\_{i+1|N}F\_{p,i}+H\_i^\*R\_{e, i}^{-1}H\_i.$$

Writing all of this as an algorithm,