[SciPy-User] deconvolution of 1-D signals

[Joe Kington]

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.

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