[SciPy-User] convolve/deconvolve

josef.pktd@gmai... josef.pktd@gmai...
Fri Feb 1 22:54:20 CST 2013


On Fri, Feb 1, 2013 at 11:02 AM, jkhilmer@chemistry.montana.edu
<jkhilmer@chemistry.montana.edu> wrote:
> I have some old code for Richardson-Lucy deconvolution, although the method
> is so simple there's no reason not to try it yourself.
>
> Jonathan
>
>
>
> On Fri, Feb 1, 2013 at 7:50 AM, <josef.pktd@gmail.com> wrote:
>>
>> On Fri, Feb 1, 2013 at 9:50 AM,  <josef.pktd@gmail.com> wrote:
>> > On Fri, Feb 1, 2013 at 9:23 AM, François Boulogne
>> > <fboulogne@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

just occured to me: leading zeros just shift the array,
lets try to drop them,
here's an almost roundtrip with "full"

>>> filter_
array([0, 1, 1, 1, 1, 0])
>>> a = np.round(np.random.randn(22), 2)
>>> res2 = signal.convolve(a, filter_, 'full')
>>> s, z = signal.deconvolve(res2, filter_[1:])   #drop leading zero in filter
>>> a
array([-0.96, -0.99,  1.3 ,  0.76, -1.1 , -0.3 , -0.51,  0.86,  1.53,
       -0.11, -0.21,  0.46, -0.48, -0.64, -0.85,  0.02,  0.21, -0.5 ,
        2.09,  1.47, -1.73,  0.01])
>>> s
array([ 0.  , -0.96, -0.99,  1.3 ,  0.76, -1.1 , -0.3 , -0.51,  0.86,
        1.53, -0.11, -0.21,  0.46, -0.48, -0.64, -0.85,  0.02,  0.21,
       -0.5 ,  2.09,  1.47, -1.73,  0.01])
>>> len(a)
22
>>> len(s)
23
>>> np.allclose(a, s[1:])
True

Josef


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


More information about the SciPy-User mailing list