I'd like to achieve a fourier series development for a x-y-dataset using numpy and scipy.
At first I want to fit my data with the first 8 cosines and plot additionally only the first harmonic. So I wrote the following two function defintions:
# fourier series defintions
tau = 0.045
def fourier8(x, a1, a2, a3, a4, a5, a6, a7, a8):
return a1 * np.cos(1 * np.pi / tau * x) + \
a2 * np.cos(2 * np.pi / tau * x) + \
a3 * np.cos(3 * np.pi / tau * x) + \
a4 * np.cos(4 * np.pi / tau * x) + \
a5 * np.cos(5 * np.pi / tau * x) + \
a6 * np.cos(6 * np.pi / tau * x) + \
a7 * np.cos(7 * np.pi / tau * x) + \
a8 * np.cos(8 * np.pi / tau * x)
def fourier1(x, a1):
return a1 * np.cos(1 * np.pi / tau * x)
Then I use them to fit my data:
# import and filename
filename = 'data.txt'
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
z, Ua = np.loadtxt(filename,delimiter=',', unpack=True)
tau = 0.045
# plot data
fig = plt.figure()
ax1 = fig.add_subplot(111)
p1, = plt.plot(z,Ua)
# fits
popt, pcov = curve_fit(fourier8, z, Ua)
# further plots
Ua_fit8 = fourier8(z,*popt)
Ua_fit1 = fourier1(z,popt[0])
p2, = plt.plot(z,Ua_fit8)
p3, = plt.plot(z,Ua_fit1)
plt.show()
which works as desired:
But know I got stuck making it generic for arbitary orders of harmonics, e.g. I want to fit my data with the first fifteen harmonics and plot the first three harmonics.
How could I achieve that without defining fourier1
, fourier2
, fourier3
... , fourier15
?
Here is the example text file data.txt
to play with it:
1.000000000000000021e-03,2.794682735905079767e+02
2.000000000000000042e-03,2.792294526290349950e+02
2.999999999999999629e-03,2.779794770690260179e+02
4.000000000000000083e-03,2.757183469104809888e+02
5.000000000000000104e-03,2.734572167519349932e+02
5.999999999999999258e-03,2.711960865933900209e+02
7.000000000000000146e-03,2.689349564348440254e+02
8.000000000000000167e-03,2.714324329982829909e+02
8.999999999999999320e-03,2.739299095617229796e+02
1.000000000000000021e-02,2.764273861251620019e+02
1.100000000000000110e-02,2.789248626886010243e+02
1.199999999999999852e-02,2.799669443683339978e+02
1.299999999999999940e-02,2.795536311643609793e+02
1.400000000000000029e-02,2.791403179603880176e+02
1.499999999999999944e-02,2.787270047564149991e+02
1.600000000000000033e-02,2.783136915524419805e+02
1.699999999999999775e-02,2.673604939531239779e+02
1.799999999999999864e-02,2.564072963538059753e+02
1.899999999999999953e-02,2.454540987544889958e+02
2.000000000000000042e-02,2.345009011551709932e+02
2.099999999999999784e-02,1.781413355804160119e+02
2.200000000000000219e-02,7.637540203022649621e+01
2.299999999999999961e-02,-2.539053151996269975e+01
2.399999999999999703e-02,-1.271564650701519952e+02
2.500000000000000139e-02,-2.289223986203409993e+02
2.599999999999999881e-02,-2.399383538664330047e+02
2.700000000000000316e-02,-2.509543091125239869e+02
2.800000000000000058e-02,-2.619702643586149975e+02
2.899999999999999800e-02,-2.729862196047059797e+02
2.999999999999999889e-02,-2.786861050144170235e+02
3.099999999999999978e-02,-2.790699205877460258e+02
3.200000000000000067e-02,-2.794537361610759945e+02
3.300000000000000155e-02,-2.798375517344049968e+02
3.399999999999999550e-02,-2.802213673077350222e+02
3.500000000000000333e-02,-2.776516459805940258e+02
3.599999999999999728e-02,-2.750819246534539957e+02
3.700000000000000511e-02,-2.725122033263129993e+02
3.799999999999999906e-02,-2.699424819991720028e+02
3.899999999999999994e-02,-2.698567311502329744e+02
4.000000000000000083e-02,-2.722549507794930150e+02
4.100000000000000172e-02,-2.746531704087540220e+02
4.199999999999999567e-02,-2.770513900380149721e+02
4.299999999999999656e-02,-2.794496096672759791e+02
4.400000000000000439e-02,-2.800761105821779893e+02
4.499999999999999833e-02,-2.800761105821779893e+02
4.599999999999999922e-02,-2.800761105821779893e+02
4.700000000000000011e-02,-2.800761105821779893e+02
4.799999999999999406e-02,-2.788333722531979788e+02
4.900000000000000189e-02,-2.763478955952380147e+02
5.000000000000000278e-02,-2.738624189372779938e+02
5.100000000000000366e-02,-2.713769422793179729e+02
5.199999999999999761e-02,-2.688914656213580088e+02
5.299999999999999850e-02,-2.715383673199499981e+02
5.400000000000000633e-02,-2.741852690185419874e+02
5.499999999999999334e-02,-2.768321707171339767e+02
5.600000000000000117e-02,-2.794790724157260229e+02
5.700000000000000205e-02,-2.804776351435970128e+02
5.799999999999999600e-02,-2.798278589007459800e+02
5.899999999999999689e-02,-2.791780826578950041e+02
5.999999999999999778e-02,-2.785283064150449945e+02
6.100000000000000561e-02,-2.778785301721940186e+02
6.199999999999999956e-02,-2.670252067497989970e+02
6.300000000000000044e-02,-2.561718833274049985e+02
6.400000000000000133e-02,-2.453185599050100052e+02
6.500000000000000222e-02,-2.344652364826150119e+02
6.600000000000000311e-02,-1.780224826854309867e+02
6.700000000000000400e-02,-7.599029851345700592e+01
6.799999999999999101e-02,2.604188565851649884e+01
6.900000000000000577e-02,1.280740698304900036e+02
7.000000000000000666e-02,2.301062540024639986e+02
7.100000000000000755e-02,2.404921248105050040e+02
7.199999999999999456e-02,2.508779956185460094e+02
7.299999999999999545e-02,2.612638664265870148e+02
7.400000000000001021e-02,2.716497372346279917e+02
7.499999999999999722e-02,2.773051723900500178e+02
7.599999999999999811e-02,2.782301718928520131e+02
7.699999999999999900e-02,2.791551713956549747e+02
7.799999999999999989e-02,2.800801708984579932e+02
7.900000000000001465e-02,2.810051704012610116e+02
8.000000000000000167e-02,2.785107135689390248e+02
8.099999999999998868e-02,2.760162567366169810e+02
8.200000000000000344e-02,2.735217999042949941e+02
8.300000000000000433e-02,2.710273430719730072e+02
8.399999999999999134e-02,2.706544464035359852e+02
8.500000000000000611e-02,2.724031098989830184e+02
8.599999999999999312e-02,2.741517733944299948e+02
8.699999999999999400e-02,2.759004368898779944e+02
8.800000000000000877e-02,2.776491003853250277e+02
8.899999999999999578e-02,2.783445666445250026e+02
8.999999999999999667e-02,2.790400329037249776e+02
The return value popt contains the best-fit values of the parameters. The return value pcov contains the covariance (error) matrix for the fit parameters. From them we can determine the standard deviations of the parameters, just as we did for linear least chi square.
The returned parameter covariance matrix pcov is based on scaling sigma by a constant factor. This constant is set by demanding that the reduced chisq for the optimal parameters popt when using the scaled sigma equals unity. In other words, sigma is scaled to match the sample variance of the residuals after the fit.
The SciPy API provides a 'curve_fit' function in its optimization library to fit the data with a given function. This method applies non-linear least squares to fit the data and extract the optimal parameters out of it.
I would approach this by defining a single function fourier
that accepts a variable number of arguments, using a loop through the cosines to evaluate the function for each of your input points:
def fourier(x, *a):
ret = a[0] * np.cos(np.pi / tau * x)
for deg in range(1, len(a)):
ret += a[deg] * np.cos((deg+1) * np.pi / tau * x)
return ret
Because this function takes a variable number of arguments, you need to provide a starting position of the correct dimensionality as a hint to the curve_fit
function so it knows the number of cosines to use. In the example you provide with 8 cosines:
popt, pcov = curve_fit(fourier, z, Ua, [1.0] * 8)
Now, the use case you're asking about (fitting with 15 harmonics and plotting the first three) can be accomplished with:
# Fit with 15 harmonics
popt, pcov = curve_fit(fourier, z, Ua, [1.0] * 15)
# Plot data, 15 harmonics, and first 3 harmonics
fig = plt.figure()
ax1 = fig.add_subplot(111)
p1, = plt.plot(z,Ua)
p2, = plt.plot(z, fourier(z, *popt))
p3, = plt.plot(z, fourier(z, popt[0], popt[1], popt[2]))
plt.show()
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With