Approximating the area of a chirp by fitting a polynomial

Once in a while we need to estimate the area of a dataset in which we are interested. This area could give us, for example, force (mass vs acceleration) or electric power (electric current vs charge).

One way to do that is fitting a curve on our data, and let's face it: this is not that easy. In this post we will work on this issue using Python and its packages. If you do not have Python installed on your system, check here how to proceed.

We used the form import package as nickname to avoid typing the entire package name. Note that we got only pyplot from matplotlib, which provides a MATLAB-like plotting framework. Then, we use np, sp and plt instead of numpy, scipy and matplotlib.pyplot respectively.

On our example, let's calculate the area of a chirp. Scipy's function signal.chirp is an easy function to generate a chirp varying in time. This is its general form:

On our example, the time interval begins on 0 and ends on 10. It has 100 samples. The following code generates this interval and stores it on the time_interval variable.

time_interval = np.linspace(0, 10, num=100)

Now, using time_interval, let's generate our chirp. Its values of f0, t1 and f1 are 0, 1 and 1, respectively. We will generate a quadratic chirp, using method='quadratic'. This code generates this chirp and stores its points on example_chirp.

Let's use plot the chirp using stem's default values. We use the following code:

plt.stem(time_interval, example_chirp)

The figure representing our discrete chirp is given on Figure 1.Figure 1. Our example chirp, generated from sp.signal.chirp and plotted using plt.stem.

Now we will fit a curve on this chirp. On this example we use Chebyshev's polynomials, that can be fitted on a dataset using numpy's polyomial.Chebyshev.fit function. There are many alternatives, such as Legendre, Laguerre and Hermite. Check these possibilities at the scipy online documentation.

The simplest form of polynomial.Chebyshev.fit follows:

np.polynomial.Chebyshev.fit(x, y, deg)

Where:

x are the x-coordinates of the sample points.

y are the y-coordinates of the sample points.

deg is the degree of the polynomial to be fitted.

For more information and arguments to this function, please check the scipy documentation. As an alternative, type np.polynomial.Chebyshev.fit? on your Python prompt.

To show the behavior of this function, we will use three different arguments for deg: 20, 45 and 115. The following code will generate these different polynomials based on our data and store them on the variables cheby_chirp20, cheby_chirp45 and cheby_chirp115.

Let's plot these results. Note that the results of polynomial.Chebyshev.fit are functions, not points; we need to provide the plot interval as an argument. We plot these three polynomials against the original data as subplots. They will be represented as green lines ('g') with width equals to 2 (lw=2). We also add a label to each of them, which will appear on the legend of each plot. Check the matplotlib examples for more ideas on subplots.

Now, let's use scipy's integrate.quad function to calculate the area of these polynomials. This function uses a technique from QUADPACK, a Fortran library. The simplest form of integrate.quad is given below:

scipy.integrate.quad(func, a, b)

Where:

func is the function to be integrated.

a is the lower limit of integration.

b is the upper limit of integration.

This function returns the integral value and the absolute error in the result. For more information on integrate.quad, its arguments and usage examples, check the scipy documentation. You can also type sp.integrate.quad? on your Python prompt.

Let's define the integration interval between zero and pi. Then, we integrate the three polynomials. Check the code below:

We need to present area and error now. This can be done using print, which provides a nice output:

print('The area of cheby20 in this interval is {0}. The absolute error is {1}'.format(cheby20_area, cheby20_err))
print('The area of cheby45 in this interval is {0}. The absolute error is {1}'.format(cheby45_area, cheby45_err))
print('The area of cheby115 in this interval is {0}. The absolute error is {1}'.format(cheby115_area, cheby115_err))

The function print presents a message. Using format we can also present the variable values on {0}, {1}, {2}, and so on. More on the print function can be seen on the Python documentation. As always, you can type print? on your Python prompt.

Let's check the results:

The area of cheby20 in this interval is 0.4637693754972824. The absolute error is 1.5334233932461724e-13
The area of cheby45 in this interval is -8.792912670367262. The absolute error is 2.9585315540800773e-12
The area of cheby115 in this interval is 151230239254.8722. The absolute error is 50.06303513050079

The two first polynomials wielded acceptable results, but there is something wrong about cheby_chirp115. We can confirm that checking its absolute error (around 50). Let's see what is going on, plotting the integration areas of these polynomials.

We define the integration interval with 50 points using linspace to plot the results. This interval will be attributed to integ_full. Again, we will use plt.subplots:

We employed plt.fill_between to plot the integration areas. For comparison, we plotted also the original points and the polynomials. Note that we added LaTeX code to the label; it is necessary to add an r before the string. The areas of integration are shown on Figure 3.

Figure 3. The integration areas highlighted on the fitted polynomials.

Check how the integration interval of cheby_chirp115 is messy! cheby_chirp45 is not good too. cheby_chirp20 is satisfactory from 0 to 1 approximately, but we need more than that on the rest of our data.

Lesson learned: be careful choosing the degree of the fitted polynomial! Lower degrees may not represent your data nicely; however, higher degrees may lead to a poor conditioned fit. The solution could be to search for a good balance.

All right, this is it! Hope you enjoyed this post. Are there anything wrong? Do you have any ideas to improve it? Please share them with us! Also, you can download the codes on GitHub.

Thank you for joining us! See you next time!

Additional notes:

jms_nh proposed another tools for fitting the chirp: "The waveforms you've shown really aren't very polynomial-like and would be better handled by curve-fitting using cubic splines or Savitsky-Golay filtering". He has a point: according to the SciPy documentation, "Fits using Chebyshev series are usually better conditioned than fits using power series, but much can depend on the distribution of the sample points and the smoothness of the data. If the quality of the fit is inadequate splines may be a good alternative".

Comments:

The waveforms you've shown really aren't very polynomial-like and would be better handled by curve-fitting using cubic splines or Savitsky-Golay filtering.

Also, if you do have a Chebyshev polynomial P(x), the way to integrate it is not numerically using something like scipy.integrate.quad, but rather algebraically by manipulating the components to get the Chebyshev coefficients of the indefinite integral of P(x), and evaluating the integral directly.

[ - ]

Comment by alexandrejaguar●November 16, 2015

Hey jms_nh, thanks again for your reply!1. sure thing. Actually, the SciPy documentation says that too:"Fits using Chebyshev series are usually better conditioned than fits using power series, but much can depend on the distribution of the sample points and the smoothness of the data. If the quality of the fit is inadequate splines may be a good alternative."I will write this consideration as an Additional Note to the post, ok?Again, keep in mind that this post is only a thought on the balance of a fit and its consequences.

2. sure thing². I could also use the Simpson's rule, composite trapezoidal, ... but this is the subject of a next post ;)"Thanks for your thoughts! Have a nice one!

[ - ]

Comment by jms_nh●November 16, 2015

Why not just get the continuous-time signal intended by the chirp function, and integrate that? (numerically, with more sample points than shown here, or symbolically)

[ - ]

Comment by alexandrejaguar●November 16, 2015

Hi jms_nh, thank you for your reply!Keep in mind that I used the discrete chirp points as an example, you can use any discrete data points you want. The purpose here is to show how to act on a discrete sample; sorry if this is not clear on the post. I'll think about it and try to improve the motivation, ok? Thanks!

[ - ]

Comment by lagus●November 16, 2015

Educative sir Linus

[ - ]

Comment by alexandrejaguar●November 16, 2015

Hi Linus, thank you for your reply!I really appreciate it! Feel free to ask a question or point a doubt if you'd like.Thanks again!

[ - ]

Comment by groger●February 20, 2016

Do you know if the python/scipy library can be compiled for use in an embedded application (i.e., by GCC) for use in a compiler like IAR?

thanks!

[ - ]

Comment by alexandrejaguar●February 20, 2016

Hi groger,I don't know how/if one could do it. I know there exists py_compile [1] and compileall [2], however. I just don't know how to use them.Check it out if you'd like:[1] https://docs.python.org/3.5/library/py_compile.html#module-py_compile[2] https://docs.python.org/3.5/library/compileall.html#module-compileallThanks! Have a nice one!

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.