In [1]:

    

%matplotlib inline
from pylab import *
from numpy import *
from numpy.random import poisson, normal


    

In [2]:

    

xmin = 0
xmax = 10
xlen = xmax-xmin
x = linspace(xmin, xmax, 200)
signal = zeros_like(x)
I = (x>(xmin+xlen*0.4)) & (x<(xmin+xlen*0.6))
signal[I] = x[I]
plot(x, signal)


    

    
Out[2]:





[<matplotlib.lines.Line2D at 0x7ff4dd12e5c0>]






    

In [3]:

    

from scipy.stats import norm
rv = norm(loc=xlen/2, scale=0.2)
psf = rv.pdf(x)
psf /= sum(psf)
plot(x, psf)


    

    
Out[3]:





[<matplotlib.lines.Line2D at 0x7ff4ce5e4710>]






    

In [4]:

    

csignal = convolve(psf, signal, mode="same")
plot(x, csignal)


    

    
Out[4]:





[<matplotlib.lines.Line2D at 0x7ff4ce54e828>]






    

In [5]:

    

from numpy import fft as F
from numpy.fft import fft, ifft, ifftshift
H = fft(psf)
S = fft(csignal)
DS = ifftshift(ifft(S/H))
plot(x, DS)
plot(x, signal)


    

    




/usr/local/lib/python3.4/dist-packages/numpy/core/numeric.py:462: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)






    
Out[5]:





[<matplotlib.lines.Line2D at 0x7ff4ce55bc18>]






    

In [6]:

    

semilogy(x, abs(fftshift(H)))


    

    
Out[6]:





[<matplotlib.lines.Line2D at 0x7ff4ce5a2240>]






    

In [7]:

    

semilogy(x, abs(fftshift(S/H)))


    

    
Out[7]:





[<matplotlib.lines.Line2D at 0x7ff4ce2639e8>]






    

In [8]:

    

semilogy(x, abs(fftshift(fft(signal))))


    

    
Out[8]:





[<matplotlib.lines.Line2D at 0x7ff4ce3737f0>]






    

In [9]:

    

# naive regularization
DS = ifftshift(ifft(S/(H+1e-11)))
plot(x, DS)
plot(x, signal)


    

    
Out[9]:





[<matplotlib.lines.Line2D at 0x7ff4ce16c748>]






    

In [11]:

    

# Wiener deconvolution
lamb = 1e-12
WDS = ifftshift(ifft(S*conj(H)/(H*conj(H) + lamb**2)))
plot(x, WDS)


    

    
Out[11]:





[<matplotlib.lines.Line2D at 0x7ff4ce0059e8>]






    

In [12]:

    

noise = normal(loc=0.0, scale=0.1, size=len(x))
ncsignal = csignal+noise
NS = fft(ncsignal)
plot(x, ncsignal)


    

    
Out[12]:





[<matplotlib.lines.Line2D at 0x7ff4ce187e80>]






    

In [13]:

    

# naive regularization
DS = ifftshift(ifft(NS/(H+.2)))
plot(x, DS)
plot(x, signal)


    

    
Out[13]:





[<matplotlib.lines.Line2D at 0x7ff4ce00a710>]






    

In [14]:

    

# Wiener deconvolution
lamb = 0.005
WDS = ifftshift(ifft(NS*conj(H)/(H*conj(H) + lamb)))
plot(x, WDS)
plot(x, signal)


    

    
Out[14]:





[<matplotlib.lines.Line2D at 0x7ff4ce18a240>]






    

In [15]:

    

def deconvolve_iterative(data, kernel, niter):
    # http://dx.doi.org/10.1086/111605
    from scipy.signal import convolve
    P = kernel
    I = data
    O = convolve(I, P, mode="same")
    eps = 1e-10 # this is to avoid division by zero
    for i in range(niter):
        denom = convolve(O, P, mode="same") + eps
        fact = convolve(I/denom, P[::-1], mode="same")
        O = fact*O
        O[O<0] = 0
    return O


    

In [16]:

    

LRDS = deconvolve_iterative(ncsignal, psf, 8)
plot(x, LRDS)
plot(x, signal)


    

    




/usr/local/lib/python3.4/dist-packages/numpy/core/fromnumeric.py:2507: VisibleDeprecationWarning: `rank` is deprecated; use the `ndim` attribute or function instead. To find the rank of a matrix see `numpy.linalg.matrix_rank`.
  VisibleDeprecationWarning)






    
Out[16]:





[<matplotlib.lines.Line2D at 0x7ff4cdfcbf60>]






    

In [17]:

    

# more realistic example:
xmin = 0
xmax = 10
xlen = xmax-xmin
x = linspace(xmin, xmax, 200)
signal = exp(-(x-xmax*0.4)**2/(0.2))*20
signal += exp(-(x-xmax*0.55)**2/(0.4))*30
plot(x, signal)

#plot(x, psf)


    

    
Out[17]:





[<matplotlib.lines.Line2D at 0x7ff4cc6070b8>]






    

In [18]:

    

rv = norm(loc=xlen/2, scale=0.5)
psf = rv.pdf(x)
psf /= sum(psf)
csignal = convolve(psf, signal, mode="same")
plot(x, csignal)


    

    
Out[18]:





[<matplotlib.lines.Line2D at 0x7ff4cc5def28>]






    

In [19]:

    

from numpy.random import poisson
seed(12)
ncsignal = poisson(csignal)
plot(x, ncsignal)


    

    
Out[19]:





[<matplotlib.lines.Line2D at 0x7ff4cc53b438>]






    

In [22]:

    

H = fft(psf)
NS = fft(ncsignal)
lamb = 0.001
WDS = ifftshift(ifft(NS*conj(H)/(H*conj(H) + lamb)))
plot(x, WDS)
plot(x, signal)


    

    
Out[22]:





[<matplotlib.lines.Line2D at 0x7ff4cc4d8240>]






    

In [29]:

    

LRDS = deconvolve_iterative(ncsignal, psf, 10)
plot(x, LRDS)
plot(x, signal)


    

    
Out[29]:





[<matplotlib.lines.Line2D at 0x7ff4cc1ce7f0>]






    

In [ ]: