Which IIR filters approximate a Gaussian filter?

Question

So it recently dawned on me that Bessel filters, despite being listed along with the other common types, are really an oddball that belongs in a different "class", and I'm trying to learn more about it.

The rectangular magnitude response represents the ideal frequency domain response, for the transition band is zero and the stopband has infinite attenuation. The Gaussian magnitude response, on the other hand, represents the ideal time-domain response, in that no overshoots occur in the impulse response and the step response. Many of the responses attained in practice are approximations to these ideal ones ^source

So a brickwall filter is convolution with a sinc function, and has these frequency domain properties:

• Flat passband
• Zero stopband
• Infinite roll-off rate/no transition band

It's non-causal and unrealizable because of the infinite tails in both directions. It is approximated by these IIR filters, with the approximation improving as order increases:

• Butterworth (maximally flat passband)
• Chebyshev (maximum roll-off rate with stopband or passband ripple)
• Elliptic (maximum roll-off rate with stopband and passband ripple)
• Legendre (maximum roll-off rate with monotonic passband)

---

The Gaussian filter is convolution with a Gaussian function, and has these time domain properties:

• Zero overshoot
• Minimal rise and fall time
• Minimal group delay

It's unrealizable for the same reasons as the sinc function, and can be approximated by these IIR filters, more closely as order increases:

• Bessel (maximally flat group delay) ^{according to 1 and 2}

Here's Bessel filters of increasing order along with a Gaussian dashed line I picked merely because it seemed to fit the trend ($e^{-{1 \over 2}(\pi \omega)^2}$):

So my questions are:

Is everything right so far? If so, are there other IIR filters that approximate the Gaussian? What are they optimized for? Maybe one that minimizes overshoot?

If you search for "IIR Gaussian" you can find a few things (Deriche? van Vliet?), but I don't know if they're really the same as a Bessel or if they optimize for some other property, etc.

Answer 1

The Deriche and van Vliet filters are heuristics. In both cases they choose the locations of poles and zeros to minimize either the RMS difference or the maximum difference of the filter's impulse response from a Gaussian.

Both filters are causal-anti-causal pairs. So I think they have no phase error or group delay, but you need to be able to run them backwards on the data as well as forwards. This makes them popular in image processing, but perhaps restricts their applicability elsewhere.

That they are heuristics is evidenced by the wealth of papers that tweak them. For example a google search (while I was looking for the link to the Deriche paper) turned up this one that tries to solve the problem that the Deriche derivative-of-gaussian filter doesn't have an exactly 0 DC response. There are also some interesting issues about initializing the boundary conditions correctly.

I have found the following overview a good resource: Dave Hale, Recursive Gaussian filters, Colorado School of Mines Center for Wave Phenomena report CWP-546.

Answer 2

I think you gave a nice summary of existing analytic solutions for discrete-time IIR filters. But I would also add Bessel filters to the list of filters approximating ideal frequency-selective filter characteristics. Its magnitude response does not show as sharp a transition as the other filter types of the same order, but this is the price you pay for an almost linear phase in the passband. So the Bessel filter is a compromise between a frequency-selective magnitude response and a good phase response.

For approximating a Gaussian filter with IIR filters, I do not know of any analytic solutions, apart from the Bessel filter you mentioned. But note that the Bessel filter was not meant to approximate a Gaussian filter, so I'm not sure how good it really is in approximating such a filter. If you really want an IIR filter for this purpose, I would suggest you go for numerical approximation of the Gaussian filter. There are several options how to do this.

You could try to approximate the Gaussian filter in the frequency domain. The problem is that you have to make some decision concerning the desired phase response. A pure magnitude approximation with minimum phase response will very likely result in very poor time-domain properties. If you specify a linear desired phase, then you get a complex approximation problem (because you approximate the complex frequency response with magnitude and phase). Even though such an approximation problem can be quite tough to solve, there exist methods in the literature.

A simpler and probably better approach is to approximate the Gaussian filter in the time domain. Prony's method would be a good starting point.

Please note that these are just my thoughts on the topic. I have not tried to design an IIR Gaussian filter myself. I would actually go for an FIR implementation unless there are very good reasons against it.

EDIT: just a few more remarks concerning the question whether a Bessel filter approximates a Gaussian or not. I do not know of any meaningful error criterion which the Bessel filter minimizes in approximating a Gaussian filter. I would be glad to learn about it though. People may claim that the impulse response of a Bessel filter looks similar to Gaussian, or that its frequency response resembles a Gaussian, but I have not yet seen any proof that Bessel filters approximate a Gaussian in any sense, and that the approximation error goes to zero as the filter order increases. I do not deny that it is more similar to a Gaussian than the other standard filters (Butterworth, Chebyshev, etc.), but this is not important for the question.

See below four plots of impulse responses of Bessel filters (orders 5, 10, 15, 20), designed in Octave (function besself). As you can see, the ringing in the tail does not decrease with increasing filter order, and I do not see how these filters approximate a Gaussian, and if so, according to which optimality criterion. However, if anybody can enlighten me about this, I would be more than happy.

Answer 3

I might be a tad late to this, but I'll only reply to the part about the "similarities" between the (analog) Bessel and Gaussian.

They are not the same. The Bessel filter is meant to approximate an ideal delay:

$$B(s)=\mathrm{e}^{-s\tau}\tag{1}$$

while the Gaussian filter tries to approximate a Gaussian bell:

$$|G(j\omega)|^2=\mathrm{e}^{-\alpha\omega^2}\tag{2} \\ \alpha=\dfrac{\ln{2}}{\omega_p}$$

In making the filters, the Bessel filter is built starting from the premise that what's needed is a flat group delay, so the derivation starts from there and reaches the Bessel polynomials (see this on ee.se, for example).

The Gaussian, on the other hand, starts from a Taylor expansion from the generating function:

$$\mathrm{e}^{\alpha\omega^2}\approx 1+\sum_{k=1}^N{\dfrac{\alpha^k}{k!}\omega^{2k}}$$

If the errors seem to get smaller and smaller with increasing order, that's because the approximations gain higher and higher terms, so the errors are really the residuals from very high polynomials. But their characteristics differ quite a bit. And for high orders, when searching for differences, they need to be sought at a logarithmic scale.

Here are comparisons for two 8th order filters, where the Bessel filter has a frequency scaling applied to match the -3 dB frequency of the Gaussian. Normally, that messes up the group delay (possibly related on ee.se), but here is just for comparison:

The Bessel filter has a flat(ter) group delay, while the Gaussian has slower phase roll-off and, thus, better temporal delay. Between the two, the Bessel filter has better higher attenuation. Also, the Bessel filter will always have a little overshoot, even the Gaussian for small orders, but the no overshoot criterion applies to the Gaussian, not the Bessel.

In short, both are for temporal filtering, but the Bessel is for flat group delay, while the Gaussian is for no overshoot -- not minimal delay, that can be achieved with a very poor filter, but minimal delay with the fastest rising time and no overshoot.

Answer 4

_{^{I hesitated whether to reuse the already existent answer but, this deals with the OP, directly, whereas the other one deals with the differences between the two filters.}}

I see that the proposed solutions involve some hefty algorithms and whatnot. The way I see it, even those algorithms will not provide an exact solution for low orders but, then, low orders can't really be "exact". So, I was wondering, why not use the BLT on the already existent analog Gaussian (approximation)? After all, what's needed is the magnitude since, Gaussian filters are meant to have an $\text{e}^{-t^2}$ type of response. And, since the Fourier transform is the same (or, should be, given that it's an approximation), it follows that the bilinear transform can be applied. Plus, dirty deeds, done dirt cheap.

I tested this and, as long as the bandwidth is at least half the Nyquist or less, the results are quite impressive. This is related to lowpass Gaussian, whose derivation is described in the other answer. So, for a 2nd order, the lowpass prototype is:

$$\begin{align} H_2(s)&=\dfrac{2.04028}{s^2+2.63931s+2.04028} \tag{1} \\ f_s&=2 \\ G_2(z)&=H\left(2f_s\dfrac{1-z^{-1}}{1+z^{-1}}\right) \\ {}&=\dfrac{0.0713446+0.142689z^{-1}+0.0713446z^{-2}}{1-0.643056z^{-1}+0.261667z^{-2}} \tag{2} \\ h_{_2}(t)&=1.86628\cdot\text{e}^{-0.659828t}\sin(0.273309t) \tag{3} \\ g_{_2}[n]&= \begin{cases} n=0\quad 0.0713446 \\ n>0\quad 2.04124\cdot 0.511534^n\cos(0.303572n-1.66958) \end{cases} \tag{4} \end{align}$$

The impulse responses may not be exactly "on point" but, the step response is a different matter (I've used $(t-1)$ for the analog impulse response and $(t-1/2)$ for the step response to align the "continuous" response with the points, whose index starts at 1):

And this for $f_s=2$. For 3 Hz it looks like this:

And this works even better for higher orders (5th order, $f_s=2$):

$$H_5(s)=\dfrac{27.3859}{(s+1.77368)(s^2+2.89154s+4.52647)(s^2+3.40318s+3.41107)} \tag{5}$$

(I could type the IIR but, really, it's nothing special, just your average bilinear transform)

You can use $f_s<2$ but, at some point the poles will get too affected by the conversion, and the result will start gathering oscillations/overshoots. However, this is to be expected.

---

Here's a comparison of the magnitudes for the 2nd order, where the analog is blue with H(%i*x*f_s), and the IIR is G(exp(%i*x)):