[SciPy-User] deconvolution of 1-D signals

[Anne Archibald]

I realize this discussion has gone rather far afield from efficient 1D
deconvolution, but we do a funny thing in radio interferometry, and
I'm curious whether this is normal for other kinds of deconvolution as
well.

In radio interferometry we obtain our images convolved with the
so-called "dirty beam", a convolution kernel that has a nice narrow
peak but usually a chaos of monstrous sidelobes often only marginally
smaller than the main lobe. We use a different regularization
condition to do our deconvolution: we treat the underlying image as a
modest collection of point sources. (One can see why this appeals to
astronomers.) Through an iterative process (the "CLEAN" algorithm and
its many descendants) we obtain an estimate of this underlying image.
But we very rarely actually work with this image directly. We normally
convolve it with a sort of idealized version of our kernel without all
the sidelobes. This then gives an image one might have obtained from a
normal telescope the size of the interferometer array. (Apart from all
the CLEAN artifacts.)

What I'm wondering is, is this final step of convolving with an
idealized version of the kernel standard practice elsewhere?

>From one point of view it could just be parochiality, astronomers
being so accustomed to smudgy images that we have to convert anything
else to this format. But I think that at the least it softens the
effect of the rather strict regularization assumption behind CLEAN -
which amounts to "no extended sources". It probably makes us less
sensitive to shortcuts in CLEAN implementations. I think, though, that
this trick may be useful for many applications of deconvolution.
Rather than try to translate the image from the observed kernel to
some ideal Dirac-delta kernel, this tries to convert it from the
observed kernel to a similar but simpler kernel; one would expect the
impact of a deconvolution artifact to be related to the magnitude of
the difference between kernels.

In terms of 1D Fourier deconvolution, this is saying, after
deconvolution, that we don't really need all those high frequencies
amplified so much anyway, and smoothing them back down with a nice
clean easy-to-understand kernel. In these terms, in fact, it makes
perfect sense to use a wider kernel than necessary for this smoothing
if one is interested in larger-scale features.

Anne

On 31 July 2011 20:51, David Baddeley <david_baddeley at yahoo.com.au> wrote:
> Hi Ralf,
> I do a reasonable amount of (2 & 3D) deconvolution of microscopy images and
> the method I use depends quite a lot on the exact properties of the signal.
> You can usually get away with fft based convolutions even if your signal is
> not periodic as long as your kernel is significantly smaller than the signal
> extent.
> As Joe mentioned, for a noisy signal convolving with the inverse or
> performing fourier domain division doesn't work as you end up amplifying
> high frequency noise components. You thus need some form of regularisation.
> The thresholding of fourier components that Joe suggests does this, but you
> might also want to explore more sophisticated options, the simplest of which
> is probably Wiener filtering
> (http://en.wikipedia.org/wiki/Wiener_deconvolution).
> If you've got a signal which is constrained to be positive, it's often
> useful to introduce a positivity constraint on the deconvolution result
> which generally means you need an iterative algorithm. The choice of
> algorithm should also depend on the type of noise that is present in your
> signal -  my image data is constrained to be +ve and typically has either
> Poisson or a mixture of Poisson and Gaussian noise and I use either the
> Richardson-Lucy or a weighted version of ICTM (Iterative Constrained
> Tikhonov-Miller) algorithm. I can provide more details of these if required.
> cheers,
> David
>
>
>
> ________________________________
> From: Ralf Gommers <ralf.gommers at googlemail.com>
> To: SciPy Users List <scipy-user at scipy.org>
> Sent: Mon, 1 August, 2011 5:56:49 AM
> Subject: [SciPy-User] deconvolution of 1-D signals
>
> Hi,
>
> For a measured signal that is the convolution of a real signal with a
> response function, plus measurement noise on top, I want to recover the real
> signal. Since I know what the response function is and the noise is
> high-frequency compared to the real signal, a straightforward approach is to
> smooth the measured signal (or fit a spline to it), then remove the response
> function by deconvolution. See example code below.
>
> Can anyone point me towards code that does the deconvolution efficiently?
> Perhaps signal.deconvolve would do the trick, but I can't seem to make it
> work (except for directly on the output of np.convolve(y, window,
> mode='valid')).
>
> Thanks,
> Ralf
>
>
> import numpy as np
> from scipy import interpolate, signal
> import matplotlib.pyplot as plt
>
> # Real signal
> x = np.linspace(0, 10, num=201)
> y = np.sin(x + np.pi/5)
>
> # Noisy signal
> mode = 'valid'
> window_len = 11.
> window = np.ones(window_len) / window_len
> y_meas = np.convolve(y, window, mode=mode)
> y_meas += 0.2 * np.random.rand(y_meas.size) - 0.1
> if mode == 'full':
>     xstep = x[1] - x[0]
>     x_meas = np.concatenate([ \
>         np.linspace(x[0] - window_len//2 * xstep, x[0] - xstep,
> num=window_len//2),
>         x,
>         np.linspace(x[-1] + xstep, x[-1] + window_len//2 * xstep,
> num=window_len//2)])
> elif mode == 'valid':
>     x_meas = x[window_len//2:-window_len//2+1]
> elif mode == 'same':
>     x_meas = x
>
> # Approximating spline
> xs = np.linspace(0, 10, num=500)
> knots = np.array([1, 3, 5, 7, 9])
> tck = interpolate.splrep(x_meas, y_meas, s=0, k=3, t=knots, task=-1)
> ys = interpolate.splev(xs, tck, der=0)
>
> # Find (low-frequency part of) original signal by deconvolution of smoothed
> # measured signal and known window.
> y_deconv = signal.deconvolve(ys, window)[0]  #FIXME
>
> # Plot all signals
> fig = plt.figure()
> ax = fig.add_subplot(111)
>
> ax.plot(x, y, 'b-', label="Original signal")
> ax.plot(x_meas, y_meas, 'r-', label="Measured, noisy signal")
> ax.plot(xs, ys, 'g-', label="Approximating spline")
> ax.plot(xs[window.size//2-1:-window.size//2], y_deconv, 'k-',
>         label="signal.deconvolve result")
> ax.set_ylim([-1.2, 2])
> ax.legend()
>
> plt.show()
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>