diff --git a/dat_ident.py b/dat_ident.py new file mode 100644 index 0000000..8cfd02c --- /dev/null +++ b/dat_ident.py @@ -0,0 +1,70 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 16 10:22:54 2016 + +@author: Baptiste Marechal +""" + +import csv, time, glob, datetime, os +import matplotlib.pyplot as plt +import matplotlib.dates as md + +#from scipy import signal +#from scipy.optimize import curve_fit + +tic = time.time() + +os.chdir('/home/user/sicav_data/Manip/2016/2016-03/') + +def getColumn(filename, column): + results = [] + for dat_file in sorted(glob.glob(filename)): + file_result = csv.reader(open(dat_file), delimiter='\t') + results = results + map(float,[result[column] for result in file_result]) + return results + +t = getColumn('*-lakeshore.dat', 0) +T1 = getColumn('*-lakeshore.dat', 2) +T2 = getColumn('*-lakeshore.dat', 3) +T3 = getColumn('*-lakeshore.dat', 4) +T4 = getColumn('*-lakeshore.dat', 5) + +""" +n = 400 +t = [t[n*i] for i in range(len(t)/n)] +T1 = [T1[n*i] for i in range(len(T1)/n)] +T2 = [T2[n*i] for i in range(len(T2)/n)] +T3 = [T3[n*i] for i in range(len(T3)/n)] +T4 = [T4[n*i] for i in range(len(T4)/n)] + +def func(U, a0, b1, b2, y0): + sys = signal.lti([a0],[b2, b1, 1]) + y = sys.output(U, t, y0) + return y[1] + +print('Fitting...') +popt, cov = curve_fit(func, T1, T4) +print(popt) +Yfit = func(T1, *popt) +""" +timetamps = [datetime.datetime.fromtimestamp(ti) for ti in t] +datenums=md.date2num(timetamps) + +plt.subplots_adjust(bottom=0.35) +plt.xticks(rotation=90) +ax=plt.gca() +xfmt = md.DateFormatter('%Y-%m-%d %H:%M:%S') +ax.xaxis.set_major_formatter(xfmt) + +plt.plot(datenums, T1, label = 'Table') +plt.plot(datenums, T2, label = 'Link st.') +plt.plot(datenums, T3, label = 'PT2') +plt.plot(datenums, T4, label = 'Reg. st.') +#plt.plot(datenums, Yfit, label = 'Fit') + +plt.legend() +plt.grid() +plt.show() + +toc = time.time() - tic +print(toc) diff --git a/first_order_ident.py b/first_order_ident.py new file mode 100644 index 0000000..aeb6c42 --- /dev/null +++ b/first_order_ident.py @@ -0,0 +1,49 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 11 11:04:35 2016 + +@author: Baptiste Marechal +""" + +from scipy import signal, linspace, pi, randn, ones +from scipy.optimize import curve_fit +import matplotlib.pyplot as plt +import time + +tic = time.time() + +'''function to optimize''' +def func(u, K, tau, y0): + sys = signal.lti(K, [tau, 1]) + y = sys.output(u, t, y0) + return y[1] + +'''input square signal''' +global t +t = linspace(0, 20, 1000) +Umes = -0.5*(signal.square(2*pi*0.1*t)+ones(len(t)))+ones(len(t))+randn(len(t))/50 + +'''noisy output signal''' +p = [1, 1, 0] +Ymes = func(Umes, *p)+randn(len(Umes))/50 + +'''input and output signals plot''' +plt.plot(t, Umes, label = 'Umes') +plt.plot(t, Ymes, label = 'Ymes') + + +'''optimization with non-linear least squares method''' +popt, cov = curve_fit(func, Umes, Ymes) +print(p) +print(popt) + +'''estimated response plot''' +Yfit = func(Umes, *popt) +plt.plot(t, Yfit, label ='Yfit') + +plt.legend() +plt.grid() +plt.show() + +toc = time.time() - tic +print(toc) \ No newline at end of file