Shlomo Engelberg
Moving average filters output the average of N samples, and it is easy to see (and to prove) that they are low-pass filters. A simple, causal moving average filter satisfies \[{y}_{n} = \frac{1}{N}\mathop{\sum}\limits_{{k} = {0}}^{{N}{-}{1}}{{x}_{{n}{-}{k}}}{.} \tag{1} \]
Because of their simplicity and intuitive appeal, they are often preferred to more complicated low-pass filters when one would like to remove high-frequency noise and the demands on the filter are not too great.
In this column, we explain how implementing a moving average filter using a recursive formulation (where the current output is a function of the current and previous inputs and previous outputs) is possible. We then show why it can be problematic and how to deal with that problem (and this is the “tip”). Finally, we describe circumstances under which the problem does not exist even though one might have thought that it should (and this is the “trick”).
Let ${x}_{k} = {x}{(}{kT}_{s}{)}$, ${T}_{s}\gt{0}$ be samples of x(t). One can define a moving average filter by using (1), and this can be used as a template for implementing a moving average filter; looking at (1), one would be inclined to calculate each sample of the output by summing the last N values sampled and dividing the sum by N.
One can, however, express the output of a moving average filter, yn, as \begin{align*}{y}_{n} & = \frac{1}{N}\mathop{\sum}\limits_{{k} = {0}}^{{N}{-}{1}}{{x}_{{n}{-}{k}}} \\ & = {x}_{n}{/}{N} + \frac{1}{N}\mathop{\sum}\limits_{{k} = {1}}^{{N}{-}{1}}{{x}_{{n}{-}{k}}} \\ & = {x}_{n}{/}{N} + \frac{1}{N}\mathop{\sum}\limits_{{k} = {0}}^{{N}{-}{2}}{{x}_{{(}{n}{-}{1}{)}{-}{k}}} \\ & = {x}_{n}{/}{N} + \frac{1}{N}\mathop{\sum}\limits_{{k} = {0}}^{{N}{-}{1}}{{x}_{{(}{n}{-}{1}{)}{-}{k}}}{-}{x}_{{n}{-}{N}}{/}{N}\end{align*} which is equivalent to \[{y}_{n} = {x}_{n}{/}{N} + {y}_{{n}{-}{1}}{-}{x}_{{n}{-}{N}}{/}{N} \tag{2} \] and this points to a different possible implementation. One can implement the filter by setting the current output, yn, equal to the sum of the previous output value, ${y}_{{n}{-}{1}}$, and the current sample of the input divided by ${N},{x}_{n}{/}{N}$, less the oldest sample of the input that was “part of” the previous value of the output divided by ${N},{x}_{{n}{-}{N}}{/}{N}$. If one would like to minimize the number of operations needed to calculate each value of the output, yn (if an efficient implementation of the moving average filter is important), the method based on (2) is to be preferred. (See, for example, [2].)
Is a moving average filter stable? The short answer is that it must be. Considering (1), it is clear that the filter is a finite-impulse-response (FIR) filter, and all FIR filters are stable [1].
Considering the (two-sided) Z-transform of (2), we find that \[{Y}{(}{z}{)} = {X}{(}{z}{)/}{N} + {z}^{{-}{1}}{Y}{(}{z}{)}{-}{z}^{{-}{N}}{X}{(}{z}{)/}{N}{.}\]
The transfer function of the filter would seem to be \[{H}{(}{z}{)} = \frac{1}{N}\frac{{1}{-}{z}^{{-}{N}}}{{1}{-}{z}^{{-}{1}}},{z}\ne{0},{1}{.}\]
Note, however, that the apparent singularity at ${z} = {1}$ is removable, and in fact, it is not properly speaking a singular point of the transfer function. Dividing the numerator of H(z) by the denominator, one finds that \[{H}{(}{z}{)} = \left({{1} + {z}^{{-}{1}} + \cdots{z}^{{-}{(}{N}{-}{1}{)}}}\right){/}{N},{z}\ne{0}\] and this is the version of the transfer function one arrives at if one starts from (1) and proceeds naïvely. Here it is clear that H(0) is perfectly well defined and is equal to 1. If one would like to write the transfer function in closed form, the correct way to do so is to write \begin{align*}{H}{(}{z}{)} = \begin{cases}{\begin{array}{cc}{\frac{1}{N}\frac{{1}{-}{z}^{{-}{N}}}{{1}{-}{z}^{{-}{1}}},}&{{z}\ne{0},{1}} \\ {1}&{{z} = {1}}\end{array}}\end{cases}{.} \tag{3} \end{align*}
When written this way, it is clear that 1 is contained in the transfer function’s region of convergence, and this is yet another way to see that the filter is stable [1]. As we find in the section “An Implementation Issue with the Recursive Implementation,” when rounding errors are added to the picture, the system’s stability becomes a somewhat more involved question.
Floating point numbers such as C’s float data type are versatile, but when adding a small floating point number to a large one, one generally loses some accuracy. Floating point numbers generally have a certain number of bits allocated to storing a “number,” called the mantissa, and a certain number of bits allocated to storing an exponent, the power of two by which the mantissa, a binary fraction, is multiplied. If num is the number to be stored, we find that ${\text{num}} = {\text{mantissa}}\,{\times}\,{2}^{\text{exponent}}$. The absolute value of the mantissa is often required to be greater than or equal to one and less than two, so that its first digit is always a one (which makes it possible not to store that binary digit).
In order to fix ideas, consider a simple, nonpractical, example. Suppose one has four bits dedicated to the mantissa and two to the exponent. Consider the sum ${1}{.}{01}{\text{b}}\,{\times}\,{2}^{{11}{\text{b}}} + {1}{.}{01}{\text{b}}\,{\times}\,{2}^{{00}{\text{b}}}$. Writing the summands out in binary, we find that we are adding $1010\text{b}$ to $1.010\text{b}.$ As the mantissa is limited to four bits, the sum is ${1011}{\text{b}} = {1}{.}{011}{\text{b}}\,{\times}\,{2}^{{11}\hspace{1pt}{\text{b}}}$. Because of the limited number of bits the mantissa has, the sum is not accurate. It is off by $0.01\text{b}$.
Suppose that one uses the recursive implementation based on (2), but at each stage, there is a small amount of noise caused by rounding errors. Then (2) becomes \[{y}_{n} = {x}_{n}{/}{N} + {y}_{{n}{-}{1}}{-}{x}_{{n}{-}{N}}{/}{N} + {r}_{n} \tag{4} \] where rn is the noise term.
The Z-transform of (4) is \begin{align*}{Y}{(}{z}{)} & = {X}{(}{z}{)/}{N} + {z}^{{-}{1}}{Y}{(}{z}{)} \\ & \quad{-}{z}^{{-}{N}}{X}{(}{z}{)/}{N} + {R}{(}{z}{)}\end{align*} or \begin{align*}{N}{(}{z}{-}{1}{)}{Y}{(}{z}{)} & = {(}{z}{-}{z}^{{-}{(}{N}{-}{1}{)}}{)}{X}{(}{z}{)} \\ & \quad + {NzR}{(}{z}{).}\end{align*}
Thus, we find that \begin{align*}{Y}{(}{z}{)} & = \frac{1}{N}\frac{{z}{-}{z}^{{-}{(}{N}{-}{1}{)}}}{{z}{-}{1}}{X}{(}{z}{)} \\ & \quad + \frac{z}{{z}{-}{1}}{R}{(}{z}{),}{z}\ne{0},{1}{.}\end{align*}
It is clear that the first term on the right-hand side is the same as the transfer function in (3), and it is the transfer function of a stable system. The second term is, however, another story. This term truly has a single pole at z = 1, and this is the sign of a marginally stable system. In fact, \[\frac{z}{{z}{-}{1}}{,|}{z}{|}{>}{1}\] is the transfer function of a summer. Thus, the output of a moving average filter that is implemented recursively and that starts operating at time n = 0 is composed of the desired output (related to the input signal, xn) and the sum of all the previous rounding errors; it is \[{y}_{n} = \frac{1}{N}\mathop{\sum}\limits_{{k} = {0}}^{{N}{-}{1}}{{x}_{{n}{-}{k}}} + \mathop{\sum}\limits_{{k} = {0}}^{n}{{r}_{k}}{.}\]
If there are no rounding errors, the last term does not cause any problems. If there are rounding errors, as one would expect when adding and subtracting floating point numbers, for example, then one expects to see the sum of these (generally very small errors) affect the output of the system by causing a small change to the output, and one expects the maximum size of this change to increase with time.
As yn is the sum of N terms of the form ${x}_{n}/N$ and of noise-related terms, yn should be substantially larger than any given ${x}_{n}/N$. As we saw in the section “Rounding Errors in Floating Point Calculations,” adding small floating point numbers to large ones can lead to rounding errors. Thus, (4), which includes a noise term, seems to be an appropriate model for a moving average filter implemented using floating point numbers, the moving average filter is marginally stable with respect to rounding errors, and we expect the output voltage to tend to drift slightly with time as the sum of the (very, very small) rounding errors changes (very slowly).
If one stores all of one’s numbers as fixed-point numbers or integers, then as long as there is no overflow or malfunction when adding and subtracting, there is no “rounding error,” (2) is a complete description of the filter, and the filter is stable. Thus, if one works entirely in integer or fixed point arithmetic, the filter is properly characterized by (2), and one can implement the filter recursively without any stability issues.
Thus, the tip is to use some form of fixed point or integer arithmetic when implementing a moving average filter via recursion. When working in the C programming language, for example, use ints rather than floats in the moving average filter. (See [2], for example, for a somewhat different presentation of the problem and this solution.)
When implementing a moving average filter in C, it is often actually convenient to work with C ints rather than C floats. On the microcontroller we used (the ADuC841, a member of the 8051 family), when working with ints, it is convenient to store the (integer) values read from the analog to digital converter (ADC) in an array and not convert the values to voltages. When the time comes to output values to the digital to analog converter (DAC), the values are already properly scaled, as the ADC and DAC can be set to scale values in the same way. Working with ints makes a lot of sense from the point of view of effective and efficient programming, even if it leads to there being no point at which the int-based program “knows” the numerical value of the voltage of the input or the output (in volts) and to the engineer needing to work with unnormalized quantities and not volts.
We now consider a case where the tip turns out to be unnecessary, and that brings us to the trick. We wanted to write a program to demonstrate the problem with using floating point numbers when implementing a moving average filter recursively by implementing such a filter on an ADuC841, which has a 12-bit ADC. To calculate the voltage corresponding to the 12-bit value returned by the ADC, one must multiply the 12-bit value (when considered as an integer) by the voltage of the ADuC841’s reference voltage, 2.5 V, and divide by 212. In order to demonstrate the problem with using floats, we stored the measured voltages and used them (and not the 12-bit values returned by the ADC) in our calculations. We confidently expected to find that the output of the filter would be the expected output and a small “dc drift” because of the buildup of the rounding error. We did not find such a drift.
The reason is actually fairly clear and has to do with the way floating point operations are implemented. In the C we used, a floating point number is stored as a single sign-bit, (what is effectively) a 24-bit mantissa (where the mantissa is actually stored in 23 bits, and the first bit of the mantissa, the “unstored bit,” is always one), and an 8-bit exponent [3]. Multiplication by 2.5 is the same as adding twice a number to half the number, and for binary numbers, this is the same as shifting the number left by one and adding it to the number shifted right by one. Multiplying a floating point number by 2.5 lengthens its mantissa by two bits. As the number read from the ADC uses not more than 12 out of the 24 bits of the mantissa that it is “entitled to,” this is not a problem. Similarly, dividing by ${4},{096} = {2}^{12}$ is simply changing the exponent (reducing it by 12), and again, no precision is lost. In total, at most 14 bits of the mantissa are in use, so even after adding 32 such numbers to one another (as we do in our moving average filter), no more than 19 of the 24 bits are in use. Adding the new measurement does not lead to any imprecision, and the subtractions also come off without a hitch.
If you want to see the effects of rounding error, divide by 3,000 (and multiply by it when necessary) instead of by 4,096. Then you should see a very slow drift in the “dc” value of the signal.
We wrote a program that implements a moving average filter that averages the last 32 samples using C floats, and we examined the output of the filter when its input was a sine wave “riding on” 1 V (dc). The formula used to convert samples to floating point was
sample_val = (ADCDATAH*256
+ ADCDATAL)* 2.5
/ CONVERSION_FACTOR;
where ADCDATAL contains the lower eight bits of the ADC reading, and ADCDATAH contains the upper four bits of the ADC reading (and in the program we wrote, its four most significant bits are always zero). At compile time, the user could choose to make CONVERSION_FACTOR either 4,096.0 (the correct conversion factor) or 3,000.0 (as though the correct conversion factor was 3,000.0) by making use of the following preprocessor commands and either leaving in or commenting out the first line (that #defines EXACT).
#define EXACT
#ifdef EXACT
#define CONVERSION_FACTOR 4096.0
#else
#define CONVERSION_FACTOR 3000.0
#endif
When CONVERSION_FACTOR was set to 4,096.0, we expected to see the dc offset remain constant. When CONVERSION_FACTOR was set to 3,000.0, we expected to see a very slow drift because of the information in the least significant bits that was ignored.
We measured the output of the system using the PicoScope 2204A. When setting CONVERSION_FACTOR to 4,096.0 and inputting a 10-Hz sine wave for a period of just over 2 h, we found a change in the dc value of about 1.0 mV (and from time to time, the dc value returned to the value from which it started). When using a CONVERSION_FACTOR of 3,000.0, we found that after a bit more than 1 h and 40 min, the dc value had increased by almost 17 mV.
In this column, we provide a description of a problem that moving average filters can suffer from, and we explain why that problem should not actually rear its head if the values read in from that A/D are stored as integers but that we expect trouble from floating point numbers. The tip, which is not new, is to use ints rather than floats. We then explained why under certain conditions, one can store the numbers as floats and not suffer any ill effects. In particular, when using the ADuC841, the most natural formula for conversion from the integer value of the A/D reading to the floating point value of the measured voltage leads to no loss of accuracy, to no rounding errors, and even adding thirty-two such measurements does not bring us to the point where there will be any rounding errors. To receive a copy of the code used to implement the moving average filter, please send an e-mail to shlomoe@g.jct.ac.il.
Shlomo Engelberg (shlomoe@g.jct.ac.il) received his bachelor’s and master’s degrees in engineering from The Cooper Union, New York, and his Ph.D. degree in mathematics from New York University’s Courant Institute. He is a member of the Department of Electrical and Electronics Engineering, Jerusalem College of Technology, Jerusalem 9116001, Israel. His research interests include applied mathematics, instrumentation and measurement, signal processing, coding theory, and control theory. He is a Senior Member of IEEE.
[1] S. Engelberg, Digital Signal Processing: An Experimental Approach. London, U.K.: Springer-Verlag, 2008.
[2] S. W. Smith, The Scientist and Engineer’s Guide to Digital Signal Processing, 2nd ed. San Diego, CA, USA: California Technical Publishing, 1999. Accessed: Dec. 18, 2022. [Online] . Available: https://www.analog.com/en/education/education-library/scientist_engineers_guide.html
[3] “Floating-point numbers.” ArmKeil. Accessed: Dec. 19, 2022. [Online] . Available: https://www.keil.com/support/man/docs/c51/c51_ap_floatingpt.htm
Digital Object Identifier 10.1109/MSP.2023.3294721