File:Detection pics savitzky golay.svg

#!/usr/bin/python3
# coding: utf-8

"""nom : creationDonnees.py
auteur : User:cdang
date de création : 2023-04-19
dates de modification : 
— 2023-04-20 : ajout de la courbe lissée et des cordes sur la sortie graphique
----------------------------------------------------------------------------
version de Python : 3
module requis : NumPy, Pandas
----------------------------------------------------------------------------
Objectif : crée un fichier de données, nuage de points (x, y) pour tester le programme de recherche de pics.

Entrées
-------
    Paramètres de la fonction générant le nuage de points (codé en dur).

Sorties
-------
    Un fichier CSV de deux colonnes correspondant aux coordonnées x et y des point d'une courbe
    (nombres réels convertis en chaînes de caractères).
"""

import numpy as np
import pandas as pd

# ****************
# ****************
# ** Paramètres **
# ****************
# ****************

facteurSigma = 0.03/np.sqrt(2*np.pi) # rapport entre la hauteur et la largeur d'un pic

x1 = 0.3 ; h1 = 0.5 # position et hauteur des pics
x2 = 0.75 ; h2 = 0.1
x3 = 0.8 ; h3 = 1

epsilon = 0.01 # amplitude du bruit

n = 200 # nombre de points

# chemin d'accès et nom du fichier à créer
chemin = "~/dummypath/"
nomFichier = "donnees.csv"

# **************
# **************
# ** Fonction **
# **************
# **************

def f(x, mu, sigma):
    """Fonction pour générer un pic.
    Entrées :
    — x : vecteur de réels (abscisses) ;
    — mu : réel, position du pic ;
    — sigma : réel, écart type de la fonction gaussienne.
    Sortie : y, vecteur de réels (ordonnées) de même dimension que x.
    y = 1/(σ × √(2π)) × exp((x - μ)²/(2σ²))."""
    
    unSurSigma = 1/sigma
    y = (unSurSigma*facteurSigma)*np.exp(-0.5*unSurSigma*unSurSigma*(x - mu)*(x - mu))
    
    return y

# *************************
# *************************
# ** Programme principal **
# *************************
# *************************

X = np.linspace(0, 1, n)
Y = f(X, x1, facteurSigma/h1) + f(X, x2, facteurSigma/h2) + f(X, x3, facteurSigma/h3) + epsilon*np.random.normal(X)

resultat = pd.DataFrame(np.concatenate((X.reshape(-1, 1), Y.reshape(-1, 1)), axis=1))

resultat.to_csv(chemin+nomFichier)

#!/usr/bin/python3
# coding: utf-8

"""nom : detcectionPics.py
auteur : User:cdang
date de création : 2023-04-19
dates de modification : 
— 2023-04-20 : ajout de la courbe lissée et des cordes sur la sortie graphique
— 2023-05-02 : ajout d'un filtrage sur la dérivée seconde
— 2023-05-10 : factorisation des fonctions de Savitzky-Golay ; filtrage des cordes
----------------------------------------------------------------------------
version de Python : 3
module requis : NumPy, Matplotlib, Pandas
----------------------------------------------------------------------------
Objectif : détecte les pics à partir d'une courbe, sous la forme d'un nuage de poins (x, y).

Entrées
-------
    Fichier CSV, paramètres de filtrage (codés en dur).

Sorties
-------
    Une liste de position et un graphique montrant les positions trouvées par rapport au nuage de points.
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# ****************
# ****************
# ** Constantes **
# ****************
# ****************

chemin = "~/dummypath/"
nomFichier = "donnees.csv"

fenetre = 9 # nombre de points de la fenêtre glissante pour le lissage
largeurMin = 5 # largeur minimale de la corde d'un pic, en nombre de points, pour le filtrage
hauteurMin = 0 # hauteur minimale d'un pic par rapport à sa corde, pour le filtrage
profondeurSeconde = -150 # profondeur minimale de la dérivée seconde, pour filtrage

# ***************
# ***************
# ** Fonctions **
# ***************
# ***************

# ******************
# * Savitzky-Golay *
# ******************

def lissage(y, fenetre):
    """Effectue un lissage selon l'algotithme de Savitzky-Golay, en supposant que les abscisses sont équidistantes
    (pas constant).
    Entrées :
    — y : vecteur de réels, ordonnée des points ;
    — fenetre : entier, largeur de la fenêtre glissante (nombre de points).
    Sortie : yliss : vecteur de réels de même dimension que y, ordonnées des points lissés."""
    
    if fenetre == 5:
        yliss = 1/35*(-3*(y[:-4] + y[4:]) + 12*(y[1:-3] + y[3:-1]) + 17*y[2:-2]) # 5 points
        
    elif fenetre == 7:
        yliss = 1/21*(-2*(y[:-6] + y[6:]) + 3*(y[1:-5] + y[5:-1]) + 6*(y[2:-4] + y[4:-2]) + 7*y[3:-3]) # 7 points
        
    elif fenetre == 9:
        yliss = 1/231*(-21*(y[:-8] + y[8:]) + 14*(y[1:-7] + y[7:-1]) + 39*(y[2:-6] + y[6:-2]) + 54*(y[3:-5] + y[5:-3]) +
                    59*y[4:-4]) # 9 points
    return yliss

def derivee(y, h, fenetre):
    """Effectue la dérivée selon l'algotithme de Savitzky-Golay, en supposant que les abscisses sont équidistantes
    (pas constant).
    Entrées :
    — y : vecteur de réels, ordonnée des points ;
    — h : réel, distance entre l'abscisse de deux points consécutifs (pas)
    — fenetre : entier, largeur de la fenêtre glissante (nombre de points).
    Sortie : yprime : vecteur de réels de même dimension que y, ordonnées des dérivées."""
    
    if fenetre == 5:
        yprime = 1/(10*h)*(2*(y[4:] - y[:-4]) + y[3:-1] - y[1:-3]) # 5 points
        
    elif fenetre == 7:
        yprime = 1/(28*h)*(3*(y[6:] - y[:-6]) + 2*(y[5:-1] - y[1:-5]) + y[4:-2] - y[2:-4]) # 7 points
        
    elif fenetre == 9:
        yprime = 1/(60*h)*(4*(y[8:] - y[:-8]) + 3*(y[7:-1] - y[1:-7]) + 2*(y[6:-2]- y[2:-6]) +
                           y[5:-3] - y[3:-5]) # 9 points
    return yprime

def deriveeSeconde(y, h, fenetre):
    """Effectue la dérivée seconde selon l'algotithme de Savitzky-Golay, en supposant que les abscisses sont équidistantes
    (pas constant).
    Entrées :
    — y : vecteur de réels, ordonnée des points ;
    — h : réel, distance entre l'abscisse de deux points consécutifs (pas)
    — fenetre : entier, largeur de la fenêtre glissante (nombre de points).
    Sortie : yseconde : vecteur de réels de même dimension que y, ordonnées des dérivées secondes."""
    
    if fenetre == 5:
        yseconde = 1/(7*h*h)*(2*(y[:-4] + y[4:]) - (y[1:-3] + y[3:-1]) - 2*y[2:-2]) # 5 points
        
    elif fenetre == 7:
        yseconde = 1/(42*h*h)*(5*(y[:-6] + y[6:]) - 3*(y[2:-4] + y[5:-1]) - 4*y[3:-3]) # 7 points
        
    elif fenetre == 9:
        yseconde = 1/(462*h*h)*(28*(y[:-8] + y[8:]) + 7*(y[1:-7] + y[7:-1]) - 8*(y[2:-6] + y[6:-2]) -
                                17*(y[3:-5] + y[5:-3]) - 20*y[4:-4]) # 9 points
    return yseconde

# ******************
# Détection de pic *
# ******************

def detectionPics(x, y, fenetre, largeurMin, hauteurMin):
    """Détecte les pics dans une courbe (nuage de points).
    Entrées :
    — x : vecteur de réels, abscisse des points ;
    — y : vecteur de réels de même dimension que x, ordonnée des points ;
    — fenetre : entier, largeur de la fenêtre de lissage (nombre de poins d'abscisse) ;
    — largeurMin : entier, largeur minimale de la corde d'un pic (pour éliminer les parasites) ;
    — hauteurMin : réel, hauteur minimale d'un pic par rapport à la corde (pour éliminer les parasites)
    La corde est le segment reliant les points d'inflexions sur les flancs d'un pic.
    Sorties :
    — sommets : vecteur d'entiers, index de positions des sommets dans le vecteur x ;
    — pics : vecteur d'entiers, index de positions des cordes dans le vecteur x ;
    — yliss : vecteur de réels de même dimension que x, ordonnées du nuage de point lissé ;
    — yprime : vecteur de réels de même dimension que x, dérivées ;
    — yseconde : vecteur de réels de même dimension que x, dérivées secondes."""
    
    h = np.mean(np.diff(x)) # pas d'échantillonnage
        
    # bords de la fenêtre de lissage
    decalage = int(0.5*(fenetre - 1)) # décalage dû à la largeur de la fenêtre
    
    # Préparation des données
    # lissage, dérivations
    yliss = lissage(y, fenetre) # lissage de Savitzky-Golay
    yprime = derivee(y, h, fenetre)
    yseconde = deriveeSeconde(y, h, fenetre)
    
    chgtsgn = np.sign(yseconde[:-1]*yseconde[1:]) # détecte les changements de signe de la dérivée seconde
    booleen = chgtsgn <= 0
    
    pics = [[0, 0]] # initialisation de la matrice de pics 
    
    cont = True
    n = len(booleen)
    i = 0
    
    # Détection des pics et création des cordes
    # (segments reliant les points d'inflexion)
    while cont:
        if booleen[i]: # première rencontre d'un changement de signe : ouverture d'un pic
            j = i + 1
            cont2 = (j < n)
            #print("j = "+str(j))
            while cont2: # ouvre la recherche de la fin du pic
                if booleen[j]:
                    if (j - i) >= largeurMin: # si les événements sont suffisamment distants (filtrage en largeur de pic)
                        cont2 = False # arrête la boucle
                        pics = np.append(pics, [[i, j]], axis=0) # créee un pic
                    i = j # fait avancer l'indice 
                    cont2 = False # ce n'était pas un pic, arrêt de la boucle
                else: # tant qu'il n'y a pas de nouvelle inversion de signe
                    j = j + 1 # on avance 
                cont2 = cont2 and (j < n)
        i = i + 1 # avance l'indice
        cont = (i < n) # condition d'arrêt
    
    pics = np.delete(pics, 0, 0) # retrait de la première ligne
    
    m = len(pics) # nombre de pics identifiés
    
    # Détermination de la position des sommets
    sommets = np.array([]) # initialisation
    cordes = np.array([]).reshape(0,2)
    
    for i in range(m): # pour chaque pic identifié
        posmin = np.argmin(yseconde[pics[i, 0]:pics[i, 1]]) #position du pic au minimum de la dérivée seconde
        j = pics[i, 0] + posmin # indice du sommet
        ygauche = yliss[pics[i, 0]] # ordonnée de la corde à gauche
        ydroite = yliss[pics[i, 1]] # ordonnée de la corde à droite
        ycentre = yliss[j] # ordonnée du sommet
        if ((ycentre - 0.5*(ygauche + ydroite)) >= hauteurMin) and (yseconde[j] <= profondeurSeconde):
            # filtrage en hauteur de pic par rapport au milieu de la corde, et en profondeur de la dérivée seconde
            sommets = np.append(sommets, j)
            cordes = np.append(cordes, pics[i, :].reshape((1, 2)), axis = 0)
    
    return (sommets.astype(np.int64)+decalage, cordes.astype(np.int64)+decalage, yliss, yprime, yseconde)

# *************************
# *************************
# ** Programme principal **
# *************************
# *************************

# ***********************
# * lecture des données *
# ***********************

donnees = pd.read_csv(chemin+nomFichier).to_numpy()

x = donnees[:, 1]
y = donnees[:, 2]

# **********************
# * Détection des pics *
# **********************

(sommets, cordes, yliss, yprime, ysecond) = detectionPics(x, y, fenetre, largeurMin, hauteurMin, profondeurSeconde)

# *******************
# * Sortie des pics *
# *******************

listePics = np.concatenate((x[sommets].reshape(-1, 1), yliss[sommets].reshape(-1, 1)), axis=1)
print("Liste des pics :\n")
print(listePics)

# *********
# * Tracé *
# *********

decalage = int(0.5*(fenetre - 1)) # décalage dû à la largeur de la fenêtre

fig = plt.figure(figsize = [10, 6])

for i in sommets:
    plt.plot(x[i]*np.array([1, 1]), [0, 1], color="lightgreen", linestyle="-", linewidth=0.5)

for i in range(len(cordes)):
    i1 = cordes[i, 0]
    i2 = cordes[i, 1]
    plt.plot([x[i1], x[i2]], [yliss[i1-decalage], yliss[i2-decalage]], color="lightgreen", linestyle="-", linewidth=0.5)

plt.plot(x[decalage:-decalage], yliss,
         color="black", linewidth = "1", label = "courbe lissée")
plt.plot(x, y,
         color="blue", linewidth = "0.5", label = "données brutes")

plt.xlabel("x (u.a.)")
plt.ylabel("y (u.a.)")
plt.legend()

plt.savefig(chemin+"detection_pics_savitzky_golay.svg", format="svg")