fitv2.py

import numpy as np
from scipy.optimize import minimize
from scipy.signal import lfilter, butter

Fs = 50
Ts = 1 / Fs  # sec
tmax = 10  # sec
noise_level = 0.05
system_order = 1

# time
t = np.arange(0, tmax, Ts)

# Input data
input = (t > tmax / 10).astype(np.uint)

# The system is an exponential filter
system = np.array(butter(1, 0.2 / (Fs / 2)))
y = lfilter(*system, input)
noise = np.random.normal(0, noise_level, y.shape)
y_noisy = y + noise

# Fit the system
x0 = np.zeros(system.shape)
x0[:, 0] = 1  # the starting system is output=input
# The cost function
cost = lambda x: np.linalg.norm(lfilter(*x.reshape(system.shape), input) - y_noisy)
found_system = minimize(cost, x0.ravel()).x.reshape(system.shape)

# Plot fitting
import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 1)
ax.plot(t, input, label="Input")
ax.plot(t, y_noisy, label="System response")
ax.plot(t, lfilter(*found_system, input), label="Found system")
ax.legend()
ax.grid()
ax.set_xlabel("time")
ax.set_ylabel("amplitude")
plt.show()