Determining the mean of of a sampled signal (improvements to time averaging method)

Question

I am sampling a signal at 1kHz, which is essentially a single-frequency sinusoid with a DC offset (and a small amount of noise). I don't know the precise frequency of the sinusoid (it is at most 10Hz), and I want to record the mean value of the signal.

Currently I take samples for an arbitrary number of seconds and find the average value (with a longer sample time corresponding to a lower error value, clearly):

$$ \langle x \rangle = \frac{1}{N} \sum_{i=0}^N x_i $$

I feel like there should be a better way to do this, however, given that I know something about the form of the signal (namely, it consists of a single-frequency sine wave). If I could determine the frequency, then averaging over a complete number of periods would be an improvement.

Using a FFT to find the frequency seems to be tricky though, since (as far as I can tell) the discrete fourier transform will yield values that are periodic in the frequency domain (ie. the discrete frequencies returned are $ f =i \times f_s \mathbin{/} N$, where $ f_s $ is the sampling frequency; $ f $ is thus always periodic within the $ N $ samples I've measured).

Can anyone suggest a good way to do this - either via FFT, or a better method of mean-determination?

Answer 1

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}$.

Answer 2

tl;dr: The last couple of paragraphs actually give possible answers to the question, including windowing.

With any situation where you try to improve upon a method, it is worthwhile looking at how much improvement is available. That way you can decide how valuable the improvement is, to decide whether it is worth pursuing.

The error you are trying to avoid is basically that introduced by having partial cycles of the sinusoid in the signal. If there are a complete number of periods, then taking the time average is the best you are going to get, since the sinusoid disappears in the sum (I'm assuming zero-mean noise here).

The plot below shows the amount of error that will be introduced by incomplete cycles in the signal. This shows the mean of a sinusoidal signal (unit amplitude) over a 1 second, with $f_s = 1$kHz. If you sample for more than one second, the error will reduce (as you note), and the error will depend on the amplitude of the sinusoid. For reference, with two seconds of sampling, you end up with this much error:

You note correctly that you won't necessarily be able to measure the exact frequency of your signal with a DFT. One method that springs to mind for resolving this could be to calculate the FFT of your signal after chopping various numbers of samples from the end. Do this to minimise the amount of spectral leakage (skirting) in the spectrum, then you will have found the closest approximation to a complete number of cycles of the sinusoid, and can assume the signal mean is about the best you're going to get.

You can detect minimisation of the leakage by identifying and removing the peak value from the spectrum, then finding the mean of the surrounding few values in the spectrum. When that mean is minimal (it starts rising again as you continue to remove samples from the signal), you have your sweet spot.

EDIT: forgot to mention, another simple method would be to window your signal (for example, with a Hamming window). Then find the average of the windowed signal, and scale it back up to counter the effect of the window. The scaling value should be $$ \frac{N}{\sum_{n=1}^N w[n]} $$ where $N$ is the number of samples in your signal. Experiment with different windows to see what works best.

Answer 3

Your sampled signal looks like

$$x(n) = A\sin\left (\frac{2\pi fn}{f_s} + \phi\right ) + B + e(n)$$

where $A$ is the amplitude, $\phi$ is the phase, $f$ is the frequency of the sinusoidal signal, $f_s$ is the sampling frequency, $B$ is the DC offset, and $e(n)$ is some additive noise, hopefully with zero mean. $A$, $f$, $\phi$, $B$, and $e(n)$ are unknown, and you would like to know $B$.

Your idea of averaging the signal is a perfectly valid choice for estimating $B$ (if the additive noise $e(n)$ has indeed a mean value of zero). The problem you have noted is of course that you never know if there is an integer number of periods in your sampled data. To improve on your solution you could try the following. Take a large block of data. Your block length $N$ must satisfy $$N\gg \frac{f_s}{f}\tag{1}$$ Then create sub-blocks of sufficient length (still satisfying (1)). To get your sub-blocks simply cut out overlapping blocks from your original block of data. In this way you introduce a random phase of your signal with respect to the beginning and end of the data blocks, and if you average over each sub-block, and then average the results of all sub-blocks, you'll get a better total estimate. Say your estimate based on the $i^{th}$ sub-block is

$$\hat{B_{i}}=\frac{1}{M}\sum_{n=n_i}^{n_i+M-1}x(n)$$

with $M$ the number of data samples in the sub-blocks, then your final estimate of the DC offset would be

$$\hat{B} = \frac{1}{L}\sum_{i=0}^{L-1}\hat{B_{i}}$$

where $L$ is the number of sub-blocks.