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.
-
1$\begingroup$ 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. $\endgroup$– Matt L.Aug 8, 2018 at 18:47
-
1$\begingroup$ Not sure why the downvote ... (?) $\endgroup$– Matt L.Aug 8, 2018 at 19:02
-
$\begingroup$ 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. $\endgroup$– burnpanckAug 8, 2018 at 20:01
-
$\begingroup$ Matt L's answer below appears to be the best answer given your specific question; you should select it as such unless you disagree. As far as implementation I would recommend seriously consider proceeding with a decimating CIC structure; after such a long moving average the higher sampling rate is not necessary and this would result in by far the simplest solution with performance matching the moving average filter you describe. See this post for further details on that dsp.stackexchange.com/questions/31548/… $\endgroup$– Dan BoschenAug 22, 2018 at 15:26
-
$\begingroup$ @DanBoschen: Indeed, there have been many good suggestions here how to solve this apparently surprisingly difficult question. I guess you are right though that a decimated CIC is probably the best given the additional constraint of small number of taps (which I didn't explicitly specify in the question). $\endgroup$– burnpanckAug 23, 2018 at 7:09
6 Answers
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.
-
$\begingroup$ What $L_2$ error values did you get from these? $\endgroup$ Aug 13, 2018 at 13:13
-
1$\begingroup$ @OlliNiemitalo: Between 3e-4 (1st order) and 7e-5 (5th order) when defined as $\sum_n|h_{FIR}[n]-h_{IIR}[n]|^2$. $\endgroup$– Matt L.Aug 13, 2018 at 14:09
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.
(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.
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])$$
-
1$\begingroup$ 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. $\endgroup$ Aug 8, 2018 at 21:32
-
1$\begingroup$ 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. $\endgroup$ Aug 8, 2018 at 21:33
-
1$\begingroup$ and see "smoothing, forecasting and prediction of discrete time series" by robert brown for details on the derivation of the smoothing constant $\lambda$. $\endgroup$ Aug 8, 2018 at 21:34
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} $$
-
$\begingroup$ 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. $\endgroup$ Aug 8, 2018 at 21:35
-
$\begingroup$ 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 $\endgroup$– user28715Aug 8, 2018 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?