Détection du signal de crête dans les données de séries temporelles en temps réel


Mise à jour: L’algorithme le plus performant à ce jour est celui-ci .


Cette question explore des algorithmes robustes permettant de détecter les pics soudains dans les données de séries temporelles en temps réel.

Considérons le jeu de données suivant:

p = [1 1 1.1 1 0.9 1 1 1.1 1 0.9 1 1.1 1 1 0.9 1 1 1.1 1 1 1 1 1.1 0.9 1 1.1 1 1 0.9 1, ... 1.1 1 1 1.1 1 0.8 0.9 1 1.2 0.9 1 1 1.1 1.2 1 1.5 1 3 2 5 3 2 1 1 1 0.9 1 1 3, ... 2.6 4 3 3.2 2 1 1 0.8 4 4 2 2.5 1 1 1]; 

(Format Matlab mais il ne s’agit pas de la langue mais de l’algorithme)

Tracé de données

Vous pouvez clairement voir qu’il y a trois grands pics et quelques petits pics. Cet dataset est un exemple spécifique de la classe de jeux de données de la série de temps dont il s’agit. Cette classe de jeux de données présente deux caractéristiques générales:

  1. Il y a du bruit de base avec une moyenne générale
  2. Il existe de grands « pics » ou « points de données supérieurs » qui s’écartent considérablement du bruit.

Supposons également ce qui suit:

  • la largeur des pics ne peut pas être déterminée au préalable
  • la hauteur des pics s’écarte clairement des autres valeurs
  • l’algorithme utilisé doit calculer en temps réel (changez donc avec chaque nouveau sharepoint données)

Pour une telle situation, il faut construire une valeur limite qui déclenche des signaux. Cependant, la valeur limite ne peut pas être statique et doit être déterminée en temps réel sur la base d’un algorithme.


Ma question: qu’est-ce qu’un bon algorithme pour calculer ces seuils en temps réel? Existe-t-il des algorithmes spécifiques pour de telles situations? Quels sont les algorithmes les plus connus?


Des algorithmes robustes ou des informations utiles sont tous très appréciés. (peut répondre dans n’importe quelle langue: il s’agit de l’algorithme)

Algo z-score lissé (détection de pic avec seuil robuste)

J’ai construit un algorithme qui fonctionne très bien pour ces types de jeux de données. Il est basé sur le principe de la dispersion : si un nouveau sharepoint donnée est un nombre x donné d’écart-types par rapport à une moyenne mobile, l’algorithme signale (également appelé z-score ). L’algorithme est très robuste car il construit une moyenne mobile et une déviation distinctes , de sorte que les signaux ne corrompent pas le seuil. Les signaux futurs sont donc identifiés avec approximativement la même précision, indépendamment de la quantité de signaux précédents. L’algorithme prend 3 entrées: lag = the lag of the moving window , threshold = the z-score at which the algorithm signals et influence = the influence (between 0 and 1) of new signals on the mean and standard deviation . Par exemple, un lag de 5 utilisera les 5 dernières observations pour lisser les données. Un threshold de 3,5 signalera si un sharepoint données est à 3,5 écarts-types de la moyenne mobile. Et une influence de 0,5 donne aux signaux la moitié de l’influence des points de données normaux. De même, une influence de 0 ignore complètement les signaux pour recalculer le nouveau seuil. Une influence de 0 est donc l’option la plus robuste (mais suppose une stationnarité ); mettre l’option d’influence à 1 est moins robuste. Pour les données non stationnaires, l’option d’influence doit donc être placée entre 0 et 1.

Cela fonctionne comme suit:

Pseudocode

 # Let y be a vector of timeseries data of at least length lag+2 # Let mean() be a function that calculates the mean # Let std() be a function that calculates the standard deviaton # Let absolute() be the absolute value function # Settings (the ones below are examples: choose what is best for your data) set lag to 5; # lag 5 for the smoothing functions set threshold to 3.5; # 3.5 standard deviations for signal set influence to 0.5; # between 0 and 1, where 1 is normal influence, 0.5 is half # Initialise variables set signals to vector 0,...,0 of length of y; # Initialise signal results set filteredY to y(1),...,y(lag) # Initialise filtered series set avgFilter to null; # Initialise average filter set stdFilter to null; # Initialise std. filter set avgFilter(lag) to mean(y(1),...,y(lag)); # Initialise first value set stdFilter(lag) to std(y(1),...,y(lag)); # Initialise first value for i=lag+1,...,t do if absolute(y(i) - avgFilter(i-1)) > threshold*stdFilter(i-1) then if y(i) > avgFilter(i-1) then set signals(i) to +1; # Positive signal else set signals(i) to -1; # Negative signal end # Make influence lower set filteredY(i) to influence*y(i) + (1-influence)*filteredY(i-1); else set signals(i) to 0; # No signal set filteredY(i) to y(i); end # Adjust the filters set avgFilter(i) to mean(filteredY(i-lag),...,filteredY(i)); set stdFilter(i) to std(filteredY(i-lag),...,filteredY(i)); end 

Les règles empiriques pour la sélection de bons parameters pour vos données se trouvent à l’ annexe 3 (ci-dessous).


Démo

Démonstration de l'algorithme de seuillage robuste

Le code Matlab pour cette démo se trouve à la fin de cette réponse. Pour utiliser la démo, lancez-la simplement et créez vous-même une série chronologique en cliquant sur le graphique supérieur. L’algorithme commence à fonctionner après avoir dessiné le nombre de lag des observations.


Annexe 1: Code Matlab et R pour l’algorithme

Code matlab

 function [signals,avgFilter,stdFilter] = ThresholdingAlgo(y,lag,threshold,influence) % Initialise signal results signals = zeros(length(y),1); % Initialise filtered series filteredY = y(1:lag+1); % Initialise filters avgFilter(lag+1,1) = mean(y(1:lag+1)); stdFilter(lag+1,1) = std(y(1:lag+1)); % Loop over all datapoints y(lag+2),...,y(t) for i=lag+2:length(y) % If new value is a specified number of deviations away if abs(y(i)-avgFilter(i-1)) > threshold*stdFilter(i-1) if y(i) > avgFilter(i-1) % Positive signal signals(i) = 1; else % Negative signal signals(i) = -1; end % Make influence lower filteredY(i) = influence*y(i)+(1-influence)*filteredY(i-1); else % No signal signals(i) = 0; filteredY(i) = y(i); end % Adjust the filters avgFilter(i) = mean(filteredY(i-lag:i)); stdFilter(i) = std(filteredY(i-lag:i)); end % Done, now return results end 

Exemple:

 % Data y = [1 1 1.1 1 0.9 1 1 1.1 1 0.9 1 1.1 1 1 0.9 1 1 1.1 1 1,... 1 1 1.1 0.9 1 1.1 1 1 0.9 1 1.1 1 1 1.1 1 0.8 0.9 1 1.2 0.9 1,... 1 1.1 1.2 1 1.5 1 3 2 5 3 2 1 1 1 0.9 1,... 1 3 2.6 4 3 3.2 2 1 1 0.8 4 4 2 2.5 1 1 1]; % Settings lag = 30; threshold = 5; influence = 0; % Get results [signals,avg,dev] = ThresholdingAlgo(y,lag,threshold,influence); figure; subplot(2,1,1); hold on; x = 1:length(y); ix = lag+1:length(y); area(x(ix),avg(ix)+threshold*dev(ix),'FaceColor',[0.9 0.9 0.9],'EdgeColor','none'); area(x(ix),avg(ix)-threshold*dev(ix),'FaceColor',[1 1 1],'EdgeColor','none'); plot(x(ix),avg(ix),'LineWidth',1,'Color','cyan','LineWidth',1.5); plot(x(ix),avg(ix)+threshold*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5); plot(x(ix),avg(ix)-threshold*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5); plot(1:length(y),y,'b'); subplot(2,1,2); stairs(signals,'r','LineWidth',1.5); ylim([-1.5 1.5]); 

Code R

 ThresholdingAlgo <- function(y,lag,threshold,influence) { signals <- rep(0,length(y)) filteredY <- y[0:lag] avgFilter <- NULL stdFilter <- NULL avgFilter[lag] <- mean(y[0:lag]) stdFilter[lag] <- sd(y[0:lag]) for (i in (lag+1):length(y)){ if (abs(y[i]-avgFilter[i-1]) > threshold*stdFilter[i-1]) { if (y[i] > avgFilter[i-1]) { signals[i] <- 1; } else { signals[i] <- -1; } filteredY[i] <- influence*y[i]+(1-influence)*filteredY[i-1] } else { signals[i] <- 0 filteredY[i] <- y[i] } avgFilter[i] <- mean(filteredY[(i-lag):i]) stdFilter[i] <- sd(filteredY[(i-lag):i]) } return(list("signals"=signals,"avgFilter"=avgFilter,"stdFilter"=stdFilter)) } 

Exemple:

 # Data y <- c(1,1,1.1,1,0.9,1,1,1.1,1,0.9,1,1.1,1,1,0.9,1,1,1.1,1,1,1,1,1.1,0.9,1,1.1,1,1,0.9, 1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3, 2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1) lag <- 30 threshold <- 5 influence <- 0 # Run algo with lag = 30, threshold = 5, influence = 0 result <- ThresholdingAlgo(y,lag,threshold,influence) # Plot result par(mfrow = c(2,1),oma = c(2,2,0,0) + 0.1,mar = c(0,0,2,1) + 0.2) plot(1:length(y),y,type="l",ylab="",xlab="") lines(1:length(y),result$avgFilter,type="l",col="cyan",lwd=2) lines(1:length(y),result$avgFilter+threshold*result$stdFilter,type="l",col="green",lwd=2) lines(1:length(y),result$avgFilter-threshold*result$stdFilter,type="l",col="green",lwd=2) plot(result$signals,type="S",col="red",ylab="",xlab="",ylim=c(-1.5,1.5),lwd=2) 

Ce code (les deux langues) donnera le résultat suivant pour les données de la question originale:

Exemple de seuil du code Matlab


Implémentations dans d'autres langues:

  • Golang (Xeoncross)

  • Python (R Kiselev)

  • Swift (moi)

  • Groovy (JoshuaCWebDeveloper)

  • C ++ (brad)

  • C ++ (Animesh Pandey)

  • rust (swizard)

  • Scala (Mike Roberts)

  • Kotlin (leoderprofi)

  • Ruby (Kimmo Lehto)

  • Fortran [pour la détection de résonance] (THo)


Annexe 2: Code de démonstration Matlab (cliquez pour créer des données)

 function [] = RobustThresholdingDemo() %% SPECIFICATIONS lag = 5; % lag for the smoothing threshold = 3.5; % number of st.dev. away from the mean to signal influence = 0.3; % when signal: how much influence for new data? (between 0 and 1) % 1 is normal influence, 0.5 is half %% START DEMO DemoScreen(30,lag,threshold,influence); end function [signals,avgFilter,stdFilter] = ThresholdingAlgo(y,lag,threshold,influence) signals = zeros(length(y),1); filteredY = y(1:lag+1); avgFilter(lag+1,1) = mean(y(1:lag+1)); stdFilter(lag+1,1) = std(y(1:lag+1)); for i=lag+2:length(y) if abs(y(i)-avgFilter(i-1)) > threshold*stdFilter(i-1) if y(i) > avgFilter(i-1) signals(i) = 1; else signals(i) = -1; end filteredY(i) = influence*y(i)+(1-influence)*filteredY(i-1); else signals(i) = 0; filteredY(i) = y(i); end avgFilter(i) = mean(filteredY(i-lag:i)); stdFilter(i) = std(filteredY(i-lag:i)); end end % Demo screen function function [] = DemoScreen(n,lag,threshold,influence) figure('Position',[200 100,1000,500]); subplot(2,1,1); title(sprintf(['Draw data points (%.0f max) [settings: lag = %.0f, '... 'threshold = %.2f, influence = %.2f]'],n,lag,threshold,influence)); ylim([0 5]); xlim([0 50]); H = gca; subplot(2,1,1); set(H, 'YLimMode', 'manual'); set(H, 'XLimMode', 'manual'); set(H, 'YLim', get(H,'YLim')); set(H, 'XLim', get(H,'XLim')); xg = []; yg = []; for i=1:n try [xi,yi] = ginput(1); catch return; end xg = [xg xi]; yg = [yg yi]; if i == 1 subplot(2,1,1); hold on; plot(H, xg(i),yg(i),'r.'); text(xg(i),yg(i),num2str(i),'FontSize',7); end if length(xg) > lag [signals,avg,dev] = ... ThresholdingAlgo(yg,lag,threshold,influence); area(xg(lag+1:end),avg(lag+1:end)+threshold*dev(lag+1:end),... 'FaceColor',[0.9 0.9 0.9],'EdgeColor','none'); area(xg(lag+1:end),avg(lag+1:end)-threshold*dev(lag+1:end),... 'FaceColor',[1 1 1],'EdgeColor','none'); plot(xg(lag+1:end),avg(lag+1:end),'LineWidth',1,'Color','cyan'); plot(xg(lag+1:end),avg(lag+1:end)+threshold*dev(lag+1:end),... 'LineWidth',1,'Color','green'); plot(xg(lag+1:end),avg(lag+1:end)-threshold*dev(lag+1:end),... 'LineWidth',1,'Color','green'); subplot(2,1,2); hold on; title('Signal output'); stairs(xg(lag+1:end),signals(lag+1:end),'LineWidth',2,'Color','blue'); ylim([-2 2]); xlim([0 50]); hold off; end subplot(2,1,1); hold on; for j=2:i plot(xg([j-1:j]),yg([j-1:j]),'r'); plot(H,xg(j),yg(j),'r.'); text(xg(j),yg(j),num2str(j),'FontSize',7); end end end 

Annexe 3: Règles de base pour la configuration de l'algorithme

lag : le paramètre lag détermine le degré de lissage de vos données et l'adaptation de l'algorithme aux modifications de la moyenne à long terme des données. Plus vos données sont stationnaires , plus vous devez inclure des décalages (cela devrait améliorer la robustesse de l'algorithme). Si vos données contiennent des tendances variables dans le temps, vous devez tenir compte de la rapidité avec laquelle vous souhaitez que l'algorithme s'adapte à ces tendances. Par exemple, si vous placez un lag à 10, cela prend 10 «périodes» avant que le seuil de l'algorithme ne soit ajusté à des changements systématiques de la moyenne à long terme. Choisissez donc le paramètre de lag fonction du comportement de tendance de vos données et de l'adaptation de l'algorithme.

influence : ce paramètre détermine l'influence des signaux sur le seuil de détection de l'algorithme. Si mis à 0, les signaux n'ont aucune influence sur le seuil, de sorte que les signaux futurs sont détectés sur la base d'un seuil calculé avec une moyenne et un écart type non influencé par les signaux passés. Une autre façon de penser à cela est que si vous mettez l'influence à 0, vous supposez implicitement la stationnarité (c'est-à-dire que peu importe le nombre de signaux, la série chronologique retourne toujours à la même moyenne à long terme). Si ce n'est pas le cas, vous devez placer le paramètre d'influence entre 0 et 1, selon la mesure dans laquelle les signaux peuvent influencer systématiquement la tendance variant dans le temps des données. Par exemple, si les signaux entraînent une rupture structurelle de la moyenne à long terme des séries chronologiques, le paramètre d’influence doit être élevé (proche de 1) pour que le seuil puisse s’ajuster rapidement à ces changements.

threshold : le paramètre du seuil est le nombre d'écarts types par rapport à la moyenne mobile au-delà duquel l'algorithme classera un nouveau sharepoint données comme étant un signal. Par exemple, si un nouveau sharepoint données a un écart type de 4,0 au-dessus de la moyenne mobile et que le paramètre seuil est défini sur 3,5, l'algorithme identifiera le sharepoint données en tant que signal. Ce paramètre doit être défini en fonction du nombre de signaux attendus. Par exemple, si vos données sont normalement dissortingbuées, un seuil (ou: z-score) de 3,5 correspond à une probabilité de signalisation de 0,00047 (à partir de ce tableau ), ce qui implique un signal une fois tous les 2128 points de données (1 / 0.00047) . Le seuil influence donc directement la sensibilité de l'algorithme et donc aussi la fréquence à laquelle il signale. Examinez vos propres données et déterminez un seuil sensible qui déclenche le signal de l'algorithme lorsque vous le souhaitez (des essais et erreurs peuvent être nécessaires ici pour atteindre un seuil adéquat).


AVERTISSEMENT: Le code ci-dessus passe toujours en boucle sur tous les points de données chaque fois qu'il est exécuté. Lors de la mise en œuvre de ce code, veillez à diviser le calcul du signal en une fonction distincte (sans la boucle). Ensuite, lorsqu'un nouveau sharepoint données arrive, mettez à jour avgFilter , avgFilter et stdFilter une fois. Ne recalculez pas les signaux pour toutes les données à chaque fois qu'il y a un nouveau sharepoint données (comme dans l'exemple ci-dessus), ce qui serait extrêmement inefficace et lent!

D'autres manières de modifier l'algorithme (pour des améliorations potentielles) sont les suivantes:

  1. Utiliser la médiane au lieu de la moyenne
  2. Utiliser une mesure d’échelle robuste , telle que le MAD, au lieu de l’écart type
  3. Utilisez une marge de signalisation pour que le signal ne change pas trop souvent
  4. Changer le fonctionnement du paramètre d'influence
  5. Traiter différemment les signaux de haut en bas (traitement asymésortingque)

Citations académiques (connues) de cette réponse:

  • MC Catalbas, T. Cegovnik, J. Sodnik et A. Gulten (2017). Détection de la fatigue du conducteur basée sur des mouvements oculaires saccadés , 10ème Conférence internationale sur le génie élecsortingque et électronique (ELECO), pp. 913-917.

  • P. Willems (2017). Ambiances affectives contrôlées pour les personnes âgées , Mémoire de maîsortingse, Université de Twente.

  • O. Lo, WJ Buchanan, P. Griffiths et R. Macfarlane (2018), Méthodes de mesure de distance pour l'amélioration des réseaux de détection , de sécurité et de communication des menaces internes , vol. 2018, article ID 5906368.


Si vous utilisez cette fonction quelque part, merci de me créditer ou de répondre à cette question. Si vous avez des questions concernant cet algorithme, publiez-les dans les commentaires ci-dessous ou contactez-moi sur LinkedIn .


Voici l’implémentation Python / numpy de l’algorithme z-score lissé (voir la réponse ci-dessus ). Vous pouvez trouver l’ essentiel ici .

 #!/usr/bin/env python # Implementation of algorithm from https://stackoverflow.com/a/22640362/6029703 import numpy as np import pylab def thresholding_algo(y, lag, threshold, influence): signals = np.zeros(len(y)) filteredY = np.array(y) avgFilter = [0]*len(y) stdFilter = [0]*len(y) avgFilter[lag - 1] = np.mean(y[0:lag]) stdFilter[lag - 1] = np.std(y[0:lag]) for i in range(lag, len(y)): if abs(y[i] - avgFilter[i-1]) > threshold * stdFilter [i-1]: if y[i] > avgFilter[i-1]: signals[i] = 1 else: signals[i] = -1 filteredY[i] = influence * y[i] + (1 - influence) * filteredY[i-1] avgFilter[i] = np.mean(filteredY[(i-lag):i]) stdFilter[i] = np.std(filteredY[(i-lag):i]) else: signals[i] = 0 filteredY[i] = y[i] avgFilter[i] = np.mean(filteredY[(i-lag):i]) stdFilter[i] = np.std(filteredY[(i-lag):i]) return dict(signals = np.asarray(signals), avgFilter = np.asarray(avgFilter), stdFilter = np.asarray(stdFilter)) 

Vous trouverez ci-dessous le test sur le même jeu de données qui donne le même tracé que dans la réponse originale pour R / Matlab

 # Data y = np.array([1,1,1.1,1,0.9,1,1,1.1,1,0.9,1,1.1,1,1,0.9,1,1,1.1,1,1,1,1,1.1,0.9,1,1.1,1,1,0.9, 1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3, 2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1]) # Settings: lag = 30, threshold = 5, influence = 0 lag = 30 threshold = 5 influence = 0 # Run algo with settings from above result = thresholding_algo(y, lag=lag, threshold=threshold, influence=influence) # Plot result pylab.subplot(211) pylab.plot(np.arange(1, len(y)+1), y) pylab.plot(np.arange(1, len(y)+1), result["avgFilter"], color="cyan", lw=2) pylab.plot(np.arange(1, len(y)+1), result["avgFilter"] + threshold * result["stdFilter"], color="green", lw=2) pylab.plot(np.arange(1, len(y)+1), result["avgFilter"] - threshold * result["stdFilter"], color="green", lw=2) pylab.subplot(212) pylab.step(np.arange(1, len(y)+1), result["signals"], color="red", lw=2) pylab.ylim(-1.5, 1.5) 

Une approche consiste à détecter les pics basés sur l’observation suivante:

  • Le temps t est un pic si (y (t)> y (t-1)) && (y (t)> y (t + 1))

Il évite les faux positifs en attendant la fin de la tendance haussière. Ce n’est pas exactement “en temps réel” dans le sens où il va manquer le pic d’un dt. la sensibilité peut être contrôlée en exigeant une marge de comparaison. Il y a un compromis entre la détection bruyante et le délai de détection. Vous pouvez enrichir le modèle en ajoutant plus de parameters:

  • pic si (y (t) – y (t-dt)> m) && (y (t) – y (t + dt)> m)

dt et m sont des parameters pour contrôler la sensibilité en fonction du temps

Voici ce que vous obtenez avec l’algorithme mentionné: entrer la description de l'image ici

voici le code pour reproduire le tracé en python:

 import numpy as np import matplotlib.pyplot as plt input = np.array([ 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1.1, 1. , 0.8, 0.9, 1. , 1.2, 0.9, 1. , 1. , 1.1, 1.2, 1. , 1.5, 1. , 3. , 2. , 5. , 3. , 2. , 1. , 1. , 1. , 0.9, 1. , 1. , 3. , 2.6, 4. , 3. , 3.2, 2. , 1. , 1. , 1. , 1. , 1. ]) signal = (input > np.roll(input,1)) & (input > np.roll(input,-1)) plt.plot(input) plt.plot(signal.nonzero()[0], input[signal], 'ro') plt.show() 

En définissant m = 0.5 , vous pouvez obtenir un signal plus propre avec un seul faux positif: entrer la description de l'image ici

Dans le traitement du signal, la détection des pics se fait souvent par transformée en ondelettes. Vous faites essentiellement une transformée en ondelettes discrète sur vos données de séries chronologiques. Les passages à zéro dans les coefficients de détail renvoyés correspondront aux pics du signal de série temporelle. Vous obtenez différentes amplitudes de crête détectées à différents niveaux de coefficient de détail, ce qui vous donne une résolution à plusieurs niveaux.

GH Palshikar a trouvé un autre algorithme dans des algorithmes simples pour la détection de pics dans les séries temporelles .

L’algorithme se présente comme suit:

 algorithm peak1 // one peak detection algorithms that uses peak function S1 input T = x1, x2, …, xN, N // input time-series of N points input k // window size around the peak input h // typically 1 <= h <= 3 output O // set of peaks detected in T begin O = empty set // initially empty for (i = 1; i < n; i++) do // compute peak function value for each of the N points in T a[i] = S1(k,i,xi,T); end for Compute the mean m' and standard deviation s' of all positive values in array a; for (i = 1; i < n; i++) do // remove local peaks which are “small” in global context if (a[i] > 0 && (a[i] – m') >( h * s')) then O = O + {xi}; end if end for Order peaks in O in terms of increasing index in T // retain only one peak out of any set of peaks within distance k of each other for every adjacent pair of peaks xi and xj in O do if |j – i| <= k then remove the smaller value of {xi, xj} from O end if end for end 

Avantages

  • Le papier fournit 5 algorithmes différents pour la détection de pointe
  • Les algorithmes fonctionnent sur les données brutes de séries chronologiques (aucun lissage nécessaire)

Désavantages

  • Difficile de déterminer préalablement k et h
  • Les pics ne peuvent pas être plats (comme le troisième pic de mes données de test)

Exemple:

entrer la description de l'image ici

Voici une implémentation de l’algorithme Smoothed z-score (ci-dessus) dans Golang. Il suppose une tranche de []int16 (échantillons PCM 16 bits). Vous pouvez trouver une idée ici .

 /* Settings (the ones below are examples: choose what is best for your data) set lag to 5; # lag 5 for the smoothing functions set threshold to 3.5; # 3.5 standard deviations for signal set influence to 0.5; # between 0 and 1, where 1 is normal influence, 0.5 is half */ // ZScore on 16bit WAV samples func ZScore(samples []int16, lag int, threshold float64, influence float64) (signals []int16) { //lag := 20 //threshold := 3.5 //influence := 0.5 signals = make([]int16, len(samples)) filteredY := make([]int16, len(samples)) for i, sample := range samples[0:lag] { filteredY[i] = sample } avgFilter := make([]int16, len(samples)) stdFilter := make([]int16, len(samples)) avgFilter[lag] = Average(samples[0:lag]) stdFilter[lag] = Std(samples[0:lag]) for i := lag + 1; i < len(samples); i++ { f := float64(samples[i]) if float64(Abs(samples[i]-avgFilter[i-1])) > threshold*float64(stdFilter[i-1]) { if samples[i] > avgFilter[i-1] { signals[i] = 1 } else { signals[i] = -1 } filteredY[i] = int16(influence*f + (1-influence)*float64(filteredY[i-1])) avgFilter[i] = Average(filteredY[(i - lag):i]) stdFilter[i] = Std(filteredY[(i - lag):i]) } else { signals[i] = 0 filteredY[i] = samples[i] avgFilter[i] = Average(filteredY[(i - lag):i]) stdFilter[i] = Std(filteredY[(i - lag):i]) } } return } // Average a chunk of values func Average(chunk []int16) (avg int16) { var sum int64 for _, sample := range chunk { if sample < 0 { sample *= -1 } sum += int64(sample) } return int16(sum / int64(len(chunk))) } 

Ce problème ressemble à celui que j’ai rencontré dans un cours sur les systèmes hybrides / embarqués, mais cela était lié à la détection des défauts lorsque l’entrée d’un capteur était bruyante. Nous avons utilisé un filtre de Kalman pour estimer / prédire l’état caché du système, puis nous avons utilisé une parsing statistique pour déterminer la probabilité qu’un défaut se soit produit . Nous travaillions avec des systèmes linéaires, mais des variantes non linéaires existent. Je me souviens que l’approche était étonnamment adaptative, mais elle exigeait un modèle de la dynamic du système.

Dans la topologie computationnelle, l’idée d’une homologie persistante conduit à une solution efficace, rapide comme un sorting. Il ne détecte pas seulement les pics, il quantifie la “signification” des pics de manière naturelle, ce qui vous permet de sélectionner les pics les plus significatifs pour vous.

Résumé de l’algorithme. Dans un paramètre à 1 dimension (série temporelle, signal à valeur réelle), l’algorithme peut être facilement décrit par la figure suivante:

Les pics les plus persistants

Considérez le graphe de fonction (ou son ensemble de sous-niveaux) comme un paysage et considérez un niveau d’eau décroissant commençant au niveau infini (ou 1,8 dans cette image). Alors que le niveau diminue, des maxima locaux apparaissent sur les îles. Aux minima locaux, ces îles fusionnent. Un détail de cette idée est que l’île qui est apparue plus tard dans le temps est intégrée à l’île la plus ancienne. La “persistance” d’une île est son heure de naissance moins son temps de mort. Les longueurs des barres bleues représentent la persistance, qui est la “signification” mentionnée ci-dessus d’un pic.

Efficacité. Il n’est pas trop difficile de trouver une implémentation qui fonctionne en temps linéaire – en fait, c’est une simple boucle simple – après le sorting des valeurs de la fonction. Cette mise en œuvre doit donc être rapide dans la pratique et facilement mise en œuvre également.

Les références. Une description de l’histoire entière et des références à la motivation de l’homologie persistante (un domaine de la topologie algébrique informatique) peut être trouvée ici: https://www.sthu.org/blog/13-perstopology-peakdetection/index.html

Voici une implémentation C ++ de l’algorithme z-score lissé de cette réponse

 std::vector smoothedZScore(std::vector input) { //lag 5 for the smoothing functions int lag = 5; //3.5 standard deviations for signal float threshold = 3.5; //between 0 and 1, where 1 is normal influence, 0.5 is half float influence = .5; if (input.size() <= lag + 2) { std::vector emptyVec; return emptyVec; } //Initialise variables std::vector signals(input.size(), 0.0); std::vector filteredY(input.size(), 0.0); std::vector avgFilter(input.size(), 0.0); std::vector stdFilter(input.size(), 0.0); std::vector subVecStart(input.begin(), input.begin() + lag); avgFilter[lag] = mean(subVecStart); stdFilter[lag] = stdDev(subVecStart); for (size_t i = lag + 1; i < input.size(); i++) { if (std::abs(input[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1]) { if (input[i] > avgFilter[i - 1]) { signals[i] = 1; //# Positive signal } else { signals[i] = -1; //# Negative signal } //Make influence lower filteredY[i] = influence* input[i] + (1 - influence) * filteredY[i - 1]; } else { signals[i] = 0; //# No signal filteredY[i] = input[i]; } //Adjust the filters std::vector subVec(filteredY.begin() + i - lag, filteredY.begin() + i); avgFilter[i] = mean(subVec); stdFilter[i] = stdDev(subVec); } return signals; } 

C++ Implementation

 #include  #include  #include  #include  #include  #include  #include  using namespace std; typedef long double ld; typedef unsigned int uint; typedef std::vector::iterator vec_iter_ld; /** * Overriding the ostream operator for pretty printing vectors. */ template std::ostream &operator<<(std::ostream &os, std::vector vec) { os << "["; if (vec.size() != 0) { std::copy(vec.begin(), vec.end() - 1, std::ostream_iterator(os, " ")); os << vec.back(); } os << "]"; return os; } /** * This class calculates mean and standard deviation of a subvector. * This is basically stats computation of a subvector of a window size qual to "lag". */ class VectorStats { public: /** * Constructor for VectorStats class. * * @param start - This is the iterator position of the start of the window, * @param end - This is the iterator position of the end of the window, */ VectorStats(vec_iter_ld start, vec_iter_ld end) { this->start = start; this->end = end; this->compute(); } /** * This method calculates the mean and standard deviation using STL function. * This is the Two-Pass implementation of the Mean & Variance calculation. */ void compute() { ld sum = std::accumulate(start, end, 0.0); uint slice_size = std::distance(start, end); ld mean = sum / slice_size; std::vector diff(slice_size); std::transform(start, end, diff.begin(), [mean](ld x) { return x - mean; }); ld sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0); ld std_dev = std::sqrt(sq_sum / slice_size); this->m1 = mean; this->m2 = std_dev; } ld mean() { return m1; } ld standard_deviation() { return m2; } private: vec_iter_ld start; vec_iter_ld end; ld m1; ld m2; }; /** * This is the implementation of the Smoothed Z-Score Algorithm. * This is direction translation of https://stackoverflow.com/a/22640362/1461896. * * @param input - input signal * @param lag - the lag of the moving window * @param threshold - the z-score at which the algorithm signals * @param influence - the influence (between 0 and 1) of new signals on the mean and standard deviation * @return a hashmap containing the filtered signal and corresponding mean and standard deviation. */ unordered_map> z_score_thresholding(vector input, int lag, ld threshold, ld influence) { unordered_map> output; uint n = (uint) input.size(); vector signals(input.size()); vector filtered_input(input.begin(), input.end()); vector filtered_mean(input.size()); vector filtered_stddev(input.size()); VectorStats lag_subvector_stats(input.begin(), input.begin() + lag); filtered_mean[lag - 1] = lag_subvector_stats.mean(); filtered_stddev[lag - 1] = lag_subvector_stats.standard_deviation(); for (int i = lag; i < n; i++) { if (abs(input[i] - filtered_mean[i - 1]) > threshold * filtered_stddev[i - 1]) { signals[i] = (input[i] > filtered_mean[i - 1]) ? 1.0 : -1.0; filtered_input[i] = influence * input[i] + (1 - influence) * filtered_input[i - 1]; } else { signals[i] = 0.0; filtered_input[i] = input[i]; } VectorStats lag_subvector_stats(filtered_input.begin() + (i - lag), filtered_input.begin() + i); filtered_mean[i] = lag_subvector_stats.mean(); filtered_stddev[i] = lag_subvector_stats.standard_deviation(); } output["signals"] = signals; output["filtered_mean"] = filtered_mean; output["filtered_stddev"] = filtered_stddev; return output; }; int main() { vector input = {1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0, 1.0, 1.0, 1.0, 1.1, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 1.1, 1.0, 0.8, 0.9, 1.0, 1.2, 0.9, 1.0, 1.0, 1.1, 1.2, 1.0, 1.5, 1.0, 3.0, 2.0, 5.0, 3.0, 2.0, 1.0, 1.0, 1.0, 0.9, 1.0, 1.0, 3.0, 2.6, 4.0, 3.0, 3.2, 2.0, 1.0, 1.0, 0.8, 4.0, 4.0, 2.0, 2.5, 1.0, 1.0, 1.0}; int lag = 30; ld threshold = 5.0; ld influence = 0.0; unordered_map> output = z_score_thresholding(input, lag, threshold, influence); cout << output["signals"] << endl; } 

Here is a Groovy (Java) implementation of the smoothed z-score algorithm ( see answer above ).

 /** * "Smoothed zero-score alogrithm" shamelessly copied from https://stackoverflow.com/a/22640362/6029703 * Uses a rolling mean and a rolling deviation (separate) to identify peaks in a vector * * @param y - The input vector to analyze * @param lag - The lag of the moving window (ie how big the window is) * @param threshold - The z-score at which the algorithm signals (ie how many standard deviations away from the moving mean a peak (or signal) is) * @param influence - The influence (between 0 and 1) of new signals on the mean and standard deviation (how much a peak (or signal) should affect other values near it) * @return - The calculated averages (avgFilter) and deviations (stdFilter), and the signals (signals) */ public HashMap> thresholdingAlgo(List y, Long lag, Double threshold, Double influence) { //init stats instance SummaryStatistics stats = new SummaryStatistics() //the results (peaks, 1 or -1) of our algorithm List signals = new ArrayList(Collections.nCopies(y.size(), 0)) //filter out the signals (peaks) from our original list (using influence arg) List filteredY = new ArrayList(y) //the current average of the rolling window List avgFilter = new ArrayList(Collections.nCopies(y.size(), 0.0d)) //the current standard deviation of the rolling window List stdFilter = new ArrayList(Collections.nCopies(y.size(), 0.0d)) //init avgFilter and stdFilter (0..lag-1).each { stats.addValue(y[it as int]) } avgFilter[lag - 1 as int] = stats.getMean() stdFilter[lag - 1 as int] = Math.sqrt(stats.getPopulationVariance()) //getStandardDeviation() uses sample variance (not what we want) stats.clear() //loop input starting at end of rolling window (lag..y.size()-1).each { i -> //if the distance between the current value and average is enough standard deviations (threshold) away if (Math.abs((y[i as int] - avgFilter[i - 1 as int]) as Double) > threshold * stdFilter[i - 1 as int]) { //this is a signal (ie peak), determine if it is a positive or negative signal signals[i as int] = (y[i as int] > avgFilter[i - 1 as int]) ? 1 : -1 //filter this signal out using influence filteredY[i as int] = (influence * y[i as int]) + ((1-influence) * filteredY[i - 1 as int]) } else { //ensure this signal remains a zero signals[i as int] = 0 //ensure this value is not filtered filteredY[i as int] = y[i as int] } //update rolling average and deviation (i - lag..i-1).each { stats.addValue(filteredY[it as int] as Double) } avgFilter[i as int] = stats.getMean() stdFilter[i as int] = Math.sqrt(stats.getPopulationVariance()) //getStandardDeviation() uses sample variance (not what we want) stats.clear() } return [ signals : signals, avgFilter: avgFilter, stdFilter: stdFilter ] } 

Below is a test on the same dataset that yields the same results as the above Python / numpy implementation .

  // Data def y = [1d, 1d, 1.1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 0.9d, 1d, 1.1d, 1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 1d, 1d, 1d, 1.1d, 0.9d, 1d, 1.1d, 1d, 1d, 0.9d, 1d, 1.1d, 1d, 1d, 1.1d, 1d, 0.8d, 0.9d, 1d, 1.2d, 0.9d, 1d, 1d, 1.1d, 1.2d, 1d, 1.5d, 1d, 3d, 2d, 5d, 3d, 2d, 1d, 1d, 1d, 0.9d, 1d, 1d, 3d, 2.6d, 4d, 3d, 3.2d, 2d, 1d, 1d, 0.8d, 4d, 4d, 2d, 2.5d, 1d, 1d, 1d] // Settings def lag = 30 def threshold = 5 def influence = 0 def thresholdingResults = thresholdingAlgo((List) y, (Long) lag, (Double) threshold, (Double) influence) println y.size() println thresholdingResults.signals.size() println thresholdingResults.signals thresholdingResults.signals.eachWithIndex { x, idx -> if (x) { println y[idx] } } 

Here is a (non-idiomatic) Scala version of the smoothed z-score algorithm :

 /** * Smoothed zero-score alogrithm shamelessly copied from https://stackoverflow.com/a/22640362/6029703 * Uses a rolling mean and a rolling deviation (separate) to identify peaks in a vector * * @param y - The input vector to analyze * @param lag - The lag of the moving window (ie how big the window is) * @param threshold - The z-score at which the algorithm signals (ie how many standard deviations away from the moving mean a peak (or signal) is) * @param influence - The influence (between 0 and 1) of new signals on the mean and standard deviation (how much a peak (or signal) should affect other values near it) * @return - The calculated averages (avgFilter) and deviations (stdFilter), and the signals (signals) */ private def smoothedZScore(y: Seq[Double], lag: Int, threshold: Double, influence: Double): Seq[Int] = { val stats = new SummaryStatistics() // the results (peaks, 1 or -1) of our algorithm val signals = mutable.ArrayBuffer.fill(y.length)(0) // filter out the signals (peaks) from our original list (using influence arg) val filteredY = y.to[mutable.ArrayBuffer] // the current average of the rolling window val avgFilter = mutable.ArrayBuffer.fill(y.length)(0d) // the current standard deviation of the rolling window val stdFilter = mutable.ArrayBuffer.fill(y.length)(0d) // init avgFilter and stdFilter y.take(lag).foreach(s => stats.addValue(s)) avgFilter(lag - 1) = stats.getMean stdFilter(lag - 1) = Math.sqrt(stats.getPopulationVariance) // getStandardDeviation() uses sample variance (not what we want) // loop input starting at end of rolling window y.zipWithIndex.slice(lag, y.length - 1).foreach { case (s: Double, i: Int) => // if the distance between the current value and average is enough standard deviations (threshold) away if (Math.abs(s - avgFilter(i - 1)) > threshold * stdFilter(i - 1)) { // this is a signal (ie peak), determine if it is a positive or negative signal signals(i) = if (s > avgFilter(i - 1)) 1 else -1 // filter this signal out using influence filteredY(i) = (influence * s) + ((1 - influence) * filteredY(i - 1)) } else { // ensure this signal remains a zero signals(i) = 0 // ensure this value is not filtered filteredY(i) = s } // update rolling average and deviation stats.clear() filteredY.slice(i - lag, i).foreach(s => stats.addValue(s)) avgFilter(i) = stats.getMean stdFilter(i) = Math.sqrt(stats.getPopulationVariance) // getStandardDeviation() uses sample variance (not what we want) } println(y.length) println(signals.length) println(signals) signals.zipWithIndex.foreach { case(x: Int, idx: Int) => if (x == 1) { println(idx + " " + y(idx)) } } val data = y.zipWithIndex.map { case (s: Double, i: Int) => Map("x" -> i, "y" -> s, "name" -> "y", "row" -> "data") } ++ avgFilter.zipWithIndex.map { case (s: Double, i: Int) => Map("x" -> i, "y" -> s, "name" -> "avgFilter", "row" -> "data") } ++ avgFilter.zipWithIndex.map { case (s: Double, i: Int) => Map("x" -> i, "y" -> (s - threshold * stdFilter(i)), "name" -> "lower", "row" -> "data") } ++ avgFilter.zipWithIndex.map { case (s: Double, i: Int) => Map("x" -> i, "y" -> (s + threshold * stdFilter(i)), "name" -> "upper", "row" -> "data") } ++ signals.zipWithIndex.map { case (s: Int, i: Int) => Map("x" -> i, "y" -> s, "name" -> "signal", "row" -> "signal") } Vegas("Smoothed Z") .withData(data) .mark(Line) .encodeX("x", Quant) .encodeY("y", Quant) .encodeColor( field="name", dataType=Nominal ) .encodeRow("row", Ordinal) .show return signals } 

Here’s a test that returns the same results as the Python and Groovy versions:

 val y = List(1d, 1d, 1.1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 0.9d, 1d, 1.1d, 1d, 1d, 0.9d, 1d, 1d, 1.1d, 1d, 1d, 1d, 1d, 1.1d, 0.9d, 1d, 1.1d, 1d, 1d, 0.9d, 1d, 1.1d, 1d, 1d, 1.1d, 1d, 0.8d, 0.9d, 1d, 1.2d, 0.9d, 1d, 1d, 1.1d, 1.2d, 1d, 1.5d, 1d, 3d, 2d, 5d, 3d, 2d, 1d, 1d, 1d, 0.9d, 1d, 1d, 3d, 2.6d, 4d, 3d, 3.2d, 2d, 1d, 1d, 0.8d, 4d, 4d, 2d, 2.5d, 1d, 1d, 1d) val lag = 30 val threshold = 5d val influence = 0d smoothedZScore(y, lag, threshold, influence) 

vegas chart of result

Gist here

I needed something like this in my android project. Thought I might give back Kotlin implementation.

 /** * Smoothed zero-score alogrithm shamelessly copied from https://stackoverflow.com/a/22640362/6029703 * Uses a rolling mean and a rolling deviation (separate) to identify peaks in a vector * * @param y - The input vector to analyze * @param lag - The lag of the moving window (ie how big the window is) * @param threshold - The z-score at which the algorithm signals (ie how many standard deviations away from the moving mean a peak (or signal) is) * @param influence - The influence (between 0 and 1) of new signals on the mean and standard deviation (how much a peak (or signal) should affect other values near it) * @return - The calculated averages (avgFilter) and deviations (stdFilter), and the signals (signals) */ fun smoothedZScore(y: List, lag: Int, threshold: Double, influence: Double): Triple, List, List> { val stats = SummaryStatistics() // the results (peaks, 1 or -1) of our algorithm val signals = MutableList(y.size, { 0 }) // filter out the signals (peaks) from our original list (using influence arg) val filteredY = ArrayList(y) // the current average of the rolling window val avgFilter = MutableList(y.size, { 0.0 }) // the current standard deviation of the rolling window val stdFilter = MutableList(y.size, { 0.0 }) // init avgFilter and stdFilter y.take(lag).forEach { s -> stats.addValue(s) } avgFilter[lag - 1] = stats.mean stdFilter[lag - 1] = Math.sqrt(stats.populationVariance) // getStandardDeviation() uses sample variance (not what we want) stats.clear() //loop input starting at end of rolling window (lag..y.size - 1).forEach { i -> //if the distance between the current value and average is enough standard deviations (threshold) away if (Math.abs(y[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1]) { //this is a signal (ie peak), determine if it is a positive or negative signal signals[i] = if (y[i] > avgFilter[i - 1]) 1 else -1 //filter this signal out using influence filteredY[i] = (influence * y[i]) + ((1 - influence) * filteredY[i - 1]) } else { //ensure this signal remains a zero signals[i] = 0 //ensure this value is not filtered filteredY[i] = y[i] } //update rolling average and deviation (i - lag..i - 1).forEach { stats.addValue(filteredY[it]) } avgFilter[i] = stats.getMean() stdFilter[i] = Math.sqrt(stats.getPopulationVariance()) //getStandardDeviation() uses sample variance (not what we want) stats.clear() } return Triple(signals, avgFilter, stdFilter) } 

sample project with verification graphs can be found at github .

entrer la description de l'image ici

Here is my attempt at creating a Ruby solution for the “Smoothed z-score algo” from the accepted answer:

 module ThresholdingAlgoMixin def mean(array) array.reduce(&:+) / array.size.to_f end def stddev(array) array_mean = mean(array) Math.sqrt(array.reduce(0.0) { |a, b| a.to_f + ((b.to_f - array_mean) ** 2) } / array.size.to_f) end def thresholding_algo(lag: 5, threshold: 3.5, influence: 0.5) return nil if size < lag * 2 Array.new(size, 0).tap do |signals| filtered = Array.new(self) initial_slice = take(lag) avg_filter = Array.new(lag - 1, 0.0) + [mean(initial_slice)] std_filter = Array.new(lag - 1, 0.0) + [stddev(initial_slice)] (lag..size-1).each do |idx| prev = idx - 1 if (fetch(idx) - avg_filter[prev]).abs > threshold * std_filter[prev] signals[idx] = fetch(idx) > avg_filter[prev] ? 1 : -1 filtered[idx] = (influence * fetch(idx)) + ((1-influence) * filtered[prev]) end filtered_slice = filtered[idx-lag..prev] avg_filter[idx] = mean(filtered_slice) std_filter[idx] = stddev(filtered_slice) end end end end 

And example usage:

 test_data = [ 1, 1, 1.1, 1, 0.9, 1, 1, 1.1, 1, 0.9, 1, 1.1, 1, 1, 0.9, 1, 1, 1.1, 1, 1, 1, 1, 1.1, 0.9, 1, 1.1, 1, 1, 0.9, 1, 1.1, 1, 1, 1.1, 1, 0.8, 0.9, 1, 1.2, 0.9, 1, 1, 1.1, 1.2, 1, 1.5, 1, 3, 2, 5, 3, 2, 1, 1, 1, 0.9, 1, 1, 3, 2.6, 4, 3, 3.2, 2, 1, 1, 0.8, 4, 4, 2, 2.5, 1, 1, 1 ].extend(ThresholdingAlgoMixin) puts test_data.thresholding_algo.inspect # Output: [ # 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, # 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, # 1, 1, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0 # ] 

Here is an altered Fortran version of the z-score algorithm . It is altered specifically for peak (resonance) detection in transfer functions in frequency space (Each change has a small comment in code).

The first modification gives a warning to the user if there is a resonance near the lower bound of the input vector, indicated by a standard deviation higher than a certain threshold (10% in this case). This simply means the signal is not flat enough for the detection initializing the filters properly.

The second modification is that only the highest value of a peak is added to the found peaks. This is reached by comparing each found peak value to the magnitude of its (lag) predecessors and its (lag) successors.

The third change is to respect that resonance peaks usually show some form of symmetry around the resonance frequency. So it is natural to calculate the mean and std symmesortingcally around the current data point (rather than just for the predecessors). This results in a better peak detection behavior.

The modifications have the effect that the whole signal has to be known to the function beforehand which is the usual case for resonance detection (something like the Matlab Example of Jean-Paul where the data points are generated on the fly won’t work).

 function PeakDetect(y,lag,threshold, influence) implicit none ! Declaring part real, dimension(:), intent(in) :: y integer, dimension(size(y)) :: PeakDetect real, dimension(size(y)) :: filteredY, avgFilter, stdFilter integer :: lag, ii real :: threshold, influence ! Executing part PeakDetect = 0 filteredY = 0.0 filteredY(1:lag+1) = y(1:lag+1) avgFilter = 0.0 avgFilter(lag+1) = mean(y(1:2*lag+1)) stdFilter = 0.0 stdFilter(lag+1) = std(y(1:2*lag+1)) if (stdFilter(lag+1)/avgFilter(lag+1)>0.1) then ! If the coefficient of variation exceeds 10%, the signal is too uneven at the start, possibly because of a peak. write(unit=*,fmt=1001) 1001 format(1X,'Warning: Peak detection might have failed, as there may be a peak at the edge of the frequency range.',/) end if do ii = lag+2, size(y) if (abs(y(ii) - avgFilter(ii-1)) > threshold * stdFilter(ii-1)) then ! Find only the largest outstanding value which is only the one greater than its predecessor and its successor if (y(ii) > avgFilter(ii-1) .AND. y(ii) > y(ii-1) .AND. y(ii) > y(ii+1)) then PeakDetect(ii) = 1 end if filteredY(ii) = influence * y(ii) + (1 - influence) * filteredY(ii-1) else filteredY(ii) = y(ii) end if ! Modified with respect to the original code. Mean and standard deviation are calculted symmesortingcally around the current point avgFilter(ii) = mean(filteredY(ii-lag:ii+lag)) stdFilter(ii) = std(filteredY(ii-lag:ii+lag)) end do end function PeakDetect real function mean(y) !> @brief Calculates the mean of vector y implicit none ! Declaring part real, dimension(:), intent(in) :: y integer :: N ! Executing part N = max(1,size(y)) mean = sum(y)/N end function mean real function std(y) !> @brief Calculates the standard deviation of vector y implicit none ! Declaring part real, dimension(:), intent(in) :: y integer :: N ! Executing part N = max(1,size(y)) std = sqrt((N*dot_product(y,y) - sum(y)**2) / (N*(N-1))) end function std 

For my application the algorithm works like a charm! entrer la description de l'image ici

If the boundary value or other criteria depends on future values, then the only solution (without a time-machine, or other knowledge of future values) is to delay any decision until one has sufficient future values. If you want a level above a mean that spans, for example, 20 points, then you have to wait until you have at least 19 points ahead of any peak decision, or else the next new point could completely throw off your threshold 19 points ago.

Your current plot does not have any peaks… unless you somehow know in advance that the very next point isn’t 1e99, which after rescaling your plot’s Y dimension, would be flat up until that point.

Instead of comparing a maxima to the mean, one can also compare the maxima to adjacent minima where the minima are only defined above a noise threshold. If the local maximum is > 3 times (or other confidence factor) either adjacent minima, then that maxima is a peak. The peak determination is more accurate with wider moving windows. The above uses a calculation centered on the middle of the window, by the way, rather than a calculation at the end of the window (== lag).

Note that a maxima has to be seen as an increase in signal before and a decrease after.

It is unclear why in the answer

Peak signal detection in realtime timeseries data

we set

 set filteredY(i) to influence*y(i) + (1-influence)*filteredY(i-1); 

only in peak mode

 if absolute(y(i) - avgFilter(i-1)) > threshold*stdFilter(i-1) 

but else we just

 set filteredY(i) to y(i); 

Why we do not smooth signal in no-peak mode?