In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import HTML
from numpy import linalg as LA
import matplotlib.animation as animation
In [1]:
def kalman(n,Nx):
    eta = np.random.normal(0, 0.1, n+1)
    t = np.linspace(0,n,n+1)
    yn  = np.sin(0.1*t + 4*np.sin(0.01*t)) + eta
    yground  = np.sin(0.1*t + 4*np.sin(0.01*t))
    err = np.vstack((eta,np.roll(eta,-1)))
    Ny = np.dot(err,err.T)/n




    M = np.array([(1,0),(1,-1)])
    A = np.eye(2)
    x = np.zeros([1,n])
    x[0,0:2] = np.ones(2)
    xest = np.array([(1),(1)])
    #xhist = np.array([(1,1),(1,1)])
    xhist = np.ones(n)
    B = np.eye(2)
    E = np.eye(2)

    for i in range(2,n):    
        ### Estimate new observable
        yest = np.sin(np.dot(B,xest))

        ### Get true y
        ytrue = yn[i-1:i+1]

        ### Update Kalman matrix
        temp = LA.inv(np.dot(B,np.dot(E,B.T))+Ny)
        K = np.dot(np.dot(E,B.T),temp) 

        ### estimate new state
        xestnew = xest + np.dot(K,(ytrue - yest))

        ###
        Eestnew = np.dot((np.eye(2) - np.dot(K,B)),E)

        #xtemp[:,0] = xest
        #xtemp[:,1:3] = xhist[:,0:2]
        #xhist = xtemp
        ###
        xest = np.dot(A,xestnew)
        #xest = xest + (xhist[:,1]-xhist[:,2])
        xhist[i] = xest[0]
        ###
        E = Eestnew + Nx
    return xhist,yn,yground
In [3]:
g = 10**-1
Nx = np.array([(g**2,0),(0,g**2)])
n = 700
xhist,yn,yground = kalman(n,Nx)
fig=plt.figure(figsize=(8, 5), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()

ax.plot(yground,linewidth = 2,color = 'k',label ='ground truth')
ax.plot(yn,linewidth = 2,color = 'blue',alpha = 0.5,label ='noisy data')
ax.plot(xhist,linewidth = 2,color = 'red',alpha = 0.5,label ='kalman')



handles, labels = ax.get_legend_handles_labels()
plt.legend(handles, labels)
plt.legend(loc=0)
plt.ylabel('Signal')
plt.xlabel('Time')
plt.title(' Noise = 10e-1')
Out[3]:
<matplotlib.text.Text at 0x11070def0>
In [4]:
g = 10**-3
Nx = np.array([(g**2,0),(0,g**2)])
n = 700
xhist,yn,yground = kalman(n,Nx)
fig=plt.figure(figsize=(8, 5), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()

ax.plot(yground,linewidth = 2,color = 'k',label ='ground truth')
ax.plot(yn,linewidth = 2,color = 'blue',alpha = 0.5,label ='noisy data')
ax.plot(xhist,linewidth = 2,color = 'red',alpha = 0.5,label ='kalman')



handles, labels = ax.get_legend_handles_labels()
plt.legend(handles, labels)
plt.legend(loc=0)
plt.ylabel('Signal')
plt.xlabel('Time')
plt.title(' Noise = 10e-3')
Out[4]:
<matplotlib.text.Text at 0x111136780>
In [5]:
g = 10**-5
Nx = np.array([(g**2,0),(0,g**2)])
n = 700
xhist,yn,yground = kalman(n,Nx)
fig=plt.figure(figsize=(8, 5), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()

ax.plot(yground,linewidth = 2,color = 'k',label ='ground truth')
ax.plot(yn,linewidth = 2,color = 'blue',alpha = 0.5,label ='noisy data')
ax.plot(xhist,linewidth = 2,color = 'red',alpha = 0.5,label ='kalman')



handles, labels = ax.get_legend_handles_labels()
plt.legend(handles, labels)
plt.legend(loc=0)
plt.ylabel('Signal')
plt.xlabel('Time')
plt.title(' Noise = 10e-5')
Out[5]:
<matplotlib.text.Text at 0x1118bcc50>