%matplotlib inline
from pylab import *
from numpy import *
from numpy.random import poisson, normal
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)
[<matplotlib.lines.Line2D at 0x7ff4dd12e5c0>]
from scipy.stats import norm
rv = norm(loc=xlen/2, scale=0.2)
psf = rv.pdf(x)
psf /= sum(psf)
plot(x, psf)
[<matplotlib.lines.Line2D at 0x7ff4ce5e4710>]
csignal = convolve(psf, signal, mode="same")
plot(x, csignal)
[<matplotlib.lines.Line2D at 0x7ff4ce54e828>]
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)
[<matplotlib.lines.Line2D at 0x7ff4ce55bc18>]
semilogy(x, abs(fftshift(H)))
[<matplotlib.lines.Line2D at 0x7ff4ce5a2240>]
semilogy(x, abs(fftshift(S/H)))
[<matplotlib.lines.Line2D at 0x7ff4ce2639e8>]
semilogy(x, abs(fftshift(fft(signal))))
[<matplotlib.lines.Line2D at 0x7ff4ce3737f0>]
# naive regularization
DS = ifftshift(ifft(S/(H+1e-11)))
plot(x, DS)
plot(x, signal)
[<matplotlib.lines.Line2D at 0x7ff4ce16c748>]
# Wiener deconvolution
lamb = 1e-12
WDS = ifftshift(ifft(S*conj(H)/(H*conj(H) + lamb**2)))
plot(x, WDS)
[<matplotlib.lines.Line2D at 0x7ff4ce0059e8>]
noise = normal(loc=0.0, scale=0.1, size=len(x))
ncsignal = csignal+noise
NS = fft(ncsignal)
plot(x, ncsignal)
[<matplotlib.lines.Line2D at 0x7ff4ce187e80>]
# naive regularization
DS = ifftshift(ifft(NS/(H+.2)))
plot(x, DS)
plot(x, signal)
[<matplotlib.lines.Line2D at 0x7ff4ce00a710>]
# Wiener deconvolution
lamb = 0.005
WDS = ifftshift(ifft(NS*conj(H)/(H*conj(H) + lamb)))
plot(x, WDS)
plot(x, signal)
[<matplotlib.lines.Line2D at 0x7ff4ce18a240>]
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
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)
[<matplotlib.lines.Line2D at 0x7ff4cdfcbf60>]
# 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)
[<matplotlib.lines.Line2D at 0x7ff4cc6070b8>]
rv = norm(loc=xlen/2, scale=0.5)
psf = rv.pdf(x)
psf /= sum(psf)
csignal = convolve(psf, signal, mode="same")
plot(x, csignal)
[<matplotlib.lines.Line2D at 0x7ff4cc5def28>]
from numpy.random import poisson
seed(12)
ncsignal = poisson(csignal)
plot(x, ncsignal)
[<matplotlib.lines.Line2D at 0x7ff4cc53b438>]
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)
[<matplotlib.lines.Line2D at 0x7ff4cc4d8240>]
LRDS = deconvolve_iterative(ncsignal, psf, 10)
plot(x, LRDS)
plot(x, signal)
[<matplotlib.lines.Line2D at 0x7ff4cc1ce7f0>]