3
$\begingroup$

I am attempting to create a Kalman filter to track a sine wave (I am using a linear Kalman filter example assuming I already know the frequency of the sine wave) - the example I am using is derived on pages 194-196 of "Fundamentals of Kalman Filtering: A Practical Approach" 2nd edition by Paul Zarchan and Howard Musoff.

It is working to track the AC part of the signal, however the offset of the sine wave from the $x$-axis is not correct, it seems to be tied to my value of $R$, I'm guessing it is because in the derivation of the model no offset is included but I cannot see where an offset term should be inserted, how should I go about this?

The matrices describing my system are like so:

$$\mathbf \Phi_k = \begin{bmatrix} \cos(\omega T_s) & \dfrac{\sin(\omega T_s)}{\omega} \\ -\omega \sin(\omega T_s) & \cos(\omega T_s) \\\end{bmatrix}$$

Where the state matrix is:

$$\mathbf X = \begin{bmatrix} x \\ 0 \end{bmatrix}$$

And I have set $\mathbf Q$ and $\mathbf P$ like so:

$$\mathbf Q = \begin{bmatrix} \dfrac{T_S^3}{3} & \dfrac{T_s^2}{2} \\ \dfrac{T_s^2}{2} & T_s \end{bmatrix}, \quad\mathbf P = \begin{bmatrix} 9999999999999 & 0 \\ 0 & 9999999999999 \end{bmatrix} $$

I am now attemping to use this to track a generated signal of a noisy sine wave of amplitude and frequncy $1$ with an offset of $300$.

If $R$ is set to $0.1$ I get the following output of my filter:

Kalman filter output

Which is offset by $300/20$, as can be seen from this plot where I add an offset to the Kalman filter output:

Kalman filter output with offset

Changing $R$ changes this offset (as one might expect as $R$ denotes the uncertainty in your data).

How can I add this offset into my model such that my Kalman filter can tracking the sine wave correctly regardless of offset?

EDIT: I've been thinking more about my question and realised that in a real signal with noise this offset could be considered a separate constant signal with frequency $0$ which could be extracted by using a Kalman filter fitting a constant value (or simply a moving average filter) and that the offset is not inherently part of the sine wave.

  • In which case how might I remove the offset in the Kalman filter output entirely? (e.g. if I used a bandpass filter no offset would be present).
  • Since my model does not include an offset I might expect no offset at all but I do see an offset of $300/20$, as discussed. Why does this occur if there is no such offset parameter in my model?
$\endgroup$
  • 1
    $\begingroup$ It seems you need to add a term to follow the DC Level. Since it is additive, just think how would you build a Kalman Filter for DC. Actually I think the solution will match very very intuitive solution of Running Mean (Auto Regressive Form). $\endgroup$ – Royi Dec 7 '16 at 18:00
6
$\begingroup$

Well, in continuous time, a sinusoid with a bias can be seen as the output of the linear system \begin{align*} \begin{bmatrix}\dot x_1\\\dot x_2\\\dot x_3\end{bmatrix} &= \begin{bmatrix}0 & 1 & 0\\-\omega^2&0&0\\0 &0 &0\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}\\ y &= \begin{bmatrix}1&0&1\end{bmatrix} \begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix} \end{align*} As the initial condition $(x_1(0),x_2(0),x_3(0))$ vary in $\mathbb R^3$, the system generates all the possible $y(t)$ of the form $$ y(t) = a\cos(\omega t+\phi) + b $$ The system is observable for all $\omega\ne 0$,therefore a Kalman Filter is possible. With $x:=col(x_1,x_2,x_3)$ and $$ A=\begin{bmatrix}0 & 1 & 0\\-\omega^2&0&0\\0 &0 &0\end{bmatrix},\qquad C=\begin{bmatrix}1&0&1\end{bmatrix} $$ The filter has the form \begin{align*} \dot{\hat x} = A\hat x -L(C\hat x-y) \end{align*} Where $L$ is chosen according to the Kalman theory, on the basis of the covariance matrices of the noises and the solution of the Riccati Equation. For discrete time systems the procedure is the same.
EDIT:-----
If you only use a model for the oscillator you are actually restricting the $A$ map to the $x_1$-$x_2$ subspace. That's equivalent to take $C=[1\;0\;0]$ in the above system. In this latter case you loose observability (actually you do not have neither detectability) and therefore you cannot observe the bias.

$\endgroup$
  • $\begingroup$ In this case what is x1, x2 and x3? $\endgroup$ – SomeRandomPhysicist Dec 8 '16 at 14:07
  • 1
    $\begingroup$ Sorry I updated the equation. $x$ is a state space representation of the model you use to represent the process that generates the biased sinusoid. The Kalman filter is model based, thus you ll need a model $\endgroup$ – LJSilver Dec 8 '16 at 14:21
  • $\begingroup$ ahh, yes, in my model without a bias I simply had $A = \begin{bmatrix} 0 & 1 \\ -\omega^2 & 0 \end{bmatrix}$ with the state space representation $x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}$ where $x_1 = x$ and $x_2 = \dot{x}$ should I simply add a term to this as you have done so that $x_3 =$ Offset? $\endgroup$ – SomeRandomPhysicist Dec 8 '16 at 14:31
  • $\begingroup$ Currently I have the model $\begin{bmatrix} \dot{x} \\ \ddot{x} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -\omega^2 & 0 \end{bmatrix} \begin{bmatrix} x \\ \dot{x} \end{bmatrix}$ $\endgroup$ – SomeRandomPhysicist Dec 8 '16 at 14:33
  • 1
    $\begingroup$ yes, you should redo everything with the new model. Essentially your model only captures the oscillating term (in fact it is equivalent to take $C=[1\;0\;0]$ in my answer. If this is the case then you loose observability! and the bias is invisible). The third state is $\dot x_3=0$, that means that $x_3(t)=x_3(0)=bias$. With the model I suggest you are going to observe also the bias, no matter how big is the initial error. $\endgroup$ – LJSilver Dec 8 '16 at 17:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.