[SciPy-User] deconvolution of 1-D signals

[Ralf Gommers]

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