I want to get the best fit of multiple files at a time. I have more than 700 files with two different feature some of them can be fitted with a single frequency equation [(amp * sin(x2piomega+phase)) * exp(-xdecay)] and the remaining files can be with the (first harmonics) equation [(amp * np.sin(x2np.piomega+phase)+ ampbnp.sin(x4np.piomega+phaseb)) * np.exp(-xdecay)]. So I want to pass all files to the first equation and other files that cannot be nicely fitted with the first equation(this means the frequency is very high or low compared to a neighbor) automatically pass to the second equation to get the best fit. My aim is to get a very small standard deviation of the mean frequency. Here I have applied the if-else condition but I could not get the expected result. So I think there is some error with the conditional command. Here is my code (I want to fit those green peaks(in figure)with second formula)
import numpy as np
from numpy import loadtxt
import lmfit
import glob
import os
def resid(params, x, ydata):
decay = params['decay'].value
phase = params['phase'].value
omega = params['omega'].value
amp = params['amp'].value
ampb=params['ampb'].value
phaseb=params['phaseb'].value
if params['omega']> 9.70 and params['omega'].value< 10.40:
y_model = (amp * np.sin(x*2*np.pi*omega+phase)) * np.exp(-x*decay)
return y_model - ydata
else:
y_model = (amp * np.sin(x*2*np.pi*omega+phase)+ ampb*np.sin(x*4*np.pi*omega+phaseb)) * np.exp(-x*decay)
return y_model - ydata
filelist = glob.glob("C:/Users/USER/Desktop/more_data/Lower_frequencies/*.txt")
#print(filelist)
for file in filelist:
data = np.loadtxt(file)
x = data[:,0]
y = data[:,1]
file_name=file
temp = os.path.splitext(file_name)[0]
params = lmfit.Parameters()
params.add('phase', 0.0, min=-np.pi, max=np.pi)
params.add('omega', 10.32, min=9, max=12)
params.add('amp', 0.2, min=0, max=5.0)
params.add('decay', 4, min=0, max=10.0)
params.add('ampb', 0.03, min=0.0, max=5)
params.add('phaseb', 0.0, min=-np.pi, max=np.pi)
fit = lmfit.minimize(resid, params, args=(x, y), method='least_squares')
params = {}
for name, param in fit.params.items():
elif name == 'omega':
params['omega'] = param.value
params['omega_stderr'] = param.stderr
e
f = open('c:/Users/USER/Desktop/outputA2.txt','a')
print(
os.path.basename(temp),
params['omega'],
params['decay'],
fit.success,
fit.nfev,
sep="\t",
file=f
)

Aucun commentaire:
Enregistrer un commentaire