Is there an algorithm for second order filtering of data where the filter frequency is not known a priori?

Question

It is well known that a moving average algorithm done in the time domain is equivalent to a filter with frequency response $\mathrm{sinc}(\omega\tau)$ where $\tau$ is the averaging time. (see this related answer)

This has the following beneficial property: you are streaming a time series of data ${x_n}$, and the average at any point ($a_n$) is just: $$ a_n = a_{n-1} \frac{n-1}{n} + \frac{x_n}{n} .$$

Thus you may apply the above recursive algorithm for an arbitrary amount of time ($\tau$), and when you stop, the value you have is filtered by $\mathrm{sinc}(\omega\tau)$, and has a correspondingly reduced variance. Now the $\mathrm{sinc}$ function is a first order low pass, modulated by a $\sin$ envelope. So in effect you have done a first order low pass where the characteristic low pass time constant $\tau$ is equal to the length of the data stream, and $\tau$ was not necessarily known before you started.

My question is: is there some analogous procedure which allows for an (approximate) second order low pass where the time constant is not known a priori?

A possibility is to "average the averages" but that requires keeping all the averages in memory. Is there some law preventing such a procedure with small memory requirements?

Answer 1

The formula $$ a_n = a_{n-1} \frac{n-1}{n} + \frac{x_n}{n} .$$ is a way to compute the average of $n$ samples, and update this average whenever a new sample is obtained, without storing all of the samples individually.

Ignoring the part about online computation, I think your question is equivalent to "is there a way to estimate the mean of N identically distributed random variables that has lower standard deviation than the arithmetic average?" If so, the answer is: no, there isn't.

I would also worry about the numerical behavior of this formula for large $n$.

Answer 2

You can "average the averages" the same way as you average your input signal. It can be done by the same recursive procedure without storing "all the averages". The only thing you need to do is store two numbers instead of one.

Let $x_n$ be the data to be averaged and let $y_n$ be the output of the first averaging procedure:

$$y_n=\alpha y_{n-1}+(1-\alpha)x_n,\quad 0<\alpha <1$$

Applying the same type of recursion again (just with a possibly different time constant) results in the final output $z_n$:

$$z_n=\beta z_{n-1}+(1-\beta)y_n,\quad 0<\beta <1$$

You can also write the total procedure as a single second order recursion (eliminating $y_n$):

$$z_n=(\alpha+\beta)z_{n-1}+\alpha\beta z_{n-2}+(1-\alpha)(1-\beta)x_n$$

So you have second order recursive filter which only needs to store two past output values. If you want a second-order system, this is the minimum storage possible.

Answer 3

This answer would not have been possible without Matt L.'s answer, as well as some out of band communication with nibot.

Let's see one way to derive the formula for computing the average that is given in the question.

Starting from a set of numbers $\{x_n\}$, we have the definition of the average up to the $n$th sample: $$ a_n = \frac{1}{n}\sum_{j=1}^n x_j = \frac{1}{n}s_n,$$ and $s_n$ is the sum of all samples up to $n$. Now, $s_n$ can be defined recursively: $$ s_n=s_{n-1}+x_n,$$ and given that $na_n=s_n$, we have: $$ a_n=\frac{n-1}{n}a_{n-1}+\frac{x_n}{n}.$$ And we have the averaging formula from the question.

Now we want to basically perform this averaging operation again on the $a_n$ samples. So we just repeat the same formula, but now for the averages of $a_n$.

$$d_n = \frac{n-1}{n}d_{n-1}+\frac{a_n}{n},$$ $$d_n = \frac{n-1}{n}d_{n-1}+\frac{n-1}{n^2}a_{n-1}+\frac{x_n}{n^2}, $$

but we can replace $a_{n-1}$ in terms of $d_{n-1}$ and $d_{n-2}$.

$$d_n = \frac{n-1}{n}d_{n-1}+\frac{n-1}{n^2}\left[(n-1)d_{n-1}-(n-2)d_{n-2}\right]+\frac{x_n}{n^2}, $$

and finally after simplification

$$d_n =\frac{(2n-1)(n-1)}{n^2}d_{n-1}-\frac{(n-1)(n-2)}{n^2}d_{n-2}+\frac{x_n}{n^2}.$$

Now this set of numbers is equivalent to "averaging the averages" and only requires two stored values!

Below I plot a signal which is random noise where the RMS is 20 times the mean value. I also show the first and second order averages.

As one can see, the second order average takes longer to approach the true mean value, but it has smaller fluctuations relative to the mean.

The fluctuations get smaller as more and more samples are recorded, so it has the added benefit that the time scale of the effective low-pass filter is always increasing. If this were a simple low pass filter with a fixed pole frequency, then at some point we would be throwing away information from very old samples. This filter uses information from all samples, regardless of how old they are.

Finally, I think this recipe can be repeated and the average can be done to any order.

Answer 4

Yes, you can do a second-order low-pass filter without using lots of memory. The key is to use the fact that convolution is a linear operation. You want to do the following: $$ y(t) = (x(t)*f_1(t))*f_2(t) $$ where $f_1(t)$ and $f_2(t)$ are your two moving average filters of unknown a priori width. If we use the associative property of linearity we can do the following: $$ y(t) = (x(t)*f_1(t))*f_2(t) = x(t)*(f_1(t)*f_2(t)) $$ You create a new filter by convolving the two averaging filters, and then using that composite filter to filter your data.