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

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?

• Also perhaps relevant: dsp.stackexchange.com/a/892/64 – nibot Mar 26 '14 at 21:28
• @MattL. : The formula is correct, and does depend on the time index $n$. It's a way to calculate the mean of a signal that you don't know the length of when you start sampling. – Peter K. Mar 27 '14 at 13:43
• @PeterK. OK, but then we're talking about a time-varying system, and we cannot think in terms of frequency responses anymore. BTW, also with an LTI recursion (i.e. with fixed multipliers) the length of the input signal doesn't really matter. – Matt L. Mar 27 '14 at 16:02

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

• Unless the OP has a different interpretation of what is there some analogous procedure which allows for an (approximate) second order low pass where the time constant is not known a priori? means, I believe this is correct (+1). – Peter K. Mar 27 '14 at 13:49
• I think this answer is correct if one assumes that the measurement samples are independent, or equivalently, that the noise is white. However, the problem I was trying to solve had a noise power spectral density rising like $\omega^2$, and thus the measurements are correlated. – achoo5000 Oct 1 '14 at 2:23

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.

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.

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.

• "Associative property of convolution" I suppose. – Matt L. Mar 26 '14 at 21:06
• @MattL. It is my understanding that linearity implies associativity. Is this not the case? – Jim Clay Mar 26 '14 at 21:50
• When I read your answer, I was sure that you actually meant to say "associative property of convolution", because it is always some type of binary operation that is either associative or not, and you used the associativity of convolution. I think we cannot talk about the 'associative property of linearity', because 'linearity' is no binary operation. I didn't mean to be nit-picky, but maybe I was ... But anyway, your question is interesting (as to the relation between linearity and associativity) and I must admit that I have no satisfactory answer to it. – Matt L. Mar 27 '14 at 11:33
• Linearity is for maps between vector spaces (or more abstractly, modules). The answer turns out to be yes, and you can find details in the general case in a book on category theory, probably like The Joy of Cats by Adamek et al. – Batman Mar 27 '14 at 19:29