# -*- coding: utf-8 -*-
"""
Created on Wed Aug 17 00:37:27 2016

@author: david
"""

from numpy import *
from matplotlib.pyplot import *
from pylab import plot, show, title, xlabel, ylabel, subplot
from scipy import fft, arange
import serial


#http://stackoverflow.com/questions/11686720/is-there-a-numpy-builtin-to-reject-outliers-from-a-list
def reject_outliers(data, m = 2.):
    d = np.abs(data - np.median(data))
    mdev = np.median(d)
    s = d/mdev if mdev else 0.
    return data[s<m]




#http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html

whichmode=0

if(whichmode==0): # read data from serial device
    ser = serial.Serial('/dev/ttyACM0', 115200)
    size=1000000
    arr=empty(size, dtype=float)
    i=0
    fails=0
#    while i in range(0,size):
    for x in np.nditer(arr, op_flags=['readwrite']):
        a=ser.readline()
        a=a.strip()
        if a.isdigit():
            x[...]=float(a)
            i=i+1
        else:
            fails=fails+1
            print a
	if(mod(i,1000)==0):
		print i

#    if(len(a)>2 and i>5): 
#            arr[i]=int(a)
#    print ser.readline()
    print fails 
    savetxt("mydata",arr)

    
    
elif(whichmode==1): # load LT Spice output

    threshold=1.7

    transition=empty(0)

    data=loadtxt(open("/home/david/Documents/python/put04bb50.txt","rb"),skiprows=1)

    print data.shape[0]
    print data.shape[1]

    for i in range(0,data.shape[0]-2):
        if(data[i,1]>=threshold and data[i+1,1]<threshold):
            transition=append(transition,data[i,0])

    print transition.shape[0]
        
    arr=empty(0)
    
    for i in range(0,transition.shape[0]-1):
        delta=transition[i+1]-transition[i]
#        if(delta>=0.005 and delta<0.015):
        arr=append(arr,delta)
        
    arr=arr*16*1000000
    
else: # load saved data

    arr=loadtxt(open("/home/david/Documents/python/PutData/4V/mydata","rb"))



arr=reject_outliers(arr,m=7)

arr=arr/16

median=np.median(arr)
mean=np.mean(arr)
sd=np.std(arr)

xstr="median %d mean %d sd %d delta %f" % (median,mean,sd,sd/mean)

figure()

suptitle(xstr)

hist(arr, bins=50)


farr=arr/1000000
farr=1/farr


figure()

hist(farr, bins=50)



x=arr[1::1]
x=append(x,x[0])
y=arr[0::1]


figure()
scatter(x, y)

figure()

hist2d(x,y,100)

#http://stackoverflow.com/questions/20105364/how-can-i-make-a-scatter-plot-colored-by-density-in-matplotlib


figure()


bins=20

wave=np.zeros(bins)
wavecount=np.zeros(bins)

now=0
periodicity=20000
periodperbin=periodicity/bins

for i in range(0,arr.shape[0]):
    nowp=now#+period[i]/2
    now=now+arr[i]
    wave[int(mod(nowp,periodicity)/periodperbin)]+=arr[i]
    wavecount[int(mod(nowp,periodicity)/periodperbin)]+=1


for i in range(0,bins):
    wave[i]/=wavecount[i]


plot(wave)

#http://glowingpython.blogspot.co.uk/2011/08/how-to-plot-frequency-spectrum-with.html        


figure()

y=arr-median

Fs=1000000/median

n = len(y) # length of the signal
k = arange(n)
T = n/Fs
frq = k/T # two sides frequency range
frq = frq[range(n/2)] # one side frequency range

Y = fft(y)/n # fft computing and normalization
Y = Y[range(n/2)]

frq=frq[10:]
Y=Y[10:]
 
plot(frq,abs(Y)) # plotting the spectrum