Optimal $ n $ -th Order IIR /AR Approximation of a Moving Average Filter

Question

I would like to approximate a moving average filter with an IIR filter of much lower order than the tap-length of the moving average filter. Optimality shall refer to the $L_2$ norm of the impulse response (up to changes in the total gain). Ideally, I'd also like to constrain the impulse response to be monotonic (i.e. no ripple), but I fear that this requirement along will already fully constrain the filter. If that is the case, I might be willing to take some ripple in turn for a sharper roll-off.

Answer 1

I can show you some low order IIR approximations to an FIR moving average filter. In the figure below you see $3$ (infinite) impulse responses that approximate a moving average of length $N=600$. The filter orders are $1$, $2$, and $5$, respectively, and they all approximate the desired response in a least squares sense. I used the equation error method to design those filters. Only the filter with order $1$ is monotonic, but it's also the worst filter in terms of the $L_2$ error.

It's hard to find good solutions for higher order IIR filters because the problem is highly non-linear and you tend to get stuck in a local minimum. I haven't been able to find any useful solutions for filter orders higher than $5$.

I should add that I just used the equation error method and Prony's method. With Prony's method (not shown in the plot), if you increase the order of the numerator $M$ (and leave the denominator order small), you can make the first $M+1$ samples of the designed impulse response equal the first $M+1$ samples of the desired impulse response, after which the response decreases exponentially. But these filters have a higher $L_2$ approximation error than the ones shown in the figure.

Answer 2

This is not a full answer. It explores some basic approximations.

For a moving average as long as 600 samples it is informative to look at impulse responses of continuous-time filters as approximations of those of the desired discrete-time filters. Let's have as the ideal desired impulse response a boxcar function normalized to a width of 1 spanning $x = 0\ldots1$ and a height of 1. The square of the $L_2$ norm of the difference between the actual and the ideal impulse response can be calculated identically from the impulse response or the frequency response by squaring the error and integrating over the domain.

The impulse response of a 1st order filter is a right-sided exponential function with three parameters: start, decay rate, and height. The limiting impulse response of a cascade of an infinite number of such filters is a Gaussian function, by the central limit theorem. A Gaussian function has also three parameters: the center, width, and height. Guessing the start and center parameters and then minimizing numerically the square of the time domain $L_2$ norm gives impulse responses:

$$\left\{\begin{array}{ll}1.4807e^{- 1.1862x}&\text{if }x \ge 0,\\ 0 &\text{otherwise}\end{array}\right.$$

for the exponential function with start at $x = 0$ and a square of the $L_2$ norm of $0.16203$ and

$$0.44504e^{3.9200x - 3.9200x^2}$$

for the Gaussian function with center at $x = \frac{1}{2}$ and a better square of the $L_2$ norm of $0.10990$.

Figure 1. 1st order (red) and the limiting shape of the infinite order filter (blue) impulse response and the boxcar function (black dotted).

In a discrete-time implementation, the left-side tail of the Gaussian impulse response could be removed without computational penalty, removing its contribution to the total error.

Answer 3

(Update: I just realized that the first part of this covering CIC structures is basically what Hilmar has already answered-- I'll leave this up since it offers more graphics and details in case that his helpful to anyone but it is indeed the same answer)

This may not be optimum and although a feedback structure is involved it is strictly not IIR, but want to point out the efficiency of the cascade integrator-comb structures for implementing boxcar moving average filters- and importantly how the results match exactly and not an approximation. So a 600 tap moving average filter can be implemented with simply one subtraction after a 601 sample delay followed by an accumulator.

These structures are quite popular when resampling is involved (often after a moving average the higher sampling rate of the input is no longer required and a lower output sampling rate would be even more efficient).

The equivalency is easily proven from the Z-transform for each of these structures. Specific to the example structures in the graphic:

The Z-transform for the accumulator is:

$$\frac{1}{1-z^{-1}}$$

The Z-transform for the comb (delay and subtract) filter is:

$$1-z^{-4}$$

And the Z-Transform for the moving average filter is:

$$1+z^{-1}+z^{-2}+z^{-3}$$

Combining confirms the equivalence:

$$1+z^{-1}+z^{-2}+z^{-3} = \frac{1}{1-z^{-1}}(1-z^{-4})$$ $$(1-z^{-1})(1+z^{-1}+z^{-2}+z^{-3})= (1-z^{-4})$$

Another favorite of mine as a true IIR alternative to a moving average filter is the exponentially weighted moving average filter, where assuming sufficient bit precision, very tight (in other words effectively long delay) filtering can be achieved by using an alpha that approaches 1. This is demonstrated in the following graphics:

In fixed precision structures the order of the two operations shown in the block diagram of the exponentially weighted moving average filter is very important! In most cases when dealing with fixed precision it is best to allow digital filters to grow the signal first and do the scaling down after the accumulation. This also explains the motivation for extended precision accumulators- the accumulator shown needs to be designed with sufficient precision to prevent any overflow for the given $\alpha$ value chosen (the Infinite Summation identity $\sum_{n=0}^\infty \alpha^n =\frac{1}{1-\alpha}$ is useful in this case to determine worst case maximum value; resulting in $gain = \frac{1}{1-0.99}$ or 100 in the example shown, so in this case we would need to add $log2(100) = 7$ bits in the accumulator). Also clear from the summation identity is how the scaling term itself was determined such that the gain of the structure is 1. The order of these operations is no different than the reason to put an amplifier prior to an attenuator in an analog system to get the best noise figure, while paying attention to clipping and maximum input signal, for those familiar with that.

Answer 4

Depends a bit on your application. A moving average filter is a low pass filter and one with many lobes and pretty poor stop band rejection at that. Depending on what specific requirements you have, you may be better off doing a simple butterworth lowpass filter.

If you are worried about execution cycles, there is a very efficient way to implement it, if you use a rectangular windows. Assuming your window length is L you can simply do

$$y[n] = y[n-1]+\frac{1}{L} \cdot (x[n]-x[n-L])$$

Answer 5

Assuming you mean by a moving average is a boxcar filter. The suggestion by @Hilmar will work but has problems with roundoff errors as N gets large.

Instead I suggest you try factoring the boxcar in terms of a cascade of sparse FIR filters, look at:

Mitra, Sanjit K., et al. "General polynomial factorization-based design of sparse periodic linear arrays." IEEE transactions on ultrasonics, ferroelectrics, and frequency control 57.9 (2010): 1952-1966.

The paper is focused on Transmit/Receive Linear Arrays but is equally applicable to cascades of FIR Filters. As an example: $$ (1 + z^{-1})(1 + z^{-2})(1 + z^{-4}) = 1 + z^{-1} + z^{-2} + z^{-3} + z^{-4} + z^{-5}+ z^{-6} + z^{-7} $$

Answer 6

i would think that if your system has no complex poles, only real poles, then you could make the impulse response to be monotonic. the impulse response would be the sum of decaying exponential functions. you can always normalize its area to 1 and that would make it a moving weighted average.

how do you prefer your moving average to be weighted?