Improve an exponential curve fit of nearly exponential data

Question

I have recorded some data about the color of an LED that varies with the 8bit signal sent to the LED driver, the signal can vary between 0 and 255.

Exponential curve fitting seems to work very well to represent the LED's behavior. I have had good results with the following formula:

x * signal ** ex 
y * signal ** ey 
z * signal ** ez

In Python, I use the following function:

from scipy.optimize import curve_fit

def fit_func_xae(x, a, e):
    # Curve fitting function
    return a * x**e

# X, Y, Z are real colorimetric values that are measured by a physical instrument

(aX, eX), cov = curve_fit(fit_func3xa, signal, X)
(aY, eY), cov = curve_fit(fit_func3xa, signal, Y)
(aZ, eZ), cov = curve_fit(fit_func3xa, signal, Z)

Note:

In colorimetry, we represent the color of the LED in the CIE XYZ color space, which is a linear space that works in a similar way as a linear RGB color space. Even if it is an aproximation, you can think of XYZ as a synonym of (linear) RGB. So a color can be represented as a triplet of linear values X, Y, Z.

here is the data behind the curves. for each 8bit parameter sent to the LED driver, there are 3 measures

Signal 
[  3.   3.   3.   5.   5.   5.   7.   7.   7.  10.  10.  10.  15.  15.
  15.  20.  20.  20.  30.  30.  30.  40.  40.  40.  80.  80.  80. 160.
 160. 160. 240. 240. 240. 255. 255. 255.]

X, Y, Z
[[9.93295448e-05 8.88955748e-04 6.34978556e-04]
 [9.66399391e-05 8.86031926e-04 6.24680520e-04]
 [1.06108685e-04 8.99010175e-04 6.41577838e-04]
 [1.96407653e-04 1.70210146e-03 1.27178991e-03]
 [1.84965943e-04 1.67927596e-03 1.24985475e-03]
 [1.83770476e-04 1.67905297e-03 1.24855580e-03]
 [3.28537613e-04 2.75382195e-03 2.14639821e-03]
 [3.17804246e-04 2.74152647e-03 2.11730825e-03]
 [3.19167905e-04 2.74977632e-03 2.11142769e-03]
 [5.43770342e-04 4.09314433e-03 3.33793380e-03]
 [5.02493149e-04 4.04392581e-03 3.24784452e-03]
 [5.00712102e-04 4.03456071e-03 3.26803716e-03]
 [1.48001671e-03 1.09367632e-02 9.59283037e-03]
 [1.52082180e-03 1.09920985e-02 9.63624777e-03]
 [1.50153844e-03 1.09623592e-02 9.61724422e-03]
 [3.66206564e-03 2.74730946e-02 2.51982924e-02]
 [3.64074861e-03 2.74283157e-02 2.52187294e-02]
 [3.68719991e-03 2.75033778e-02 2.51691331e-02]
 [1.50905917e-02 1.06056566e-01 1.06534373e-01]
 [1.51370269e-02 1.06091182e-01 1.06790424e-01]
 [1.51654172e-02 1.06109863e-01 1.06943957e-01]
 [3.42912601e-02 2.30854413e-01 2.43427207e-01]
 [3.42217124e-02 2.30565972e-01 2.43529454e-01]
 [3.41486993e-02 2.30807320e-01 2.43591644e-01]
 [1.95905112e-01 1.27409867e+00 1.37490536e+00]
 [1.94923951e-01 1.26934278e+00 1.37751808e+00]
 [1.95242984e-01 1.26805844e+00 1.37565458e+00]
 [1.07931878e+00 6.97822521e+00 7.49602715e+00]
 [1.08944832e+00 7.03128378e+00 7.54296884e+00]
 [1.07994964e+00 6.96864302e+00 7.44011991e+00]
 [2.95296087e+00 1.90746191e+01 1.99164655e+01]
 [2.94254973e+00 1.89524517e+01 1.98158118e+01]
 [2.95753358e+00 1.90200667e+01 1.98885050e+01]
 [3.44049055e+00 2.21221159e+01 2.29667049e+01]
 [3.43817829e+00 2.21225393e+01 2.29363833e+01]
 [3.43077583e+00 2.21158929e+01 2.29399652e+01]]

_

The Problem:

Here's a scatter plot of some of my LED's XYZ values, together with a plot of the exponential curve fitting obtained with the code above:

It all seems good... until you zoom a bit:

On this zoom we can also see that the curve is fitted on multiple measurements:

At high values, Z values (blue dots) are always higher than Y values (green dots). But at low values, Y values are higher than Z values.

The meaning of this is that the LED changes in color depending on the PWM applied, for some reason (maybe because the temperature rises when more power is applied).

This behavior cannot be represented mathematically by the formula that I have used for the curve fit, however, the curve fit is so good for high values that I am searching for a way to improve it in a simple and elegant way.

Do you have any idea how this could be done? I have tried unsuccessfully to add more parameters, for example I tried to use:

x1 * signal ** ex + x2 * signal ** fx

instead of:

x * signal ** ex

but that causes scipy to overflow.

My idea was that by adding two such elements I could still have a funtion that equals 0 when signal = 0, but that increases faster at low values than a simple exponential.

Do you get a better fit if you exclude the highest values from the raw data? — Vin, Nov 15 '22 at 04:14
measurements without any error bars are meaningless. So get your error bars, then add them to your fit. Also are you aware of overfitting concept? Getting your fit to pass through all points by just tampering with your model without any sound theory to back it up is generally a bad idea. — Julien, Nov 15 '22 at 04:19
Maybe fitting in the log-space is better here. If it is linear for large values but drops faster for smaller ones you may test something like `x^( p+1 ) / ( x0 + x^p )`. Posting some data would be nice. — mikuszefski, Nov 15 '22 at 07:10
@Julien, sounds like a good idea, but how do I add error bars to the fit? Isn't fitting on multiple measurements of the same values enough? I edited the question to reflect this — adrienlucca.net, Nov 15 '22 at 23:03
Hi, thanks for the data. Not so sure about by formula now, but you get a good idea of your data on a log-log-scale — mikuszefski, Nov 16 '22 at 06:46
If the data cannot be fitted properly with a classic GOG (Gain-Offset-Gamma) model, the next step is to create a 3x1D LUT or even a 3D LUT. This is typically what is used in hardware and software calibration. — Kel Solaar, Nov 17 '22 at 17:53
@KelSolaar Is there a way to linearize a system to go from a set of CIE XYZ coordinates to a PWM and vice-versa using a LUT? How? — adrienlucca.net, Nov 18 '22 at 00:21
@Vin not really, that just translates the problem, when low values are better, high values are bad — adrienlucca.net, Nov 18 '22 at 03:04
@adrienlucca.net: This would be a good start: https://www.researchgate.net/publication/225223121_An_Accurate_Characterization_of_CRT_Monitor_II_Proposal_for_an_Extension_to_CIE_Method_and_Its_Verification The GOG model is very simple, e.g. f(x) = (ax + b)^c — Kel Solaar, Nov 19 '22 at 21:57
@KelSolaar thanks, (nice to talk to a color expert by the way!) however the GOG model works under the assumption that "A property of all luminescent material is that their normalized spectral radiance is invariant.", which is only true for phosphor-converted LEDs. In my case, the LED is monochromatic thus the ratio X/Y and Z/Y change when Y increases... any idea how to take this into account? — adrienlucca.net, Nov 20 '22 at 15:11
Do you mean that that its chromaticity coordinates are changing as a function of the PWM? — Kel Solaar, Nov 22 '22 at 08:06
@KelSolaar sorry I did not see your last question. YES, this is exactly what happens. Chromaticity changes at low PWMs. IMO it is due to a change in temperature inside the LED. See for example: https://www.researchgate.net/publication/267266688_Temperature-Dependent_Internal_Quantum_Efficiency_of_Blue_High-Brightness_Light-Emitting_Diodes/download — adrienlucca.net, Dec 10 '22 at 23:44
@KelSolaar actually this reference is not the best because 450nm blue LEDs do not change that much. The change is more spectacular with amber, orange, red and photo-red LEDs. For example, a photo-red LED can appear less bright at 1A than at 350mA because the peak wavelength moves towards IR, thus it emits more energy but this energy is less visible. — adrienlucca.net, Dec 10 '22 at 23:49

mikuszefski · Answer 1 · 2022-11-18T09:32:24.003

The data shows two steps in the log-log plot so I used an approach already used here.

Code is as follows:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
signal = np.array( [
    3.0,   3.0,   3.0,
    5.0,   5.0,   5.0,
    7.0,   7.0,   7.0,
    10.0,  10., 10.,
    15.0,  15., 15.,
    20.0,  20.,  20.,
    30.0,  30.,  30.,
    40.0,  40.,  40.,
    80.0,  80., 80.,
    160.0, 160., 160.,
    240.0, 240., 240.,
    255.0, 255., 255.
] )

data = np.array( [
    [9.93295448e-05, 8.88955748e-04, 6.34978556e-04],
    [9.66399391e-05, 8.86031926e-04, 6.24680520e-04],
    [1.06108685e-04, 8.99010175e-04, 6.41577838e-04],
    [1.96407653e-04, 1.70210146e-03, 1.27178991e-03],
    [1.84965943e-04, 1.67927596e-03, 1.24985475e-03],
    [1.83770476e-04, 1.67905297e-03, 1.24855580e-03],
    [3.28537613e-04, 2.75382195e-03, 2.14639821e-03],
    [3.17804246e-04, 2.74152647e-03, 2.11730825e-03],
    [3.19167905e-04, 2.74977632e-03, 2.11142769e-03],
    [5.43770342e-04, 4.09314433e-03, 3.33793380e-03],
    [5.02493149e-04, 4.04392581e-03, 3.24784452e-03],
    [5.00712102e-04, 4.03456071e-03, 3.26803716e-03],
    [1.48001671e-03, 1.09367632e-02, 9.59283037e-03],
    [1.52082180e-03, 1.09920985e-02, 9.63624777e-03],
    [1.50153844e-03, 1.09623592e-02, 9.61724422e-03],
    [3.66206564e-03, 2.74730946e-02, 2.51982924e-02],
    [3.64074861e-03, 2.74283157e-02, 2.52187294e-02],
    [3.68719991e-03, 2.75033778e-02, 2.51691331e-02],
    [1.50905917e-02, 1.06056566e-01, 1.06534373e-01],
    [1.51370269e-02, 1.06091182e-01, 1.06790424e-01],
    [1.51654172e-02, 1.06109863e-01, 1.06943957e-01],
    [3.42912601e-02, 2.30854413e-01, 2.43427207e-01],
    [3.42217124e-02, 2.30565972e-01, 2.43529454e-01],
    [3.41486993e-02, 2.30807320e-01, 2.43591644e-01],
    [1.95905112e-01, 1.27409867e+00, 1.37490536e+00],
    [1.94923951e-01, 1.26934278e+00, 1.37751808e+00],
    [1.95242984e-01, 1.26805844e+00, 1.37565458e+00],
    [1.07931878e+00, 6.97822521e+00, 7.49602715e+00],
    [1.08944832e+00, 7.03128378e+00, 7.54296884e+00],
    [1.07994964e+00, 6.96864302e+00, 7.44011991e+00],
    [2.95296087e+00, 1.90746191e+01, 1.99164655e+01],
    [2.94254973e+00, 1.89524517e+01, 1.98158118e+01],
    [2.95753358e+00, 1.90200667e+01, 1.98885050e+01],
    [3.44049055e+00, 2.21221159e+01, 2.29667049e+01],
    [3.43817829e+00, 2.21225393e+01, 2.29363833e+01],
    [3.43077583e+00, 2.21158929e+01, 2.29399652e+01]
] )

def determine_start_parameters(  x , y, edge=9 ):
    logx = np.log( x )
    logy = np.log( y )
    xx = logx[ :edge ]
    yy = logy[ :edge ]
    (ar1, br1), _ = curve_fit( lambda x, slope, off: slope * x + off, xx , yy )
    xx = logx[ edge : -edge ]
    yy = logy[ edge : -edge]
    (ar2, br2), _ = curve_fit( lambda x, slope, off: slope * x + off, xx , yy )
    xx = logx[ -edge : ]
    yy = logy[ -edge : ]
    (ar3, br3), _ = curve_fit( lambda x, slope, off: slope * x + off, xx , yy )
    cross1r = ( br2 - br1 ) / ( ar1 - ar2 )
    mr = ar1 * cross1r + br1
    cross2r = ( br3 - br2 ) / ( ar2 - ar3 )
    xx0r = [ mr, ar1, ar2 , ar3, cross1r, cross2r, 1 ]
    return xx0r

def func(
    x, b,
    m1, m2, m3,
    a1, a2,
    p
):
    """
    continuous approxiation for a two-step function
    used to fit the log-log data
    p is a sharpness parameter for the transition
    """
    out = b - np.log(
    ( 1 + np.exp( -m1 * ( x - a1 ) )**abs( p ) )
    ) / p + np.log(
    ( 1 + np.exp( m2 * ( x - a1 ) )**abs( p ) )
    ) / p - np.log(
    ( 1 + np.exp( m3 * ( x - a2 ) )**abs( p ) )
    ) / abs( p )
    return out

def expfunc(
    x, b,
    m1, m2, m3,
    a1, a2,
    p
):
    """
    remapping to the original data
    """
    xi = np.log( x )
    eta = func(
    xi, b,
    m1, m2, m3,
    a1, a2,
    p)
    return np.exp(eta)

def expfunc2(
    x, b,
    m1, m2, m3,
    a1, a2,
    p
):
    """
    simplified remapping
    """
    aa1 = np.exp( a1 )
    aa2 = np.exp( a2 )
    return (
        np.exp( b ) * (
        ( 1 + ( x / aa1 )**( m2 * p ) ) / 
        ( 1 + ( x / aa2 )**( m3 * p ) ) /
        ( 1 + ( aa1 / x )**( m1 * p ) )
        )**( 1 / p )
    )

logsig = np.log( signal )
logred =  np.log( data[:,0] )
loggreen =  np.log( data[:,1] )
logblue =  np.log( data[:,2] )

### getting initial parameters
### red

xx0r = determine_start_parameters( signal, data[ :, 0 ] )
xx0g = determine_start_parameters( signal, data[ :, 1 ] )
xx0b = determine_start_parameters( signal, data[ :, 2 ] )
print( xx0r )
print( xx0g )
print( xx0b )

xl = np.linspace( 1, 6, 150 )
tl = np.linspace( 1, 260, 150 )


solred = curve_fit( func, logsig, logred, p0=xx0r )[0]
solgreen = curve_fit( func, logsig, loggreen, p0=xx0g )[0]
solblue = curve_fit( func, logsig, logblue, p0=xx0b )[0]

print( solred )
print( solgreen )
print( solblue )

fig = plt.figure()
ax = fig.add_subplot( 2, 1, 1 )
bx = fig.add_subplot( 2, 1, 2 )

ax.scatter(  np.log( signal ), np.log( data[:,0] ),  color = 'r' )
ax.scatter(  np.log( signal ), np.log( data[:,1] ),  color = 'g' )
ax.scatter(  np.log( signal ), np.log( data[:,2] ),  color = 'b' )

ax.plot( xl, func( xl, *solred ),  color = 'r' )
ax.plot( xl, func( xl, *solgreen ),  color = 'g' )
ax.plot( xl, func( xl, *solblue ),  color = 'b' )

bx.scatter( signal, data[:,0],  color = 'r' )
bx.scatter( signal, data[:,1],  color = 'g' )
bx.scatter( signal, data[:,2],  color = 'b' )
bx.plot( tl, expfunc2( tl, *solred), color = 'r'  )
bx.plot( tl, expfunc2( tl, *solgreen), color = 'g'  )
bx.plot( tl, expfunc2( tl, *solblue), color = 'b'  )


plt.show()

Which results in

How did you find the values: "xx0=[ -6, 1, 2.5,1.5, 2.5, 4, 10]" ? — adrienlucca.net, Nov 17 '22 at 22:15
your solution seems to work with this data despite many warnings such as "RuntimeWarning: overflow encountered in exp out = b - np.log((1 + np.exp(-m1 * (x - a1)) ** p)) / p + np.log((1. + np.exp(m2 *(x - a1)) ** p)) / p - np.log((1. + np.exp(m3 * (x - a2)) ** p)) / p". However, it does not work with other data from other LEDs :/ Maybe if I understood the way you created the xx0 parameters... I could figure out a solution — adrienlucca.net, Nov 17 '22 at 22:36
@adrienlucca.net hi, you may remove the warnings by changing `p` to `abs( p )` — mikuszefski, Nov 18 '22 at 07:46
@adrienlucca.net concerning the start parameters, they are actually found manually. You may automatize by: (1) double-log the data. (2) check initial slope from first two or three points. Same for final slope. Intermediate slope maybe via slope from first to last. (3) Offset is near third `y` value. (4) `a1,2` are somewhere in the middle. (5) `p` probably smaller to start softer. Maybe `p=1` — mikuszefski, Nov 18 '22 at 08:01
@adrienlucca.net edited the code to calc the initial values. Maybe this works — mikuszefski, Nov 18 '22 at 09:33

Improve an exponential curve fit of nearly exponential data

1 Answers1