Approximating inverse of unstable difference of Gaussians filter

I am trying to invert a difference of Gaussians (DoG) filter. The inverse is not stable and so I am trying to find an approximation applied to a specific input. The DoG filter increases contrast at edges, and I'm trying to find an inverse that decreases contrast at the edges.

The approximate inverse is not immediately clear and my attempts have gotten me close but not quite there.

All code is available in this Google Colab notebook.

Input

My input is a 1D series of ascending steps. It's created as such:

luminances = np.repeat(np.linspace(0, 1, 5), 100)

And looks like this: DoG Filter

My DoG filter is made with the following code:

filter_len = 81
excitatory_gaussian = scipy.signal.gaussian(filter_len, 1.4)
inhibitory_gaussian = scipy.signal.gaussian(filter_len, 16) * -0.043
filter = excitatory_gaussian + inhibitory_gaussian
filter /= sum(filter)
plt.plot(filter)

It looks like this: Convolution

When I convolve my input with the DoG filter, I see higher contrast at the edges, as expected: I can create the same output with either convolve or lfilter:

deconvolved_luminances = signal.convolve(luminances, filter)
deconvolved_luminances = signal.lfilter(filter, , luminances)

Inverting Attempts

I want to find the inverse of the DoG filter such that, when convolving this input with the inverted filter, the output is ascending steps with decreased contrast at the edges and can be convolved with my original DoG filter to output an approximate of the original input.

signal.deconvolve

Seemingly because the inverse of the DoG filter isn't stable, using signal.deconvolve results in a bad output.

deconvolved_luminances, remainder = scipy.signal.deconvolve(luminances, filter) As expected, the output is the same when using lfilter instead of deconvolve:

deconvolved_luminances = scipy.signal.lfilter(, filter, luminances)

Inverting half of filter

It appears that the right half of the filter is stable. If I take this right half (half_filter_r = filter[len(filter)//2:]) and deconvolve it with my signal (s_deconvolved_r = signal.lfilter(, half_filter_r, s), I get a sensical output, which looks like this: Taking the left half of the filter (half_filter_l = filter[:len(filter)//2 + 1]) and doing the same, however, produces an output that suggests the left half of the filter is unstable: I can assume the appropriate output of the left half is simply the inverse of the right half's output (s_deconvolved_l = max(s_deconvolved_r) - s_deconvolved_r[::-1]) and can add this supposed left half output together with the actual right half output (s_deconvolved_r + s_deconvolved_l) and I end up with something that looks correct: However, assuming the left half's output based on the right half seems suspicious and seems that there should be a better way. Why would the left half be unstable while the right half is stable? Also, when I convolve my filter with the above output, I get a result that looks similar to the original filter but not quite right (hit my limit on images I'm allowed to post, but it's in the colab notebook).

Thanks for any help!

If you have a filter $$\sum a_n z^n$$, the deconvolution filter $$1/ \sum a_n z^n$$ attempts to cancel the earliest sample of the input by adding a scaled version of the filter, then proceed to the next sample (like a gaussian elimination process).

If the first sample in the filter is small, the gain applied to the filter must be large, and that will cause the subsequent samples to become even higher.

Yes, the original filter has roots outside the unit circle: so when it's used inverted it will be unstable.

One approach might be to form the minimum phase equivalent:

minphase_filter = signal.minimum_phase(filter)

that will guarantee roots inside the unit circle.

The resulting impulse response is: this yields the filtered signal: and this allows s to be deconvolved stably: This is a tricky problem. Roughly speaking, the root cause here is the original filter "destroys" information which cannot be recovered, so you have to make some trade offs

The best way to design the inverse is to formulate the inversion as a least square error problem where the error criteria is used to dial in the trade offs for your specific application.

Let's say your transfer function is $$H(z)$$. Let's say we want to design "sort of inverse" filter, $$G(z)$$. So ideally we would have

$$H(z) \cdot G(z) = 1$$ or $$G(z) = \frac{1}{H(z)}$$

However that only works if $$H(z)$$ is reasonably well behaved, which yours isn't. I n particular it's not minimum phase and it has "dense comb" at high frequencies which if inverted, will result in narrow band peaks with gains exceeding 100 dB.

You can try estimating the frequency range over which you think you can reasonable invert the filter. This is a function of your filter and also the expected noise floor (physical or numerical) and then define a target function $$T(z)$$ . A good choice of $$T(z)$$ might be a low pass filter with a cutoff similar to your original filter and it should also include enough bulk delay to allow the inverse to be causal. I would probably start with twice the filter length. Then you can solve for $$G(z)$$ by minimizing

$$|G(z) \cdot H(z) - T(z)|^2 = min$$

That's pretty straight forward for an FIR filter but difficult for an IIR, which requires some intelligent search algorithm.

You may also have to throw in some weighting function and potential some non-linear constraints for maximum coefficients value and/or gain.

I would start with an FIR and start fiddling with the target function, weights, constrains, filter length, bulk delay and see if that's "good enough" for what you to do.