# Recursive Implementation of the Gaussian Filter (1D & 2D)

I'm trying to implement an IIR form to approximate the Gaussian Blur Filter.
I'm working with the article "Recursive Implementation of the Gaussian Filter" by Ian T. Young and Lucas J. van Vliet.

They suggest a form and way to calculate the coefficients as given by:

I'm trying to reproduce their example for $q = 5.0$.

There is the MATLAB code I wrote:

qFactor = 5;

b0Coeff = 1.57825 + (2.44413 * qFactor) + (1.4281 * qFactor * qFactor) + (0.422205 * qFactor * qFactor * qFactor);
b1Coeff = (2.44413 * qFactor) + (2.85619 * qFactor * qFactor) + (1.26661 * qFactor * qFactor * qFactor);
b2Coeff = (-1.4281 * qFactor * qFactor) + (-1.26661 * qFactor * qFactor * qFactor);
b3Coeff = 0.422205 * qFactor * qFactor * qFactor;

normalizationCoeff = 1 - ((b1Coeff + b2Coeff + b3Coeff) / b0Coeff);

vDenCoeff = [b0Coeff, b1Coeff, b2Coeff, b3Coeff] / b0Coeff;

vXSignal = zeros(61, 1);
vXSignal(31) = 10;

vYSignal = filter(normalizationCoeff, vDenCoeff, vXSignal);
vYSignal = filter(normalizationCoeff, vDenCoeff, vYSignal(end:-1:1));

figure();
plot(vYSignal);


Now, the result I get is this:

Namely, I get a filter which isn't stable.
Yet it seems I get the same coefficients as they get in their example.

What am I missing?
Has anyone ever implemented this method?

Thank You.

• i've never had the occasion of implementing a gaussain filter other than as an FIR on a big table of data. normally, the only time i am implementing an IIR filter is for a real-time application and i have never done a gaussian filter in that context. but i have read about them (at least in an analog filter textbook). i dunno how to do an all-pole gaussian, but i imagine that there is a way to approximate it for analog filter. then it depends on if the time-domain or frequency-domain performance is more important to you whether you choose impulse-invariant or bilinear transform to convert it. Commented Mar 15, 2015 at 1:14
• Related to dsp.stackexchange.com/questions/49583.
– Royi
Commented Jul 23, 2021 at 19:19
• Related to dsp.stackexchange.com/questions/50576.
– Royi
Commented Jul 23, 2021 at 19:20
• Here is Lucas van Vliet’s original C code, written for that paper. It has been modified over the years to fit into the DIPlib library, and then to work in a C++ compiler, but the core of the code is mostly the same. github.com/DIPlib/diplib/blob/master/src/linear/gaussiir.cpp Commented Jan 12, 2022 at 15:20

The answer was simple, the article uses the coefficients value on one hand where the MATLAB implementation on the other.
Namely, a minus sign should be added.

Here's the correct code:

qFactor = 5;

b0Coeff = 1.57825 + (2.44413 * qFactor) + (1.4281 * qFactor * qFactor) + (0.422205 * qFactor * qFactor * qFactor);
b1Coeff = (2.44413 * qFactor) + (2.85619 * qFactor * qFactor) + (1.26661 * qFactor * qFactor * qFactor);
b2Coeff = (-1.4281 * qFactor * qFactor) + (-1.26661 * qFactor * qFactor * qFactor);
b3Coeff = 0.422205 * qFactor * qFactor * qFactor;

normalizationCoeff = 1 - ((b1Coeff + b2Coeff + b3Coeff) / b0Coeff);

vDenCoeff = [b0Coeff, -b1Coeff, -b2Coeff, -b3Coeff] / b0Coeff;

vXSignal = zeros(61, 1);
vXSignal(31) = 10;

vYSignal = filter(normalizationCoeff, vDenCoeff, vXSignal);
vYSignal = filter(normalizationCoeff, vDenCoeff, vYSignal(end:-1:1));

figure();
plot(vYSignal);

• It is really nice of you to post a question and get back to share you self developed answer.
– Mark
Commented Dec 28, 2020 at 19:42