Estimating a Low Frequency Signal Corrupted by High Frequency Noise

Question

This is a follow up to the question How to Pick a Journal for Publishing a Research Work? question. Suppose we take discrete samples of a low frequency signal corrupted by high frequency noise. We know the signal has no frequencies higher than 20% of the Nyquist rate. We also know all the noise is at or above 50% of the Nyquist rate, and below the Nyquist rate. I will use the following as a specific example to experiment with.

\begin{align}&\textrm{measurement}(t) = \textrm{signal}(t) + \textrm{noise}(t)\\ &\textrm{signal}(t) = 4.4+\sin(0.06\pi t-0.1) + \sin(0.14\pi t+0.07) + \sin(0.2\pi t+0.4)\\ &\textrm{noise}(t) = \sin(0.06+0.5\pi t) + \sin(0.1-0.68\pi t) - \sin(0.04-0.74\pi t) + \sin(0.03-0.93\pi t) \end{align}

The measurement is sampled at integer values of ($t$). A plot of the measurement samples and the continuous signal is below.

How well can we estimate $\textrm{signal}(t)$ at integer values of ($t$) using the current and past measurements? I was able to filter it using the Kalman filter below.

A plot of the signal and Kalman filter output is in the next plot. The final plot shows the error in the Kalman filter output which has a maximum magnitude of 1.865 and the RMS error is 0.7895. Can we make a Kalman filter that does better than the Kalman filter above? What known filter method would give the best estimates, and how good would the estimates be?

I learned about discrete lowpass filters 25 years ago, and I always assumed a problem such as the one above was a common application for them. At that time I didn't know this problem is called estimation, and I didn't know about the Kalman filter. Then I recently took this course, and another course Georgia Tech provided that goes into the details of estimation algorithms. However, the Georgia Tech courses mention nothing about using a discrete lowpass filter for estimation. I check several books on estimation, and found nothing about discrete lowpass filters. I work with a computer scientist who has been doing research in data fusion for 20 years, and he knows nothing about z-transforms, transfer functions, etc. Why is it that the estimation community ignores the use of a discrete lowpass filter for estimation? As far as I can tell it's the best approach to estimate the signal above.

Answer 1

To answer your final question:

Why is it that the estimation community ignores the use of a discrete lowpass filter for estimation? As far as I can tell it's the best approach to estimate the signal above.

That's because you're feeding them the wrong signal model.

TL;DR: Use the right tool for the job!

Gory Details

Any time you start in with the Kalman Filter, you are assuming that the signal (measurement in your equations above) was generated as: $$ x[k+1] = \mathbf{A} x[k] + u[k]\\ \mbox{measurement}[k] = \mathbf{H} x[k] + v[k] $$ where $$ \mathbf{A} = \left [ \begin{array}{cc} 1 & 1 \\ 0 & 1 \end{array} \right]\\ \mathbf{H} = \left [ \begin{array}{cc} 1 & 0 \end{array} \right]\\ $$ and $u$ and $v$ are zero mean, Gaussian distributed, independent random variables with covariances $$ \mathbf{Q} = \left [ \begin{array}{cc} 1.6 & 2.4 \\ 2.4 & 4.8 \end{array} \right] \times 10^{-9}\\ \mathbf{R} = 2.35 \times 10^{-6}\\ $$ At least for the equations you cite. More general application of the KF is possible, but usually not done.

So, one thing that might be done to improve your estimates is to ensure that $\mathbf{Q}$ and $\mathbf{R}$ are closer to reality. From the look of things, your process noise ($u[k]$) has about the same variance as your measurement noise ($v[k]$). Make them closer.

More importantly, your signal is correlated with your noise because, rather than stochastic noise, you have harmonic noise. Even worse: you have synchronized harmonic noise (meaning the noise moves in lockstep with your signal).

So: the problem is that by using the Kalman Filter your signal model is all wrong.

A better model for the signal is in the frequency domain: the signal is low frequency. The noise is high frequency. So, a low pass filter will do a much better job of improving things that the Kalman Filter can hope to do.

If I apply the Kalman filter to your problem (with some modifications to your equations to take account of the difference between $\mathbf{X}[k|k-1]$ and $\mathbf{X}[k|k]$) then I get the picture below.

The black curve with dots is the measurement. The red curve is signal. The blue curve is the output of the Kalman filter. The green curve is the output of a 5 point moving average filter (a simple low pass filter).

The sum of squared errors for the Kalman filter is 118.63722. The same figure for the low pass filter is 14.18046. Clearly, the signal model of a low pass filter fits better because the error is smaller.

R code to implement this is below.

---

# 26489
T <- 128;
t <- 0:(T-1)/T*70

signal <- 4.4 + sin(0.06*pi*t - 0.1) + sin(0.14*pi*t + 0.07) + sin(0.2*pi*t + 0.4)
noise <- sin(0.06 + 0.5*pi*t) + sin(0.1 - 0.68*pi*t) - sin(0.04 - 0.74*pi*t) + sin(0.03 - 0.93*pi*t)

measurement <- signal + noise


xkm1km1 <- matrix(c(4.4, 0),2,1)
Pkm1km1 <- matrix(c(1,0,0,1),2,2)
H <- matrix(c(1,0),1,2) 
A <- matrix(c(1,1,0,1),2,2)
Q <- 10^-6 * matrix(c(1.6,2.4,2.4,4.8),2,2)
R <- 10^-6 * matrix(c(2.35),1,1)

library("MASS") # For pseudo inverse ginv()

zhat <- t*0

for (k in 1:T)
{
  xkkm1 <- A %*% xkm1km1
  Pkkm1 <- A %*% Pkm1km1 %*% t(A) + Q
  K <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  z <- matrix(c(measurement[k]), 1, 1)
  xkm1km1 <- xkkm1 + K %*% (z - H %*% xkkm1)
  Pkm1km1 <- (matrix(c(1,0,0,1),2,2) - K %*% H) %*% Pkkm1  
  zhat[k] <- as.numeric(H %*% xkkm1)
}

plot(t, measurement, type="b", pch=19)
lines(t, signal, col="red",  lwd=10)
lines(t, zhat, col="blue", lwd=5)

lpf <-  filter(measurement, c(0.2,0.2,0.2,0.2,0.2), circular = TRUE)
lines(t,lpf, col="green", lwd=4)

errs <- c(sum((signal - zhat)^2), sum((signal - lpf)^2))
print(errs)