[SciPy-User] deconvolution of 1-D signals

[David Baddeley]

Hi Ralf,

5-15 times smaller would probably be fine, although you might want to watch the 
edges in the reconstruction - if they're at different dc levels you'll get edge 
artifacts (within ~ 1 kernel width of the edges). I'd tend to avoid spline 
filtering (or any form of noise reduction) before deconvolution, as this will 
also transform the data in a way not explained by the model you're using to 
deconvolve with. 

Weiner filtering is a 2 liner -> 

H = fft(kernel)
deconvolved = ifftshift(ifft(fft(signal)*np.conj(H)/(H*np.conj(H) + lambda**2)))

where lambda is your regularisation parameter, and white noise is assumed. There 
are various methods for choosing lambda optimally, but most people tend to use 
trial and error.

Iterative methods are typically ~1-2 orders of magnitude slower than a Weiner 
filter, but with fast fft libraries and modern computers still quite reasonable 
for modest data sizes (a 3D image stack of ~ 512x512x50 pixels will tend to be 
done in a bit under a minute, can't really comment on 1D data, but unless your 
signal is very long I'd expect it to be significantly quicker). Ffts scale with 
O(nlogn) so will generally dramatically outperform things based on a simple 
convolution or filtering approaches (O(n**2)) for large n. This might make an 
iterative approach using ffts faster than something like scipy.signal.deconvolve 
if your kernel is large.

cheers,
David



________________________________
From: Ralf Gommers <ralf.gommers at googlemail.com>
To: David Baddeley <david_baddeley at yahoo.com.au>; SciPy Users List 
<scipy-user at scipy.org>
Sent: Mon, 1 August, 2011 6:03:10 PM
Subject: Re: [SciPy-User] deconvolution of 1-D signals




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/f81b43ee/attachment.html>