0

I would like to find a fit for my data points which are periodic but non-trivial (they're not following a single sin, cos or tan function). These data points follow a periodic rectangle-ish/pulse-ish shape. The Fourier series should be able to help as it can be crafted to fit to any desired shape.

I tried using the Fourier series up to different degrees of precision (4/8/12 cosines) but it never really matched the desired shape.

I also tried working with linear combinations of an arbitrary amount of polynomial functions or using arbitrary Fourier series approximation but I didn't really understand these methods so it didn't work at all.

This is the code I'm using to plot my data and to fit it. (Here with a Fourier series fit with 4 cosines)

# Define the time scale of the plot
# In the end, the function should stretch from -5e-3 to 5e3 but the given points are all within the range t ϵ (0,4e-3)
t = np.linspace(-5e-3,5e-3,1000)

# I used estimation by eye to find points of a graph displayed on an oscilloscope
# These points represent a periodic function of pulses (kinda) or the voltage being charged and discharged
# Every 0th component of a tuple represents the t (time) coordinate, every 1st component the U (voltage) coordinate
coord1 = [(0,-0.3),(0.185e-3,-0.29775),(0.21e-3,-0.2785),(0.3e-3,-0.27),(0.4e-3,-0.2645),(0.5e-3,-0.26),
         (0.57e-3,-0.259),(0.595e-3,-0.27),(0.605e-3,-0.28),(0.623e-3,-0.3015),(0.64e-3,-0.316),(0.675e-3,-0.31),
         (0.7e-3,-0.305),(0.75e-3,-0.3),(0.8e-3,-0.298),(0.9e-3,-0.3),(1e-3,-0.3),(1.2e-3,-0.3),(1.24e-3,-0.279),
         (1.3e-3,-0.27),(1.4e-3,-0.266),(1.485e-3,-0.259),(1.58e-3,-0.2545),(1.62e-3,-0.269),(1.625e-3,-0.28),(1.65e-3,-0.29),
         (1.67e-3,-0.3125),(1.7e-3,-0.31),(1.745e-3,-0.301),(1.8e-3,-0.299),(1.985e-3,-0.2985),(2.21e-3,-0.297),
         (2.25e-3,-0.28),(2.3e-3,-0.2675),(2.4e-3,-0.26),(2.5e-3,-0.2575),(2.6e-3,-0.255),(2.65e-3,-0.26),
         (2.66e-3,-0.269),(2.665e-3,-0.279),(2.68e-3,-0.3),(2.72e-3,-pi/10),(2.75e-3,-0.305),(2.8e-3,-0.3),
         (3e-3,-0.3),(3.22e-3,-0.299),(3.3e-3,-0.298),(3.33e-3,-0.279),(3.4e-3,-0.268),(3.45e-3,-0.265),(3.5e-3,-0.262),
         (3.55e-3,-0.26),(3.62e-3,-0.2575),(3.73e-3,-0.2585),(3.745e-3,-0.269),(3.75e-3,-0.28),(3.755e-3,-0.3),
         (3.76e-3,-0.314),(3.795e-3,-0.309),(3.85e-3,-0.305),(3.92e-3,-0.3),(4e-3,-0.295)]

# Extract t and U coordinates separately
t1 = np.array([point[0] for point in coord1])
U1 = np.array([point[1] for point in coord1])

# I tried using the Fourier series. I tried using 4, 8 or 12 components but it never really matched the given shape
# Here, a Fourier series of 4 cosine components is used
def fou4(x,tau,a1,a2,a3,a4,U0):
    return a1*np.cos(np.pi/tau *x) + a2*np.cos(2*np.pi/tau *x) + a3*np.cos(3*np.pi/tau *x) + a4*np.cos(4*np.pi/tau *x) + U0

# The U0 shift is needed for the correct positioning of the function along the y axis

# Use curve_fit to determine optimal parameters for fitting the Fourier series to the given data set
para1,pcov1 = curve_fit(fou4,t1,U1,p0=[0.0005,1,1,1,1,-0.3])

# Plotting
# The plot should contain the given data points and a periodic fit to them
fig,ax = plt.subplots(figsize=(10,10))
# Plot data points
ax.plot(t1,U1,marker="x",ls="--",alpha=0.5,markersize=6,color="blue",label="Channel 2 points, eye estimation")
# Plot fit
ax.plot(t,fou4(t,*para1),color="#00DDDD",alpha=0.85,ls="-",label=r"Fit for Ch.2")

# The x and y limits are set to zoom in on the given data points so that one can easily compare the fit with the real data
ax.set(title="Plot for image of 2023/02/09 15:24",xlabel="Time [s]",ylabel="Voltage [V]",xlim=(-0.0005,0.005),ylim=(-0.35,-0.25))
ax.grid()
ax.legend()

The result is a fit shape that does not match the points and that especially completely ignores the downward peaks below -0.3.

I thought that maybe using only 4 cosines is just not precise enough but using 8 or 12 cosines the fit becomes even worse. It might work for better p0 and bounds but I can't figure these out.

my output

Jay
  • 9
  • 1
  • Does this post answer your question? https://stackoverflow.com/questions/32590720/create-scipy-curve-fitting-definitions-for-fourier-series-dynamically?rq=2 – Suraj Shourie Aug 07 '23 at 15:29
  • I've tried these methods already and either it didn't work for my specific problem or I just did not do it right. – Jay Aug 07 '23 at 20:30
  • This looks like charging, switch-off with oscillation. So maybe a proper model is available. – mikuszefski Sep 01 '23 at 13:56

1 Answers1

0

I would like to suggest using a Gaussian Process (GP). I have recently answered a question using GPs, take a look here to understand better what they are and how they work. In that question I have used a noise-free GP, but for this case we definetely need to model the noise.

Since your data show a periodic behavior, an Exp-Sine-Squared kernel is likely to be what you need.

Here is the code:

import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ExpSineSquared

coordinates = [
    (0, -0.3),
    (0.185e-3, -0.29775),
    (0.21e-3, -0.2785),
    (0.3e-3, -0.27),
    (0.4e-3, -0.2645),
    (0.5e-3, -0.26),
    (0.57e-3, -0.259),
    (0.595e-3, -0.27),
    (0.605e-3, -0.28),
    (0.623e-3, -0.3015),
    (0.64e-3, -0.316),
    (0.675e-3, -0.31),
    (0.7e-3, -0.305),
    (0.75e-3, -0.3),
    (0.8e-3, -0.298),
    (0.9e-3, -0.3),
    (1e-3, -0.3),
    (1.2e-3, -0.3),
    (1.24e-3, -0.279),
    (1.3e-3, -0.27),
    (1.4e-3, -0.266),
    (1.485e-3, -0.259),
    (1.58e-3, -0.2545),
    (1.62e-3, -0.269),
    (1.625e-3, -0.28),
    (1.65e-3, -0.29),
    (1.67e-3, -0.3125),
    (1.7e-3, -0.31),
    (1.745e-3, -0.301),
    (1.8e-3, -0.299),
    (1.985e-3, -0.2985),
    (2.21e-3, -0.297),
    (2.25e-3, -0.28),
    (2.3e-3, -0.2675),
    (2.4e-3, -0.26),
    (2.5e-3, -0.2575),
    (2.6e-3, -0.255),
    (2.65e-3, -0.26),
    (2.66e-3, -0.269),
    (2.665e-3, -0.279),
    (2.68e-3, -0.3),
    (2.72e-3, -np.pi / 10),
    (2.75e-3, -0.305),
    (2.8e-3, -0.3),
    (3e-3, -0.3),
    (3.22e-3, -0.299),
    (3.3e-3, -0.298),
    (3.33e-3, -0.279),
    (3.4e-3, -0.268),
    (3.45e-3, -0.265),
    (3.5e-3, -0.262),
    (3.55e-3, -0.26),
    (3.62e-3, -0.2575),
    (3.73e-3, -0.2585),
    (3.745e-3, -0.269),
    (3.75e-3, -0.28),
    (3.755e-3, -0.3),
    (3.76e-3, -0.314),
    (3.795e-3, -0.309),
    (3.85e-3, -0.305),
    (3.92e-3, -0.3),
    (4e-3, -0.295),
]


def main() -> None:
    t_ = np.array([point[0] for point in coordinates])
    u_ = np.array([point[1] for point in coordinates])

    t = np.array(t_).reshape(-1, 1)
    u = np.array(u_).reshape(-1, 1)

    gp_regressor = train_gaussian_process(t, u)
    x = np.linspace(0, 4e-3, 400).reshape(-1, 1)
    mean = gp_regressor.predict(x)

    plt.figure(figsize=(14, 6))
    plt.scatter(t, u, label="Observations")
    plt.plot(x, mean, linewidth=1.2)
    plt.show()


def train_gaussian_process(x: np.ndarray, y: np.ndarray) -> GaussianProcessRegressor:
    """
    Train Gaussian Process.

    Parameters
    ----------
    x : np.ndarray
    y : np.ndarray

    Returns
    -------
    GaussianProcessRegressor

    """
    kernel = ExpSineSquared(length_scale=1.0e-4, periodicity=1.0e-3)
    gaussian_process = GaussianProcessRegressor(
        kernel=kernel, alpha=1e-4, n_restarts_optimizer=9, random_state=42
    )
    gaussian_process.fit(x, y)
    return gaussian_process


if __name__ == "__main__":
    main()

and here is the plot: enter image description here

as you can see we are now getting a better fit and the model is also able to capture the downward peak. Fine-tune the hyperparameters to get exactly what you need.

blunova
  • 2,122
  • 3
  • 9
  • 21