# Optimal n-th order IIR approximation to a moving average filter

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.

• For a first order filter, the answers to this question should be sufficient. For higher order filters you could try Prony's method for approximating a given impulse response. – Matt L. Aug 8 at 18:47
• Not sure why the downvote ... (?) – Matt L. Aug 8 at 19:02
• Thanks Matt! I have seen the question you mentioned, but first-order is very restrictive. I was hoping to approximate a 600-tap moving average with order of magnitude 10 taps IIR. – burnpanck Aug 8 at 20:01

This is not a full answer, but I want to bring another point of view, that the time lag of the moving average filter matters.

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

So when you create your IIR approximation, it may be beneficial to allow the impulse response rise a bit, like the Gaussian function in Fig. 1, already before the moving average proper switches on.

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])$$

• Well, I'm worried about storage space for the taps, not execution cycles. As you mention, the efficient boxcar implementation solves the latter not the former. One should add that these need to be implemented in integer/fixed-point arithmetic, or be prepared for brownian walk-away due to round-off error. – burnpanck Aug 8 at 21:32
• Along the lines of what Hilmar said, if you use exponential smoothing with $\lambda = 2/(n+1)$ and $y_{t} = (1- \lambda) y_{t-1} + \lambda x_{t}$ you get a very similar response. but " similar" is defined in time domain as equal average age of input observations. I don't know what it implies in frequency domain. – mark leeds Aug 8 at 21:33
• and see "smoothing, forecasting and prediction of discrete time series" by robert brown for details on the derivation of the smoothing constant $\lambda$. – mark leeds Aug 8 at 21:34

(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:

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

• Yes, I'm trying to approximate a boxcar filter, but with the constraint of not having access to delays large enough for the time window. – burnpanck Aug 8 at 21:35
• Those goals tend to act in opposition. This seems like work. You could try a Pade approximation. Without knowing why you need an effective 600 tap boxcar, a trade off is Ill posed – Stanley Pawlukiewicz Aug 8 at 23:18

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?

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.

• What $L_2$ error values did you get from these? – Olli Niemitalo Aug 13 at 13:13
• @OlliNiemitalo: Between 3e-4 (1st order) and 7e-5 (5th order) when defined as $\sum_n|h_{FIR}[n]-h_{IIR}[n]|^2$. – Matt L. Aug 13 at 14:09