Cree dinámicamente definiciones de ajuste de curvas para series de Fourier dinámicamente.

Me gustaría lograr un desarrollo de la serie de Fourier para un conjunto de datos xy con números y scipy .

Al principio quiero ajustar mis datos con los primeros 8 cosenos y trazar, además, solo el primer armónico . Así que escribí las siguientes dos definiciones de funciones:

# 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) 

Luego los uso para ajustar mis datos:

 # 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() 

que funciona como se desea:

introduzca la descripción de la imagen aquí

Pero sepa que me quedé atascado haciéndolo genérico para órdenes arbitrarias de armónicos, por ejemplo, quiero ajustar mis datos con los primeros quince armónicos y trazar los primeros tres armónicos.

¿Cómo podría lograr eso sin definir fourier1 , fourier2 , fourier3 …, fourier15 ?


Aquí está el ejemplo de archivo de texto data.txt para jugar con él:

 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 

Me acercaría a esto definiendo una única función de fourier que acepta un número variable de argumentos, utilizando un bucle a través de los cosenos para evaluar la función para cada uno de sus puntos de entrada:

 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 

Debido a que esta función toma un número variable de argumentos, debe proporcionar una posición inicial de la dimensionalidad correcta como una sugerencia a la función curve_fit para que sepa la cantidad de cosenos que debe usar. En el ejemplo usted proporciona 8 cosenos:

 popt, pcov = curve_fit(fourier, z, Ua, [1.0] * 8) 

Ahora, el caso de uso que está preguntando (encajando con 15 armónicos y trazando los tres primeros) se puede lograr con:

 # 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() 

introduzca la descripción de la imagen aquí

basado en este hilo.

¿Pasando argumentos adicionales usando scipy.optimize.curve_fit?

 # 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 def make_fourier(na, nb): def fourier(x, *a): ret = 0.0 for deg in range(0, na): ret += a[deg] * np.cos((deg+1) * np.pi / tau * x) for deg in range(na, na+nb): ret += a[deg] * np.sin((deg+1) * np.pi / tau * x) return ret return fourier popt, pcov = curve_fit(make_fourier(15,15), z, Ua, [0.0]*30) # 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, (make_fourier(15,15))(z, *popt)) #p3, = plt.plot(z, (make_fourier(8,8))(z, popt[0], popt[1], popt[2])) plt.show()