bmarechal / py-ident

Download as
Browse Code ยป

Commit c21bab35d9e337135a3c2e5b9a6181fb2d633b91

Authored by bmarechal
1 parent 8fe481d742
Exists in master

replace 4 spaces by tabs

Showing 2 changed files with 11 additions and 11 deletions Inline Diff

# -*- coding: utf-8 -*- 1 1 # -*- coding: utf-8 -*-
2 2
import csv, time, glob, datetime, os 3 3 import csv, time, glob, datetime, os
import matplotlib.pyplot as plt 4 4 import matplotlib.pyplot as plt
import matplotlib.dates as md 5 5 import matplotlib.dates as md
6 6
#from scipy import signal 7 7 #from scipy import signal
#from scipy.optimize import curve_fit 8 8 #from scipy.optimize import curve_fit
9 9
tic = time.time() 10 10 tic = time.time()
11 11
os.chdir('/home/user/sicav_data/Manip/2016/2016-03/') 12 12 os.chdir('/home/user/sicav_data/Manip/2016/2016-03/')
13 13
def getColumn(filename, column): 14 14 def getColumn(filename, column):
results = [] 15 15 results = []
for dat_file in sorted(glob.glob(filename)): 16 16 for dat_file in sorted(glob.glob(filename)):
file_result = csv.reader(open(dat_file), delimiter='\t') 17 17 file_result = csv.reader(open(dat_file), delimiter='\t')
results = results + map(float,[result[column] for result in file_result]) 18 18 results = results + map(float,[result[column] for result in file_result])
return results 19 19 return results
20 20
t = getColumn('*-lakeshore.dat', 0) 21 21 t = getColumn('*-lakeshore.dat', 0)
T1 = getColumn('*-lakeshore.dat', 2) 22 22 T1 = getColumn('*-lakeshore.dat', 2)
T2 = getColumn('*-lakeshore.dat', 3) 23 23 T2 = getColumn('*-lakeshore.dat', 3)
T3 = getColumn('*-lakeshore.dat', 4) 24 24 T3 = getColumn('*-lakeshore.dat', 4)
T4 = getColumn('*-lakeshore.dat', 5) 25 25 T4 = getColumn('*-lakeshore.dat', 5)
26 26
""" 27 27 """
n = 400 28 28 n = 400
t = [t[n*i] for i in range(len(t)/n)] 29 29 t = [t[n*i] for i in range(len(t)/n)]
T1 = [T1[n*i] for i in range(len(T1)/n)] 30 30 T1 = [T1[n*i] for i in range(len(T1)/n)]
T2 = [T2[n*i] for i in range(len(T2)/n)] 31 31 T2 = [T2[n*i] for i in range(len(T2)/n)]
T3 = [T3[n*i] for i in range(len(T3)/n)] 32 32 T3 = [T3[n*i] for i in range(len(T3)/n)]
T4 = [T4[n*i] for i in range(len(T4)/n)] 33 33 T4 = [T4[n*i] for i in range(len(T4)/n)]
34 34
def func(U, a0, b1, b2, y0): 35 35 def func(U, a0, b1, b2, y0):
sys = signal.lti([a0],[b2, b1, 1]) 36 36 sys = signal.lti([a0],[b2, b1, 1])
y = sys.output(U, t, y0) 37 37 y = sys.output(U, t, y0)
return y[1] 38 38 return y[1]
39 39
print('Fitting...') 40 40 print('Fitting...')
popt, cov = curve_fit(func, T1, T4) 41 41 popt, cov = curve_fit(func, T1, T4)
print(popt) 42 42 print(popt)
Yfit = func(T1, *popt) 43 43 Yfit = func(T1, *popt)
""" 44 44 """
timetamps = [datetime.datetime.fromtimestamp(ti) for ti in t] 45 45 timetamps = [datetime.datetime.fromtimestamp(ti) for ti in t]
datenums=md.date2num(timetamps) 46 46 datenums=md.date2num(timetamps)
47 47
plt.subplots_adjust(bottom=0.35) 48 48 plt.subplots_adjust(bottom=0.35)
plt.xticks(rotation=90) 49 49 plt.xticks(rotation=90)
ax=plt.gca() 50 50 ax=plt.gca()
xfmt = md.DateFormatter('%Y-%m-%d %H:%M:%S') 51 51 xfmt = md.DateFormatter('%Y-%m-%d %H:%M:%S')
ax.xaxis.set_major_formatter(xfmt) 52 52 ax.xaxis.set_major_formatter(xfmt)
53 53
plt.plot(datenums, T1, label = 'Table') 54 54 plt.plot(datenums, T1, label = 'Table')
plt.plot(datenums, T2, label = 'Link st.') 55 55 plt.plot(datenums, T2, label = 'Link st.')
plt.plot(datenums, T3, label = 'PT2') 56 56 plt.plot(datenums, T3, label = 'PT2')
plt.plot(datenums, T4, label = 'Reg. st.') 57 57 plt.plot(datenums, T4, label = 'Reg. st.')
#plt.plot(datenums, Yfit, label = 'Fit') 58 58 #plt.plot(datenums, Yfit, label = 'Fit')
59 59
plt.legend() 60 60 plt.legend()
first_order_ident.py
Diff comments   View file @ c21bab3
# -*- coding: utf-8 -*- 1 1 # -*- coding: utf-8 -*-
2 2
from scipy import signal, linspace, pi, randn, ones 3 3 from scipy import signal, linspace, pi, randn, ones
from scipy.optimize import curve_fit 4 4 from scipy.optimize import curve_fit
import matplotlib.pyplot as plt 5 5 import matplotlib.pyplot as plt
import time 6 6 import time
7 7
tic = time.time() 8 8 tic = time.time()
9 9
'''function to optimize''' 10 10 '''function to optimize'''
def func(u, K, tau, y0): 11 11 def func(u, K, tau, y0):
sys = signal.lti(K, [tau, 1]) 12 12 sys = signal.lti(K, [tau, 1])
y = sys.output(u, t, y0) 13 13 y = sys.output(u, t, y0)
return y[1] 14 14 return y[1]
15 15
'''input square signal''' 16 16 '''input square signal'''
global t 17 17 global t
t = linspace(0, 20, 1000) 18 18 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 19 19 Umes = -0.5*(signal.square(2*pi*0.1*t)+ones(len(t)))+ones(len(t))+randn(len(t))/50
20 20
'''noisy output signal''' 21 21 '''noisy output signal'''
p = [1, 1, 0] 22 22 p = [1, 1, 0]
Ymes = func(Umes, *p)+randn(len(Umes))/50 23 23 Ymes = func(Umes, *p)+randn(len(Umes))/50
24 24
'''input and output signals plot''' 25 25 '''input and output signals plot'''
plt.plot(t, Umes, label = 'Umes') 26 26 plt.plot(t, Umes, label = 'Umes')
plt.plot(t, Ymes, label = 'Ymes') 27 27 plt.plot(t, Ymes, label = 'Ymes')
28 28
29 29
'''optimization with non-linear least squares method''' 30 30 '''optimization with non-linear least squares method'''
popt, cov = curve_fit(func, Umes, Ymes) 31 31 popt, cov = curve_fit(func, Umes, Ymes)
print(p) 32 32 print(p)
print(popt) 33 33 print(popt)
34 34
'''estimated response plot''' 35 35 '''estimated response plot'''
Yfit = func(Umes, *popt) 36 36 Yfit = func(Umes, *popt)