在 python 中拟合周期图：scipy.interpolate.splrep 的参数，曲线方程？

Question

[原问题]

根据下面的数据，我需要一个随着时间的流逝而无限增加的曲线方程。如何获得？

[问题更新]

我需要为scipy.interpolate.splrep 指定正确的参数。有人可以帮忙吗？

另外，有没有办法从 b 样条的系数中得到方程？

[替代问题]

如何使用傅里叶级数分解信号获得拟合？

这个图似乎是线性图的组合，一个周期函数 pf1 增加了四倍，一个更大的周期函数导致 pf1 无限期地一次又一次地发生。剧情的难度就是这个问题被问到的原因。

数据：

Time elapsed in sec.    TX + RX Packets
(0,0)
(10,2422)
(20,2902)
(30,2945)
(40,3059)
(50,3097)
(60,4332)
(70,4622)
(80,4708)
(90,4808)
(100,4841)
(110,6081)
(120,6333)
(130,6461)
(140,6561)
(150,6585)
(160,7673)
(170,8091)
(180,8210)
(190,8291)
(200,8338)
(210,8357)
(220,8357)
(230,8414)
(240,8414)
(250,8414)
(260,8414)
(270,8414)
(280,8414)
(290,8471)
(300,8471)
(310,8471)
(320,8471)
(330,8471)
(340,8471)
(350,8471)
(360,8471)
(370,8471)
(380,8471)
(390,8471)
(400,8471)
(410,8471)
(420,8528)
(430,8528)
(440,8528)
(450,8528)
(460,8528)
(470,8528)
(480,8528)
(490,8528)
(500,8528)
(510,9858)
(520,10029)
(530,10129)
(540,10224)
(550,10267)
(560,11440)
(570,11773)
(580,11868)
(590,11968)
(600,12039)
(610,13141)

我的代码：

import numpy as np
import matplotlib.pyplot as plt

points = np.array(
    [(0,0), (10,2422), (20,2902), (30,2945), (40,3059), (50,3097), (60,4332), (70,4622), (80,4708), (90,4808), (100,4841), (110,6081), (120,6333), (130,6461), (140,6561), (150,6585), (160,7673), (170,8091), (180,8210), (190,8291), (200,8338), (210,8357), (220,8357), (230,8414), (240,8414), (250,8414), (260,8414), (270,8414), (280,8414), (290,8471), (300,8471), (310,8471), (320,8471), (330,8471), (340,8471), (350,8471), (360,8471), (370,8471), (380,8471), (390,8471), (400,8471), (410,8471), (420,8528), (430,8528), (440,8528), (450,8528), (460,8528), (470,8528), (480,8528), (490,8528), (500,8528), (510,9858), (520,10029), (530,10129), (540,10224), (550,10267), (560,11440), (570,11773), (580,11868), (590,11968), (600,12039), (610,13141)]
    )
# get x and y vectors
x = points[:,0]
y = points[:,1]

# calculate polynomial
z = np.polyfit(x, y, 3)
print z
f = np.poly1d(z)

# calculate new x's and y's
x_new = np.linspace(x[0], x[-1], 50)
y_new = f(x_new)

plt.plot(x,y,'o', x_new, y_new)
plt.xlim([x[0]-1, x[-1] + 1 ])
plt.show()

我的输出：

我的代码 2：

import numpy as N
from scipy.interpolate import splprep, splev

x = N.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, 400, 410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530, 540, 550, 560, 570, 580, 590, 600, 610])
y = N.array([0, 2422, 2902, 2945, 3059, 3097, 4332, 4622, 4708, 4808, 4841, 6081, 6333, 6461, 6561, 6585, 7673, 8091, 8210, 8291, 8338, 8357, 8357, 8414, 8414, 8414, 8414, 8414, 8414, 8471, 8471, 8471, 8471, 8471, 8471, 8471, 8471, 8471, 8471, 8471, 8471, 8471, 8528, 8528, 8528, 8528, 8528, 8528, 8528, 8528, 8528, 9858, 10029, 10129, 10224, 10267, 11440, 11773, 11868, 11968, 12039, 13141])

# spline parameters
s=1.0 # smoothness parameter
k=3 # spline order
nest=-1 # estimate of number of knots needed (-1 = maximal)

# find the knot points
tckp,u = splprep([x,y],s=s,k=k,nest=nest,quiet=True,per=1)

# evaluate spline, including interpolated points
xnew,ynew = splev(N.linspace(0,1,400),tckp)

import pylab as P
data,=P.plot(x,y,'bo-',label='data')
fit,=P.plot(xnew,ynew,'r-',label='fit')
P.legend()
P.xlabel('x')
P.ylabel('y')

P.show()

我的输出 2：

Answer 1

看起来你有反应动力学：

#%%
import numpy as np
from scipy.integrate import odeint
from scipy import optimize
from matplotlib import pyplot as plt
#%%
data = []
with open('data.txt', 'r') as f:
    for line in f:
        data.append(line.strip(' \n ()').split(','))
data = np.array(data,dtype=float)
data = data[0:-1].T

#%%
slope = np.diff(data[1])
index = np.where(slope>1000)
index = np.append(index, len(data[0]) -1 )
plt.plot(data[0],data[1],'.')
plt.plot(data[0,index],data[1,index],'ro')
plt.plot(data[0,1:],np.diff(data[1]))

从这里开始，我假设反应从每个标记点（红色）开始。我确信代码可以写得更干净，但这是第一次快速而肮脏的黑客攻击。您可以使用 scipy curvefit 或类似方法来拟合反应常数 k

#%%
time = data[0,index]

def model(y,t,k):
    dydt = np.zeros(2)
    dydt[0] = -k*y[0]
    dydt[1] = k*y[0]

    return dydt


def res(k):
    y_hat = []
    t_hat = []
    for i in xrange(len(index) -1):
        '''
        I assume that at every timepoint the reaction is initiated by
        adding y[i + 1] - y[i] new datapackages. Over time they are 
        converted to sent packages. All packages which do not react,
        do not take part in the next cycle.
        '''
        y0 = [data[1, index[i+1]] - data[1, index[i]], 0]
        t0 = data[0, index[i]:index[i+1]]
        y_int,info = odeint(model, y0, t0, args=(k,), full_output = 1 )
        # I am not very happy about the construct below, but could
        # not find a better solution.
        y_hat.append(list(y_int[:,1]))
        t_hat.append(list(t0))
    return y_hat,t_hat

k = 2e-1
y,t = res(k)
''' It may be possible to play with y0[1] in the model in order 
to avoid the construct below. But since I started every reaction at y[1]=0
I have to add the last y value from the last step. This is a bit of a hack,
since data[1, index[i]] is not necessarily the corresponding solution. But hey, It seems to work.
'''
y_hat = [subitem + data[1, index[i]] for i,item in enumerate(y) for subitem in item]
t_hat = [subitem for item in t for subitem in item]
y_hat = np.array(y_hat,dtype=float)
t_hat = np.array(t_hat,dtype=float)


#%%
plt.plot(data[0],data[1],'.')
plt.plot(data[0,index],data[1,index],'ro')
plt.plot(t_hat,y_hat,'go')

另一种方法可能是（可能在物理上更正确）在每个时间点添加一个高斯峰的 CDF。