# Fitting a rectangle-function to a signal in Python

I have a measured signal (current of a motor, turning on and off again) to which I want to fit a rectangular function in python. I came up with a reasonable rect(...) function and the scipy.optimize.curve_fit is running, but the results are far from perfect (see code below).

import numpy as np
from scipy.optimize import curve_fit

def rect(x ,a=1, b=1, x0=0, y0=0):
return a*np.where(np.abs(np.asarray(x)-x0)<=b/2, 1., 0.)+y0

# time, current
t,y = np.array([-0.20, -0.061726,
-0.15, -0.061358,
-0.10, -0.061036,
-0.05, -0.060445,
0.00, -0.044412,
0.05,  0.257997,
0.10,  0.277212,
0.15,  0.287594,
0.20,  0.257026,
0.25,  0.269109,
0.30,  0.278451,
0.35,  0.262615,
0.40,  0.291011,
0.45,  0.269262,
0.50,  0.271740,
0.55,  0.269469,
0.60,  0.278860,
0.65,  0.280009,
0.70,  0.309264,
0.75,  0.307818,
0.80,  0.311503,
0.85,  0.311120,
0.90,  0.319188,
0.95,  0.306017,
1.00,  0.324442,
1.05,  0.308988,
1.10,  0.287356,
1.15,  0.306467,
1.20,  0.830262,
1.25, -0.061024,
1.30, -0.061441,
1.35, -0.061425,
1.40, -0.061275,
1.45, -0.060958,
1.50, -0.060957,
1.55, -0.060451,
1.60, -0.061107,
1.65, -0.060990,
1.70, -0.061227,
1.75, -0.061214,
1.80, -0.061618]).reshape((-1,2)).T

p0 = np.array([.5, .5, .5, 0.1]) # initial point
bounds =([-1., .1, -.2, -1.],
[ 1., 2., 1.5,  1.])
popt, pcov = curve_fit(rect,t,y,p0=p0,bounds=bounds,check_finite=True)

pmanual = np.array([.35, 1.2, .65, -.06]) # manual set of parameters, which is close to optimal
# plot
import matplotlib.pyplot as plt
plt.figure(1); plt.clf()
plt.plot(t,y,label='original data')
plt.plot(t,rect(t,*p0),label='initial point')
plt.plot(t,rect(t,*popt),label='curve_fit result')
plt.plot(t,rect(t,*pmanual),label='expected result',ls='--')
plt.legend()


The curve_fit does not seem to converge to a reasonable set of parameters. Changing the values in p0 gives different results but none of these seem to be close to the optimum. How can I improve the code/data to robustly find a good set of parameters?

• Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
– Community Bot
Aug 23 at 19:27
• Welcome to Scicomp! Maybe you don't need to fit. Sweep the array from front to back with a reasonable threshold, use the outermost threshold passes as the bounds for your rect, use the average of points in between for your amplitude. Aug 24 at 6:17

Using the following script, attached to your data, the result obtained is plotted in cyan

from scipy.optimize import minimize
import matplotlib.pyplot as plt

def bump(x,x1,h0,dh,dx,a):
if x <= x1:
return h0
x2 = x1 + dx
if x > x1 and x <= x2:
return h0 + (dh/dx)*(x - x1)
x3 = x2 + a
if x > x2 and x <= x3:
return h0 + dh
x4 = x3 + dx
if x > x3 and x <= x4:
return h0 + dh - (dh/dx)*(x - x3)
if x > x4:
return h0

def objective(x):
(x1,h0,dh,dx,a) = x
error = 0
for k in range(len(y)):
error += abs(y[k] - bump(t[k],x1,h0,dh,dx,a))
return error

x0 = [0.,0.,1.,0.1,1.1]
bds = ((-0.5, None),(-0.5, None),(0, None),(0, None),(0.8, None))

res = minimize(objective, x0, method = 'SLSQP', bounds=bds)
(x1,h0,dh,dx,a) = list(res.x)
print(res.fun)
print('h0 = ', h0,' dh = ', dh,' dx = ',dx,' a = ', a,' x1 = ', x1)

val = []
for k in range(len(y)):
val.append(bump(t[k],x1,h0,dh,dx,a))

plt.plot(t,val)


• The code works as expected. Could you please explain a bit, what's the idea behind the redefined rect-function, the objective function and the chosen optimizer?
– n4pK
Aug 30 at 18:29
• First, it appears that a trapezoid shape with rise and fall flanks is more appropriate. Second, the curve-fitting procedure uses as well a minimization process which is unknown to the user. I prefer an explicit constrained optimization method choice. Aug 30 at 21:05