我试着用正弦函数来确定一个点的读数. 你可以在图像上看到结果并不令人满意. 我不能设法调整curve_fit的各种参数来很好地适应函数.我能做些什么来提高我的成绩?

from scipy.optimize import curve_fit as cf
import numpy as np
from matplotlib import pyplot as plt

xdata = np.array([0.  , 0.02, 0.04, 0.06, 0.08, 0.1 , 0.12, 0.14, 0.16, 0.18, 0.2 ,
       0.22, 0.24, 0.26, 0.28, 0.3 , 0.32, 0.34, 0.36, 0.38, 0.4 , 0.42,
       0.44, 0.46, 0.48, 0.5 , 0.52, 0.54, 0.56, 0.58, 0.6 , 0.62, 0.64,
       0.66, 0.68, 0.7 , 0.72, 0.74, 0.76, 0.78, 0.8 , 0.82, 0.84, 0.86,
       0.88, 0.9 , 0.92, 0.94, 0.96, 0.98, 1.  , 1.02, 1.04, 1.06, 1.08,
       1.1 , 1.12])
ydata = np.array([37.5, 36.4, 37.1, 37.3, 38.2, 38.4, 38.1, 36.7, 34.3, 32.2, 33.1,
       31.8, 33.4, 35.7, 37.8, 38.3, 38.1, 37. , 34.9, 32.5, 31.6, 31.7,
       33.5, 35.5, 37.8, 38.4, 38.3, 36.9, 34.8, 32.2, 33.9, 31.6, 33.3,
       35.5, 37.6, 38.3, 38.2, 36.7, 34.8, 32.4, 31.5, 32.1, 33.3, 35.5,
       37.2, 38.3, 38.3, 36.8, 34.5, 32.3, 31.6, 31.8, 33.3, 35.7, 37.8,
       38.4, 38.2])


def sin_fun(x,a,b,c,d):
    return a*np.sin(b*x+c)+d

p_opt,p_cov=cf(sin_fun,xdata,ydata )
print(p_opt)
 
plt.plot(xdata,sin_fun(xdata,*p_opt))
plt.plot(xdata,ydata, 'r')
plt.show()

enter image description here

推荐答案

总而言之,@jared的答案是正确的.对于特定的正弦(或通常是任何周期函数),FFT是一种更可靠的方法,特别是当它与第二遍峰值内插相结合时:

import matplotlib.pyplot as plt
import numpy as np

amplitude = np.array((
    37.5, 36.4, 37.1, 37.3, 38.2, 38.4, 38.1, 36.7, 34.3, 32.2, 33.1,
    31.8, 33.4, 35.7, 37.8, 38.3, 38.1, 37. , 34.9, 32.5, 31.6, 31.7,
    33.5, 35.5, 37.8, 38.4, 38.3, 36.9, 34.8, 32.2, 33.9, 31.6, 33.3,
    35.5, 37.6, 38.3, 38.2, 36.7, 34.8, 32.4, 31.5, 32.1, 33.3, 35.5,
    37.2, 38.3, 38.3, 36.8, 34.5, 32.3, 31.6, 31.8, 33.3, 35.7, 37.8,
    38.4, 38.2))

timestep = 0.02
time = np.arange(0, timestep*(len(amplitude) - 0.5), timestep)

# first pass
spectrum = np.fft.rfft(amplitude, norm='forward')
spectral_amplitude = np.abs(spectrum)
freqs = np.fft.rfftfreq(n=time.size, d=timestep)
peak_index = 1 + spectral_amplitude[1:].argmax()

# Instead of peak_freq = freqs[peak_index]
# use parabolic peak interpolation from http://www.ericjacobsen.org/fe2/quadterp.txt
a, b, c = spectrum[peak_index-1: peak_index+2]
norm_frac = np.real((a - c)/(2*b - a - c))
peak_freq_interpolated = freqs[peak_index]*(1 + norm_frac) - freqs[peak_index - 1]*norm_frac

# second pass
# shorten the timeseries length so that the peak lands in the next-left bin
excerpt = int(np.round(time.size * freqs[peak_index-1] / peak_freq_interpolated))
spectrum = np.fft.rfft(amplitude[-excerpt:], norm='forward')
spectral_amplitude = np.abs(spectrum)
freqs = np.fft.rfftfreq(n=excerpt, d=timestep)
peak_index = 1 + spectral_amplitude[1:].argmax()
peak_freq = freqs[peak_index]
phase_correction = (1 - (len(time) - excerpt)*peak_freq*timestep)*2*np.pi
peak_phase = np.angle(spectrum[peak_index]) + phase_correction
peak_amplitude = 2*spectral_amplitude[peak_index]

print(f'Amplitude: {peak_amplitude:.3f}')
print(f'Frequency (Hz):', peak_freq)
print(f'Phase (deg): {np.rad2deg(peak_phase):.2f}')

reconstructed = peak_amplitude * np.cos(
    2*np.pi*peak_freq*time + peak_phase
) + spectral_amplitude[0]

ax_top: plt.Axes
ax_bottom: plt.Axes
fig, (ax_top, ax_bottom) = plt.subplots(nrows=2)

ax_top.plot(time, amplitude, label='experiment')
ax_top.plot(time, reconstructed, label='reconstructed')
ax_top.set_xlabel('Interpolated time (s)')
ax_top.set_ylabel('Amplitude')
ax_top.legend(loc='lower right')

ax_bottom.semilogy(freqs, spectral_amplitude)
ax_bottom.set_xlabel('Frequency (Hz)')
ax_bottom.set_ylabel('Amplitude')

plt.show()
Amplitude: 3.425
Frequency (Hz): 5.0
Phase (deg): 169.70

fit

Python相关问答推荐

提取两行之间的标题的常规表达

将DF中的名称与另一DF拆分并匹配并返回匹配的公司

根据另一列中的nan重置值后重新加权Pandas列

对整个 pyramid 进行分组与对 pyramid 列子集进行分组

如何标记Spacy中不包含特定符号的单词?

如何在Windows上用Python提取名称中带有逗号的文件?

按列分区,按另一列排序

' osmnx.shortest_track '返回有效源 node 和目标 node 的'无'

两个pandas的平均值按元素的结果串接元素.为什么?

如何获得每个组的时间戳差异?

从Windows Python脚本在WSL上运行Linux应用程序

在单次扫描中创建列表

如何使用使用来自其他列的值的公式更新一个rabrame列?

当单元测试失败时,是否有一个惯例会抛出许多类似的错误消息?

将一个双框爆炸到另一个双框的范围内

Python pint将1/华氏度转换为1/摄氏度°°

pandas fill和bfill基于另一列中的条件

修改.pdb文件中的值并另存为新的

用由数据帧的相应元素形成的列表的函数来替换列的行中的值

为什么在更新Pandas 2.x中的列时,数据类型不会更改,而在Pandas 1.x中会更改?