Recursive Implementation of the Gaussian Filter

Hello,

I'm trying to implement the article "Recursive Implementation of the Gaussian Filter".

This article suggest an IIR Filter as an approximation of the Gaussian Blur. This is the suggested method:

Namely it is an order 4 IIR Filter.

I tried to reproduce the results for q = 5 as given in the article (See "Example").

Here is my 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);

I get the correct number for all coefficients, yet the result is:

What am I missing?

Has anyone managed to make it work?

Thank You.

Answer 1

[Royi Avital]

he 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);