Astronomie

Comprendre le périodogramme lomb-scargle

Comprendre le périodogramme lomb-scargle

J'essaie d'utiliser le périodogramme pour dire quand un signal est périodique ou non en suivant le didacticiel pour le périodogramme d'astropie Lomb-scargle ici.

http://docs.astropy.org/en/stable/stats/lombscargle.html

J'ai simulé certaines données, une qui est une sinusoïde (période = 200) et une qui est une gaussienne asymétrique (c'est-à-dire un événement transitoire unique). L'espoir était que le périodogramme sélectionne la période de l'objet périodique et donne une période pour le transitoire qui impliquerait qu'un seul événement s'est produit dans la fenêtre.

Malheureusement, les résultats n'ont aucun sens. J'ai joint le code à la fin et les chiffres générés ci-dessous. En raison du bruit aléatoire dans ma simulation, chaque résultat est différent et je fournis deux exemples des résultats ci-dessous. J'utilise la même méthode décrite dans le lien ci-dessus où nous utilisons la fonctionLombScargle.model()sur la meilleure fréquence d'ajustement defréquence[argmax[puissance]]

La ligne rouge est la vraie fonction à partir de laquelle j'ai simulé les données. Le vert correspond le mieux au périodogramme. Les tracés de droite sont le PSD du périodogramme.

Exemple 1

Ici, nous pouvons voir que la meilleure fréquence d'ajustement pour la sinusoïde (graphique en haut à droite) choisie est de 0,105 (c'est-à-dire une période de 9 à 10 jours), ce qui n'est pas proche d'une fréquence de 0,05 à laquelle je m'attendrais pour quelque chose avec une période de 200 jours , pourtant, lorsque je transmets cette fréquence de meilleur ajustement de 0,105 à l'ajusteur du modèle lomb-scargle, une courbe périodique bien adaptée sort avec une période de 200 jours ?

Cela n'a pas de sens.

Exemple 2

Ici, j'ai réexécuté le code et cette fois, les résultats sont inversés ? Il correspondait à une très grande période du transitoire, de sorte que je peux dire avec confiance qu'il s'agit d'un événement transitoire unique, mais l'ajustement sinusoïdal est terrible. La meilleure fréquence est toujours de 0,105 (période = 10), mais l'installateur du modèle lomb-scargle superpose quelque chose qui semble avoir une période de 60 jours, ce qui est faux ?

Pourrais-je obtenir des éclaircissements si je fais quelque chose de mal? On m'a dit que le périodogramme est l'outil de facto pour les données échantillonnées de manière inégale comme celle-ci encore… les résultats semblent horribles la moitié du temps.

Pour clarifier mes questions

  1. Peut-on expliquer comment, dans le premier tracé, la fréquence la mieux ajustée de 0,105 que j'alimente dans l'ajusteur du modèle d'astropie lomb-scargle crée en quelque sorte une sinusoïde correctement adaptée avec une fréquence de 0,05 ? Quelle est l'explication ?

  2. Pourquoi y a-t-il 5 pics forts dans le graphique en haut à droite pour le périodogramme de l'exemple 1 alors que je n'en attends qu'un ? Les deux du milieu sont proches de la valeur réelle de 0,05 (à 0,045 et 0,055)

Voici le code court utilisé pour simuler et tracer les données

importer numpy en tant que np importer matplotlib.pyplot en tant que plt importer scipy.stats en tant que sc de astropy.stats importer LombScargle importer des mathématiques #simuler les paramètres data_range = [i for i in range (1 1001)] number_of_samples = 50 gauss_skew = sc.skewnorm. pdf skew = -10 période = 200 location = data_range[int(len(data_range)/2)] y1= [(2 * (1. + np.sin(2. * np. pi * x/period)) + np .random.normal(loc =0.0, scale = 0.5)) for x in data_range] error1 = [np.random.normal(loc = 0.0, scale = 1) for x in data_range] y2 = [(1000* (gauss_skew( x,skew,loc=location ,scale = 50)) + np.random.normal(loc =0.0, scale = 1)) for x in data_range] error2 = [np.random.normal(loc = 0.0, scale = 1 ) for x in data_range] sample_rate = int(len(data_range)/number_of_samples)# Pour éclaircir un peu les données y1 = y1[::sample_rate] y2 = y2[::sample_rate] error1 = error1[::sample_rate] error2 = erreurs2[::sample_rate] x1= x2 =data_range[::sample_rate] vérité1 = [2* (1. + np.sin(2. * np. pi * x/period)) pour x dans data_range] vérité2 = [1000 * (gau ss_skew(x,-10,loc=location ,scale = 50)) pour x dans data_range] vérités = [vérité1,vérité2] x= [x1,x2] y=[y1,y2] erreurs = [erreurs1,erreurs2] fig ,ax = plt.subplots(nrows=2,ncols=2) ax[0][0].errorbar(x1,y1,yerr=errors1,fmt='o') ax[0][0].set_xlabel(' time') ax[1][0].errorbar(x2,y2,yerr=errors2,fmt='o') ax[1][0].set_xlabel('time') ax[0][1].set_xlabel ('fréquence') ax[1][1].set_xlabel('fréquence') pour i dans la plage (0,2) : fréquence, puissance = LombScargle(x[i],y[i],errors[i]) .autopower() #Obtenir la meilleure fréquence d'ajustement comme dans le tutoriel best_fréquence = fréquence[np.argmax(power)] print('meilleure fréquence:',best_fréquence) t_fit = np.linspace(x[i][0], math .floor(x[i][-1]),num =50) #Ajuster la fréquence de meilleur ajustement #tracer le meilleur meilleur modèle basé sur le meilleur ajustement y_fit = LombScargle(x[i], y[i], erreurs[ i]).model(t_fit, meilleure_fréquence) ax[i][0].plot(t_fit,y_fit,'g') ax[i][0].plot(data_range,truths[i],'r') # Tracer le PSD ax[i][1].plot(fréquence, puissance) ax[i][1].axvline(x=meilleure_fréquence,color='noir', ls="--") plt.show()

Je pense que je vois ce qui se passe. Votre meilleur ajustement est en fait une fréquence de battement entre la fréquence réelle (0.005) et le double de la fréquence d'échantillonnage (0.05). Cela produit en effet un modèle qui passe par vos points de données, mais comme vous n'avez tracé que 50 points dans le modèle, vous n'avez pas pu le voir.

Si vous modifiez cette ligne t_fit = np.linspace(x[i][0], math.floor(x[i][-1]),num =1000)

de sorte qu'il montre le modèle à plus de points, vous pouvez voir qu'il est en effet à une fréquence beaucoup plus élevée.

L'astuce est de connaître (approximativement) la réponse avant de commencer.

Si vous limitez votre plage de fréquences dans la commande Lomb-Scargle en dessous de la fréquence de Nyquist :

fréquence, puissance = LombScargle(x[i],y[i],errors[i]).autopower(nyquist_factor=1.5)

Vous trouverez le comportement que vous attendiez. par exemple. (NB j'ai réduit la taille de la barre d'erreur pour me concentrer sur votre problème, mais j'ai également fixé la taille de la barre d'erreur à un seul positif valeur. Votre code introduisait des erreurs très petites et même négatives et je ne sais pas comment le code les traiterait).

Dans des données inégalement espacées, les choses ne sont pas si simples car la fréquence de Nyquist n'est pas bien définie.


Remarque : ce tracé contient une erreur¶

Après la mise sous presse du livre, un lecteur a souligné que cette intrigue contenait une faute de frappe. Au lieu de transmettre les données bruitées à la routine Lomb-Scargle, nous avions transmis les données sous-jacentes non bruitées. Cela a provoqué une surestimation de la puissance de Lomb-Scargle.

Pour cette raison, nous ajoutons deux tracés supplémentaires à ce script : le premier reproduit le tracé actuel sans la faute de frappe. On y voit que pour les données bruitées, la période n'est détectée que pour

30 points en dix périodes.

Dans le deuxième graphique supplémentaire, nous augmentons la ligne de base et le nombre de points par un facteur de dix. Avec cette configuration, le pic est détecté et les aspects qualitatifs de la discussion ci-dessus restent vrais.


Programme AFNI : 3dLombScargle

Faire un périodogramme ou un spectre d'amplitude d'une série temporelle qui a un
taux d'échantillonnage non constant. Les spectres émis par ce programme sont
'unilatérale', de sorte qu'elles représentent la demi-amplitude ou la puissance
associés à une fréquence, et ils nécessiteraient un facteur de 2 à
compte à la fois des solutions de fréquence de déplacement à droite et à gauche
de la transformée de Fourier (voir ci-dessous 'SORTIE' et 'NOTE').

L'application de cette fonctionnalité à
séries temporelles d'état de repos qui peuvent avoir été censurées. La théorie derrière
les mathématiques et les algorithmes de ceci sont dus à des groupes séparés, principalement
dans le domaine des applications astrophysiques : Vaníček (1969, 1971),
Lomb (1976), Scargle (1982) et Press & Rybicki (1989). Bravo à eux.

Cette implémentation particulière est due à Press & Rybicki (1989), par
traduire essentiellement leur implémentation Fortran publiée en C,
en utilisant GSL pour la FFT, au lieu de realft() de NR, et en faisant
plusieurs ajustements basés sur cela.

L'adaptation Lomb-Scargle a été faite avec des changements assez minimes ici par
PA Taylor (v1.4, juin 2016).

+ UTILISATION :
Saisir une série temporelle volumétrique 4D (ensemble de données BRIK/HEAD ou NIFTI)
ainsi qu'un fichier 1D optionnel de 0 et de 1 qui définit quels points
censurer (c'est-à-dire que chaque 0 représente un point/volume à censurer)
si aucun fichier 1D n'est entré, le programme vérifiera les volumes qui sont
uniformément zéro et considérer ceux qui sont censurés.

La sortie est un périodogramme LS, décrivant les magnitudes spectrales
jusqu'à une certaine "fréquence maximale" - le maximum par défaut est ce que
la fréquence de Nyquist de la série temporelle *aurait été* sans
toute censure. (Il est intéressant de noter que cette analyse peut en fait être
légitimement appliqué dans les cas pour estimer le contenu fréquentiel >Nyquist.
Wow!)

Le spectre de fréquence sera dans la plage [df, f_N], où :
df = 1/T, et T est la durée totale de la série chronologique non censurée
f_N = 1/dt, et dt est le temps d'échantillonnage (c'est-à-dire TR)
et l'intervalle de fréquences est également df.
Ces plages et tailles de pas doivent être *indépendantes* de la censure
qui est une belle propriété du Lomb-Scargle-iness.

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ SORTIE :
1) PREFIX_time.1D : un fichier 1D des points temporels échantillonnés (en unités de
secondes) des données analysées (et éventuellement censurées)
base de données.
2) PREFIX_freq.1D : un fichier 1D des points d'échantillonnage fréquentiel (en unités
de 1/seconde) du périodogramme/spectre de sortie
base de données.
3) PREFIX_amp+orig :ensemble de données volumétriques contenant un dérivé LS
ou un spectre d'amplitude (par défaut, nommé 'amp') ou un
Spectre de puissance PREFIX_pow+orig (voir '-out_pow_spec', nommé 'pow')
un par voxel.
Veuillez noter que l'amplitude et la puissance de sortie
les spectres sont «unilatéraux», pour représenter le
*demi* amplitude ou puissance d'une fréquence donnée
(voir la remarque suivante).

+ UNE NOTE SUR Fourier+Parseval (veuillez pardonner le maladroit
mise en page):
Dans la formulation utilisée ici, pour une série temporelle x[n] de longueur N,
la valeur du périodogramme S[k] est liée à la valeur d'amplitude |X[k]| :
(1) S[k] = (|X[k]|)**2,
pour chaque k-ième harmonique.

Le théorème de Parseval relie les fluctuations temporelles aux amplitudes spectrales,
indiquant que (pour les séries en temps réel avec une moyenne nulle) :
(2) sum_n < x[n]**2 >= (1/N) * sum_k< |X[k]|**2 >,
= (1/N) * somme_k< S[k] >,
où n=0,1. N-1 et k=0,1. N-1 (NB : A[0]=0, pour une moyenne nulle
séries). Le LHS est essentiellement la variance de la série chronologique
(fois N-1). Ce qui précède est dérivé des mathématiques de la transformée de Fourier, et
les spectres de Lomb-Scargle sont des approximations de Fourier, donc ce qui précède
devrait tenir approximativement, si tout se passe bien.

Un autre résultat lié à Fourier est que pour les séries temporelles réelles et discrètes,
les valeurs spectrales d'amplitudes/puissances sont symétriques et périodiques en N.
Par conséquent, |X[k]| = |X[-k]| = |X[N-k-1]| (dans un tableau de base zéro
compte)
la distinction entre les fréquences indexées positives et négatives
peut être considérée comme signifiant des ondes se déplaçant à droite et à gauche, qui
les deux contribuent à la puissance totale d'une fréquence spécifique.
Le résultat est que l'on pourrait écrire la formule de Parseval comme :
(3) sum_n < x[n]**2 >= (2/N) * sum_l< |X[l]|**2 >,
= (2/N) * somme_l< S[l] >,
où n=0,1. N-1 et l=0,1. (N/2)-1 (notez maintenant le facteur 2
apparaissant sur les relations RHS). Ces symétries/considérations
sont la raison pour laquelle

Les valeurs de fréquence N/2 sont sorties ici (nous supposons
que seules les séries temporelles à valeur réelle sont entrées), sans aucune perte de
informations.

De plus, en vue d'exprimer l'amplitude globale
ou la puissance d'une fréquence donnée, que de nombreuses personnes pourraient vouloir utiliser pour
estimer les paramètres spectraux de « connectivité fonctionnelle » tels que ALFF,
fALFF, RSFA, etc. (en utilisant par exemple 3dAmptoRSFC), on
notez que l'amplitude ou la puissance *totale* d'une fréquence donnée serait
être:
A[k] = 2*|X[k]|
P[k] = 2*S[k] = 2*|X[k]|**2 = 0.5*A[k]**2
au lieu de simplement celui de la partie mobile gauche/droite. Ces types de
les grandeurs (A et P) sont également appelées spectres « bifaces ». le
La relation de Parseval résultante pourrait alors s'écrire :
(4) somme_n < x[n]**2 >= (1/(2N)) * somme_l< A[l]**2 >,
= (1/N) * somme_l< P[l] >,
où n=0,1. N-1 et l=0,1. (N/2)-1. D'une certaine manière, il semble juste plus facile
pour sortir les valeurs unilatérales, X et S, de sorte que le Parsevalian
les règles de sommation se ressemblent davantage.

Avec tout cela à l'esprit, les résultats 3dLombScargle sont sortis comme
suit. Pour les amplitudes, les env. Relation Parsevellian
devrait tenir entre la série temporelle « trouée » x[m] de M points et
la série de fréquences Y[l] de L

M/2 points (où <|Y[l]|> se rapproche
les amplitudes de Fourier <|X[l]|>en nombre de points censurés
diminue et M->N):
(5) sum_m < x[m]**2 >= (1/L) * sum_l< Y[l]**2 >,
où m=0,1. M-1 et l=0,1. L-1. Pour le spectre de puissance T[l]
de L

valeurs M/2, alors :
(6) sum_m < x[m]**2 >= (1/L) * sum_l < T[l] >
pour les mêmes plages de sommations.

Donc, veuillez considérer cela lorsque vous utilisez les sorties d'ici. 3dAmpVersRSFC
est préparé pour cela lors du calcul des paramètres spectraux (à partir de
amplitude).

+ COURSE :
-prefix PREFIX : nom du préfixe de sortie pour le volume de données, fichier 1D du point temporel
et le fichier de fréquence 1D.
-inset FILE :série temporelle de volumes, un ensemble de données volumétriques 4D.

-censor_1D C1D : une seule ligne ou colonne de 1 (conserver) et de 0 (censuré)
décrivant quels volumes de FILE sont conservés dans le
l'échantillonnage et qui sont censurés, respectivement. le
la longueur de la liste des nombres doit être du
même longueur que le nombre de volumes dans FILE.
S'il n'est pas entré, le programme recherchera les sous-briques
de tous les zéros et supposer que ceux-ci sont censurés.
-censor_str CSTR :chaîne de sélection de volumes de style AFNI à *conserver*
l'analyse. Tel que:
'[0..4,7,10..$]'
Pourquoi nous l'appelons une "chaîne de censure" alors qu'elle est
vraiment la liste des volumes à conserver. eh bien, ça a fait
sens à l'époque. Les futurs historiens peuvent se battre avec
encre à ce sujet.

-mask MASK :optionnel, masque de volume à analyser en plus, tout
voxel avec des valeurs uniformément nulles dans le temps
produire un spectre zéro.

-out_pow_spec : commutateur pour afficher le spectre d'amplitude des fréquences
au lieu du périodogramme. Dans la formulation utilisée
ici, pour une série temporelle de longueur N, la puissance spectrale
la valeur S est liée à la valeur d'amplitude X comme :
S = (X)**2. (Sans cette option, la sortie par défaut est
spectre d'amplitude.)

-nyq_mult N2 :Les périodogrammes L-S peuvent inclure des fréquences au-dessus de ce que
serait généralement considéré comme Nyquist (ici défini
comme:
f_N = 0.5*(nombre d'échantillons)/(intervalle de temps total)
Par défaut, la fréquence maximale sera ce que
f_N *aurait* été si aucune censure de points n'avait été
s'est produit. (Cela facilite la comparaison des spectres L-S
à travers un groupe avec le même protocole de scan, même si
il y a de légères différences dans la censure, par sujet.)
Les valeurs acceptables sont >0. (Pour ceux qui lisent le
papiers d'algorithme, cela définit le paramètre 'hifac'.)
Si vous n'avez pas de bonne raison pour changer cela,
ne le change pas !
-nifti : basculer vers la sortie du fichier de volume *.nii.gz
(le format par défaut est BRIK/HEAD).

+ EXEMPLE :
3dLombScargle -préfixe LSout -inset TimeSeries.nii.gz
-mask mask.nii.gz -censor_1D censor_list.txt


____________________________________________________________________________
Cette page a été générée automatiquement le mercredi 23 juin 15:19:40 EDT 2021


Résultats

Buffles d'Afrique

Le LSP des emplacements de Cilla était caractéristique d'une marche aléatoire apériodique (pas de pic), alors que le LSP des emplacements de Pepper suggérait un cycle quotidien (Fig. 3). En appliquant le modèle d'observations manquantes de Pepper à la série temporelle de localisation de Cilla, nous avons récupéré, comme prévu, un signal de périodicité quotidienne (Fig. 3). Cela a illustré que le modèle périodique a été créé par autocorrélation dans la manière dont les observations manquaient, et non par le comportement des animaux. le P-valeur de 0,29 à partir de 150 simulations dans le test du modèle nul appliqué aux emplacements de Pepper a confirmé que le modèle périodique observé dans l'enregistrement de suivi de Pepper était artefactuel.

La gauche: Périodogramme des données de suivi du buffle Cilla selon le programme d'échantillonnage original. Centre: Périodogramme des données de Cilla lors du rééchantillonnage pour imiter le programme d'échantillonnage de Pepper. Droite: Périodogramme des données de suivi de Pepper. Les flèches pointent vers les pics correspondant à la période d'un jour

En appliquant le LSP au record de vitesse de Cilla (vs. record de localisation précédemment), nous avons trouvé, comme prévu, une périodicité quotidienne significative (P-valeur <0.01 à partir de 150 simulations Fichier supplémentaire 2 : Figure C5). Cilla activité était donc périodique même si son utilisation de l'espace ne l'était pas (Fig. 3). Cela illustre comment les modèles périodiques d'utilisation de l'espace sont découplés de la périodicité des niveaux d'activité chez cette espèce.

Albatros agités

Les quatre individus d'albatros présentaient une périodicité significative (tous P-valeurs <0.005 à partir de 1000 simulations chacune) avec une période estimée allant de 22 à 35 jours. Chez l'individu qui a effectué les voyages les plus longs, atteignant les eaux à 1800 km de sa colonie de reproduction avec une périodicité estimée à 29 jours, l'association avec la phase lunaire était très prononcée (Fig. 4). Pendant trois cycles surveillés, l'oiseau a quitté la colonie peu de temps avant la nouvelle lune et a commencé son voyage de retour peu de temps après la pleine lune, vraisemblablement pour que lui et son compagnon puissent se nourrir pendant les périodes de clair de lune.

Suivi des données d'un albatros agité. La gauche: tracé avec une échelle de couleurs indiquant la phase de la lune. Les couleurs bleues sont proches de la nouvelle lune et les couleurs rouges de la pleine lune. Droite: périodogrammes. Le périodogramme noir provient des données de suivi et le périodogramme rouge provient du programme d'échantillonnage. La présence d'un pic dans le périodogramme noir pendant une période d'un cycle lunaire (ligne verticale), et son absence dans le périodogramme rouge, indiquent que le modèle périodique n'est pas causé par le programme d'échantillonnage

Loups à crinière

La périodicité d'un jour était statistiquement significative pour les huit individus (tous P-valeurs <0,05 à partir de 1000 simulations chacune Fig. 5, panneau de droite), même si le modèle périodique d'utilisation de l'espace n'était pas directement évident dans les données de suivi brutes (Fig. 5, panneau de gauche). En d'autres termes, ces huit loups ont montré une tendance significative à patrouiller leur domaine vital le long d'itinéraires répétés quotidiennement, mais ce schéma était obscurci par le bruit de fond important (les composantes stochastiques apériodiques des chemins).

Suivi des données des loups à crinière. La gauche: tracé pour un individu (Amadeo) selon l'heure de la journée, pour montrer la difficulté à détecter visuellement le modèle périodique en raison de la composante stochastique importante dans le processus de mouvement. Droite: périodogramme multi-individuel d'Amadeo et de 7 autres loups, montrant la périodicité d'un jour et les séries d'harmoniques associées (flèches rouges)


Observabilité des composantes spectrales au-delà de la limite de Nyquist dans les signaux échantillonnés de manière non uniforme

L'identification d'un composant de signal dont la fréquence dépasse la limite de Nyquist est un problème difficile en théorie du signal ainsi que dans certains domaines d'applications spécifiques comme l'astronomie et les biosciences. Une conséquence du théorème d'échantillonnage bien connu pour un signal échantillonné uniformément est que la composante spectrale au-dessus de la limite de Nyquist est repliée dans la plage de fréquences inférieure, ce qui rend impossible la distinction entre les composantes vraies et repliées. L'échantillonnage non uniforme, cependant, offre une possibilité de réduire les composants aliasés et de découvrir des informations au-dessus de la limite de Nyquist. Dans cet article, nous proposons une analyse théorique de la réduction des composantes aliasées dans le périodogramme non paramétrique pour deux schémas d'échantillonnage : le modèle d'échantillonnage aléatoire et le modèle d'échantillonnage généré par la modulation de fréquence d'impulsions intégrale (IPFM), ce dernier largement accepté comme la fréquence cardiaque modèle de chronométrage. Une formule générale qui relie la variance des écarts de synchronisation d'un schéma régulier avec la suppression de composant aliasé est proposée. La relation dérivée est illustrée par des périodogrammes de Lomb-Scargle appliqués sur des données simulées. Les données expérimentales présentées consistant en le signal respiratoire dérivé de l'électrocardiogramme et le signal de fréquence cardiaque soutiennent également la possibilité de détecter des fréquences au-dessus de la limite de Nyquist dans la condition connue sous le nom de repliement cardiaque.

1. Introduction

L'aliasing est un phénomène bien connu dans la théorie du signal qui se produit si un signal en temps continu est échantillonné à une fréquence

en dessous du double de la fréquence maximale

du signal échantillonné, c'est-à-dire

. La fréquence d'échantillonnage réduite de moitié,

, est appelée fréquence de Nyquist et peut être interprétée comme la fréquence maximale présumée du signal échantillonné à condition qu'un échantillonnage ait été effectué à une fréquence suffisamment élevée. Si la fréquence du signal échantillonné dépasse la limite de Nyquist, cette fréquence apparaît à une mauvaise position, dans la plage

, sauf si un filtre d'anticrénelage a été appliqué. Ces règles ne sont cependant pas pleinement applicables dans les données échantillonnées de manière non uniforme.

L'échantillonnage non uniforme peut être exploité intentionnellement dans une conception de système, ou il peut refléter des conditions restrictives spécifiques dans la mesure et l'exploration des données dans certains domaines scientifiques, tels que l'astronomie, la sismologie, la génétique, les mesures biophysiques et la physiologie cardiaque. Les données ne peuvent pas être acquises à des instants prescrits, elles sont manquantes ou corrompues. Dans d'autres situations, des temps d'échantillonnage essentiellement irréguliers sont dictés par la nature, comme c'est le cas dans l'analyse de la variabilité de la fréquence cardiaque. L'échantillonnage non uniforme rend l'analyse d'un signal plus compliquée, car la grande majorité des procédures d'analyse de données et de traitement de signal numérique ont été développées pour des signaux échantillonnés de manière uniforme. D'un autre côté, les signaux échantillonnés de manière non uniforme sont plus informatifs et des composants au-delà de la limite de Nyquist peuvent être détectés.

La notion de fréquence de Nyquist a besoin d'être clarifiée lorsque l'on considère des signaux échantillonnés de manière non uniforme. Il n'y a pas de définition bien acceptée si la fréquence d'échantillonnage n'est pas constante [1]. Certains auteurs définissent la fréquence de Nyquist dans le cas non uniforme comme

est le plus petit intervalle de temps entre les données, tandis que d'autres déduisent la fréquence de Nyquist à partir des temps d'échantillonnage moyens

. Une manière plus générale de définir la fréquence de Nyquist pour des échantillons non uniformes utilise un concept de fenêtre spectrale, c'est-à-dire la transformée de Fourier au carré de l'amplitude (FT) du motif d'échantillonnage. La fenêtre spectrale atteint son maximum à la fréquence zéro et un maximum secondaire proche du maximum principal à la fréquence que l'on dit être la fréquence de Nyquist doublée. Dans cet article, nous considérerons des temps d'échantillonnage qui s'écartent des instants réguliers de manière homogène, c'est-à-dire sans grappes marquées ou sans longs écarts. Dans une telle situation, les concepts de temps d'échantillonnage moyen et de fenêtre spectrale donnent des résultats pratiquement identiques.

L'analyse spectrale des données échantillonnées de manière non uniforme a une histoire relativement ancienne en raison de l'importance théorique et pratique de ce problème. Diverses approches ont été proposées pour traiter les échantillons non uniformes. Une revue complète de l'état de l'art peut être trouvée dans [1]. Nous limiterons notre considération aux méthodes couramment utilisées dans la recherche des données de fréquence cardiaque qui sont utilisées dans cet article comme données démonstratives pour illustrer les phénomènes analysés. La méthode la plus simple pour traiter des données échantillonnées de manière non uniforme consiste à ignorer l'irrégularité d'échantillonnage, à déplacer les échantillons dans une grille régulière et à les traiter comme un signal échantillonné uniformément. Une telle approche entraîne une réduction des amplitudes spectrales et l'irrégularité d'échantillonnage se traduit par un bruit spectral. Les composantes au-dessus de la fréquence de Nyquist ne peuvent pas être distinguées des fréquences repliées. Les méthodes basées sur l'interpolation atténuent les fréquences plus élevées, ce qui les rend impropres à l'analyse des composants au-dessus de la fréquence de Nyquist. Une méthode largement utilisée qui ne suppose pas un rééchantillonnage du signal est connue sous le nom de périodogramme de Lomb-Scargle, qui est essentiellement un ajustement par les moindres carrés des données mesurées par des sinusoïdes [2-6].

Le périodogramme de Lomb-Scargle (LS) [2] est une méthode spécialement conçue pour les échantillons acquis de manière non uniforme. Cette méthode a été proposée à l'origine pour l'analyse de séries temporelles astronomiques et a également trouvé une application réussie dans l'analyse de données cardiovasculaires, telles que celles des analyses de la variabilité de la fréquence cardiaque [7, 8], des rythmes circadiens [9] et de la séquence génétique. Analyse. Les résultats de LP sont dans de nombreuses situations presque identiques au périodogramme classique (Schuster) généralisé aux échantillons non uniformes [10] et ils peuvent donc être déduits du FT de données échantillonnées de manière non uniforme. Une telle approche est utilisée dans [6] où plusieurs propriétés de la FT échantillonnée de manière non uniforme sont présentées, bien qu'elles soient appelées de manière imprécise les propriétés de la densité spectrale de puissance de Lomb.

Contrairement à notre article, les analyses [6] ont été effectuées sans une attention particulière pour quantifier un effet du repliement, c'est-à-dire la capacité à découvrir des composants au-dessus de la limite de Nyquist. Cet article étudie l'utilité d'un échantillonnage non uniforme pour réduire à la fois les composantes aliasées et révéler les composantes supérieures au taux de Nyquist. Nous fournirons une analyse théorique d'une situation de modèle, lorsqu'un signal sinusoïdal est échantillonné selon un modèle d'échantillonnage non uniforme prédéfini. Une formule générale qui relie la variance des écarts de synchronisation d'un schéma régulier avec la suppression de composant aliasé est proposée. Une observabilité théoriquement prouvée des composants au-delà de la limite de Nyquist est illustrée par des exemples de simulation et des données cardio-respiratoires obtenues à partir d'une expérience spécialement aménagée à cet effet.

Le plan du document est le suivant. La section 2 résume les définitions des périodogrammes utilisés dans l'analyse spectrale des données non uniformes. Les sections suivantes présentent une analyse spectrale théorique de la sinusoïde échantillonnée avec deux schémas d'échantillonnage. Dans la section 3 de l'article, nous dériverons les formules pour les densités spectrales de puissance (PSD) de signaux sinusoïdaux échantillonnés de manière non uniforme à condition que les temps d'échantillonnage obéissent à un modèle aléatoire. La section 4 fournit un traitement du spectre sinusoïdal lorsque le modèle d'échantillonnage a été généré par la modulation de fréquence d'impulsions intégrale (IPFM). Une formule quantitative décrivant une réduction de composant aliasé est dérivée dans la section 5. Des analyses démonstratives de données simulées sont présentées dans la section 6. Les analyses spectrales des données expérimentales - la respiration dérivée de l'électrocardiogramme (ECG) et les signaux de variabilité de la fréquence cardiaque (HRV) - illustrent performance de la méthode LS dans le cas de ce qu'on appelle le repliement cardiaque, lorsque des composantes au-dessus de la fréquence de Nyquist ont été provoquées par l'arythmie sinusale respiratoire.

2. Périodogrammes pour les signaux échantillonnés de manière non uniforme

Le périodogramme conventionnel basé sur la transformée de Fourier au carré de l'amplitude peut être généralisé également pour les signaux échantillonnés de manière non uniforme sous la forme [4] :

sont des échantillons de données prélevés aux instants

. Notez que les racines d'une analyse de type Fourier avec des exponentielles complexes résident dans un ajustement par les moindres carrés des données à la sinusoïde complexe

. Lomb [2] a formulé le périodogramme au moyen de l'ajustement par les moindres carrés des données observées à valeur réelle au modèle de la forme :

La variable de temporisation redondante

a été introduit dans le modèle (2) afin d'orthogonaliser les termes sinus et cosinus, qui apparaissent dans une formulation des équations normales du problème d'ajustement. Le PSD LS est ensuite calculé comme [4]

À condition que les termes sommés de sinus et de cosinus au carré dans (4) soient à la fois égaux et s'approchent de

, le LS (4) se rapproche du périodogramme classique (1). Généralement, (1) et (4) ne sont pas équivalents. Le LS est préféré en raison de propriétés statistiques utiles pour les tests de signification des pics spectraux [4, 11]. D'autre part, (1) est mathématiquement plus traitable.

3. Modèle d'échantillonnage aléatoire

Le signal analysé prendra la forme d'une sinusoïde complexe d'amplitude

, échantillonné aux instants

Les instants d'échantillonnage sont supposés être les nombres aléatoires obéissant au modèle suivant :

est la période d'échantillonnage moyenne,

est une séquence de variables aléatoires indépendantes et distribuées de manière identique. Le comportement moyen de (1), approximativement valable également pour le spectre LS, peut être déduit en appliquant l'espérance sur la transformée de Fourier de (5) :

Après substitution (5), (6) en (7) et réarrangement, on obtient

Dans (8), nous pouvons reconnaître FT d'une sinusoïde uniformément échantillonnée exprimable par le noyau de Dirichlet périodique :

La moyenne FT (8) peut être réécrite en fonction de la fonction caractéristique décalée par la fréquence

où la fonction caractéristique

de l'irrégularité de synchronisation

, caractérisé par la fonction de densité de probabilité

Par conséquent, la fonction caractéristique atténue les pics périodiques—alias ou images spectrales dans (11), alors qu'elle conserve le vrai pic à la fréquence

L'espérance résultante du FT (11) ne tient pas complètement compte d'un effet aléatoire de l'irrégularité d'échantillonnage. Par conséquent, l'espérance d'une magnitude au carré FT, c'est-à-dire PSD par sa définition, est évaluée pour avoir une vision plus complète sur le spectre d'un signal échantillonné irrégulièrement :

Le terme de sommation dans (13) après substitution (5) est

, on a du fait de l'indépendance présumée :

où nous avons supposé une fonction de probabilité symétrique dans (12) et donc une fonction caractéristique à valeur réelle.

, l'espérance (15) est 1 et elle peut être formellement développée comme

dans (14) et en tenant compte de (15), (16), la sommation (14) se simplifie en

Le terme de sommation dans (17) peut être reconnu comme le FT d'une sinusoïde uniformément échantillonnée multipliée par la fenêtre triangulaire. Afin de faciliter la lisibilité de la formule, ce FT sera noté

La formule finale devient alors

Contrairement à (11), l'expression PSD (19), outre les pics liés aux sinusoïdes, comprend une portion de spectre lisse liée au caractère aléatoire des instants de synchronisation. Depuis

, cette portion lisse présente une vallée proche de la fréquence vraie de la sinusoïde complexe. La composition du PSD à partir des termes présents dans (19) est clarifiée dans la figure 1.


4. Modèle d'échantillonnage généré par le modèle IPFM

Le modèle de modulation de fréquence d'impulsions intégrale (IPFM) est largement accepté comme modèle de synchronisation du nœud sino-auriculaire dans la modélisation de la variabilité de la fréquence cardiaque. C'est un modèle d'intégration et de feu, également utilisé pour décrire une activité neuronale. Le modèle IPFM suppose qu'un signal modulant est intégré et qu'une impulsion de déclenchement de battement est générée lorsque le signal intégré atteint un seuil (Figure 2). Les temps d'occurrence des battements (représentant les temps d'échantillonnage)

généré par le modèle IPFM peut être défini comme [12]

est la période cardiaque moyenne, et

est un signal modulant sans dimension interprété comme le changement fractionnaire de la fréquence cardiaque instantanée par rapport à


Le modèle d'échantillonnage peut être écrit sous la forme d'une séquence d'impulsions de Dirac :

La transformée de Fourier de (21) est appelée spectre des comptes (SPC) [12] :

Le spectre d'un signal échantillonné

est lié au spectre vrai

par la convolution fréquentielle :

Le spectre des dénombrements a été largement analysé dans plusieurs travaux [12–14]. Nous allons résumer ici les faits nécessaires pour évaluer l'effet du repliement spectral. Le SPC défini par (22) peut être réécrit en termes de signal de modulation comme

Fait intéressant, le spectre du signal modulant

est directement inclus dans (24). A côté de l'impulsion de Dirac située à l'origine de la fréquence, un nombre infini de termes modulés en fréquence (FM) convolués avec

sont présents. Les porteuses de ces termes modulés FM produisant des pics en (24) sont situées à des multiples de la fréquence d'échantillonnage moyenne, ce qui explique une origine des composantes repliées dans le spectre (23). Pour le signal modulant sinusoïdal :

est la fréquence de modulation, le spectre (24) peut être étendu à l'aide des fonctions de Bessel du premier type [15]. Les amplitudes spectrales des porteuses suivent la fonction de Bessel d'ordre 0

, the argument of which is given by the gradually increased index of the frequency modulation:

is the mean sampling frequency corresponding to the first carrier

5. Aliasing Reduction Index

For quantification of the aliased components amplitudes, we will restrict our considerations to the situation when the input frequency is in the range

. Unlike the derivation presented in Section 3, a real-valued sinusoid, comprising of positive and negative frequencies, will be considered in this treatment. It is well known that the spectrum of such a uniformly sampled signal is mirrored around the Nyquist frequency. If the frequency of the sampled real-valued signal exceeds the Nyquist limit, this frequency is folded over the Nyquist frequency and appears at a wrong position, in the range

. Specifically, if an input signal has frequency

, the aliased component appears at the frequency

. If the spectrum is plotted in a range above

, have equal amplitudes and there is no possibility to distinguish them. For nonuniformly sampled signal, this statement does not hold, and differences can be observed. A component with higher amplitude is considered to be the true component and other component at mirrored frequency as its alias. We introduce a quantity to describe the observed difference between the true and aliased components and name this quantity as the aliasing reduction index (ARI):

is the spectral amplitude at the true (input) frequency

is the spectral amplitude at the corresponding aliased frequency

. Next, we will be dealing with an evaluation of the proposed ARI for both considered sampling models—the random sampling pattern and the IPFM-generated sampling pattern.

The reduction of an aliased component in the situation with the random sampling according to (6) is deduced directly from (11). For a true component, the argument in (11) is

, while the component aliased into

is originated from the negative input frequency

. The ARI can be then written in terms of the characteristic function as

The result (28) depends on a particular probability distribution function which describes the random sampling jitter. For small-to-moderate values of ARI, a more general version of (28) can be derived using power series expansion of the characteristic function. Taking into account the presumed symmetry of the probability density function and neglecting the high-order terms in the expansion, the approximated characteristic function is expressed by the variance of the timing deviations:

in (28), we obtain an approximate formula for the ARI in term of the relative timing standard deviation

which is independent of a particular probability distribution type.

The proposed index can be evaluated also for the IPFM-generated sampling pattern. The spectrum (23) of a real-valued sinusoid, nonuniformly sampled in IPFM-generated time instants, is composed from SPC shifted in the frequency direction by

. A dominant aliased term is associated with the 1st FM carrier, that appears, after

. Its amplitude is determined by the 0th-order Bessel function of the first kind for

. The resulting formula for the ARI estimation, analogical to (28), becomes

ARI can be approximated by

In order to facilitate a comparison of the obtained formula (33) with the ARI derived for the random sampling model, we will rewrite (33) by using the standard deviation of the timing instants. Integration of the IPFM model formula (20) gives

The deviations of the sampling times from a regular scheme is thus sinusoidal, with amplitude

, and the standard deviation estimated as the RMS value of a continuous-time sinusoid is

Combining (33) with (35) and (26) results in the formula:

Surprisingly, (36) has exactly the same form as (30), despite the essentially different sampling pattern.

6. Simulation Example

In order to illustrate a mechanism of the aliasing in nonuniformly sampled signal and show an extent to which an aliased component reduction can be achieved, we have performed a spectral analysis of sequences synthesized according to both described sampling models. Since this work is inspired by analysis of the cardiac aliasing effect, we have arranged the signal parameters—duration, frequency—as to fit the experimental data presented in subsequent section. Each of the simulated signals consists of 120 samples of real-valued, unit-amplitude sinusoid. The frequencies are expressed as the fractions of the mean sampling rate

. The mean sampling rate is considered as the reciprocal value of the mean sampling interval. An input frequency exceeding the Nyquist limit was set to 0.65

which is aliased to the frequency 0.35

. Conversely, an input frequency 0.35

below the Nyquist limit 0.5

and an informationless peak in the spectrum can be observed when a periodogram is plotted above the Nyquist limit. Both spectral peaks appear as identical when the sampling is uniform. The nonuniform sampling, however, introduces differences between the true component and the unwanted—aliased or imaged components. This effect can be observed in the spectra plotted in the range

instead of a conventionally used range

Figures 3 and 4 show LS periodograms for both scenarios (below and above Nyquist limit) for the sampling times randomly deviated from regular positions. Spectral amplitudes are plotted as the square roots of the LS periodograms and scaled to directly reflect amplitudes of the generated sinusoids. The sampling irregularity follows the Gaussian distribution with the variance set according to a required irregularity level. Three levels were used in simulation: low, medium, high. The low irregularity is represented by sampling jitter with the standard deviation 5% of the mean sampling time

. Such a degree of the sampling nonuniformity caused only a minor spectral amplitude difference (about 5%). The medium level is represented by two-fold reduction of alias power (3 dB), that is, ARI

0.707. The formula derived in Section 5 allowed to estimate

required in simulation as

. A maximum deviation of a sampling instant is theoretically unlimited therefore, we incorporate a practical limit as the irregularity level which ensures a monotonic increase of the sampling times. This practical maximum is used as the high level of the irregularity. Since Gaussian variables are not confident to finite interval, determination of the high level deviations must be made in probabilistic sense. Le choix

ensures that randomly chosen proximate sampling times

, fulfill the monotonicity condition

with probability 0.99. Presented figures show a gradual destruction of the aliased component (Figure 4) and spectral images (Figure 3) as the sampling irregularity level increases. In the high-level irregularity, the aliasing is effectively destroyed—the aliased peak is buried in a noise caused by sampling irregularity. Notice that a noise-like spectrum increase is evident as the sampling irregularity increases. A valley near the true frequency in noise PSD, which may be expected by examination of (19), is not evident, because it is overlapped by the smooth noise spectrum coming from the negative frequency.


(a)
(b)
(c)
(a)
(b)
(c) Spectral amplitudes of the LS spectra—nonuniformly sampled signal with the random sampling pattern. The input frequency 0.35

(a)
(b)
(c)
(a)
(b)
(c) Spectral amplitudes of the LS spectra—nonuniformly sampled signal with the random sampling pattern. The input frequency 0.65

A similar set of simulations (Figures 5 and 6) was performed for the sinusoid sampled by the pattern generated by the IPFM model driven by a sinusoidal modulating signal. The sampling time instants were obtained by means of numerical solution of the IPFM equation (34) using MATLAB function fzero. The frequency of the sampled sinusoid was set equal to the frequency of the modulating signal, a condition occurred in heart rate modulation by the respiration signal. The low and medium irregularity levels were established in the same way as in the random sampling pattern. The maximum irregularity level is considered when the normalized modulation signal amplitude

reaches the value1. Due to this restriction, the achieved attenuation cannot be as high as in the random sampling pattern scenarios. Notice that numbers of spurious peaks can be found in the spectrum which are explained by a complex structure of the SPC (24).


(a)
(b)
(c)
(a)
(b)
(c) Spectral amplitudes of periodograms—nonuniformly sampled signal with IPFM generated sampling pattern. Input frequency 0.35

(a)
(b)
(c)
(a)
(b)
(c) Spectral amplitudes of periodograms—nonuniformly sampled signal with IPFM generated sampling pattern. Input frequency

7. Cardiorespiratory Data Analysis

Beat-to-beat variations of consecutive heart periods or the instantaneous heart rate, known as the heart rate variability (HRV), become intensively studied over the last decades. Analysis of the HRV provides a noninvasive insight into the cardiovascular neural control and a spectral analysis of HRV signals is a widely used method for assessment of the autonomous nervous system. Generally, the power spectra of HRV can be divided into the high-frequency (HF, 0.15–0.4 Hz), low-frequency (LF, 0.04–0.15 Hz), and very-low-frequency range (VLF, 0.003–0.04 Hz). The LF power is modulated by both sympathetic and parasympathetic controls, while the HF power mainly reflects the parasympathetic influence linked to the respiration. Therefore, a concurrent measurement of the respiration signal is helpful for HRV interpretation. To obtain the respiration activity signal, an additional specialized instrument is not necessary because the respiration can be estimated from modulation of an electrocardiogram amplitude. Such a procedure is referred to as ECG-derived respiration [16–18]. Notice that a raw ECG signal itself includes spectral components related to heart rate oscillations and respiration [19], but they are mixed in the spectrum in a complex way [20], and thus, appropriate signals must be derived from measured ECG. Both the HRV as well as the respiration signal derivation incorporate heart beat detection, and the beat occurrence times define sampling instants, irregularly spaced by nature. Although a heart rate (HR) is essentially a discrete-time signal, computed from beat-to-beat occurrence times, it can be considered as a hypothetically continuous signal representing the autonomous nervous system activity which continually modulates the heart rate. Unlike an artificially designed system, a physiological system does not include antialiasing filter. For example, 0.4 Hz frequency considered as the highest limit of the standardized HF range requires sampling at a mean heart rate 0.8 Hz, and thus, the mean heart rate of at least 48 bpm is required for reliable spectral analysis. Therefore, a bradycardia or an extended HF range caused by the elevated respiration rate can lead to the aliasing. For this kind of “physiological aliasing,” the term cardiac aliasing has been naturalized. High frequency components modulating the heart rate have been actually observed in neonates, transplanted heart patients, and animals [21, 22]. Improper signal processing or ignoring rarely occurred conditions can then result in incorrect interpretation of a spectral analysis.

To test an ability of different methods to detect frequencies beyond the Nyquist rate in a real-world measured data, we have performed an experiment exploiting respiratory sinus arrhythmia. The measured subject was asked to breathe at an elevated rate in synchrony with a moving pattern displayed on the computer screen at the frequency of 0.8 Hz. An electrocardiogram (ECG) was recorded and two signals were extracted: the respiration and the instantaneous heart rate (HR). The respiration signal was obtained as the amplitude in bandpass-filtered ECG curve taken at maxima of QRS complexes [23]. The HR was computed from beat-to-beat time distances of the QRS maxima. Both signals are thus sampled at a variable rate dictated by the heart rate. The mean heart rate of 74 bpm (1.23 Hz) was sufficiently low to observe an aliasing.

In Figure 7(a), the respiration signal is analyzed as a sequence of values equidistantly spaced at multiples of the mean heart period. A meaningful spectral range is thus up to 0.615 Hz provided that the mean heart rate was 1.23 Hz. The respiration frequency 0.8 Hz is aliased near 0.4 Hz and spreads over a wider spectral range due to wandering of the mean heart rate. The respiration sequence was then interpolated to 4 Hz sampling frequency by means of the cubic spline interpolation. A small peak at the true respiratory frequency can be detected but the aliased component near 0.4 Hz largely exceeds the true peak in its amplitude—Figure 7(b). The Lomb-Scargle periodogram manifests a clearly visible peak at 0.8 Hz, Figure 7(c), albeit the components at aliased frequency region do not disappear completely.


(a)
(b)
(c)
(a)
(b)
(c) Spectral analysis of the ECG-derived respiration signal: (a)equidistantly repositioned QRS maxima samples with frequency axis scaled according to the mean heart period, (b) spline interpolated unevenly sampled QRS maxima, (c) the Lomb-Scargle periodogram of unevenly sampled QRS maxima.

The results of application of different methods for HRV spectrum computation are summarized in Figure 8. Repositioning of unevenly spaced samples, similarly as it was observed in Figure 7(a), yields a symmetrical spectrum that is unable to detect components above the mean heart rate halved. The cubic spline interpolation, although frequently used in HRV analysis, noticeably attenuates high-frequency components above the mean Nyquist frequency. The simplest interpolation type, nearest neighbor method, gives surprisingly better outcome, comparable to the Berger resampling method [24]. As in the case of the respiration signal analysis, only the last two methods presented, the Lomb-Scargle periodogram and the SPC, are capable to show a higher amplitude at the true position 0.8 Hz than near the aliased frequency (about 0.4 Hz). The apparent increase of the amplitudes at higher frequencies in Figure 8(f) can be explained by low-pass filtering effect of the integrator in the IPFM model that does not affect the spectrum of counts. In the low frequency range, all methods give similar results.


(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f) HRV spectra obtained by means of different methods: (a)equidistantly repositioned HR samples with frequency axis scaled according to the mean heart period unevenly sampled heart rate sequence resampled to evenly spaced signal- (b) cubic-spline interpolation, (c) nearest neighbour interpolation, (d) Berger method spectra computed without resampling-, (e) Lomb-Scargle periodogram, (f) spectrum of counts.

8. Conclusion

The fact that a spectrum of nonuniformly sampled signal conveys useful information above the Nyquist limit has been pointed out in several works [25, 26]. In this work, we have presented a theoretical treatment which quantitatively explains this phenomenon. The analyses of simulated and experimental data illustrate a possibility to identify components beyond the Nyquist limit in real-world problems.

Although the spectrum of a nonuinformly sampled signal is not alias free, unlike a uniform sampling situation, the aliased spectral amplitudes are attenuated when compared to the true components amplitudes. Analogically, components below the Nyquist limit can be observed also as images above this limit. The analytic treatment of these unwanted components constitutes a key contribution of this paper. We have studied the aliasing effect in two sampling schemes: random deviations from regular times and sampling times generated by the IPFM model. The former was chosen as a basic model which allowed to explain a mechanism of aliasing suppression. The later is inspired by a model used to describe generation of events in biological systems, such as the sinus node or the neuronal firing.

The two essentially distinctive sampling models required different way of analysis. The random model analysis relies on terms common in probability/statistics. We have shown that the characteristic function of timing jitter plays a key role in composition of the resulting spectrum. The characteristic function modulates amplitudes of aliased or imaged spectral lines. Besides the line portion of the spectrum, a smooth noise portion shaped by a term of characteristic function is present due to a random nature of the sampling process, despite the sampled signal is deterministic. The IPFM sampling scheme uses concepts originated in communication systems theory and presented analysis is purely deterministic. Consequently, the smooth spectrum is not present on the other hand, spurious peaks frequently occur due to a complex nature of the spectrum of the modulated signal. In the both situations, the relative reduction of aliased components was evaluated by a quantity named in this paper as the aliasing reduction index. Interestingly, the derived formula approximating this index for small-to-moderate sampling irregularity is identical for the both sampling models. The formula which we have found seems to be quite universal and relates a decrease of unwanted aliased components with the standard deviation of timing irregularity measured as the fraction of the mean sampling period. Validity of the formula in more general sampling schemes, such as those with dependent sampling irregularities and additive random sampling, is intended to be studied in a future work.

Mitigation of the aliasing effect by means of nonuniform sampling not only attracts scientific interest but also has practical implications. We have demonstrated the possibility to detect frequencies above the Nyquist limit in the ECG-derived respiration signal and the heart rate signal. The condition of aliasing, so-called cardiac aliasing, was elicited by increasing the respiration rate, and signals derived from recorded electrocardiogram were analyzed by different methods. Conventionally used interpolation/resampling methods were found to be unsuitable in the condition of cardiac aliasing. The Lomb-Scarglre method, the classical periodogram used for nonuniform samples and the spectrum of counts are all capable to uncover components above the Nyquist frequency. Since real-world measurements are usually contaminated by noise, a signal to be detected needs not to be a pure sinusoid, or it can be random and even nonstationary, development of a universal method able to resolve between true and aliased spectral components, with properly assigned significance measure, is another challenging task intended for future work.

Les références

  1. P. Babu and P. Stoica, “Spectral analysis of nonuniformly sampled data𠅊 review,” Digital Signal Processing, vol. 20, no. 2, pp. 359–378, 2010. View at: Publisher Site | Google Scholar
  2. N. R. Lomb, “Least-squares frequency analysis of unequally spaced data,” Astrophysics and Space Science, vol. 39, no. 2, pp. 447–462, 1976. View at: Publisher Site | Google Scholar
  3. P. Stoica, J. Li, and H. He, “Spectral analysis of nonuniformly sampled data: a new approach versus the periodogram,” IEEE Transactions on Signal Processing, vol. 57, non. 3, pp. 843–858, 2009. View at: Publisher Site | Google Scholar
  4. J. D. Scargle, “Studies in astronomical time series analysis. II—statistical aspects of spectral analysis of unevenly spaced data,” Journal d'astrophysique, vol. 263, pp. 835–853, 1982. View at: Publisher Site | Google Scholar
  5. P. Vaníპk, “Approximate spectral analysis by least-squares fit-Successive spectral analysis,” Astrophysics and Space Science, vol. 4, no. 4, pp. 387–391, 1969. View at: Publisher Site | Google Scholar
  6. P. Laguna, G. B. Moody, and R. G. Mark, “Power spectral density of unevenly sampled data by least-square analysis: performance and application to heart rate signals,” IEEE Transactions on Biomedical Engineering, vol. 45, no. 6, pp. 698–715, 1998. View at: Publisher Site | Google Scholar
  7. “Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology,” Circulation, vol. 93, no. 5, pp. 1043–1065, 1996. View at: Google Scholar
  8. G. G. Berntson, J. T. Bigger, D. L. Eckberg et al., “Heart rate variability: origins, methods, and interpretive caveats,” Psychophysiology, vol. 34, no. 6, pp. 623–648, 1997. View at: Google Scholar
  9. M. Teplan, L. Molლn, and M. Zeman, “Spectral analysis of cardiovascular parameters of rats under irregular light-dark regime,” in Proceedings of the 8th International Conference on Measurement, pp. 343–346, Smolenice, Slovak Republic, 2011. View at: Google Scholar
  10. V. V. Vityazev, “Time series analysis of unequally spaced data: intercomparison between estimators of the power spectrum,” in Proceedings of the Astronomical Data Analysis Software and Systems VI ASP Conference Series, vol. 125, pp. 166–169, 1997. View at: Google Scholar
  11. A. Schwarzenberg-Czerny, “The distribution of empirical periodograms: Lomb-Scargle and PDM spectra,” Avis mensuels de la Royal Astronomical Society, vol. 301, no. 3, pp. 831–840, 1998. View at: Google Scholar
  12. J. Mateo and P. Laguna, “Improved heart rate variability signal analysis from the beat occurrence times according to the IPFM model,” IEEE Transactions on Biomedical Engineering, vol. 47, no. 8, pp. 985–996, 2000. View at: Publisher Site | Google Scholar
  13. M. Brennan, M. Palaniswami, and P. Kamen, “Distortion properties of the interval spectrum of IPFM generated heartbeats for heart rate variability analysis,” IEEE Transactions on Biomedical Engineering, vol. 48, no. 11, pp. 1251–1264, 2001. View at: Publisher Site | Google Scholar
  14. J. Púčik, O. Ondráპk, E. Cocherová, and A. Sultan, “Spectrum of counts computation for HRV analysis,” in Proceedings of 19th International Conference Radioelektronika (RADIOELEKTRONIKA '09), pp. 255–258, April 2009. View at: Publisher Site | Google Scholar
  15. E. J. Bayly, “Spectral analysis of pulse frequency modulation in the nervous systems,” IEEE Transactions on Biomedical Engineering, vol. 15, no. 4, pp. 257–265, 1968. View at: Google Scholar
  16. G. Moody, R. Mark, A. Zoccola, and S. Mantero, “Derivation of respiratory signals from multi-lead ECGs,” IEEE Computer Society, vol. 12, pp. 113–116, 1985. View at: Google Scholar
  17. G. D. Clifford, F. Azuaje, P. E. McSharry et al., Advanced Methods and Tools for ECG Data Analysis, Artech House, 2006.
  18. D. Widjaja, J. Taelman, S. Vandeput et al., “ECG-derived respiration: comparison and new measures for respiratory variability,” in Proceedings of the Computing in Cardiology (CinC '10), pp. 149–152, September 2010. View at: Google Scholar
  19. A. Gersten, O. Gersten, A. Ronen, and Y. Cassuto, “The RR interval spectrum, the ECG signal and aliasing,” . Eprint, http://arxiv.org/abs/physics/9911017v1. View at: Google Scholar
  20. J. Šurda, S. Lovás, J. Púčik, and M. Jus, “Spectral properties of ECG signal,” in Proceedings of the of the 17th International Conference Radioelektronika, pp. 537–541, Brno, Czech Republic, 2007. View at: Google Scholar
  21. E. Toledo, I. Pinhas, D. Aravot, and S. Akselrod, “Very high frequency oscillations in the heart rate and blood pressure of heart transplant patients,” Medical and Biological Engineering and Computing, vol. 41, no. 4, pp. 432–438, 2003. View at: Google Scholar
  22. H. A. Campbell, J. Z. Klepacki, and S. Egginton, “A new method in applying power spectral statistics to examine cardio-respiratory interactions in fish,” Journal of Theoretical Biology, vol. 241, no. 2, pp. 410–419, 2006. View at: Publisher Site | Google Scholar
  23. J. Púčik, M. Uhrík, A. Sultan, and A. Šurda, “Experimental setup for cardio-respiratory interaction study,” in Proceedings of the 8th Czech-Slovak Conference, Trends in Biomedical Engineering, pp. 126–129, Bratislava, Slovakia, 2009. View at: Google Scholar
  24. R. Berger, S. Akselrod, D. Gordon, and R. Cohen, “An efficient algorithm for spectral analysis of heart rate variability,” IEEE Transactions on Bio-Medical Engineering, vol. 33, no. 9, pp. 900–904, 1986. View at: Google Scholar
  25. W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling et al., Numerical Recipes in Fortran 77: The Art of Scientific Computing, vol. 1, Cambridge University Press, 1992.
  26. G. L. Bretthorst, “Nonuniform sampling: bandwidth and aliasing,” Concepts in Magnetic Resonance A, vol. 32, no. 6, pp. 417–435, 2008. View at: Publisher Site | Google Scholar

Droits d'auteur

Copyright © 2012 Jozef Púčik et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Note: This Plot Contains an Error¶

After the book was in press, a reader pointed out that this plot contains a typo. Instead of passing the noisy data to the Lomb-Scargle routine, we had passed the underlying, non-noisy data. This caused an over-estimate of the Lomb-Scargle power.

Because of this, we add two extra plots to this script: the first reproduces the current plot without the typo. In it, we see that for the noisy data, the period is not detected for just

30 points within ten periods.

In the second additional plot, we increase the baseline and the number of points by a factor of ten. With this configuration, the peak is detected, and the qualitative aspects of the above discussion hold true.


1. Introduction

The functioning of eukaryotic cells is controlled by accurate timing of biological cycles, such as cell cycles and circadian rhythms. These are composed of an echelon of molecular events and checkpoints. At the transcription level, these events can be quantitatively observed by measuring the concentration of messenger RNA (mRNA), which is transcribed from DNA and serves as the template for synthesizing protein. To achieve this goal, in the microarray experiments, high-throughput gene chips are exploited to measure genome-wide gene expressions sequentially at discrete time points. These time series data have three characteristics. Firstly, most data sets are of small sample size, usually not more than 50 data points. Large sample sizes are not financially affordable due to high cost of gene chips. Also the cell cultures lose their synchronization and render meaningless data after a period of time. Secondly, the data are usually evenly sampled and have many time points missing. Thirdly, most data sets are customarily corrupted by experimental noise and the produced uncertainty should be addressed in a stochastic framework.

Extensive genome-wide time course microarray experiments have been conducted on organisms such as Saccharomyces cerevisiae (budding yeast) [1], human Hela [2], and Drosophila melanogaster (fruit fly) [3]. Budding yeast in [1] has served as the predominant data source for various statistical methods in search of periodically expressed genes, mainly due to its pioneering publication and relatively larger sample size compared with its peers. By assuming the signal in the cell cycle to be a simple sinusoid, Spellman et al. [1] and Whitfield et al. [2] performed a Fourier transformation on the data sampled with different synchronization methods, while Giurcaneanu [4] explored the stochastic complexity of the detection mechanism of periodically expressed genes by means of generalized Gaussian distributions. Ahdesmäki et al. [5] implemented a robust periodicity testing procedure also based on the non-Gaussian noise assumption. Alternatively, Luan and Li [6] employed guide genes and constructed cubic B-spline-based periodic functions for modeling, while Lu et al. [7] employed up to three harmonics to fit the data and proposed a periodic normal mixture model. Power spectral density estimation schemes have also been employed. Wichert et al. [8] applied the traditional periodogram on various data sets. Bowles et al. [9] compared Capon and robust Capon methods in terms of their ability to identify a predetermined frequency using evenly sampled data sets, under the assumption of a known period. Lichtenberg et al. [10] compared [167] while proposing a new score by combining the periodicity and regulation magnitude. The majority of these works dealt with evenly sampled data. When missing data points were present, either the vacancies were filled by interpolation in time domain, or the genes were discarded if there were more than 30% data samples missing.

Biological experiments generally output unequally spaced measurements. The major reasons are experimental constraints and event-driven observation. The rate of measurement is directly proportional to the occurrence of events. Therefore, an analysis based on unevenly sampled data is practically desired and technically more challenging. While providing modern spectral estimation methods for stationary processes with complete and evenly sampled data [11], the signal processing literature has witnessed an increased interest in analyzing unevenly sampled data sets, especially in astronomy, in the last decades. The harmonics exploited in discrete Fourier transform (DFT) are no longer orthogonal for uneven sampling. However, Lomb [12] and Scargle [13] demonstrated that a phase shift suffices to make the sine and cosine terms orthogonal. The Lomb-Scargle scheme has been exploited in analyzing the budding yeast data set by Glynn et al. [14]. Schwarzenberg-Czerny [15] employed one-way analysis of variance (AoV) and formulated an AoV periodogram as a method to detect sharp periodicities. However, it relies on an infeasible biological assumption, that is, the observation duration covers many cycles. Along with this line of research, Ahdesmäki et al. [16] proposed to use robust regression techniques, while Stoica and Sandgren [17] updated the traditional Capon method to cope with the irregularly sampled data. Wang et al. [18] reported a novel technique, referred to as the missing-data amplitude and phase estimation (MAPES) approach, which estimates the missing data and spectra iteratively through the expectation maximization (EM) algorithm. In general, Capon and MAPES methods possess a better spectral resolution than Lomb-Scargle periodogram. In this paper, we propose to analyze the performance of three of the most representative spectral estimation methods: Lomb-Scargle periodogram, Capon method, and the MAPES technique in the presence of missing samples and irregularly spaced samples. The following questions are to be answered in this study: do technically more sophisticated schemes, such as MAPES, achieve a better performance on real biological data sets than on simpler schemes? Is the efficiency sacrificed in using these advanced methods justifiable?

The remainder of this paper is structured as follows. In Section 2, we introduce the three spectral analysis methods, that is, Lomb-Scargle, Capon and MAPES. Hypothesis tests for periodicity detection and the corresponding -values are also formulated. The multiple testing correction is discussed. Section 3 presents simulation results. The performances of the three schemes are compared based on published cell-cycle and noncell-cycle genes of the Saccharomyces cerevisiae (budding yeast). Then the spectral analysis for the data set of Drosophila melanogaster (fruit fly) is performed, and a list of 149 genes are presented as cycle-related genes. The synchronization effects are also considered. Concluding remarks and future works constitute the last section, and full results are provided online in the supplementary materials [19].


Time-uncertain spectral analysis

Now let’s consider age uncertainties.

The GeoChronR approach to quantifying and propagating those uncertainties is to leverage the power of ensembles. Here we will illustrate this with MTM. Let us repeat the previous analysis by looping over the 1,000 age realizations output by the tuning algorithm HMM-Match. That means computing 1000 spectra. We’ve already seen that nuspectral is too slow for an ensemble job, so let’s leave it out. Also, since REDFIT is a tapered version of Lomb-Scargle, and comes with more features, let’s focus here on comparing an REDFIT and MTM in ensemble mode.

This took a few seconds with 1000 ensemble members, and the resulting output is close to 10 Mb. Not an issue for modern computers, but you can see why people weren’t doing this in the 70’s, even if they wanted to. Now, let’s plot the result:

To this we can add the periods identified as significant owing to MTM’s harmonic ratio test. geoChronR deals with it by computing at each frequency the fraction of ensemble members that exhibit a significant peak. One simple criterion for gauging the level of support for such peaks given age uncertainties is to pick out those periodicities that are identified as significant more than 50% of the time.

One could of course, impose a more stringent criterion (e.g. 80%, 90%, etc). To do that, specify which(mtm.ens.spec$sig.freq>0.8) or which(mtm.ens.spec$sig.freq>0.9) . That relatively close peaks appear significant may be an artifact of neglecting the multiple hypothesis problem (cf Vaughan et al, 2011).

So what do we find? Consistent with previous investigations (e.g. Mudelsee et al 2009), the effect of age uncertainties is felt more at high frequencies, with the effect of broadening peaks, or lumping them together. Again, the peaks that rise above the power-law background are the

40kyr peaks, which is not surprising. There is significant energy around the precessional peaks, but it is rather spread out, and hard to tie to particular harmonics. No doubt a record with higher-resolution and/or tighter chronology would show sharper precessional peaks.

Do we get the same answer in REDFIT? Let’s find out.

again, it took a bit of time. let’s plot the result:

The result is quite a bit smoother here, perhaps because of the choice of parameters.


Making sense of the lomb-scargle periodogram - Astronomy

Context. Although timing variations in close binary systems have been studied for a long time, their underlying causes are still unclear. A possible explanation is the so-called Applegate mechanism, where a strong, variable magnetic field can periodically change the gravitational quadrupole moment of a stellar component, thus causing observable period changes. One of the systems exhibiting such strong orbital variations is the RS CVn binary HR 1099, whose activity cycle has been studied by various authors via photospheric and chromospheric activity indicators, resulting in contradicting periods.
Aims: We aim at independently determining the magnetic activity cycle of HR 1099 using archival X-ray data to allow for a comparison to orbital period variations.
Methods: Archival X-ray data from 80 different observations of HR 1099 acquired with 12 different X-ray facilities and covering almost four decades were used to determine X-ray fluxes in the energy range of 2-10 keV via spectral fitting and flux conversion. Via the Lomb-Scargle periodogram we analyze the resulting long-term X-ray light curve to search for periodicities.
Results: We do not detect any statistically significant periodicities within the X-ray data. An analysis of optical data of HR 1099 shows that the derivation of such periods is strongly dependent on the time coverage of available data, since the observed optical variations strongly deviate from a pure sine wave. We argue that this offers an explanation as to why other authors derive such a wide range of activity cycle periods based on optical data. We furthermore show that X-ray and optical variations are correlated in the sense that the star tends to be optically fainter when it is X-ray bright.
Conclusions: We conclude that our analysis constitutes, to our knowledge, the longest stellar X-ray activity light curve acquired to date, yet the still rather sparse sampling of the X-ray data, along with stochastic flaring activity, does not allow for the independent determination of an X-ray activity cycle.


Voir la vidéo: Introduction à la spectrophotométrie (Juillet 2021).