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 [ ]: