deconv.py

# -*- coding: utf-8 -*-

# see:
# http://eeweb.poly.edu/iselesni/lecture_notes/least_squares/
#     LeastSquares_SPdemos/deconvolution/html/deconv_demo.html

import numpy as np
from scipy import linalg, sparse
import matplotlib.pyplot as plt
from pyeparse.utils import pupil_kernel

plt.ion()

fs = 10.
rng = np.random.RandomState(0)

data_scale = 1e2
# """
dur = 10.
h = pupil_kernel(fs)
N = int(fs * dur)
h_t = np.arange(len(h)) / fs
x = np.abs(rng.randn(N))
x *= np.hanning(N)
x *= data_scale
"""
N = 300
n = np.arange(N, dtype=float)
w = 5.
n1 = 70.
n2 = 130.
x = 2.1 * np.exp(-0.5*((n - n1) / w) ** 2)
x -= 0.5 * np.exp(-0.5 * ((n - n2) / w) ** 2) * (n2 - n)
x *= data_scale
h = n * (0.9 ** n) * np.sin(0.2 * np.pi * n)
"""

t = np.arange(N) / fs
y = np.convolve(x, h)[:len(x)]
fig, axs = plt.subplots(2, 1)
axs[0].plot(t, x, 'k')
axs[1].plot(t, y, 'k')

# construct the convolution matrix
H = np.zeros((N, N))
for ii in range(N):
    n_samp = min(len(h), N - ii)
    H[ii:ii + n_samp, ii] = h[:n_samp]
np.testing.assert_allclose(np.dot(H, x), y)
y += 0.2 * rng.randn(N)

# diagonal loading
lambda_0 = 0.1
x_diag = linalg.solve(np.dot(H.T, H) + lambda_0 * np.eye(N), np.dot(H.T, y))
axs[0].plot(t, x_diag, 'r')

# derivative regularization
e = np.ones(N)
D = sparse.spdiags([e, -2*e, e], np.arange(3), N - 2, N)
print(D.shape)
lambda_1 = 2.
x_deriv = linalg.solve(np.dot(H.T, H) + lambda_1 * D.T.dot(D), np.dot(H.T, y))
axs[0].plot(t, x_deriv, 'g')