0
$\begingroup$

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:

enter image description here

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:

enter image description here

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – robert bristow-johnson Mar 15 '15 at 1:14
1
$\begingroup$

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);
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.