Curve fitting with nth order polynomial having sine ripples

Question

I'm modeling measurement errors in a certain measuring device. This is how the data looks: high frequency sine ripples on a low frequency polynomial. My model should capture the ripples too.

The curve that fits the error should be of the form: error(x) = a0 + a1*x + a2*x^2 + ... an*x^n + Asin(x/lambda). The order n of the polynomial is not known. My plan is to iterate n from 1-9 and select the one that has the highest F-value.

I've played with numpy.polyfit and scipy.optimize.curve_fit so far. numpy.polyfit is only for polynomials, so while I can generate the "best fit" polynomial, there's no way to determine the parameters A and lambda for the sine term. scipy.optimize.curve_fit would have worked great if I already knew the order of the polynomial for the polynomial part of error(x).

Is there a clever way to use both numpy.polyfit and scipy.optimize.curve_fit to get this done? Or another library-function perhaps?

Here's the code for how I'm using numpy.polyfit to select the best polynomial:

def GetErrorPolynomial(X, Y):

    maxFval = 0.0
    for i in range(1, 10):   # i is the order of the polynomial (max order = 9)
        error_func = np.polyfit(X, Y, i)
        error_func = np.poly1d(error_func)

        # F-test (looking for the largest F value)
        numerator = np.sum(np.square(error_func(X) - np.mean(Y))) / i
        denominator = np.sum(np.square(Y - error_func(X))) / (Y.size - i - 1)
        Fval = numerator / denominator

        if Fval > maxFval:
            maxFval = Fval
            maxFvalPolynomial = error_func

    return maxFvalPolynomial

And here's the code for how I'm using curve_fit:

def poly_sine_fit(x, a, b, c, d, l):
     return a*np.square(x) + b*x + c + d*np.sin(x/l)

param, _ = curve_fit(poly_sine_fit, x_data, y_data)

It's "hardcoded" to a quadratic function, but I want to select the "best" order as I'm doing above with np.polyfit

Answer 1

I finally found a way to model the ripples and can answer my own question. This 2006 paper does curve-fitting on ripples that resemble my dataset.

First off, I did a least squares polynomial fit and then subtracted this polynomial curve from the original data. This left me with only the ripples. Applying the Fourier transform, I picked out the dominant frequencies which let me reconstruct the sine ripples. Then I simply added these ripples to the polynomial curve I had obtained in the beginning. That did it.