10

Suppose a system follows this equation: $$ x(t)=A \cos(\omega t + \phi)+\eta$$

where: $\omega = 2\pi f $ and $\eta$ is a random error

using Extended Kalman Filter, how does estimated value $\hat{x}$ be?

Royi
  • 19,608
  • 4
  • 197
  • 238
unwantednoise
  • 103
  • 1
  • 5

2 Answers2

7

This isn't quite what you're asking, because it neglects the amplitude, $A$, but it's a relatively straightforward example of application of an extended Kalman filter to the frequency tracking problem. See section 1.2 of this PDF, that I wrote some time ago.

I'd also recommend starting with B. D. O. Anderson and J. B. Moore, Optimal Filtering, Prentice-Hall, Inc., Engle- wood Cliffs, New Jersey 07632, 1979.

EKF for frequency tracking

Peter K.
  • 25,714
  • 9
  • 46
  • 91
  • I like the reference Peter, your PDF is a nice summary – Dan Boschen Dec 06 '21 at 13:31
  • 1
    @DanBoschen For an unpublished article (in the journal or conference sense), that PDF has received more citations than the IEEE TSP paper for my PhD. Oh well. ;-) – Peter K. Dec 06 '21 at 15:53
  • 1
    I'm giving a callout to Dan Simon's "Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches", Wiley, 2006. I think it'll be clear to someone who's taken a senior-level statistics class and state-space control. The only downside is that after it was published someone came up with a formal way to determine the constellation for an Unscented Kalman filter, and I can't remember who wrote the paper or when (aside from "after 2006"). – TimWescott Dec 06 '21 at 20:24
  • @TimWescott thanks,Tim! I’ll see if I can get that. – Peter K. Dec 06 '21 at 20:47
  • @TimWescott, what do you mean by the constellation of the Unscented Kalman Filter? As the UKF itself was published by Julier and Uhlmann before 2000. – Royi Dec 07 '21 at 04:45
  • The set of sigma points. If I recall correctly there was some result that came out after the Simon book was written that found a formal result for some optimal, or maybe generally better, constellation of sigma points. And -- now I can't remember what it was. – TimWescott Dec 07 '21 at 21:51
  • 1
    @TimWescott, Could it be you're talking about Cubature Kalman Filter? If so, then it is a generalization of the UKF and actually it shows that in most cases the UKF is the optimal constellation. – Royi Dec 08 '21 at 11:24
4

I'm copying my answer to Estimate and Track the Amplitude, Frequency and Phase of a Sine Signal Using a Kalman Filter which solves a more general problem with example code:

We can build a non linear dynamic model in order to estimate the parameters of a sine signal.

Let's model the signal as $ a \sin \left( \phi \right) $ where $ \phi $ is the instantaneous phase. So the model could be also written as $ a \sin \left( \omega t + \psi \right) $.

Then the model can be:

$$ {a}_{k} \sin \left( {\omega}_{k} {t}_{k} + \psi \right) = {a}_{k} \sin \left( {\phi}_{k} \right) $$

With some math and pre processing of Kalman Filter you may derive the model with the matrices:

$$ \boldsymbol{x}_{k} = \begin{bmatrix} {a}_{k} \\ {\omega}_{k} \\ {\phi}_{k} \end{bmatrix}, F = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & \Delta t & 1 \end{bmatrix}, Q = \begin{bmatrix} \Delta t {\sigma}_{a}^{2} & 0 & 0 \\ 0 & \Delta t {\sigma}_{\omega}^{2} & \frac{ {\Delta t}^{2} {\sigma}_{\omega}^{2}}{2} \\ 0 & \frac{ {\Delta t}^{2} {\sigma}_{\omega}^{2}}{2} & \frac{ {\Delta t}^{3} {\sigma}_{\omega}^{2}}{3} \end{bmatrix} $$

Where $ {\sigma}_{a}^{2} $ is the process variance of the amplitude and $ {\sigma}_{\omega}^{2} $ is the variance of the process noise of instant angular frequency.

The measurement model is a bit more tricky. The measurement model is:

$$ {z}_{k} = h \left( \boldsymbol{x}_{k} \right) = {a}_{k} \sin \left( {\phi}_{k} \right) $$

Hence the Jacobian is given by $ \frac{\partial h \left( \boldsymbol{x}_{k} \right )}{\partial \boldsymbol{x}_{k}} = \left[ \sin \left( {\phi}_{k} \right), 0, {a}_{k} \cos \left( {\phi}_{k} \right) \right] $.

Wrapping all this into a Kalman Model will yield:

enter image description here

You may see that the model can effectively track changes in the parameters.
There are other alternatives to this dynamic model but I think this is a simple and effective one.

You may also use the Unscented Kalman Filter. I implemented it at Extended Kalman Filter (EKF) for Non Linear (Coordinate Conversion - Polar to Cartesian) Measurements and Linear Predictions.

The code is available at my StackExchange Signal Processing Q76443 GitHub Repository (Look at the SignalProcessing\Q76443 folder).

Royi
  • 19,608
  • 4
  • 197
  • 238
  • Nice answer! :-) – Peter K. Dec 06 '21 at 18:46
  • 1
    @PeterK., Thanks. There are many ways to derive the model matrix for harmonic signals. – Royi Dec 06 '21 at 18:50
  • Thanks for the intuitive answer, I have a question about dT, and paramAngFreq, how do you select the appropriate dT.? do we follow Nyquist theorm, or the sampling frequency (1/dT) is at least the same with paramAngFreq.? – unwantednoise Dec 08 '21 at 09:53
  • 1
    The parameters dT is the rate new measurement is measured. Of course it has to be at least by Nyquist. You may use dT = 1 and then everything is in normalized frequency. – Royi Dec 08 '21 at 11:23
  • Then when ω=2πf on a.sin(ωt+ψ)=a.sin(ϕ), if we have fundamental frequency f=100Hz, thus the dT become atleast (1/200=0.005).? – unwantednoise Dec 08 '21 at 12:21
  • Yep. But it must be the period of sampling. What you said is basically what's the Nyquist rate. But it could be the dT = 0.001 if that's the rate new samples are sampled. – Royi Dec 08 '21 at 12:40
  • Could you explain a little bit in detail about that.?, does it mean the period of sampling is not the same as the Nyquist sampling rate.? So when the fundamental frequency is 100, the Nyquist sampling rate is at least 0.005. and why we select dT=0.001 .? What happens if we select another fundamental frequency, for example, 200, does it become 0.0025? – unwantednoise Dec 08 '21 at 12:47
  • 1
    I am not sure what you mean. Assume you have which measures a sine signal. By Sampling Theorem you must sample it at sampling rate which is larger than 2 times the bandwidth. But in many cases we sample at much higher rate. In any case, dT must be the sampling rate in practice and not the Nyquist rate. – Royi Dec 08 '21 at 16:47
  • @Royi thanks a lot. If anyone interested, i've made a wed interactive demo of this in JS at https://observablehq.com/d/a033acc0859cc0de based on https://github.com/piercus/kalman-filter – piercus Mar 24 '23 at 17:47
  • 1
    @piercus, This is really nice! I wish I knew JS. – Royi Mar 24 '23 at 18:11