[SciPy-User] deconvolution of 1-D signals

[Joe Kington]

On Sun, Jul 31, 2011 at 2:05 PM, Ralf Gommers
<ralf.gommers at googlemail.com>wrote:

>
>
> On Sun, Jul 31, 2011 at 9:22 PM, Joe Kington <jkington at wisc.edu> wrote:
>
>> Gah, I hit send too soon!
>>
>> The default eps parameter in that function should be more like 1.0e-6
>> instead of 0.1.
>>
>> You'll generally need to adjust the eps parameter to match the
>> signal-to-noise ratio of the two signals you're deconvolving.
>>
>
> Thanks Joe. I had tried something similar, but now I have a name for the
> method and a confirmation that doing something like that makes sense. I'll
> play with this idea some more tomorrow.
>

For what it's worth, the implementation I showed there is completely wrong.
I wasn't thinking clearly when I wrote that out.  The same concept should
still work, though.

The most glaring mistake is that things should be padded before converting
to the frequency domain.  More like this:

def water_level_decon(y_meas, window, eps=0.1):
    padded = np.zeros_like(y_meas)
    padded[:window.size] = window

    yfreq = np.fft.fft(y_meas)
    winfreq = np.fft.fft(padded)

    winfreq[winfreq < eps] = eps
    newfreq = yfreq / winfreq

    return np.fft.ifft(newfreq)

Hope it's useful, anyway.
-Joe




> Cheers,
> Ralf
>
>
>
>>
>> Hope it's useful, at any rate.
>> -Joe
>>
>>
>> On Sun, Jul 31, 2011 at 11:21 AM, Joe Kington <jkington at wisc.edu> wrote:
>>
>>> I'm not a signal processing expert by any means, but this is a standard
>>> problem in seismology.
>>>
>>> The problem is that your "window" has near-zero amplitude at high
>>> frequencies, so you're blowing up the high-frequency content of the noisy
>>> signal when you divide in the frequency domain.
>>>
>>> A "water-level" deconvolution is a very simple way around this, and often
>>> works well.  It also allows you to skip the spline fitting, as it's
>>> effectively doing a low-pass filter.
>>>
>>> Basically, you just:
>>> 1) convert to the frequency domain
>>> 2) replace any amplitudes below some threshold with that threshold in the
>>> signal you're dividing by (the window, in your case)
>>> 3) pad the lengths to match
>>> 4) divide (the actual deconvolution)
>>> 5) convert back to the time domain
>>>
>>> As a simple implementation (I've left out the various different modes of
>>> padding here... This is effectively just mode='same'.)
>>>
>>> def water_level_decon(ys, window, eps=0.1):
>>>     yfreq = np.fft.fft(ys)
>>>     max_amp = yfreq.max()
>>>
>>>     winfreq = np.fft.fft(window)
>>>     winfreq[winfreq < eps] = eps
>>>
>>>     padded = eps * np.ones_like(yfreq)
>>>     padded[:winfreq.size] = winfreq
>>>
>>>     newfreq = yfreq / padded
>>>     newfreq *= max_amp / newfreq.max()
>>>
>>>     return np.fft.ifft(newfreq)
>>>
>>>
>>> Hope that helps a bit.
>>>  -Joe
>>>
>>>
>>> In most cases, you'll need to adjust the eps parameter to match the level
>>> of noise you want to remove. In your particular case
>>> On Sun, Jul 31, 2011 at 9:56 AM, 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()
>>>>
>>>>
>>>> _______________________________________________
>>>> 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110731/3f90eea1/attachment.html>