[SciPy-User] convolve/deconvolve

[josef.pktd at gmail.com]

On Fri, Feb 1, 2013 at 9:50 AM,  <josef.pktd at gmail.com> wrote:
> On Fri, Feb 1, 2013 at 9:50 AM,  <josef.pktd at gmail.com> wrote:
>> On Fri, Feb 1, 2013 at 9:23 AM, François Boulogne <fboulogne at sciunto.org> wrote:
>>> Hi,
>>>
>>> I need to deconvolve a signal with a filter. I had a look in the
>>> documentation. The function exists but the docstring is missing and I'm
>>> not satisfied of the result I got from a "simple" example.
>>>
>>> filter = np.array([0,1,1,1,1,0])
>>> step = np.array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
>>> 1, 1, 1, 1])
>>> # I convolve both
>>> res = convolve(step, filter, 'valid')
>>> # and it returns a slope as expected
>>> array([0, 0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
>>>
>>> Now, I want to deconvolve.
>>> deconvolve(res, filter)
>>> # oops, it raises an exception
>>> ValueError: BUG: filter coefficient a[0] == 0 not supported yet
>>>
>>> # well, let's try this
>>> deconvolve(res, filter+1e-9)
>>> (array([  0.00000000e+00,   0.00000000e+00,   1.00000000e+09,
>>>         -9.99999999e+17,   9.99999999e+26,  -9.99999999e+35,
>>>          9.99999999e+44,  -9.99999999e+53,   9.99999999e+62,
>>>         -9.99999999e+71,   9.99999999e+80,  -9.99999999e+89]),
>>>  array([  0.00000000e+00,   0.00000000e+00,   1.11022302e-16,
>>>         -8.27130862e-08,  -4.42500000e+01,   5.46901335e+10,
>>>          8.27266814e+19,   7.56858250e+28,  -8.74285726e+37,
>>>          9.99419626e+46,   8.27205507e+55,  -8.26933326e+64,
>>>          9.99999999e+89,   9.99999999e+89,   9.99999999e+89,
>>>          1.00000000e+90,   9.99999999e+80]))
>>>
>>> It's better but I do not recognize my signal :)
>>> 1/ Am I misunderstanding or missing something?
>>> 2/ How can I do it correctly?
>>
>> AFAICS:
>>
>> not supported, maybe using the numpy polynomial might work for the deconvolution
>>
>> from the docstring of lfilter, which is used by deconvolve:
>>
>>    a : array_like
>>         The denominator coefficient vector in a 1-D sequence.  If ``a[0]``
>>         is not 1, then both `a` and `b` are normalized by ``a[0]``.
>>
>> with the normalization your 1e-9 blows up the calculations.
>>
>> (it's been a long time since I tried to figure out deconvolve, and I
>> always had 1 in the first position)

another problem in your case are the unit roots

>>> np.roots(filter_)
array([ -1.00000000e+00+0.j,  -7.77156117e-16+1.j,  -7.77156117e-16-1.j,
         0.00000000e+00+0.j])

I don't remember whether lfilter supports those. Or, if setting
explicit initial conditions would help.
Or, if deconvolution makes sense with a nonstationary sequence (infinite ?).

Josef


>>
>> Josef
>>
>>
>>>
>>> I also noted that no test exists for deconvolve() :(
>
> Volunteers ?
>
> Josef
>
>
>>>
>>> Cheers,
>>> François.
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user