[SciPy-User] deconvolution of 1-D signals

[Ralf Gommers]

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.

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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110801/97fd6c10/attachment.html>