I would like to fit a periodic rectangle function to my data. I have tried this with the least squares method using curve_fit from scipy.optimize. Unfortunately, the optimization does not converge.
Here is a minimal example with an excerpt of the data:
import numpy as np
from scipy import signal
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
x = np.array([0.00000e+00, 6.97000e-05, 1.38200e-04, 2.26700e-04, 3.09700e-04,
4.38200e-04, 5.05000e-04, 5.66700e-04, 6.33400e-04, 8.41900e-04,
9.43200e-04, 1.03100e-03, 1.18000e-03, 1.26450e-03, 1.34990e-03,
2.11920e-03, 2.59670e-03, 4.61800e-03, 4.99820e-03, 5.06740e-03,
5.14820e-03, 5.22640e-03, 5.31120e-03, 5.39800e-03, 5.46740e-03,
5.55040e-03, 5.63120e-03, 5.72520e-03, 5.83990e-03, 5.91070e-03,
6.00040e-03, 6.06850e-03, 6.17600e-03, 6.26270e-03, 6.32870e-03,
6.50600e-03, 6.74950e-03, 6.85150e-03, 6.93870e-03, 7.06570e-03,
7.12740e-03, 7.44270e-03, 7.59450e-03, 8.13300e-03, 9.61600e-03,
9.99620e-03, 1.00592e-02, 1.01345e-02, 1.02170e-02, 1.03185e-02,
1.03862e-02, 1.04760e-02, 1.05484e-02, 1.06410e-02, 1.07402e-02,
1.08027e-02, 1.09087e-02, 1.09982e-02, 1.11092e-02, 1.11762e-02,
1.12610e-02, 1.13267e-02, 1.14489e-02, 1.16935e-02, 1.17800e-02,
1.21155e-02, 1.23132e-02, 1.24477e-02, 1.25500e-02, 1.37562e-02,
1.41364e-02, 1.46140e-02])
y = np.array([-14.6564, -14.6564, -14.6564, -14.6795, -14.6795, -14.6938,
-14.6938, -14.7025, -14.7025, -14.714 , -14.714 , -14.7226,
-14.7284, -14.7284, -14.7284, -14.7571, 14.8909, 14.9701,
-14.6751, -14.6751, -14.6751, -14.691 , -14.691 , -14.7082,
-14.7082, -14.7082, -14.7169, -14.7169, -14.7255, -14.7255,
-14.7341, -14.7341, -14.7399, -14.7399, -14.7399, -14.7485,
-14.7528, -14.76 , -14.76 , -14.7687, -14.7687, -14.7787,
14.8837, 14.9096, 14.96 , -14.6766, -14.6766, -14.6766,
-14.6982, -14.6982, -14.7111, -14.7111, -14.7111, -14.7212,
-14.7212, -14.7284, -14.7284, -14.737 , -14.737 , -14.7428,
-14.7428, -14.7428, -14.7557, -14.7571, -14.76 , -14.7687,
-14.7715, -14.7787, -14.7787, 14.9327, 14.9442, 14.96 ])
########## define rectangle function via heaviside step function ############
def rectangle_1(x,Amp,f,phi):
return 2*Amp*np.heaviside(np.cos(2*np.pi*f*x+phi),1)-Amp
########## define rectangle function via signal package from scipy ###########
def rectangle_2(x,Amp,f,offset):
return Amp*signal.square(2*np.pi*f*(x-offset))
########## Fitting of the data by using curvefit from scipy.optimize ############
popt, pcov = curve_fit(rectangle_1, x, y, p0=[15,200,np.pi/2],bounds=([0,0,0],[50,300,10]))
################## plot data points and rectangle function ###############
x_sig = np.linspace(min(x),max(x),10000)
y_sig = rectangle_1(x_sig,15,200,np.pi/2) # for start parameters
y_sig = rectangle_1(x_sig,*popt) # for fit parameters
plt.figure()
plt.plot(x,y,color='b', marker='.',markerfacecolor='r',markeredgecolor='r',linewidth=0,label=r'detected Pulses')
plt.plot(x_sig,y_sig,color='b',linewidth=1,label='applied voltage')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='upper center', bbox_to_anchor= (0.5, 1.13), ncol=2, frameon=False, fontsize=12)
I have tried various definitions of the rectangle function, as well as fitting it using ODR.run(). Nothing worked.
Another method for determining the function parameters such as edge detection does not work here because the data points are not necessarily located at the corners.
Here below the results visualized.
Data and rectangular function with start parameters:
Data and rectangular function with fit parameters:
I'm really surprised that I haven't found a solution to this online yet. I mean rectangular signals are quite popular.