What I am trying to do

I am trying to code a Gaussian smoothing filter using the 4th order IIR filter described in Van Vliet's paper "Recursive Gaussian derivative filters". My code works but suffers from numerical instability, so I am looking into decomposing the 4th order filter into a parallel bank of second order filters using the method of Partial Fraction Expansion.

What I don't understand

What I don't understand is how the anti-causal filter can be taken into account, as all of what I read only seems to describe partial fractions on a single causal filter. While the Van Vliet filter is a cascade of a causal and an anti-causal filter. To clarify, I want a fully parallel pipeline, not decompose each of the causal and anti-causal filter separately.

What I understand so far

Van Vliet's filter is a serial interconnection of a causal and an anti causal 4th order filters.

$$ H(z) = H_+(z)\cdot H_-(z) = \prod_{i=1}^4 \frac{d_i - 1}{d_i - z^{-1}} \cdot \prod_{i=1}^4 \frac{d_i - 1}{z - d_i} $$

Where $d_i$ are its anti-causal poles. It has no zeros, only poles.

I started by looking at the causal filter, trying to put it in the form described in Equation 3.39 in the Partial Fraction chapter of the book "Discrete-time signal processing":

Equation 3.41

I did that by dividing by $d_i$ to get:

$$ H_+(z) = \prod_{i=1}^4 1 - d_i^{-1} \prod_{i=1}^4 \frac{1}{1 - d_i^{-1}z^{-1}} $$

Which is in the expected form by noting that $\frac{b_0}{a_0}$ takes the value of the first product. Then we evaluate Equation 4.41 to get the residues $A_k$:

Equation 4.41

And write it as a sum of fractions, where $d_k = d_i^{-1}$ is the causal pole:

$$ H_+(z) = \sum_{i=1}^4 \frac{A_k}{1 - d_k z^{-1}} $$

Since $d_1 = \bar{d_2}$ and $d_3 = \bar{d_4}$, we can combine each two fractions to get our second order section.

The problem

The problem is that I only considered the causal filter above, what about the anti-causal one?

What I tried

I tried to rewrite the anti-causal expression in a form that is as close as possible to Equation 3.39, I did that by also dividing by $d_i$ and factoring a negative sign to flip the denmorator:

$$ H_-(z) = -\prod_{i=1}^4 1 - d_i^{-1} \prod_{i=1}^4 \frac{1}{1 - d_i^{-1}z} $$

The only difference from the causal filter is the negative sign and the fact that z has a power of 1, as expected. I am not sure if Equation 3.41 still holds in this case. But even if I apply the partial fraction for the anti-causal filter, what then? I now have two polynomials that when multiplied gives a very long expression that doesn't seem very useful.

$$ H(z) = \sum_{i=1}^4 \frac{A_{k_+}}{1 - d_k z^{-1}} \cdot \sum_{i=1}^4 \frac{A_{k_-}}{1 - d_k z} $$

I saw some code that does something similar, but I can't make sense of its equations. I just know it seems possible.

  • $\begingroup$ How about just doing the causal part. The anti-causal part is typically implemented by time-flipping the signal and running it through the causal filter again anyway. Why try to expand the anti-causal part at all? $\endgroup$
    – Hilmar
    Commented Apr 26 at 20:48
  • 1
    $\begingroup$ @Hilmar My code runs on the GPU, so increasing parallelism by running all filters at the same time is beneficial for me. Further, it also simplifies boundary handing since I don't need to correct the right boundary when running the anti-causal filter. Regardless, I would like to know the solution out of curiosity, since I saw some code doing it and I was curious how it was done. $\endgroup$
    – Omar Emara
    Commented Apr 26 at 20:56
  • $\begingroup$ I'm confused: is this an offline or real time applications ? If it's real time, you can't do anti causal filters at all. If it's offline, you can parallelize things simply my running multiple files at the same time. Or is it just a single very long file ? $\endgroup$
    – Hilmar
    Commented Apr 26 at 23:51
  • $\begingroup$ @Hilmar It is an offline application, operating on a single image that could be very small or very large. Each compute unit gets assigned a row of the image to apply the filters, so if the number of compute units was more than the number of rows, you will have a number of compute units doing nothing waiting for the causal filter to finish. On the other hand, if you have a parallel bank of 4 filters, each row will have 4 compute units processing it at the same time, quadrupling the parallelism potential. $\endgroup$
    – Omar Emara
    Commented Apr 27 at 7:17
  • $\begingroup$ I understand that. But as far as I know the only way to implement the anti-causal part of the filter is to flip the input and run it again. So the only part of the filter that actually needs to be implemented is the causal one. And that you can turn into parallel sections. Is that where you having trouble with ? For example, I foudn that Matlab's residuez() has severe numerical problems and it's been basically useless for my application space, so I wrote my own. But that has nothing to do with anything anti-causal $\endgroup$
    – Hilmar
    Commented Apr 27 at 18:25


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.