Cublic Spline Interpolation of Phase Space Plot

I am creating a phase-space plot of first derivative of voltage against voltage:

I want to interpolate the plot so so it is smooth. So far, I have approached this by interpolating the voltage and first derivative of the voltage separately with time, then generating phase space plots.

Python Code (toy data example)

import numpy as np
import scipy.interpolate

interp_factor = 100
n = 12
time = np.linspace(0, 10, n)
voltage = np.array([0, 1, 2, 10, 30, 70, 140, 150, 140, 80, 40, 10])

voltage_diff = np.diff(voltage)
voltage = voltage[:-1]
time = time[:-1]

interp_function_voltage = scipy.interpolate.interp1d(time, voltage, kind="cubic")
interp_function_voltage_diff = scipy.interpolate.interp1d(time, voltage_diff, kind="cubic")

new_sample_num = interp_factor * (n - 1) + 1
new_time = np.linspace(np.min(time), np.max(time), new_sample_num)

interp_voltage = interp_function_voltage(new_time)
interp_voltage_diff = interp_function_voltage_diff(new_time)

I would like to ask:

a) is the method as implemented reasonable?

b) is there a better method by interpolating directly in the phase-space? e.g. interpolating with voltage as x and voltage_diff as y? I do not think this makes sense, because the voltage values are not uniformly spaced and there may be repeated voltage values. I also tried the scipy parametric interpolation methods (e.g. scipy.interpolate.splprep) but these threw input value error. I expect (it would be nice to have this clarified) because this is raw data, rather than well behaved parametric functions.

I guess more generally, I am wondering if it makes sense to somehow do the interpolation in the phase-space to make use of the direct relationship between voltage and voltage_diff for interpolating / smoothing.

Many thanks

It is reasonable, but your difference will be biased, maybe the best approximation for the difference could be (v[i+1] - v[i-1])/(2*dt)

Another approach is using Fourier transform smoothing

def smoother_phase_space(y, sps=1, T=1):
    Y = np.fft.rfft(y)
    yu = np.fft.irfft(Y, len(y)*sps).real * sps
    dyu = np.fft.irfft(Y * (2j * np.pi * np.fft.rfftfreq(len(y))), len(y)*sps).real
    k = np.arange(len(yu)+2) % len(yu)
    return yu[k], dyu[k] * sps / T

v, dv = smoother_phase_space(voltage, sps=1)
plt.plot(v, dv, '-ob')

v, dv = smoother_phase_space(voltage, sps=4)
plt.plot(v, dv, '-r')
plt.plot(v[::4], dv[::4], 'or')
v, dv = smoother_phase_space(voltage, sps=32)
plt.plot(v, dv, '-g')
plt.plot(v[::32], dv[::32], 'og')

try: # the data computed in the original post
    plt.plot(interp_voltage, interp_voltage_diff, '--')
except:
    pass