[SciPy-User] deconvolution of 1-D signals

[Ralf Gommers]

On Sun, Jul 31, 2011 at 9:10 PM, <josef.pktd at gmail.com> wrote:

> On Sun, Jul 31, 2011 at 1:56 PM, Ralf Gommers
> <ralf.gommers at googlemail.com> wrote:
> > 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()
>
> signal.deconvolve is essentially signal.lfilter, but I don't quite
> understand what it does.
>
> I changed 2 lines, partly by trial and error and by analogy to ARMA
> models. I'm not quite sure the following changes are correct, but at
> least it produces a nice graph
>
> It doesn't have artifacts, but y_meas is almost identical to the spline,
not to the real signal. Not sure what's going wrong there, but it didn't
perform a deconvolution.

instead of deconvolve use lfilter directly
>
> y_deconv = signal.lfilter(window, [1.],  ys[::-1])[::-1]
>
> and
>
> center lfiltered window:
>
> ax.plot(xs[window.size//2-1:-window.size//2], y_deconv[:-window.size+1],
> 'k-',
>
> If your signal is periodic, then I would go for the fft versions of
> convolution, and iir filtering.
>

The signal is not periodic.

Ralf

My initial guesses were that there is either something wrong (hidden
> assumption) about the starting values of signal convolve, or there are
> problems because of the non-stationarity.
>
> But maybe a signal expert knows better.
>
> Josef
>
> >
> >
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User at scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> >
> _______________________________________________
> 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/20110731/8cca2324/attachment.html>