Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #Fabio Eboli 2012
- from numpy import arange,array,ones
- import pylab
- from pylab import *
- from scipy import stats
- import time
- #pylab.ion()
- fig=pylab.figure(1)
- sp1=fig.add_subplot(221)
- sp2=fig.add_subplot(222)
- sp3=fig.add_subplot(223)
- sp4=fig.add_subplot(224)
- manager=pylab.get_current_fig_manager()
- data=[0]
- xdata=[0]
- datapos=[0]
- correctionvalues=[0]
- dataunstep=[0]
- dataundrift=[0]
- correctionvalues=[0]
- driftslope=0
- driftintercept=0
- numsamples=0
- outliers=0
- p1=sp1.plot(xdata,data,"r.")
- sp1.set_title("Raw Data, outliers removed")
- p2=sp2.plot(xdata,correctionvalues,"b")
- sp2.set_title("Step correction values")
- p3=sp3.plot(xdata,dataunstep,"g")
- p31=sp3.plot([0,len(dataunstep)],[driftintercept,driftintercept+driftslope*len(dataunstep)],"g")
- sp3.set_title("Steps removed and lin. drift est")
- p4=sp4.plot(xdata,dataundrift,"y")
- sp4.set_title("Drift corrected")
- def Calculations(arg):
- global data
- global datapos
- global xdata
- global dataunstep
- global dataundrift
- global driftslope
- global driftintercept
- global correctionvalues
- global numsamples
- global outliers
- data=[]
- xdata=[]
- datapos=[]
- correctionsteps=[]
- correctionvalues=[]
- dataundrift=[]
- dataunstep=[]
- stepvalue=1e-7
- threshold=5e-8
- correctionstep=0
- outliers=0
- logfile=open("..\\TI.txt")
- header=logfile.readline()
- for lines in logfile:
- data.append(float(lines))
- logfile.close()
- numsamples=len(data)
- xdata=arange(0,numsamples)
- print(str(numsamples)+" Samples")
- print("Removing negative values...")
- for i in data:
- if i<0 :
- datapos.append(i+1e-7)
- else :
- datapos.append(i)
- print("Estimating steps correction and removing ouliers...")
- xi=arange(1,len(datapos))
- correctionsteps.append(0)
- correctionvalues.append(0)
- dataunstep.append(datapos[0])
- for i in xi:
- olddata=datapos[i-1]
- newdata=datapos[i]
- diff=newdata-olddata
- correctionstep=0
- if diff>0:
- if diff> 3*threshold:
- datapos[i]=datapos[i-1] # filetring out values
- outliers=outliers+1
- else:
- while diff>threshold:
- correctionstep=correctionstep-1
- diff=newdata+correctionstep*stepvalue-olddata
- else :
- if diff< -3*threshold:
- datapos[i]=datapos[i-1] # filetring out values
- outliers=outliers+1
- else:
- while diff<-threshold:
- correctionstep=correctionstep+1
- diff=newdata+correctionstep*stepvalue-olddata
- correctionsteps.append(correctionstep)
- correctionvalues.append(correctionvalues[i-1]+correctionstep*stepvalue)
- dataunstep.append(datapos[i]+correctionvalues[i])
- xi=arange(0,len(dataunstep))
- print("Drift estimation:")
- driftline=stats.linregress(xi,dataunstep)
- driftslope=driftline[0]
- driftintercept=driftline[1]
- print("Slope="+str(driftslope)+"; Intercept:"+str(driftintercept))
- print("Drift correction...")
- for i in xi:
- dataundrift.append(dataunstep[i]-driftslope*i)
- def Plotter(arg):
- global datapos
- global dataunstep
- global dataundrift
- global driftslope
- global driftintercept
- global correctionvalues
- global numsamples
- global outliers
- p1[0].set_data(arange(0,len(datapos)),datapos)
- sp1.axis([0,len(datapos),min(datapos),max(datapos)])
- sp1.set_title("Raw Data: "+str(numsamples)+"samples, "+str(outliers)+" outliers removed")
- p2[0].set_data(arange(0,len(correctionvalues)),correctionvalues)
- sp2.axis([0,len(correctionvalues),min(correctionvalues),max(correctionvalues)])
- sp2.set_title("Step correction values")
- p3[0].set_data(arange(0,len(dataunstep)),dataunstep)
- p31[0].set_data([0,len(dataunstep)],[driftintercept,driftintercept+driftslope*len(dataunstep)])
- sp3.axis([0,len(dataunstep),min(dataunstep),max(dataunstep)])
- sp3.set_title("Steps removed and lin. drift est. : "+str(driftslope))
- p4[0].set_data(arange(0,len(dataundrift)),dataundrift)
- sp4.axis([0,len(dataundrift),min(dataundrift),max(dataundrift)])
- sp4.set_title("Drift Removed")
- manager.canvas.draw()
- Calculations(0)
- Plotter(0)
- #time.sleep(20)
- timer1=fig.canvas.new_timer(interval=50000)
- timer1.add_callback(Calculations,())
- timer=fig.canvas.new_timer(interval=50000)
- timer.add_callback(Plotter,())
- timer.start()
- timer1.start()
- pylab.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement