[SciPy-User] deconvolution of 1-D signals
David Baddeley
david_baddeley at yahoo.com.au
Sun Jul 31 20:51:28 EDT 2011
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()
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110731/a6155e33/attachment.html>
More information about the SciPy-User
mailing list