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

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?

• What if you look for zero crossings of the derivative (difference between samples)? This would give a guess as to where the peaks and troughs of the sine wave are, and therefore the period. Averaging over an integer multiple of the period would then give a measure of the DC component. Apr 29, 2013 at 4:39
• Computing the derivative of a noisy signal is very error prone. Apr 29, 2013 at 6:56
• @Matt, yes that is true. The phase from an appropriate quadrature filter could get around it. Apr 29, 2013 at 7:50
• Interesting question! This paper give a maximum likelihood approach, but it takes great pains to first estimate the frequency of your sinusoid.
– Peter K.
Apr 29, 2013 at 12:08

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.

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

• Thanks for this! Can you explain at a slightly higher-level exactly what the Butterworth filter is achieving (or how practically one applies this)? Is this just helping to more accurately identify the peak frequency when I carry out a Fourier transform? Apr 29, 2013 at 23:26
• The Butterworth filter is just a more accurate way to get a moving average ($\langle x \rangle$). The problem with the boxcar ($\frac{1}{N}\sum x_i$) is that its frequency response is a sinc function. It doesn't do a good job of attenuating all the high frequencies. I'm suggesting that you use a filter that gives better guarantees about attenuating frequencies > 5Hz or so. Apr 30, 2013 at 3:10

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.

• Your error plots seem to suggest that the estimation error is always positive, i.e. that mean estimation by averaging is biased. It is important to point out that this is not true. In this sense, your plots are a bit misleading. I suppose that you used a sine wave with zero phase and your averaging always starts at $t=0.$ Then you'll always get a positive error, but this is not a realistic scenario. Apr 29, 2013 at 9:59
• Can you possibly elaborate on your windowing suggestion? I understand that it will reduce spectral leakage when I take the FFT, but this in itself doesn't get around the problem of the discrete frequency values (which your suggestion of removing samples seems to achieve). Thanks! Apr 29, 2013 at 10:24
• @Matt Fair point, and quite right. The main thing I was trying to show was a figure for how much error to expect, and I think the plots give a good estimate for a maximum expected error.
– lxop
Apr 29, 2013 at 22:02
• @Edward the windowing will simply reduce the influence of the incomplete cycle on the final result, so the error from the incomplete cycle will be smaller
– lxop
Apr 29, 2013 at 22:04

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

• This intuitively makes sense - does this approach have a name? It seems like a downside is that N will have to be very large when f is quite low, so the sampling time will be high. Apr 29, 2013 at 10:49
• Not sure about the name. I'm afraid you can't get around the large N, because your frequency is very low compared to the sampling frequency. You obviously need a certain time to get a sufficient number of periods in your data window, and this time is quite large for a low frequency signal as yours. Apr 29, 2013 at 10:55
• It seems to me like this approach would be mathematically equivalent to windowing the original signal with a piecewise linear window (rising, then flat, then falling), and finding the average of the windowed signal. That could make implementation a little easier
– lxop
Apr 29, 2013 at 21:58