import numpy as np
import matplotlib.pyplot as plt
from analog_0_1 import Device,AnalogInput,AnalogOutput
from scipy.signal.windows import get_window 
from numpy.fft import fft,ifft
from scipy.signal import firwin,convolve
import time           
            
def interpol(x,ninter):
    # interpolation par FFT
    N = len(x)
    tfd = fft(x)
    N1 = N//2
    tfd2 = np.concatenate((tfd[0:N1],np.zeros(N*ninter),tfd[N1:N]))
    y = np.real(ifft(tfd2))*(ninter+1)
    return y           
            
def dephasage(t,u0,u1,ninter,freq):
    T = 1/freq
    te = t[1]-t[0]
    if ninter>0:
        v0 = interpol(u0,ninter)
        v1 = interpol(u1,ninter)
    d = int((ninter+1)*T/(4*te))
    N = len(v0)
    x = v1[d:N]*(v0[d:N]-1j*v0[0:N-d])
    phi = np.angle(x.mean())
    return phi

            
def mesure(freq,amp,duree):
    output.function(0,'sine',freq,amp)
    output.start(0)
    time.sleep(duree)
    fEchant = freq*100
    longueur = int(duree*fEchant)
    t = analog.sampling(fEchant,longueur)
    voltage = analog.record()
    output.stop(0)
    Ve = voltage[0,:]
    Vs = voltage[1,:]
    Ve_eff = Ve.std()
    Vs_eff = Vs.std()
    G = Vs_eff/Ve_eff
    ninter = 4
    phi = dephasage(t,Ve,Vs,ninter,freq)
    return G,phi        
             
f = np.logspace(2,5,100)
N = len(f)
list_G = np.zeros(N)
list_phi = np.zeros(N)
amp = 0.3
for k in range(N):
    duree = 80/f[k]
    G,phi = mesure(f[k],amp,duree)
    print(f[k],G,phi)
    list_G[k] = G
    list_phi[k] = phi
    
GdB = 20*np.log10(list_G)
np.savetxt("bode.txt",np.array([f,list_G,list_phi]).T)
plt.figure()
plt.subplot(211)
plt.plot(f,GdB,"b-")
plt.xlabel("f")
plt.ylabel("GdB")
plt.grid()
plt.xscale("log")
plt.subplot(212)
plt.plot(f,list_phi*180/np.pi,"b-")
plt.xlabel("f")
plt.ylabel("phi")
plt.grid()
plt.xscale("log")
plt.show()