from scipy import signal

def calcul_puissance(signal1,freq,fenetre) :
    
    ##signal est un tableau de valeurs contenant l'échantillonage du signal, freq est la frequence à laquelle cet echantillonnage a été effectué
    ##fenetre est la fenêtre de fréquences considérées comme pathologiques (tuple)
    
    freqs, psd = signal.welch(signal1,freq,nperseg=len(signal1))
    
    #On vient de calculer le Power Density Spectrum du signal
    
    pas=freqs[1]-freqs[0] #va nous permettre de calculer l'intégrale
    
    #On calcule maintenant la puissance totale du signal, et la puissance générée par les fréquences situées à l'intérieur de la fenêtre choisie
    
    integr_tot = 0
    integr_fenetre = 0
    
    for i in range(len(freqs)-1) :
        integr_tot+=pas*(psd[i]+psd[i+1])/2
        if freqs[i]>=fenetre[0] and freqs[i+1]<=fenetre[1] : #si on est dans la fenêtre
            integr_fenetre+=pas*(psd[i]+psd[i+1])/2
    
    #On renvoie les deux 
    
    return integr_fenetre,integr_tot

def calcul_Q_plusieurs_signaux(T1,T2,T3,freq,fenetre) :
    
    #Dans notre cas, on a 3 signaux à prendre en compte (accélération selon les 3 directions).
    #On va donc calculer les puissances selon chacune des 3 directions, puis moyenner).
    
    puissances_1=calcul_puissance(T1,freq,fenetre)
    puissances_2=calcul_puissance(T2,freq,fenetre)
    puissances_3=calcul_puissance(T3,freq,fenetre)
    
    return (puissances_1[0]+puissances_2[0]+puissances_3[0])/(puissances_1[1]+puissances_2[1]+puissances_3[1])

def calcul_puissance_deuxieme_methode(signal1,freq,fenetre) :
    
    ##signal est un tableau de valeurs contenant l'échantillonage du signal, freq est la frequence à laquelle cet echantillonnage a été effectué
    ##fenetre est la fenêtre de fréquences considérées comme pathologiques (tuple)
    
    freqs, psd = signal.welch(signal1,freq,nperseg=len(signal1))
    
    #On vient de calculer le Power Density Spectrum du signal
    #On va maintenant chercher à quelle fréquence a lieu son pic de puissance.
    
    pmax=0
    freqmax=0
    for i in range(len(freqs)) :
        if freqs[i]>=fenetre[0] and freqs[i]<=fenetre[1] :
            if psd[i]>pmax :
                pmax=psd[i]
                freqmax=freqs[i]
    
    pas=freqs[1]-freqs[0] #va nous permettre de calculer les intégrales
    
    integr_pic = 0 #puissance développée par les fréquences autour du pic
    integr_fenetre = 0 #puissance développée par les fréquences dans la bande pathologique
    
    for i in range(len(freqs)-1) :
        if freqs[i]>=fenetre[0] and freqs[i+1]<=fenetre[1] : #si on est dans la fenêtre
            integr_fenetre+=pas*(psd[i]+psd[i+1])/2
            if freqs[i]>=freqmax-1 and freqs[i+1]<=freqmax+1 :
                integr_pic+=pas*(psd[i]+psd[i+1])/2
                
    return integr_pic,integr_fenetre

def calcul_Q_plusieurs_signaux_deuxieme_methode(T1,T2,T3,freq,fenetre) :
    puissances_1=calcul_puissance_deuxieme_methode(T1,freq,fenetre)
    puissances_2=calcul_puissance_deuxieme_methode(T2,freq,fenetre)
    puissances_3=calcul_puissance_deuxieme_methode(T3,freq,fenetre)
    
    return (puissances_1[0]+puissances_2[0]+puissances_3[0])/(puissances_1[1]+puissances_2[1]+puissances_3[1])

def calcul_Qs(T1,T2,T3,freq,fenetre) :
    #On renvoie les résultats selon les deux méthodes de calcul
    
    return calcul_Q_plusieurs_signaux(T1,T2,T3,freq,fenetre),calcul_Q_plusieurs_signaux_deuxieme_methode(T1,T2,T3,freq,fenetre)