Fitting a sine curve to data

Measured waves are never clean

A sensor follows a swinging spring and records its displacement every 0.05 s. That gives 40 readings in 2 s. Each reading is a little off, because of noise. The true motion is a sine, and you want three numbers from it: the amplitude $A$, the frequency $f$ and the phase $\varphi$.

Fitting with curve_fit

scipy.optimize.curve_fit(model, t, y, p0=guess) changes the unknowns until the sum of the squared misses is as small as possible. You give it a function model(t, A, f, phi), the readings and a first guess p0. It returns two things:

• popt, the best values of the unknowns;
• pcov, a table from which the errors come: np.sqrt(np.diag(pcov)) gives one error for each unknown.

Here is a fit to clean data, so that you can see the idea in a few lines:

It prints [0.5 1. 0.5]. The fit finds the three numbers that made the data, starting from a rough guess.

The first guess matters

curve_fit works downhill from p0. For a sine, the sum of squared misses has many dips, one for each frequency that nearly fits. A first frequency far from the truth can slide into the wrong dip. You need a guess for $f$ that is already close.

A first guess from the FFT

The fast Fourier transform, FFT, splits a signal into sine waves of different frequencies.

• np.fft.rfft(y) does the splitting. np.abs(...) of it gives the size of each frequency.
• np.fft.rfftfreq(n, d) gives the frequencies, for $n$ readings that are $d$ seconds apart.
• The frequencies are $\dfrac{1}{nd}$ apart, which is one over the whole time of the record. A 2 s record gives frequencies 0.5 Hz apart.

It prints 3.0 Hz: the frequency of the tallest bar of the spectrum. This is the first guess for curve_fit.

Taking readings too slowly

One reading every $\dfrac{1}{f_s}$ seconds can show only frequencies below $\dfrac{f_s}{2}$. This is the Nyquist limit. A wave needs at least two readings in each cycle. Above the limit, the wave folds back and looks like a slower one. This is aliasing.

Here a 9 Hz wave is read 10 times each second. The dots fall on a 1 Hz wave:

It prints 1.0 Hz. The FFT cannot tell the readings from those of a 1 Hz wave. The cure is to read faster, more than twice the highest frequency in the signal.

Practise with new numbers

Try it in Python

The program makes 40 noisy readings, takes a first guess from the FFT, and fits a sine with curve_fit. It prints each unknown with its error. The picture shows what it finds: the readings and the fitted sine above, and the FFT below.

Common mistakes

• No first guess, or a poor one. The fit can settle on the wrong frequency. Take the guess from the FFT.
• Reading pcov as the errors. pcov holds squares. The errors are np.sqrt(np.diag(pcov)).
• Putting $\omega$ where the model expects $f$. If the model has 2 * np.pi * f * t, the fit returns a frequency in Hz, not in rad/s.
• Reading too slowly. Fewer than two readings in a cycle give a wrong frequency that looks real.
• Forgetting the units of the FFT axis. rfftfreq gives Hz only when d is in seconds.

Remember these

• Fit a sine with curve_fit(model, t, y, p0=[A, f, phi]).
• The errors are np.sqrt(np.diag(pcov)).
• Take the first guess from the tallest bar of np.abs(np.fft.rfft(y)).
• The FFT frequencies are $\dfrac{1}{nd}$ apart.
• The Nyquist limit is $\dfrac{f_s}{2}$. A faster wave aliases to a slower one.

Check yourself

Try each one on your own first. If you are stuck, open a hint. Then check with Show the answer.

1. What does np.sqrt(np.diag(pcov)) give?

2. A fit starts from a frequency far from the truth and returns a frequency that does not match the data. What went wrong, and what would you do?

3. 30 readings are 0.1 s apart. How far apart are the FFT frequencies?

4. A logger takes 100 readings each second. What is the highest frequency it can show?

5. A fit gives an amplitude of $0.784 \pm 0.016$ m. Is the amplitude known to better than 5 %?