Traitement du signal sous R
L'essentiel de cette page
Nouvelle page expérimentale. J'y colle mes scripts utiles à chaque fois que je fais du traitement de signal.
On y une trouve une fonction d'expansion dynamique ainsi qu'une méthodologie de lissage des courbes, d'extraction du bruit additif, du bruit multiplicatif, du signal non-spécifique pour n'en retenir que le signal spécifique discriminant.
1- Expansion dynamique : comment changer la moyenne sans saturer la valeur maximale ou changer l'ordre des données
Une expansion dynamique permet d'augmenter l'étendue d'un ensemble de valeurs.
Prenons un exemple simple à comprendre. J'ai des notes d'étudiants allant de 4 à 16. Je voudrais augmenter le contraste de ces notes en mettant la note minimale à 0 et la note maximale à 20.
Prenons un exemple plus concret : j'ai une image, je voudrais augmenter l'étendues des valeurs des pixels pour qu'ils aillent du noir au blanc.
On peut ainsi imposer une valeur de min ou de max et conserver la position relative des valeurs.
Dans un autre exemple, on peut aussi vouloir jouer sur la moyenne.
Si je reprends l'exemple d'étudiants. Les notes sont trop basse (8 de moyennes). Je peux vouloir augmenter ces notes pour passer à 11 de moyenne mais je ne veux pas me contenter d'ajouter 3 points car, la meilleure note de 18 pourrait se retrouver à 21 ou encore parce que je veux que la note min reste là où elle est.
Je vais donc faire une expansion dynamique qui me permet de fixer la moyenne, fixer le min et le max tout en conservant les positions relatives des valeurs et s'en saturer les valeurs.
1.1- Des valeurs pour l'exemple, imaginons des notes d'étudiants mais ce pourrait être tout signal dont on veut changer la moyenne et les contrastes sans plafonner le signal
x <- round(rnorm(120,8,2),0) ; x <- x[x>3&x<=8]
y <- round(rnorm(120,10,3),0) ; y <- y[y>10&y<=20]-2
notes <- c(x,y)
hist(notes) # les notes vont bien de 4 à 16 avec une moyenne de 8
1.2- Ma fonction d'expansion dynamique à copier-coller
exp_dyn <- function(liste,moyenne,max) {
moy_liste <- (mean(liste,na.rm=T)+median(liste,na.rm=T))/2
min <- min(liste,na.rm=T)
liste_temp <- liste
liste[(liste_temp<=moy_liste)&!is.na(liste)] <- (liste[liste_temp<=moy_liste&!is.na(liste)]-min)/max((liste[liste_temp<=moy_liste&!is.na(liste)]-min),na.rm=T)*(moyenne-min)+min
#liste[liste_temp<moy_liste&!is.na(liste)] <- (liste[liste_temp<moy_liste&!is.na(liste)]-min)/(moy_list-min)*(moyenne-min)+min
liste[liste_temp==moy_liste&!is.na(liste)] <- moyenne
liste[liste_temp>moy_liste&!is.na(liste)] <- (liste[liste_temp>moy_liste&!is.na(liste)]-moy_liste)/max((liste[liste_temp>=moy_liste&!is.na(liste)]-moy_liste),na.rm=T)*(max-moyenne)+moyenne
liste <- round(liste,0)
return(liste)
}
1.3- Un exemple d'application de la fonction exp_dyn()
notes_bis <- exp_dyn(notes,10,20)
library("vioplot")
vioplot(notes,notes_bis)
On voit bien sur cette figure que la moyenne est passée à 10 et le maximum à 20 sans changer la valeur minimale.
La fonction pourrait en revanche être améliorée puisque le minimum n'a pas été imposé. Il suffirait toutefois de faire une soustraction aux valeurs avant de faire l'expansion.
2- Corriger le bruit de fond additif et les variations multiplicatives du signal
Decriptif : normalisation du background, correction du signal additif et multiplicatif afin d'illustrer comment préserver la variance inter-groupe en supprimant la variance intra-groupe.
2.1- Ouverture d'un jeu de données, spectres réalisés sur des truffes.
Ouverture d'un jeu de données, spectres réalisés sur des truffes (Merci à Sullivan Renouard, maître de conférence pour m'avoir autorisé à faire usage de ses données).
Si je trace ces spectres sur un même graphique, on voit que la ligne de base n'est pas la même pour chacun et que les spectres d'une même espèces (même couleur) ne se superposent pas (certainement à cause d'une variation quantitative), ce qui engendre un bruit multiplicatif.
library(openxlsx)
data <- read.xlsx(file.choose())
data[,1]<- as.factor(data[,1])
# Une fonction pour voir les spectres
view <- function(data,ylim=c(0,0.2)) {
for (i in c(1:nrow(data))) {
temp <- as.numeric(unlist(data[i,-1]))
if (i==1) {plot(temp,type="l",col=data[i,1],ylim=ylim,xlab="lambda",ylab="Détection")
} else {points(temp,type="l",col=data[i,1])}
}
}
# Visualisation
view(data)
2.2- Mise à 0 de la ligne de base, suppression du bruit de fond ou background
On voit que la ligne de base apparaît à gauche (horitzontales), sur les premières valeurs. On va donc soustraire à chaque spectre, la valeur médiane de ces premières valeurs.
On prend la valeur médiane pour éviter les fluctuation d'une longueur d'onde à l'autre par "grésillement".
data2 <- data
for (i in c(1:nrow(data2))) {
temp <- as.numeric(unlist(data2[i,-1]))
# Médiane sur les 200 premières valeurs
temp <- temp-median(temp[2:201],na.rm=T)
data2[i,-1] <- temp
}
colnames(data2)
view(data2)
2.3- Suppression des effets de quantité (bruit multiplicatif)
Objectif : réduire les différences entre individus d'un même groupe (liées au fluctuations de la quantité) sans altérer les modification entre groupes différents
Le principe d'un tel traitement est le suivant : tester plusieurs plages de longueurs d'onde (on ne peut pas tester toutes les longueurs d'onde) dont on va récupérer la valeur médiane. Chaque spectre sera divisé par cette valeur. La meilleure plage sera cela qui n'altère pas la variabilité inter-espèce, mais qui réduit la variabilité intra-espèce.
On voit ici qu'on y est bien arrivé : toutes les courbes d'individus d'une même espèce se sont regroupées grâce au code ci-après.
Une analyse ACP fine nous permettra ainsi de passer de données où l'on ne distinguait pas les espèces à des données très discriminantes même si on utilise des échantillons plus ou moins concentrés (ci-dessous).
Code pour éliminer le bruit multiplicatif
1) Chercher la plage idéale pour faire une correction quantitative
nlambda <- ncol(data)-1;variance <- c()
decoupe=50
data3 <- data2
cat("Nombre de plages (lambdas) à traiter : ",nlambda/decoupe,"\n")
start = 0 # On peut forcer start à partir plus ou moins loin pour ne pas tester tout le spectre (calcul long)
k<-0
for (lambda in start:(round(nlambda/decoupe,0)-1)) {
k <- k+1
cat("Lambda",lambda,"pour k :",k,"\n")
data3 <- data2
for (i in c(1:nrow(data2))) {
temp <- as.numeric(unlist(data3[i,-1]))
# On divise l'ensemble de la courbe par une médiane locale de 50 (decoupe)
temp <- temp/median(temp[(lambda*decoupe+2):(lambda*decoupe+decoupe+1)],na.rm=T)
data3[i,-1] <- temp
}
variance[k] <- 0
#########
# Normaliser la variance globale inter-espèces
########
# Pour chaque espèce de 1 à 6, on va calculer une courbe moyenne
tab <- apply(data3[data3[,1]==unique(data3[,1])[1],-1],2,mean)
for (j in 2:length(unique(data3[,1]))) {
temp <- apply(data3[data3[,1]==unique(data3[,1])[j],-1],2,mean)
tab <- rbind(tab,temp)
}
# On va calculer la variance de ces 6 courbes : c'est la variance inter-espèce.
temp <- apply(tab,2,var)
# On va faire la somme de cette variance pour chaque lambda : variance global inter-espèce
temp <- sum(tab)
# L'ensemble du tableau divisé par une médiane locale sera divisé par la variance inter-espèce globale.
data3[,-1]<- data3[,-1]/temp
# On va ensuite calculer la variance cumulée pour chaque espèce (intra-espèce) et les addition pour avoir une variance intra globale
for (j in 1:6) {
temp <- apply(data3[data3[,1]==unique(data3[,1])[j],-1],2,var)
temp <- sum(temp)
variance[k] <- variance[k]+temp # Liste des variances des 6 espèces intra globale : pour chaque fenêtre lambda (k)
}
}
plot(variance,log="y",pch=16)
# On va chercher ainsi le paramètre k et en déduire la région de 50 lambdas qui minimise la variance intra de chaque espèce.
which(variance==min(variance,na.rm=T))-> mon_k ; print(mon_k)
lambda <- 80
Domaine <- (lambda *decoupe+2):(lambda *decoupe+decoupe+1)
print(Domaine)
#############################################################
# Correction quantitative
##############################################################
data3 <- data2
for (i in c(1:nrow(data2))) {
temp <- as.numeric(unlist(data3[i,-1]))
temp <- temp/median(temp[Domaine],na.rm=T)
data3[i,-1] <- temp
}
view(data3,ylim=c(0,1.7))
2.4- Suppression du profil de l'espèce pour ne retenir que le signal spécifique aux sous-groupes (sous-espèces)
Les courbes ci-dessus sont trop bruités, on peut donc appliquer une fonction de lissage.
Profil type d'une truffe
Profil type d'une truffe (après lissage)
Code à copier-coller pour obtenir ce profil de commun
data4 <- data3
apply(data4[,-1],2,median)->truffade
plot(truffade,type="l")
# Fonction de lissage
lissage <- function(temp,decoupe=15){
section = ((decoupe-decoupe%%2)/2)
temp2 = temp[1:section]
for (i in ((section+1):(length(temp)-section))) {
temp2 <- c(temp2,median(temp[(i-section):(i+section)]))
}
temp2 <- c(temp2,temp[(length(temp)-section+1):length(temp)])
return(temp2)
}
truffade <- lissage(truffade) ; plot(truffade,type="l")
On va ainsi pouvoir soustraire ce profil à celui de toutes les truffes afin d'identifier le signal propre à chaque truffe :
On a maintenant le signal propre à chaque espèce de truffe
data4 <- apply(data4[,-1],1,function(x){return(x-truffade)})
data4 <- t(data4)
data4 <- cbind(data3[,1],data4)
view(data4)
data4 <- apply(data4[,-1],1,lissage)
data4 <- t(data4)
data4 <- cbind(data3[,1],data4)
view(data4)
2.5- Identifier les domaines discriminants entre-espèces en appliquant une courbe moyenne à chaque espèce et l'intervalle de confiance
Parmi les cournes que nous avons ci-dessus, il y a le signal propre à toutes les truffes et celui propres à chaque sous-espèce de truffe. Identifions le signal propre au truffe par un lissage median :
Cette courbe a été exprimé en fonction des lambdas réels
Pour chaque espèce, nous avons la courbe médiane et la zone de l'intervalle de confiance.
Ainsi, à chaque lambda où une espèce ressort sur les autre est un lambda discriminant pour identifier cette espèce.
Avec une telle approche, on peut ainsi obtenir des données très discriminantes en ACP :
Chaque espèce se retrouve sur un droite :
Code ci-dessous pour reproduire ces résultats.
Code pour réaliser ce graphique à copier-coller
# COURBE MOYENNE LISSEE POUR CHAQUE ESPECE
data5 <- apply(data4[data4[,1]==unique(data4[,1])[1],-1],2,median)
for (j in 2:length(unique(data4[,1]))) {
temp <- apply(data4[data4[,1]==unique(data4[,1])[j],-1],2,median)
data5 <- rbind(data5,temp)
}
data5 <- data.frame("species"=unique(data4[,1]),data5)
data5 <- apply(data5[,-1],1,lissage)
data5 <- t(data5)
data5 <- data.frame("species"=unique(data4[,1]),data5)
view(data5,ylim=c(-0.11,0.14))
# INTERVALLE DE CONFIANCE DE CHAQUE ESPECE
library(KefiR)
data6 <- apply(data4[data4[,1]==unique(data4[,1])[1],-1],2,int.ech)
for (j in 2:length(unique(data4[,1]))) {
temp <- apply(data4[data4[,1]==unique(data4[,1])[j],-1],2,int.ech)
data6 <- rbind(data6,temp)
}
data6 <- data.frame("species"=unique(data4[,1]),data6)
view(data6,ylim=c(0,0.2))
# AFFICHAGE
lambdas <- as.numeric(colnames(data[,-1]))
colors = c("#F9180A","#5BF90A","#0A1EF9","#F9C20A","#E80AF9","#0AF9F4")
colorsalpha = paste0(colors,"BB")
plot(lambdas,unlist(data5[1,-1]),col=colors[1],type="l",ylim=c(-0.13,0.15))
polygon(c(lambdas,rev(lambdas)),
c( unlist(data5[1,-1])-unlist(data6[1,-1]),
rev(unlist(data5[1,-1])+unlist(data6[1,-1]))),
col=colorsalpha[1],border=NA)
for (sp in 2:length(unique(data[,1]))) {
points(lambdas,unlist(data5[sp,-1]),col=colors[sp],type="l")
polygon(c(lambdas,rev(lambdas)),
c( unlist(data5[sp,-1])-unlist(data6[sp,-1]),
rev(unlist(data5[sp,-1])+unlist(data6[sp,-1]))),
col=colorsalpha[sp],border=NA)
}
prcomp(data4)->myacp
plot(myacp$x[,1:2],col=colors[data[,1]],cex=2,pch=16)
library(rgl)
plot3d(myacp$x[,1:3],col=colors[data[,1]],cex=2,pch=16)