There’s an easier way to construct the Legendre polynomials using scipy.special.legendre
. You can use this to create a vector of polynomials (evaluated at the points you want to reconstruct), and then apply them with a dot product, like so:
def decoder(m, points_to_reconstruct=2):
output = np.zeros(number_of_points)
width = int(theta / dt)
window = np.linspace(0, 1, points_to_reconstruct)
D = np.asarray([legendre(i)(2 * window - 1) for i in range(q)])
return m.dot(D)
This decodes a two-dimensional array, where the first element is an approximation of u(t)
and the second is an approximation of u(t - theta)
. Increasing the value of points_to_reconstruct
gives you evenly spaced points in between. You can also replace np.linspace(0, 1, points_to_reconstruct)
with whatever points you want within the window (normalized to [0, 1]
).
However, one thing to note is that the reconstruction is approximating the window using polynomials up to order $q - 1$. If you have a step function, such as the one in your example, then it is impossible for a finite number of polynomials to fit such a curve, as it has discontinuous derivatives wherever there are discontinuous changes in the input (whereas, in contrast, polynomials have continuous derivatives). You do better if you increase $q$, though, in the same way that higher-order Taylor series form better approximations as you add more polynomials to the series.
In general, the reconstruction will work best for low-frequency inputs, where the frequency scales linearly with $q$. That is, if you double the frequency of your input signal, then you need twice as many polynomials to fit it as accurately.
Here’s an example taking your previous code, increasing q
to 24
, using ZOH discretization, and the above decoder
function:
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
from scipy.special import legendre
q = 24
theta = 0.4
number_of_points = 10000
dt = 1e-4
def get_a_and_b():
# A:
a = np.zeros([q, q])
for i in range(q):
for j in range(q):
if i < j:
a[i, j] = (2 * i + 1) * -1
else:
a[i, j] = (2 * i + 1) * ((-1) ** (i - j + 1))
# B:
b = np.zeros(q)
for i in range(q):
b[i] = (2 * i + 1) * ((-1) ** i)
return a/theta, b/theta
def encoder():
Acont, Bcont = get_a_and_b()
Adisc, Bdisc, _, _, _ = scipy.signal.cont2discrete(
(Acont, Bcont[:, None], np.zeros((1, q)), 0),
dt=dt,
method='zoh',
)
Bdisc = Bdisc.squeeze(axis=1)
m = np.zeros(q)
output = np.zeros([number_of_points, q])
for i in range(number_of_points):
if 1000 < i < 2000:
u = 1
else:
u = 0
# B_u = u * B
# m = m + dt * (np.dot(A, m) + B_u)
m = Adisc.dot(m) + Bdisc * u
output[i, :] = m
return output
def decoder(m, points_to_reconstruct=2):
output = np.zeros(number_of_points)
width = int(theta / dt)
window = np.linspace(0, 1, points_to_reconstruct)
D = np.asarray([legendre(i)(2 * window - 1) for i in range(q)])
return m.dot(D)
output = encoder()
output2 = decoder(output)
plt.plot(range(number_of_points), output2)
plt.show()
Note the ringing around the edges occurs from the higher-order polynomials trying to fit the instantaneous changes in the input. It would look much better with a filtered or low-frequency input.