algorithm - synchrone - principe de demodulation




Détection de pic du signal mesuré (6)

Ce problème a été étudié en détail.

Il existe un ensemble d'implémentations très récentes dans les classes TSpectrum* de ROOT (outil d'analyse de la physique nucléaire / des particules). Le code fonctionne dans des données à une ou trois dimensions.

Le code source ROOT est disponible, vous pouvez donc saisir cette implémentation si vous le souhaitez.

Dans la documentation de classe TSpectrum :

Les algorithmes utilisés dans cette classe ont été publiés dans les références suivantes:

[1] M. Morhac et al.: Méthodes d'élimination de fond pour les spectres gamma de coïncidences multidimensionnelles. Instruments et méthodes nucléaires dans la recherche en physique A 401 (1997) 113-122.

[2] M. Morhac et al.: Déconvolution efficace de l'or à une et deux dimensions et son application à la décomposition des spectres gamma. Instruments et méthodes nucléaires dans la recherche en physique A 401 (1997) 385-408.

[3] M. Morhac et al.: Identification des pics dans les spectres gamma de coïncidences multidimensionnelles. Instruments et méthodes nucléaires en physique de la recherche A 443 (2000), 108-125.

Les articles sont liés à partir de la documentation de classe pour ceux d'entre vous qui ne possèdent pas d'abonnement en ligne à NIM.

La version courte de ce qui est fait est que l'histogramme est aplati pour éliminer le bruit, puis les maxima locaux sont détectés par la force brute dans l'histogramme aplati.

Nous utilisons une carte d’acquisition de données pour prendre les lectures d’un appareil qui augmente son signal jusqu’à un pic puis redescend à une valeur proche de la valeur initiale. Pour trouver la valeur de pic, nous recherchons actuellement dans le tableau la lecture la plus élevée et utilisons l'index pour déterminer la synchronisation de la valeur de pic utilisée dans nos calculs.

Cela fonctionne bien si la valeur la plus élevée est le pic que nous recherchons, mais si le périphérique ne fonctionne pas correctement, nous pouvons voir un deuxième pic qui peut être supérieur au pic initial. Nous effectuons 10 lectures par seconde sur 16 appareils sur une période de 90 secondes.

Mes pensées initiales sont de parcourir les lectures en vérifiant si les points précédents et suivants sont inférieurs au courant pour trouver un pic et construire un tableau de pics. Peut-être devrions-nous examiner un nombre moyen de points de part et d’autre de la position actuelle pour tenir compte du bruit dans le système. Est-ce la meilleure façon de procéder ou existe-t-il de meilleures techniques?

Nous utilisons effectivement LabVIEW et j'ai consulté les forums LAVA . Il existe un certain nombre d'exemples intéressants. Cela fait partie de notre logiciel de test et nous essayons d'éviter d'utiliser trop de bibliothèques de VI non standard. J'espérais donc recevoir des informations sur le processus / les algorithmes impliqués plutôt que sur du code spécifique.


Cette méthode est essentiellement tirée du livre "Vision" de David Marr

Le flou gaussien de votre signal correspond à la largeur attendue de vos pics. cela élimine les pointes de bruit et vos données de phase ne sont pas endommagées.

Puis détection de bord (LOG fera l'affaire)

Ensuite, vos bords étaient les bords des entités (comme les pics). cherchez les pics entre les arêtes, triez les pics par taille et vous avez terminé.

J'ai utilisé des variantes à ce sujet et elles fonctionnent très bien.


Il existe de nombreuses méthodes classiques de détection des pics, chacune d’elles pouvant fonctionner. Vous devrez voir ce qui, en particulier, limite la qualité de vos données. Voici les descriptions de base:

  1. Entre deux points quelconques dans vos données, (x(0), y(0)) et (x(n), y(n)) , additionnez y(i + 1) - y(i) pour 0 <= i < n et appelez ce T ("voyage") et réglez R ("montée") sur y(n) - y(0) + k pour un y(n) - y(0) + k convenablement petit. T/R > 1 indique un pic. Cela fonctionne bien si une course importante due au bruit est peu probable ou si le bruit est distribué de manière symétrique autour d'une forme de courbe de base. Pour votre application, acceptez le pic le plus ancien avec un score supérieur à un seuil donné ou analysez la courbe des valeurs de déplacement par élévation pour obtenir des propriétés plus intéressantes.

  2. Utilisez des filtres adaptés pour marquer la similarité avec une forme de pic standard (utilisez essentiellement un produit scalaire normalisé contre une forme quelconque pour obtenir une métrique de similitude en cosinus)

  3. Déconnectez-vous par rapport à une forme de pic standard et vérifiez les valeurs élevées (bien que je trouve souvent que 2 est moins sensible au bruit pour une sortie d'instrumentation simple).

  4. Lisser les données et rechercher des triplets de points équidistants où, si x0 < x1 < x2, y1 > 0.5 * (y0 + y2) , ou vérifier les distances euclidiennes comme ceci: D((x0, y0), (x1, y1)) + D((x1, y1), (x2, y2)) > D((x0, y0),(x2, y2)) , qui repose sur l’inégalité du triangle. L'utilisation de ratios simples vous fournira à nouveau un mécanisme de notation.

  5. Ajustez un modèle de mélange 2-gaussien très simple à vos données (par exemple, Numerical Recipes a un joli bloc de code prêt à l'emploi). Prenez le pic précédent. Cela traitera correctement les pics qui se chevauchent.

  6. Trouvez la meilleure correspondance dans les données avec une simple courbe gaussienne, de Cauchy, de Poisson ou autre. Evaluez cette courbe sur une large plage et soustrayez-la d'une copie des données après avoir noté son emplacement de pointe. Répéter. Prenez le premier pic dont les paramètres de modèle (écart-type probablement, mais certaines applications peuvent se soucier de kurtosis ou d'autres caractéristiques) répondent à certains critères. Méfiez-vous des artefacts laissés lorsque les pics sont soustraits des données. La meilleure correspondance peut être déterminée par le type de notation de match suggéré au point 2 ci-dessus.

J'ai déjà fait ce que vous faisiez auparavant: trouver des pics dans les données de séquence d'ADN, des pics dans les dérivés estimés à partir des courbes mesurées et des pics dans les histogrammes.

Je vous encourage à faire attention à bien établir les bases. Le filtrage de Wiener ou un autre filtrage ou une simple analyse d'histogramme est souvent un moyen facile de se baser en présence de bruit.

Enfin, si vos données sont généralement bruitées et que vous les récupérez en tant que sortie asymétrique non référencée (ou même référencée, tout simplement différentielle), et si vous faites la moyenne de nombreuses observations sur chaque point de données, essayez de les trier. observations et jeter le premier et le dernier quartile et faire la moyenne de ce qui reste. Il existe une foule de telles tactiques d'élimination aberrantes qui peuvent être vraiment utiles.


Je ne connais pas grand chose à l'instrumentation, donc cela pourrait être totalement irréalisable, mais là encore, cela pourrait être une direction différente utile. Si vous savez comment les lectures peuvent échouer et qu'il existe un certain intervalle entre les pics étant donné ces défaillances, pourquoi ne pas effectuer une descente de gradient à chaque intervalle. Si la descente vous ramène dans une zone déjà explorée, vous pouvez l'abandonner. En fonction de la forme de la surface échantillonnée, cela peut également vous aider à trouver des pics plus rapidement que la recherche.


Je voudrais contribuer à ce fil un algorithme que j'ai moi - même développé :

Il est basé sur le principe de dispersion : si un nouveau point de donnée est un nombre x d'écarts-types donné 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 la même précision, quelle que soit 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 point de donnée 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 permettant de recalculer le nouveau seuil: une influence de 0 est donc l'option la plus robuste.

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,...,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,...,lag));       # Initialise first value
set stdFilter(lag) to std(y(1,...,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)
      set signals(i) to +1;                     # Positive signal
    else
      set signals(i) to -1;                     # Negative signal
    end
    # Adjust the filters
    set filteredY(i) to influence*y(i) + (1-influence)*filteredY(i-1);
    set avgFilter(i) to mean(filteredY(i-lag,i),lag);
    set stdFilter(i) to std(filteredY(i-lag,i),lag);
  else
    set signals(i) to 0;                        # No signal
    # Adjust the filters
    set filteredY(i) to y(i);
    set avgFilter(i) to mean(filteredY(i-lag,i),lag);
    set stdFilter(i) to std(filteredY(i-lag,i),lag);
  end
end

Démo

> Réponse originale


Vous pouvez appliquer un écart type à votre logique et prendre en compte les pics supérieurs à x%.







language-agnostic