I found this elegant A/D averaging algorithm that works very well. Simple integer math, based 2.
But I don't understand it. The author says it's "exponential". Is this right?
I'm trying to do both a fast and a slow average. I get confused if I should cascade two (fast average result goes into slow average) or use two separate filters.
If I change the divide by n (shift) it somehow always gives the right answer, regardless.
uint16_t adcAvg, adcTot;
adcTot -= adcAvg;
adcTot += readAdc(adcChannel);
adcAvg = adcTot >> 4;
Author and source:
"I normally use a different version of the exponential filter. Similar concept, just optimized for smaller micros using only 32bits of RAM.
...
then you can change how severe the filter is by the shift (or divide if you wanted to do a divide there instead).
Has its limitations, e.g. cant have more than 64 10 bit samples without going to bigger variables, but 64 samples is heaps. I tend to use 8 or 16 (>> 3 or >> 4) but overall its a simple solution to the problem that i use every time I need to use the ADC."
That is a standard Infinite Impulse Response (IIR) first order digital filter -- there are entire textbooks devoted to digital filters. From the link in the Hackaday article, it is based on this formula:
Exponential Filter
The last filter is a recursive filter. A recursive filter is just one that calculates a new, smoothed value (yn) by using the last smoothed value (yn – 1) and a new measurement (xn):
yn = w × xn + (1 – w) × yn – 1
The code you posted above uses this formula, with w = 1/16 (>>4 = divide by 16), but cleverly rearranged to preserve accuracy.
The "speed" of the average is determined by the shift factor. Here's an explanation of how it works (there's quite a bit of maths involved, but you can just use the results, without understanding the entire derivation):
Exponential Moving Average filter
The most important things to take away from all that is the difference equation (the very first equation), the shape of the impulse and step responses, the frequency response (Bode plot), and the formula for the cutoff frequency.
You can compare it to the Simple Moving Average filter if that's more familiar to you (just the sum of the last N samples divided by N).
This is the implementation of the EMA I use in my projects: https://github.com/tttapa/Control-Surface/blob/master/src/Helpers/EMA.hpp
Pieter
Now I see the subtle difference between the two different moving average filters. Thanks for your help.
SMA: Subtract the oldest sample in the array, like Arduino Tutorial on Smoothing
EMA: Subtract the previous average
I am averaging a level sensor's water tank slosh and find the exponential does better. It responds faster when tank fill/empty occurs.
jremington:
...but cleverly rearranged to preserve accuracy.
Except for the the fact that the value is truncated rather than rounded making the least-significant bit suspect.
Yes, one should add 8 before dividing by 16.
prairiemystic:
I am averaging a level sensor's water tank slosh and find the exponential does better. It responds faster when tank fill/empty occurs.
Exactly. That's where the "exponential" part is important: when the average value changes, the SMA evolves to this new average linearly, whereas the EMA evolves exponentially.
SMA with length 9, EMA with the same -3 dB frequency, sampled at 360 Hz
jremington:
Yes, one should add 8 before dividing by 16.
Thoughts on this implementation?
template <uint8_t K, class int_t>
class EMA {
public:
int_t filter(int_t input) {
filtered += input;
int_t output = (filtered + fixedPointAHalf) >> K;
filtered -= output;
return output;
}
private:
int_t filtered = 0;
constexpr static int_t fixedPointAHalf = 1 << (K - 1);
};
jremington:
Looks good to me!
Thanks!
I was interested to see what difference the rounding makes, so I did a little experiment:
As you can see, the orange curve is often one unit too low, because of the floor division.
Here's the python code of the experiment:
import numpy as np
from matplotlib import pyplot as plt
class EMA_shift_no_round:
def __init__(self, shiftfac: int):
self.shiftfac = shiftfac
self.z = 0
def __call__(self, x: int):
self.z += x
shifted = self.z >> self.shiftfac
self.z -= shifted
return shifted
class EMA_shift_round:
def __init__(self, shiftfac: int):
self.shiftfac = shiftfac
self.fixedPointOneHalf = 1 << (shiftfac - 1)
self.z = 0
def __call__(self, x: int):
self.z += x
shifted = (self.z + self.fixedPointOneHalf) >> self.shiftfac
self.z -= shifted
return shifted
class EMA_float:
def __init__(self, alpha: float):
self.alpha = alpha
self.y = 0
def __call__(self, x: int):
self.y = self.y * (1 - self.alpha) + x * self.alpha
return self.y
K = 7
EMA_no_round = EMA_shift_no_round(K)
EMA_round = EMA_shift_round(K)
EMA_float = EMA_float(2 ** -K)
t = np.linspace(0, 31+64+64, (32+64+64)/2)
c = 16 * np.cos(2*np.pi * t / 64) + 16
z = np.zeros((64,))
m = 32 * np.ones((256,))
x = np.concatenate((c, z, m))
y_no_round = np.array(list(map(lambda x: EMA_no_round(int(round(x))), x)))
y_round = np.array(list(map(lambda x: EMA_round(int(round(x))), x)))
y_float = np.array(list(map(lambda x: EMA_float(int(round(x))), x)))
plt.plot(x, label="in")
plt.plot(y_no_round, label="no round")
plt.plot(y_round, label="round")
plt.plot(y_float, label="float")
plt.legend()
plt.show()
Thanks for posting that clear demonstration of the how the rounding correction works.
It should also be interesting and informative to post a diagram showing the effect of other choices of divisor, e.g. 8 or 32.
PieterP:
I was interested to see what difference the rounding makes, so I did a little experiment:
Trying to control (e.g. PID) is when that tiny error can become a serious problem.
jremington:
It should also be interesting and informative to post a diagram showing the effect of other choices of divisor, e.g. 8 or 32.
import numpy as np
from matplotlib import pyplot as plt
class EMA_float:
def __init__(self, alpha: float):
self.alpha = alpha
self.y = 0
def __call__(self, x: int):
self.y = self.y * (1 - self.alpha) + x * self.alpha
return self.y
EMA_4 = EMA_float(2 ** -4)
EMA_5 = EMA_float(2 ** -5)
EMA_6 = EMA_float(2 ** -6)
EMA_7 = EMA_float(2 ** -7)
t = np.linspace(0, 31+64+64, (32+64+64)/2)
c = 16 * np.cos(2*np.pi * t / 64) + 16
z = np.zeros((64,))
m = 32 * np.ones((128+64,))
x = np.concatenate((c, z, m))
plt.plot(x, 'k:', label="in", linewidth=1)
plt.plot(list(map(lambda x: EMA_4(int(round(x))), x)), label="K = 4")
plt.plot(list(map(lambda x: EMA_5(int(round(x))), x)), label="K = 5")
plt.plot(list(map(lambda x: EMA_6(int(round(x))), x)), label="K = 6")
plt.plot(list(map(lambda x: EMA_7(int(round(x))), x)), label="K = 7")
plt.legend()
plt.tight_layout()
plt.show()
Here's another visualization, where you can clearly see the phase shift and the effect it has on different frequencies.
Link
You can also hear the result of the filter applied to audio.
