[SciPy-User] deconvolution of 1-D signals

[josef.pktd at gmail.com]

On Sun, Jul 31, 2011 at 3:41 PM, Ralf Gommers
<ralf.gommers at googlemail.com> wrote:
>
>
> 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.

I mixed up numerator and denominator for lfilter

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

In terms of IIR filter the way I understand it from the ARMA analogy,
your window is not invertible

>>> r = np.roots(window)
>>> r*r.conj()
array([ 1.+0.j,  1.+0.j,  1.+0.j,  1.+0.j,  1.+0.j,  1.+0.j,  1.+0.j,
        1.+0.j,  1.+0.j,  1.+0.j])

I'm not sure if this can work, but in IIR filters the way I know it,
the coefficient for the last observation is normalized to 1. Since all
roots are 1, the window cannot be inverted.
I think lfilter uses the same assumptions. lfilter has also a funny
way to determine initial conditions (zi) which I'm never quite sure
how to use. It doesn't matter much with invertible and stationary
processes, but in this case, I guess it does.

If I change the window to an invertible window

window_len = 11.
window = np.ones(window_len) / window_len
window[0] = 1.

then my current version, which should be close to your original
version, works, the deconvolved series looks similar to the original
series.

In either case, I think Joe's answer is more useful, since you can
directly manipulate the frequency response. I only used a simple IIR
frequency domain filter with an adaptation of fftconvolve.

Josef

>>
>> 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
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>