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