我正在尝试找到我从模拟输出的一些振荡的阻尼率。
为此，我首先找到了这些振荡的峰值，并以指数拟合拟合。
该方法在大多数情况下都可以正常工作，但是对于昨天我遇到的情况，我发现该方法无法正常工作，也没有明显的原因，我可以找到它不起作用。在某些情况下，我使用的高峰发现似乎错过了峰值，似乎存在一些问题。
这是查找峰的程序：
（注意：列表中的变量Tarr，Zarr，Alpharr先前取自以前的单元格运行）
from scipy.optimize import curve_fit
def exp(t, A, lbda):
return A * np.exp(-lbda * t)
def find_peaks(x, y):
peak_x = []
peak_vals = []
for i in range(len(y)):
if i == 0:
if y[i+1]<y[i]:
peak_x.append(x[i])
peak_vals.append(y[i])
else:
continue
if i == len(y) - 1:
continue
if (y[i-1] < y[i]) and (y[i+1] < y[i]):
peak_x.append(x[i])
peak_vals.append(y[i])
return np.array(peak_x), np.array(peak_vals)
Tarr=[Tplot0,Tplot1,Tplot2,Tplot3,Tplot4]
Zarr=[Zrel0,Zrel1,Zrel2,Zrel3,Zrel4]
Alpharr=[alpha0,alpha1,alpha2,alpha3,alpha4]
for x0,y0,alf in zip(Tarr,Zarr,Alpharr):
peak_x0,peak_y0=find_peaks(x0,y0)
popt, pcov = curve_fit(exp, peak_x0, peak_y0)
print(*[f"{val:.2f} +/- {err:.2f} __ " for val, err in zip(popt, np.sqrt(np.diag(pcov)))])
print("decay rate = ",popt[1])
plt.plot(x0,y0)
plt.plot(peak_x0,exp(peak_x0,*popt),'-r', label=str((popt[0]))+"*exp(-"+str((popt[1]))+")")
plt.xlabel("time")
plt.ylabel("Zrel ")
plt.title("n = "+str(n)+" ; alpha = "+str(alf))
plt.legend()
plt.show()
对比这两种情况：
1）
>
在第二种情况下，振荡在湿降至零之前会增加，而在第一种情况下，这些振荡非常小（大约2）。我本来希望我的算法几乎将第二个数字中的峰几乎连接起来（以提供最佳拟合）。但由于某种原因，似乎错过了前两个峰。
I am trying to find the damping rate of some oscillations that I output from a simulation.
For that, I first find the peaks of these oscillations and fit them with an exponential fit.
The method works fine for most cases but for a case I encountered yesterday I found that the method does not work properly and there was no visible reason I could find for it to not work. There seems to be some issue with the peak finding I used that somehow misses the peaks for some cases.
Here is the program for finding the peaks:
(Note: Variables in the list Tarr, Zarr, Alpharr are taken from a cell run previously)
from scipy.optimize import curve_fit
def exp(t, A, lbda):
return A * np.exp(-lbda * t)
def find_peaks(x, y):
peak_x = []
peak_vals = []
for i in range(len(y)):
if i == 0:
if y[i+1]<y[i]:
peak_x.append(x[i])
peak_vals.append(y[i])
else:
continue
if i == len(y) - 1:
continue
if (y[i-1] < y[i]) and (y[i+1] < y[i]):
peak_x.append(x[i])
peak_vals.append(y[i])
return np.array(peak_x), np.array(peak_vals)
Tarr=[Tplot0,Tplot1,Tplot2,Tplot3,Tplot4]
Zarr=[Zrel0,Zrel1,Zrel2,Zrel3,Zrel4]
Alpharr=[alpha0,alpha1,alpha2,alpha3,alpha4]
for x0,y0,alf in zip(Tarr,Zarr,Alpharr):
peak_x0,peak_y0=find_peaks(x0,y0)
popt, pcov = curve_fit(exp, peak_x0, peak_y0)
print(*[f"{val:.2f} +/- {err:.2f} __ " for val, err in zip(popt, np.sqrt(np.diag(pcov)))])
print("decay rate = ",popt[1])
plt.plot(x0,y0)
plt.plot(peak_x0,exp(peak_x0,*popt),'-r', label=str((popt[0]))+"*exp(-"+str((popt[1]))+")")
plt.xlabel("time")
plt.ylabel("Zrel ")
plt.title("n = "+str(n)+" ; alpha = "+str(alf))
plt.legend()
plt.show()
Contrast these two cases:
1)
2)
In the second case, the oscillations are increased before it damps to zero while in the first case these oscillations are very less (just around 2). I would have expected my algorithm to almost connect the peaks in the second figure (to give the best fit) .But it seems it misses the first two peaks for some reason.
发布评论
评论（1）
只是一个建议，您可能可以看一下Scipy的，
scipy.signal.find_peaks。可能您需要指定突出性值Just a suggestion, you probably can look at scipy for peak finding algorithms,
scipy.signal.find_peaks. Probably you need to specify the prominence value