Determining the mean of of a sampled signal (improvements to time averaging method)
Edit: After thinking about it some more I came up with something much simpler than the phase-locked loop.

The problem you are having is because you are filtering with a boxcar.  The boxcar filter has a lot of ripples in the frequency domain, so if you choose the wrong width you don't get good attenuation of your approximately 10Hz signal.

If you use a Butterworth filter you will get a frequency response that has no ripples.  Everything over (say) 5Hz will be attenuated by more than -3db.  A Butterworth filter is also cheap to calculate.

I went to http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html and asked for a 3rd order Butterworth low-pass filter with a cutoff frequency of 5Hz on a sample rate of 1000Hz and got the following recurrence relation:

y[n] = (  1 * x[n- 3])
     + (  3 * x[n- 2])
     + (  3 * x[n- 1])
     + (  1 * x[n- 0])

     + (  0.9390989403 * y[n- 3])
     + ( -2.8762997235 * y[n- 2])
     + (  2.9371707284 * y[n- 1])


You'd probably get a better result by using Matlab to design an elliptic filter (and probably one with higher order than 3).  You'll get sharper attenuation beyond the cutoff frequency.

Here is my original answer about phase-locked loops:

I would try a discrete time phase-locked loop.

So something like this:



The idea is to multiply your input signal by a sinusoid of the same frequency, followed by a low pass filter.  This shifts the sinusoidal part of your input signal to (nearly) zero frequency.  By observing the changes to the output of the low pass filter you are getting an estimate of how far off your estimated frequency is from the actual frequency.  So you feedback the adjustment.

The output of the multiplier in the picture above is:
$$
Be^{2\pi i \theta_n} + Ae^{2 \pi i (\theta_n + \omega n)} + Ae^{2 \pi i (\theta_n - \omega n)} + \varepsilon'(n).
$$
$\theta_n$ is the running sum of the estimated frequency, $\hat{\omega}_n$.  When $\hat{\omega}_n = \omega$ then $\theta_n-\omega n$ is constant, so $A e^{2 \pi i (\theta_n - \omega n)}$ is constant.

The derivative you want here is the angular derivative.  Given samples from the low pass filter, $r_ne^{2\pi i \alpha_n}$ and $r_{n+1}e^{2\pi i \alpha_{n+1}}$ you really want $\alpha_{n+1}-\alpha_{n}$.  Dividing low-passed sample $n+1$ by low-passed sample $n$ and taking the complex part gets you essentially the right thing.

The "adjust $\hat{\omega}_n$" box is also non-trivial.  The estimated phase error is likely to be a little noisy so you might want to low-pass filter the phase error and you might want to only adjust $\hat{\omega}_n$ by a fraction of the estimated phase error, rather than by the entire phase error.

$\hat{\omega}_n$ is the estimated frequency, and $1/\hat{\omega}_n$ is the low-pass filter width you are looking for in your question.  I think you can apply the same trick to the low-pass filter in the phase-locked loop.  By dynamically adjusting the phase-locked loop low-pass filter width to $1/\hat{\omega}_n$ you completely eliminate $Be^{2\pi i \theta_n}$.