scipy 输入可变数量的参数以拟合2D傅立叶级数

4zcjmb1e  于 2022-11-10  发布在  其他
关注(0)|答案(1)|浏览(119)

我正在研究一个使用二维傅立叶级数的优化问题,基本上我正在尝试实现以下公式:

我想在我的脚本中将其实现为阶数 m 的函数,在某种意义上说,只要提供 m,我就可以用 m 阶级数拟合一些数据,并获得级数的最佳系数(除了我已经知道的 C_00
无论如何,当涉及到使用Scipy的curve_fit函数与参数数量可变的模型(例如,如果 m=1 我有8个系数+ 2个角速度,如果 m=2 我有24个系数+ 2个角速度等)时,我很挣扎。
到目前为止,我已经写了以下内容:

c00=50
def fourier(x, v, w1, w2):
   f=c00
   k=0
   for i in range(1, m+1):
       f=f+v[k]*np.cos(i*w1*x)+v[k+1]*np.sin(i*w1*x)
       k+=2
   for j in range(1, m+1):
      for i in range(-m, m+1):
           f=f+v[k]*np.cos(i*w1*x+j*w2*x)+v[k+1]*np.sin(i*w1*x+j*w2*x)
           k+=2
return f

time = [] #some x data I get from a txt file
energy = [] # some y data I get from a txt file
popt, pcov = curve_fit(fourier, time, energy)

但当我运行脚本时,它引发了以下错误:
文件“/home/arcalinux/Desktop/fourier_series.py”,第63行,傅立叶格式f=f+v[k]np.cos(iw1x)+v[k+1]np.sin(iw1x)索引错误:标量变量的索引无效。
我还尝试了其他功能:

def fourier2(args):
   f=c00
   k=3
   for i in range(1, m+1):
       f=f+args[k]*np.cos(i*args[1]*args[0])+args[k+1]*np.sin(i*args[1]*args[0])
       k+=2
   for j in range(1, m+1):
       for i in range(-m, m+1):
          f=f+args[k]*np.cos(i*args[0]*x+j*args[1]*x)+
                          +args[k+1]*np.sin(i*args[0]*x+j*args[1]*x)
          k*=2
   return f

但是当我运行curve_fit时,我有:
ValueError:无法确定拟合参数的数目。
请注意,傅立叶级数函数正确工作,例如(随机数):

v=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
m=2
x=2
w1=3
w2=4
print(fourier(x,v,w1,w2))

49.136961626348295
同样的道理也适用于傅立叶2
我怎样才能解决这个问题呢?谢谢!
完整代码如下所示:


# importing libraries

import numpy as np
import math
from scipy.optimize import curve_fit
from statistics import mean
import matplotlib.pyplot as plt

m=2 #order of series
number_of_coeff=(2*m+1)**2 - 1 
N=number_of_coeff

# build series

coefficients=[]

# first part, j=0

for i in range(1,m+1):
   coefficients.append('c'+str(i)+'0')
   coefficients.append('s'+str(i)+'0')

# second part

for j in range(1, m+1):
   for i in range(-m, 0):
       coefficients.append('g'+str(-i)+str(j))
       coefficients.append('t'+str(-i)+str(j))
   for i in range(0, m+1):
       coefficients.append('c'+str(i)+str(j))
       coefficients.append('s'+str(i)+str(j))
print(len(coefficients))
v=np.ones(N)
m=2
x=2
w1=3
w2=4
arg=[2,3,4,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
arg3=[3,4,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]

c00=50
def fourier(x, v, w1, w2):
    print(type(v))
    print('v is ',v)
    f=c00
    k=0
    print(v[0])
    for i in range(1, m+1):
        f=f+v[k]*np.cos(i*w1*x)+v[k+1]*np.sin(i*w1*x)

        k+=2
    for j in range(1, m+1):
        for i in range(-m, m+1):
            f=f+v[k]*np.cos(i*w1*x+j*w2*x)+
                v[k+1]*np.sin(i*w1*x+j*w2*x)
            k+=2
return f
def fourier2(args):
    f=c00
    k=3
    for i in range(1, m+1):
        f=f+args[k]*np.cos(i*args[1]*args[0])+
            args[k+1]*np.sin(i*args[1]*args[0])
        k+=2
    for j in range(1, m+1):
        for i in range(-m, m+1):
            f=f+args[k]*np.cos(i*args[1]*args[0]+
                j*args[2]*args[0])+
                args[k+1]*np.sin(i*args[1]*args[0]+
                j*args[2]*args[0])
        k+=2

# print('number of args =',len(args))

return f
def fourier3(x,args):
    f=c00
    k=2
    for i in range(1, m+1):
        f=f+args[k]*np.cos(i*args[0]*x)+
            args[k+1]*np.sin(i*args[0]*x)
        k+=2
    for j in range(1, m+1):
        for i in range(-m, m+1):
            f=f+args[k]*np.cos(i*args[0]*x+j*args[1]*x)+
                args[k+1]*np.sin(i*args[0]*x+j*args[1]*x)
            k+=2

# print('number of args =',len(args))

return f

# print(fourier2(arg))

# print(fourier3(x,arg3))

# print(fourier(2,v,3,4))

# getting data

time = [] 
energy = [] 
filepath=''
filename=''
with open(filepath+filename, 'r') as f:
    for line in f: # iterate through the file
        if not line: continue 
        t, e = line.split() 
        time.append(float(t)),
        energy.append(float(e)) # add a value
c00=mean(energy)
time=np.array(time)
popt, pcov = curve_fit(fourier, time, energy)
print(popt)
tyg4sfes

tyg4sfes1#

根据scipy documentation,您可以为curve_fit函数提供p0参数,这是您对函数参数的初始猜测(这里是vw1w2)。
如果不这样做,函数将假定参数的初始值为1,这就是为什么v将得到标量而不是数组。
因此,在调用curve_fit之前,您可以执行以下操作:

initial_v = ...  # Replace with value (an array)
initial_w1 = ... # Replace with value (a scalar)
initial_w2 = ... # Replace with value (a scalar)

popt, pcov = curve_fit(fourier, time, energy, p0=[initial_v, initial_w1, initial_w2])

相关问题