Calculating Inverse Filter

Question

Does anyone has an idea how to calculate the Inverse of a 2-D filter?

Let's say I have a 3x3 filter:
0 1 0
1 1 1
0 1 0
I want to fins it's inverse. It's easy to do using DFT. But let's say I want to do it by convolution. Now, that's the problem, Matlab symbolic isn't my specialty. Assuming there's a 3X3 Inverse Filter it means convolution of the two will result in:
0 0 0
0 1 0
0 0 0
The problem is to create the right set of equations for that and solving it. Doing it with symbols is easy to think yet I couldn't do it.

Any ideas? Thanks.

P.S. I'm not sure there an Inverse Filter for this one as it has zeros in its DTFT.

Moreover, someone should allow Latex in this forum the way MathOverflow has.

Answer 1

Let h[n] be the finite impulse response of a 1D filter. What does that imply about its inverse filter? The inverse cannot possibly be FIR.

Proof: Let H(omega) G(omega) = 1 for all omega, where H is the DTFT of h[n], and G is the DTFT of g[n]. If h[n] is FIR, then g[n] must be IIR.

Of course, there are ways to approximate the inverse IIR filter with an FIR filter. A basic method is adaptive filtering, e.g., Least Mean Squares (LMS) algorithm. Or just truncate the IIR filter. You still need to worry about stability, though.

For practical purposes, there probably is no desirable solution to your specific question. Especially if, for example, this is in image processing and you are trying to inverse an FIR blur filter with FIR sharpening filter. The final image will not look that good, period, unless your sharpening filter is really, really large.

EDIT: Let y[n] = b0 x[n-0] + b1 x[n-1] + ... + bN x[n-N]. Let this equation characterize the forward system, where y is the output and x is the input. By definition, the impulse response is the output when the input is an impulse: h[n] = b0 d[n-0] + b1 d[n-1] + ... + bN d[n-N]. This impulse response has finite length N+1.

Now, consider the inverse system where x is the output and y is the input. Then the impulse response is described by the recurrence equation d[n] = b0 h[n] + b1 h[n-1] + ... + bN h[n-N]. Equivalently, b0 h[n] = d[n] - b1 h[n-1] - ... - bN h[n-N].

Without loss of generality, assume that b0 and bN are both nonzero. For any m, if h[m] is nonzero, then h[m+N] is also nonzero. Because this system has feedback, its impulse response is infinitely long. QED.

Causality does not matter. The inverse of a delay is an advance, and vice versa. Neither a delay or an advance alter the finiteness of an impulse response. Shift an infinite impulse response left or right; it is still infinite.

EDIT 2: For clarification, this proof is not related to my original "proof". One was in frequency domain, the other in time domain.

Answer 2

For a separable filter (i.e., one in which you can pull out a horizontal and vertical filter which can be applied in any order and give the same result as the composite 2-D filter), you could attempt to calculate the inverse of each on individually.

The inverse of a filter H(w) is G(w)=1/H(w), so one way to do it is to take the impulse response (the h[n] time-domain coeffs) and inverse-DFT them. It's not always easy to get an analytic expression for such a filter, so you could either compute it numerically (approximating to the desired precision) or do adaptive inverse filtering, as Steve suggested. See Widrow and Stearn's Adaptive Signal Processing for more info on this latter method.

Answer 3

In practice, one useful solution is Wiener deconvolution http://en.wikipedia.org/wiki/Wiener_deconvolution which basically, in Fourier space, divides by the spectrum of the given filter. Zeros are handled by adding a fudge term: instead of 1/H(w), use H(w)/( abs(H(w))^2 + c) where H(w) is the discrete Fourier transform of the filter h(x) (add a 2nd dimension as you like) and "w" is supposed to be omega. The constant is chosen based on the noise level of the signal or image.