Estimating a Low Frequency Signal Corrupted by High Frequency Noise
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 &amp; 1 \\
0 &amp; 1
\end{array} \right]\\
\mathbf{H} = \left [ \begin{array}{cc}
1 &amp; 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 &amp; 2.4 \\
2.4 &amp; 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 &lt;- 128;
t &lt;- 0:(T-1)/T*70

signal &lt;- 4.4 + sin(0.06*pi*t - 0.1) + sin(0.14*pi*t + 0.07) + sin(0.2*pi*t + 0.4)
noise &lt;- 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 &lt;- signal + noise


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

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

zhat &lt;- t*0

for (k in 1:T)
{
  xkkm1 &lt;- A %*% xkm1km1
  Pkkm1 &lt;- A %*% Pkm1km1 %*% t(A) + Q
  K &lt;- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  z &lt;- matrix(c(measurement[k]), 1, 1)
  xkm1km1 &lt;- xkkm1 + K %*% (z - H %*% xkkm1)
  Pkm1km1 &lt;- (matrix(c(1,0,0,1),2,2) - K %*% H) %*% Pkkm1  
  zhat[k] &lt;- 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 &lt;-  filter(measurement, c(0.2,0.2,0.2,0.2,0.2), circular = TRUE)
lines(t,lpf, col="green", lwd=4)

errs &lt;- c(sum((signal - zhat)^2), sum((signal - lpf)^2))
print(errs)