Estimating a Low Frequency Signal Corrupted by High Frequency Noise

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.

• Why do you use a Kalman filter at all if your noise and signal are that well separated in frequency? A simple low pass will probably work better for what you want. – Jazzmaniac Oct 18 '15 at 8:40
• @Jazzmaniac, I corrected a typo in the text. The signal and noise have less frequency separation than I fist indicated. I didn't use a discrete lowpass filter because I wanted to limit my self to algorithms that are in the books I found on estimation. See the paragraph I added at the end of my question. – Ted Ersek Oct 18 '15 at 13:35
• Your Correct equation for $\mathbf{X}(k)$ seems to be missing a $\mathbf{X}(k-1)$ term (or something like that). Your KF output appears to be delayed; I wouldn't have expected that. It is likely to be adding to your error. – Peter K. Oct 19 '15 at 18:23
• Peter K is right. You have not posited noise but a contaminating signal. Given the lack of uncertainty I could design a perfect filter; zeroing out the unwanted signals and passing the known (?) signal with perfect fidelity. In fact: since there is no uncertainty you could throw away the input and just generate the signal you want. The point is the problem needs to be better defined; then I am sure somebody here can answer with an optimum technique. – rrogers Oct 21 '15 at 14:46
• i use hampel filter if i have a set of points following a trend . it works everytime . in real time , there will be some delay depending on your window size – user18956 Jan 5 '16 at 2:29

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)