In [1]:
%matplotlib inline
from pylab import *
from numpy import *
from numpy.random import normal, seed

In [2]:
nelem = 40
noise = 0.2

x = linspace(0, 10, nelem)
f = sin(x)
seed(1234)
fn = f + normal(0.0, noise, size=nelem)

subplot(121)
plot(x, fn)

subplot(122)
plot(x[:-1], diff(f))
plot(x[:-1], diff(fn))

Out[2]:
[<matplotlib.lines.Line2D at 0x7fd4140f4358>]
In [3]:
# the antiderivative matrix
def A(n):
antideriv = ones((n, n))
antideriv[triu_indices(n, k=1)] = 0.
return antideriv

In [4]:
print(A(10))

[[ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 1.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 1.  1.  1.  0.  0.  0.  0.  0.  0.  0.]
[ 1.  1.  1.  1.  0.  0.  0.  0.  0.  0.]
[ 1.  1.  1.  1.  1.  0.  0.  0.  0.  0.]
[ 1.  1.  1.  1.  1.  1.  0.  0.  0.  0.]
[ 1.  1.  1.  1.  1.  1.  1.  0.  0.  0.]
[ 1.  1.  1.  1.  1.  1.  1.  1.  0.  0.]
[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  0.]
[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]]

In [5]:
# check that the inverse of antiderivative approximates derivative
inv(A(10))

Out[5]:
array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
[ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.]])
In [7]:
AI = A(nelem)

In [8]:
df = inv(A(nelem)).dot(f)
dfn = inv(A(nelem)).dot(fn)

In [9]:
plot(x, df)
plot(x, dfn)

Out[9]:
[<matplotlib.lines.Line2D at 0x7fd4141b90f0>]
In [10]:
#calculate the condition number
u, s, v = svd(inv(A(nelem)))
k = max(s)/min(s)
print(k)

51.5306511313

In [15]:
# ruzne mu 1. 5.
mu = 5.
xmu = inv(AI.transpose().dot(AI)+mu**2*eye(nelem)).dot(AI.transpose()).dot(fn)

In [16]:
plot(x, xmu)
plot(x, df)
plot(x, dfn)

Out[16]:
[<matplotlib.lines.Line2D at 0x7fd40f471e48>]
In [17]:
mu = 4.

In [18]:
plot(x, xmu)
plot(x, df)
plot(x, dfn)

Out[18]:
[<matplotlib.lines.Line2D at 0x7fd40f3dcbe0>]
In [19]:
from savitzky_golay import *
plot(x, xmu)
plot(x, df)
plot(x, dfn)
plot(x, savitzky_golay(f, 15, 3, 1, rate=1.))

Out[19]:
[<matplotlib.lines.Line2D at 0x7fd40f354828>]
In [ ]:
#another example:
nelem = 40
noise = 0.2

xmin = -1
xmax = 1
x = linspace(xmin, xmax, nelem)
f = abs(x) - abs(xmin)
fn = f + normal(0.0, noise, size=nelem)

AI = A(nelem)


In [ ]:
mu = 4.

In [ ]:
plot(x, savitzky_golay(f, 15, 4, 1, rate=1.))
plot(x, xmu)
plot(x, df)
#plot(x, dfn)

In [ ]: