[SciPy-User] deconvolution of 1-D signals

[Ralf Gommers]

On Mon, Aug 1, 2011 at 2:51 AM, 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.
>

The kernel is typically about 5 to 15 times smaller than the signal extent,
so I guess that may be problematic.

>
> 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).
>

I'm aware of the problems with high frequency noise. This is why I tried the
spline fitting - I figured that on a spline the deconvolution would be okay
because the spline is very smooth. This should be fine for my data because
the noise is much higher-frequency than the underlying signal, and the SNR
is high to start with. But maybe there are better ways. I looked for a
Python implementation of Wiener deconvolution but couldn't find one so
quickly. Is there a package out there that has it?

>
> 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.
>
> By constrained to be positive I'm guessing you mean monotonic? Otherwise I
could just add a constant offset, but that's probably not what you mean.
What's typically the speed penalty for an iterative method?

Ralf


>
>
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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110801/9fb42a79/attachment.html>