Astronomía

Dar sentido al periodograma de lomb-escarlata

Dar sentido al periodograma de lomb-escarlata


We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

Estoy tratando de usar el periodograma para saber cuándo una señal es periódica o no siguiendo el tutorial para el periodograma de astropía Lomb-scargle aquí.

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

Simulé algunos datos, uno que es un sinusoide (período = 200) y otro que es un sesgo gaussiano (es decir, un solo evento transitorio). La esperanza era que el periodograma seleccionara el período para el objeto periódico y le diera un período para el transitorio que implicaría que ocurrió un solo evento en la ventana.

Desafortunadamente, los resultados no tienen ningún sentido. He adjuntado el código al final y las figuras generadas a continuación. Debido al ruido aleatorio en mi simulación, cada resultado es diferente y a continuación proporciono dos ejemplos de los resultados. Utilizo el mismo método descrito en el enlace de arriba donde usamos la funciónLombScargle.model ()en la mejor frecuencia de ajuste defrecuencia [argmax [potencia]]

La línea roja es la verdadera función de la que simulé los datos. El verde es el que mejor se ajusta al periodograma. Los gráficos de la derecha son el PSD del periodograma.

Ejemplo 1

Aquí podemos ver que la frecuencia de mejor ajuste para la sinusoide (gráfico superior derecho) seleccionada es 0.105 (es decir, un período de 9-10 días) que no está cerca de una frecuencia de 0.05 que esperaría para algo con un período de 200 días , sin embargo, cuando le doy esta frecuencia de mejor ajuste de 0.105 al ajustador del modelo lomb-scargle, ¿aparece una curva periódica que coincide con un período de 200 días?

Esto no tiene sentido.

Ejemplo 2

Aquí volví a ejecutar el código y esta vez los resultados se cambiaron. Se ajusta a un período muy grande para lo transitorio, por lo que puedo decir con confianza que es un evento transitorio único, pero el ajuste sinusoidal es terrible. La mejor frecuencia sigue siendo 0,105 (período = 10), pero el ajustador del modelo lomb-scargle superpone algo que parece tener un período de 60 días, ¿qué está mal?

¿Puedo obtener una aclaración sobre si estoy haciendo algo mal? Me han dicho que el periodograma es la herramienta de facto para datos muestreados de manera desigual como esta todavía ... los resultados parecen horribles la mitad del tiempo.

Para aclarar mis dudas

  1. ¿Se puede explicar cómo en el primer gráfico la frecuencia de mejor ajuste de 0.105 que introduzco en el ajustador del modelo lomb-escarlata de astropía crea de alguna manera una sinusoide que coincide correctamente con una frecuencia de 0.05? Cual es la explicacion?

  2. ¿Por qué hay 5 picos fuertes en la gráfica superior derecha para el periodograma en el ejemplo 1 cuando solo espero 1? Los dos del medio están cerca del valor real de 0.05 (en 0.045 y 0.055)

Aquí está el código corto utilizado para simular y trazar los datos.

import numpy as np import matplotlib.pyplot as plt import scipy.stats as sc from astropy.stats import LombScargle import math #simulate parameters data_range = [i for i in range (1,1001)] number_of_samples = 50 gauss_skew = sc.skewnorm. pdf sesgo = -10 período = 200 ubicación = rango_datos [int (len (rango_datos) / 2)] y1 = [(2 * (1. + np.sin (2. * np. pi * x / período)) + np .random.normal (loc = 0.0, scale = 0.5)) for x in data_range] errors1 = [np.random.normal (loc = 0.0, scale = 1) for x in data_range] y2 = [(1000 * (gauss_skew ( x, sesgo, loc = ubicación, escala = 50)) + np.random.normal (loc = 0.0, scale = 1)) para x en data_range] errores2 = [np.random.normal (loc = 0.0, scale = 1 ) for x in data_range] sample_rate = int (len (data_range) / number_of_samples) # Para reducir los datos un poco y1 = y1 [:: sample_rate] y2 = y2 [:: sample_rate] errors1 = errors1 [:: sample_rate] errors2 = errores2 [:: tasa_muestra] x1 = x2 = rango_datos [:: tasa_muestra] verdad1 = [2 * (1. + np.sin (2. * np. pi * x / período)) para x en rango_datos] verdad2 = [1000 * (gau ss_skew (x, -10, loc = location, scale = 50)) for x in data_range] verdades = [verdad1, verdad2] x = [x1, x2] y = [y1, y2] errores = [errores1, errores2] 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 ('frecuencia') ax [1] [1] .set_xlabel ('frecuencia') para i en el rango (0,2): frecuencia, potencia = LombScargle (x [i], y [i], errores [i]) .autopower () #Obtenga la mejor frecuencia de ajuste como en el tutorial mejor_frecuencia = frecuencia [np.argmax (potencia)] print ('mejor frecuencia:', mejor_frecuencia) t_fit = np.linspace (x [i] [0], matemáticas .floor (x [i] [- 1]), num = 50) #Ajustar la mejor frecuencia de ajuste #trazar el mejor mejor modelo basado en el mejor ajuste y_fit = LombScargle (x [i], y [i], errores [ i]). model (t_fit, best_frequency) ax [i] [0] .plot (t_fit, y_fit, 'g') ax [i] [0] .plot (data_range, verdades [i], 'r') # Trace la PSD ax [i] [1] .plot (frecuencia, potencia) ax [i] [1] .axvline (x = best_frequency, color = 'black', ls = "-") plt.show ()

Creo que veo lo que está pasando. Su mejor ajuste es en realidad una frecuencia de pulsación entre la frecuencia real (0,005) y el doble de la frecuencia de muestreo (0,05). De hecho, esto produce un modelo que pasa por sus puntos de datos, pero debido a que solo ha trazado 50 puntos en el modelo, no ha podido verlo.

Si cambia esta línea t_fit = np.linspace (x [i] [0], math.floor (x [i] [- 1]), num = 1000)

para que muestre el modelo en más puntos, puede ver que, de hecho, está en una frecuencia mucho más alta.

El truco consiste en conocer (aproximadamente) la respuesta antes de empezar.

Si limita su rango de frecuencia en el comando Lomb-Scargle por debajo de la frecuencia de Nyquist:

frecuencia, potencia = LombScargle (x [i], y [i], errores [i]). autopower (nyquist_factor = 1.5)

Encontrarás el comportamiento que esperabas. p.ej. (Nota: reduje el tamaño de la barra de error para centrarme en su problema, pero también arreglé el tamaño de la barra de error a un solo positivo valor. Su código estaba introduciendo errores muy pequeños e incluso negativos y no sé cómo trataría el código con ellos).

En datos espaciados desigualmente, las cosas no son tan simples porque la frecuencia de Nyquist no está bien definida.


Nota: esta gráfica contiene un error¶

Después de que el libro estuvo en prensa, un lector señaló que esta trama contiene un error tipográfico. En lugar de pasar los datos ruidosos a la rutina Lomb-Scargle, pasamos los datos subyacentes no ruidosos. Esto provocó una sobreestimación del poder de Lomb-Scargle.

Debido a esto, agregamos dos gráficos adicionales a este script: el primero reproduce el gráfico actual sin errores tipográficos. En él, vemos que para los datos ruidosos, el período no se detecta por solo

30 puntos en diez períodos.

En el segundo gráfico adicional, aumentamos la línea de base y el número de puntos en un factor de diez. Con esta configuración, se detecta el pico y los aspectos cualitativos de la discusión anterior son verdaderos.


Programa AFNI: 3dLombScargle

Haga un periodograma o espectro de amplitud de una serie de tiempo que tenga un
tasa de muestreo no constante. Los espectros de salida de este programa son
'unilateral', de modo que representan la media amplitud o potencia
asociado con una frecuencia, y requerirían un factor de 2 para
tener en cuenta las soluciones de frecuencia de desplazamiento a la derecha y a la izquierda
de la transformada de Fourier (ver más abajo 'SALIDA' y 'NOTA').

De particular interés es la aplicación de esta funcionalidad a
series de tiempo en estado de reposo que pueden haber sido censuradas. La teoría detrás
las matemáticas y los algoritmos de esto se deben a grupos separados, principalmente
en el ámbito de las aplicaciones astrofísicas: Vaníček (1969, 1971),
Lomb (1976), Scargle (1982) y Press & amp Rybicki (1989). Grítales.

Esta implementación particular se debe a Press & amp Rybicki (1989), por
esencialmente traduciendo su implementación Fortran publicada a C,
mientras usa GSL para la FFT, en lugar de realft () de NR, y hace
varios ajustes basados ​​en eso.

La adaptación de Lomb-Scargle se realizó con cambios bastante mínimos aquí por
PA Taylor (v1.4, junio de 2016).

+ USO:
Ingrese una serie de tiempo volumétrica 4D (conjunto de datos BRIK / HEAD o NIFTI)
así como un archivo 1D opcional de 0 y 1 que define qué puntos
censurar (es decir, cada 0 representa un punto / volumen para censurar)
Si no se ingresa ningún archivo 1D, el programa buscará volúmenes que estén
uniformemente cero y considerarlos censurados.

La salida es un periodograma LS, que describe magnitudes espectrales
hasta una 'frecuencia máxima': el máximo predeterminado aquí es lo que
la frecuencia de Nyquist de la serie de tiempo * habría sido * sin
cualquier censura. (Curiosamente, este análisis puede ser
Aplicado legítimamente en casos para estimar contenido de frecuencia & gtNyquist.
¡Guau!)

El espectro de frecuencia estará en el rango [df, f_N], donde:
df = 1 / T, y T es la duración total de la serie de tiempo sin censura
f_N = 1 / dt, y dt es el tiempo de muestreo (es decir, TR)
y el intervalo de frecuencias también es df.
Estos rangos y tamaños de paso deben ser * independientes * de la censura
que es una bonita propiedad de Lomb-Scargle-iness.

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ SALIDA:
1) PREFIX_time.1D: un archivo 1D de los puntos de tiempo muestreados (en unidades de
segundos) del analizado (y posiblemente censurado)
conjunto de datos.
2) PREFIX_freq.1D: un archivo 1D de los puntos de muestreo de frecuencia (en unidades
de 1 / segundo) del periodograma / espectro de salida
conjunto de datos.
3) PREFIX_amp + orig: conjunto de datos volumétricos que contiene un derivado de LS
o espectro de amplitud (por defecto, llamado 'amp') o un
PREFIX_pow + espectro de potencia orig (ver '-out_pow_spec', llamado 'pow')
uno por vóxel.
Tenga en cuenta que la amplitud y potencia de salida
espectros son 'unilaterales', para representar el
* mitad * amplitud o potencia de una frecuencia determinada
(ver la siguiente nota).

+ UNA NOTA SOBRE Fourier + Parseval importa (por favor, perdone los
formateo):
En la formulación utilizada aquí, para una serie de tiempo x [n] de longitud N,
el valor del periodograma S [k] está relacionado con el valor de amplitud | X [k] |:
(1) S [k] = (| X [k] |) ** 2,
para cada k-ésimo armónico.

El teorema de Parseval relaciona las fluctuaciones del tiempo con las amplitudes espectrales,
indicando que (para series en tiempo real con media cero):
(2) suma_n = (1 / N) * suma_k <| X [k] | ** 2>,
= (1 / N) * suma_k ,
donde n = 0,1. N-1 y k = 0,1. N-1 (NB: A [0] = 0, para media cero
serie). El LHS es esencialmente la varianza de la serie de tiempo
(veces N-1). Lo anterior se deriva de las matemáticas de la transformada de Fourier, y
los espectros Lomb-Scargle son aproximaciones a Fourier, por lo que lo anterior
se puede esperar que se mantenga aproximadamente, si todo va bien.

Otro resultado relacionado con Fourier es que para series de tiempo reales y discretas,
las amplitudes espectrales / valores de potencia son simétricas y periódicas en N.
Por tanto, | X [k] | = | X [-k] | = | X [N-k-1] | (en matriz de base cero
contando)
la distinción entre frecuencias indexadas positivas y negativas
puede pensarse que significa ondas que viajan a la derecha e izquierda, que
ambos contribuyen a la potencia total de una frecuencia específica.
El resultado es que se podría escribir la fórmula de Parseval como:
(3) suma_n = (2 / N) * suma_l <| X [l] | ** 2>,
= (2 / N) * suma_l ,
donde n = 0,1. N-1 y l = 0,1. (N / 2) -1 (observe el factor de 2 ahora
que aparece en las relaciones de RHS). Estas simetrías / consideraciones
son la razón por la que

Los valores de frecuencia N / 2 se emiten aquí (asumimos
que solo se ingresan series de tiempo de valor real), sin ninguna pérdida de
información.

Además, con miras a expresar la amplitud general
o potencia de una frecuencia determinada, que muchas personas pueden querer utilizar para
estimar parámetros espectrales de 'conectividad funcional' como ALFF,
fALFF, RSFA, etc. (usando, por ejemplo, 3dAmptoRSFC), por lo tanto
tenga en cuenta que la amplitud * total * o la potencia de una frecuencia dada
ser:
A [k] = 2 * | X [k] |
P [k] = 2 * S [k] = 2 * | X [k] | ** 2 = 0.5 * A [k] ** 2
en lugar de solo el de la parte de desplazamiento izquierda / derecha. Estos tipos de
las cantidades (A y P) también se denominan espectros de "dos caras". La
La relación de Parseval resultante podría entonces escribirse:
(4) suma_n = (1 / (2N)) * suma_l ,
= (1 / N) * suma_l

,
donde n = 0,1. N-1 y l = 0,1. (N / 2) -1. De alguna manera, parece más fácil
para generar los valores unilaterales, X y S, de modo que el parsevaliano
las reglas de suma se parecen más.

Con todo eso en mente, los resultados de 3dLombScargle se generan como
sigue. Para amplitudes, los siguientes aprox. Relación pasevelliana
debe mantenerse entre la serie de tiempo 'holey' x [m] de M puntos y
la serie de frecuencias Y [l] de L

M / 2 puntos (donde <| Y [l] |> se acerca
las amplitudes de Fourier <| X [l] |> como el número de puntos censurados
disminuye y M- & gtN):
(5) suma_m = (1 / L) * suma_l ,
donde m = 0,1. M-1 y l = 0,1. L-1. Para el espectro de potencia T [l]
de L

Valores M / 2, entonces:
(6) suma_m = (1 / L) * suma_l
para los mismos rangos de sumas.

Por lo tanto, tenga en cuenta eso cuando utilice las salidas de aquí. 3dAmpToRSFC
está preparado para esto al calcular los parámetros espectrales (de
amplitudes).

+ CORRER:
-prefix PREFIX: nombre del prefijo de salida para el volumen de datos, archivo 1D de punto de tiempo
y archivo 1D de frecuencia.
-Inset FILE: serie temporal de volúmenes, un conjunto de datos volumétricos 4D.

-censor_1D C1D: una sola fila o columna de 1s (mantener) y 0s (censurado)
describiendo qué volúmenes de ARCHIVO se guardan en el
muestreo y que están censurados, respectivamente. La
La longitud de la lista de números debe ser del
la misma longitud que el número de volúmenes en FILE.
Si no se ingresa, el programa buscará subladrillos
de todos ceros y supongamos que están censurados.
-censor_str CSTR: Cadena selectora de volúmenes al estilo AFNI para * mantener * adentro
el analisis. Como:
'[0..4,7,10..$]'
¿Por qué nos referimos a ella como una 'cadena de censura' cuando es
realmente la lista de volúmenes a mantener. bueno, hizo
sentido en ese momento. Los historiadores del futuro pueden batirse en duelo con
tinta al respecto.

-mask MASK: opcional, máscara de volumen para analizar adicionalmente, cualquier
vóxel con valores uniformemente cero a lo largo del tiempo
producir un espectro cero.

-out_pow_spec: cambia para generar el espectro de amplitud de las frecuencias
en lugar del periodograma. En la formulación utilizada
aquí, para una serie de tiempo de longitud N, la potencia espectral
El valor S está relacionado con el valor de amplitud X como:
S = (X) ** 2. (Sin esta opción, la salida predeterminada es
espectro de amplitud.)

-nyq_mult N2: Los periodogramas L-S pueden incluir frecuencias por encima de
normalmente se consideraría Nyquist (aquí definido
como:
f_N = 0.5 * (número de muestras) / (intervalo de tiempo total)
Por defecto, la frecuencia máxima será la
f_N * habría * sido si no hubiera censura de puntos
Ocurrió. (Esto facilita la comparación de espectros L-S
en un grupo con el mismo protocolo de escaneo, incluso si
hay ligeras diferencias en la censura, por tema).
Los valores aceptables son & gt0. (Para aquellos que lean el
documentos de algoritmo, esto establece el parámetro 'hifac').
Si no tiene una buena razón para cambiar esto,
¡No lo cambie!
-nifti: cambia al archivo de volumen de salida * .nii.gz
(el formato predeterminado es BRIK / HEAD).

+ EJEMPLO:
3dLombScargle -prefijo LSout -inset TimeSeries.nii.gz
-mask mask.nii.gz -censor_1D censor_list.txt


____________________________________________________________________________
Esta página se generó automáticamente el miércoles 23 de junio a las 15:19:40 EDT de 2021


Resultados

Búfalos africanos

El LSP de las ubicaciones de Cilla era característico de una caminata aleatoria aperiódica (sin pico), mientras que el LSP de las ubicaciones de Pepper sugería un ciclo diario (Fig. 3). Al aplicar el patrón de observaciones faltantes de Pepper a la serie temporal de ubicación de Cilla, recuperamos, como se predijo, una señal de periodicidad diaria (Fig. 3). Esto ilustró que el patrón periódico se creó por autocorrelación en la forma en que faltaban las observaciones, no el comportamiento de los animales. La PAGUn valor de 0,29 de 150 simulaciones en la prueba del modelo nulo aplicada a las ubicaciones de Pepper confirmó que el patrón periódico observado en el registro de seguimiento de Pepper era un artefacto.

Izquierda: Periodograma de los datos de seguimiento de Buffalo Cilla según el programa de muestreo original. Centrar: Periodograma de los datos de Cilla cuando se remuestrean para imitar el programa de muestreo de Pepper. Derecha: Periodograma de los datos de seguimiento de Pepper. Las flechas apuntan a los picos correspondientes al período de un día

Al aplicar el LSP al registro de velocidad de Cilla (frente al registro de ubicación anterior), encontramos, como se predijo, una periodicidad diaria significativa (PAG-valor & lt0.01 de 150 simulaciones Archivo adicional 2: Figura C5). De Cilla actividad fue así periódica a pesar de que su uso del espacio no era (Fig. 3). Esto ilustra cómo los patrones periódicos de uso del espacio se desacoplan de la periodicidad en los niveles de actividad en esta especie.

Albatros de olas

Los cuatro individuos de albatros exhibieron una periodicidad significativa (todos PAG-valores & lt0.005 de 1000 simulaciones cada una) con un período estimado de 22 a 35 días. En el individuo que realizó los viajes más largos, llegando a aguas a 1800 km de su colonia reproductora con una periodicidad estimada de 29 días, la asociación con la fase lunar fue muy pronunciada (Fig. 4). Durante tres ciclos monitoreados, el ave abandonó la colonia poco antes de la luna nueva y comenzó su viaje de regreso poco después de la luna llena, presumiblemente para que tanto ella como su pareja pudieran alimentarse durante los períodos con algo de luz de la luna.

Seguimiento de datos de un albatros ondulado. Izquierda: trazado con una escala de colores que indica la fase lunar. Los colores azules están cerca de la luna nueva y los colores rojos de la luna llena. Derecha: periodogramas. El periodograma negro es de los datos de seguimiento y el periodograma rojo es del programa de muestreo. La presencia de un pico en el periodograma negro durante un período de un ciclo lunar (línea vertical), y su ausencia en el periodograma rojo, indican que el patrón periódico no es causado por el programa de muestreo.

Lobos de crin

La periodicidad de un día fue estadísticamente significativa para los ocho individuos (todos PAG-valores & lt0.05 de 1000 simulaciones cada una (Fig. 5, panel derecho), incluso si el patrón periódico de uso del espacio no era directamente evidente en los datos de seguimiento sin procesar (Fig. 5, panel izquierdo). En otras palabras, esos ocho lobos mostraron una tendencia significativa a patrullar su área de distribución a lo largo de rutas repetidas diariamente, pero este patrón fue oscurecido por el importante ruido de fondo (los componentes estocásticos aperiódicos de los caminos).

Seguimiento de datos de lobos de crin. Izquierda: trazado para un individuo (Amadeo) según la hora del día, para mostrar la dificultad de detectar visualmente el patrón periódico debido al importante componente estocástico en el proceso de movimiento. Derecha: periodograma multiindividual de Amadeo y otros 7 lobos, que muestra la periodicidad de un día y la serie de armónicos asociados (flechas rojas)


Observabilidad de componentes espectrales más allá del límite de Nyquist en señales muestreadas no uniformemente

La identificación de un componente de señal con una frecuencia que excede el límite de Nyquist es un problema desafiante en la teoría de señales, así como en algunas áreas de aplicaciones específicas como la astronomía y las biociencias. Una consecuencia del bien conocido teorema de muestreo para una señal muestreada uniformemente es que el componente espectral por encima del límite de Nyquist se alias en un rango de frecuencia más bajo, lo que hace imposible la distinción entre componentes verdaderos y con alias. Sin embargo, el muestreo no uniforme ofrece la posibilidad de reducir los componentes con alias y descubrir información por encima del límite de Nyquist. En este artículo, proporcionamos un análisis teórico de la reducción de componentes con alias en el periodograma no paramétrico para dos esquemas de muestreo: el patrón de muestreo aleatorio y el patrón de muestreo generado por la modulación de frecuencia de pulso integral (IPFM), este último ampliamente aceptado como frecuencia cardíaca. modelo de cronometraje. Se propone una fórmula general que relaciona la varianza de las desviaciones de tiempo de un esquema regular con la supresión de componentes con alias. La relación derivada se ilustra mediante periodogramas Lomb-Scargle aplicados sobre datos simulados. Los datos experimentales presentados que consisten en la señal de respiración derivada del electrocardiograma y la señal de frecuencia cardíaca también respaldan la posibilidad de detectar frecuencias por encima del límite de Nyquist en la condición conocida como aliasing cardíaco.

1. Introducción

El aliasing es un fenómeno bien conocido en la teoría de señales que ocurre si se muestrea una señal de tiempo continuo a una frecuencia

por debajo del doble de la frecuencia máxima

de la señal muestreada, es decir,

. La frecuencia de muestreo se redujo a la mitad,

, se denomina frecuencia de Nyquist y se puede interpretar como la frecuencia máxima supuesta de la señal muestreada siempre que se haya realizado un muestreo a una velocidad suficientemente alta. Si la frecuencia de la señal muestreada excede el límite de Nyquist, esta frecuencia aparece en una posición incorrecta, en el rango

, a menos que se haya aplicado un filtro de suavizado. Sin embargo, estas reglas no son completamente aplicables en datos muestreados de manera no uniforme.

El muestreo no uniforme se puede explotar intencionalmente en el diseño de un sistema, o puede reflejar condiciones restrictivas específicas en la medición y exploración de datos en algunos campos científicos, como astronomía, sismología, genética, mediciones biofísicas y fisiología cardíaca. Los datos no se pueden adquirir en los instantes de tiempo prescritos, faltan o están dañados. En otras situaciones, los tiempos de muestreo esencialmente irregulares vienen dictados por la naturaleza, como ocurre en el análisis de la variabilidad de la frecuencia cardíaca. El muestreo no uniforme hace que el análisis de una señal sea más complicado, porque la gran mayoría de los procedimientos de análisis de datos y procesamiento de señales digitales se han desarrollado para señales muestreadas uniformemente. Por otro lado, las señales muestreadas de manera no uniforme son más informativas y se pueden detectar componentes más allá del límite de Nyquist.

La noción de la frecuencia de Nyquist necesita cierta aclaración cuando se consideran señales muestreadas de manera no uniforme. No existe una definición bien aceptada si la frecuencia de muestreo no es constante [1]. Algunos autores definen la frecuencia de Nyquist en el caso no uniforme como

es el intervalo de tiempo más pequeño entre datos, mientras que otros infieren la frecuencia de Nyquist a partir de tiempos de muestreo promedio

. Una forma más general de definir la frecuencia de Nyquist para muestras no uniformes utiliza un concepto de ventana espectral, es decir, transformada de Fourier (FT) de magnitud al cuadrado del patrón de muestreo. La ventana espectral alcanza su máximo en la frecuencia cero y un máximo secundario cercano al máximo principal en la frecuencia que se establece que es la frecuencia de Nyquist duplicada. En este trabajo, consideraremos tiempos de muestreo que se desvían de los instantes regulares de forma homogénea, es decir, sin conglomerados marcados ni brechas largas. En tal situación, tanto el tiempo de muestreo promedio como los conceptos de ventana espectral dan resultados prácticamente idénticos.

El análisis espectral de datos muestreados de manera no uniforme tiene una historia relativamente antigua debido a la importancia teórica y práctica de este problema. Se han propuesto varios enfoques para tratar con muestras no uniformes. Se puede encontrar una revisión completa del estado de la técnica en [1]. Restringiremos nuestra consideración a los métodos comúnmente utilizados en la investigación de los datos de frecuencia cardíaca que se utilizan en este artículo como datos demostrativos para ilustrar los fenómenos analizados. El método más simple para tratar con datos muestreados de manera no uniforme es ignorar la irregularidad del muestreo, reubicar las muestras en una cuadrícula regular y tratarlas como una señal muestreada uniformemente. Tal enfoque da como resultado una reducción de las amplitudes espectrales y la irregularidad del muestreo se traduce en un ruido espectral. Los componentes por encima de la frecuencia de Nyquist no se pueden distinguir de las frecuencias con alias. Los métodos basados ​​en interpolación atenúan las frecuencias más altas, haciéndolos inadecuados para el análisis de componentes por encima de la frecuencia de Nyquist. Un método ampliamente utilizado que no asume un remuestreo de señal se conoce como periodograma de Lomb-Scargle, que es básicamente un ajuste por mínimos cuadrados de datos medidos por sinusoides [2-6].

El periodograma de Lomb-Scargle (LS) [2] es un método especialmente diseñado para muestras adquiridas de forma no uniforme. Este método se propuso originalmente para el análisis de series de tiempo astronómicas y también encontró una aplicación exitosa en el análisis de datos cardiovasculares, como los de los análisis de la variabilidad de la frecuencia cardíaca [7, 8], los ritmos circadianos [9] y la secuencia genética. análisis. Los resultados de LP son en muchas situaciones casi idénticos al periodograma clásico (Schuster) generalizado a muestras no uniformes [10] y, por lo tanto, pueden inferirse de la FT de datos muestreados no uniformemente. Este enfoque se utiliza en [6] donde se presentan varias propiedades de la FT muestreada de manera no uniforme, aunque se refieren de manera imprecisa como propiedades de la densidad espectral de potencia Lomb.

A diferencia de nuestro artículo, los análisis [6] se han realizado sin una atención especial para cuantificar un efecto del aliasing, es decir, la capacidad de descubrir componentes por encima del límite de Nyquist. Este artículo estudia la utilidad del muestreo no uniforme tanto para reducir componentes con alias como para revelar componentes por encima de la tasa de Nyquist. Proporcionaremos un análisis teórico de una situación de modelo, cuando se muestrea una señal sinusoidal de acuerdo con un patrón de muestreo no uniforme preespecificado. Se propone una fórmula general que relaciona la varianza de las desviaciones de tiempo de un esquema regular con la supresión de componentes con alias. Una observabilidad teóricamente probada de componentes más allá del límite de Nyquist se ilustra mediante ejemplos de simulación y datos cardiorrespiratorios obtenidos de un experimento especialmente diseñado para este propósito.

El esquema del documento es el siguiente. La sección 2 resume las definiciones de periodogramas utilizadas en el análisis espectral de datos no uniformes. Las secciones siguientes presentan un análisis espectral teórico de la sinusoide muestreada con dos esquemas de muestreo. En la Sección 3 del artículo, derivaremos las fórmulas para las densidades espectrales de potencia (PSD) de señales sinusoidales muestreadas de manera no uniforme, siempre que los tiempos de muestreo obedezcan a un modelo aleatorio. La Sección 4 proporciona un tratamiento del espectro sinusoide cuando el patrón de muestreo fue generado por la modulación de frecuencia de pulso integral (IPFM). En la Sección 5 se deriva una fórmula cuantitativa que describe la reducción de un componente con alias. En la Sección 6 se presentan análisis demostrativos de datos simulados. rendimiento del método LS en el caso del llamado aliasing cardíaco, cuando los componentes por encima de la frecuencia de Nyquist fueron provocados por la arritmia sinusal respiratoria.

2. Periodogramas para señales muestreadas no uniformemente

El periodograma convencional basado en la transformada de Fourier de magnitud cuadrada se puede generalizar también para señales muestreadas no uniformemente en la forma [4]: ​​& # 13

son muestras de datos tomadas en los instantes de tiempo

. Observe que las raíces de un análisis tipo Fourier con exponenciales complejos residen en un ajuste de datos por mínimos cuadrados a la sinusoide compleja

. Lomb [2] formuló el periodograma mediante el ajuste por mínimos cuadrados de los datos observados de valor real al modelo de la forma: & # 13

& # 13 La variable de retardo de tiempo redundante

& # 13 se introdujo en el modelo (2) para ortogonalizar los términos seno y coseno, que surgen en una formulación de las ecuaciones normales del problema de ajuste. El PSD LS se calcula entonces como [4] & # 13

& # 13 Siempre que los términos sumados de seno y coseno al cuadrado en (4) sean iguales y se acerquen a

, las aproximaciones LS (4) al periodograma clásico (1). Generalmente, (1) y (4) no son equivalentes. Se prefiere el LS debido a las propiedades estadísticas útiles para las pruebas de significación de los picos espectrales [4, 11]. Por otro lado, (1) es matemáticamente más manejable.

3. Patrón de muestreo aleatorio

La señal analizada tomará una forma de sinusoide compleja con la amplitud

, muestreado en los instantes de tiempo

& # 13 Se supone que los instantes de muestreo son los números aleatorios que obedecen al siguiente modelo: & # 13

es el período de muestreo medio,

es una secuencia de variables aleatorias independientes distribuidas de forma idéntica. El comportamiento promedio de (1), aproximadamente válido también para el espectro LS, se puede deducir aplicando la expectativa sobre la transformada de Fourier de (5): & # 13

& # 13 Después de la sustitución (5), (6) en (7) y reordenamiento, obtenemos & # 13

& # 13 En (8), podemos reconocer FT de una sinusoide muestreada uniformemente expresable por el núcleo de Dirichlet periódico: & # 13

& # 13 La media FT (8) se puede reescribir en términos de la función característica desplazada por la frecuencia

& # 13 donde la función característica

de la irregularidad del tiempo

, caracterizado por la función de densidad de probabilidad

& # 13 Por lo tanto, la función característica atenúa los picos periódicos (alias o imágenes espectrales en (11), mientras que retiene el pico verdadero en la frecuencia

La expectativa resultante del FT (11) no tiene en cuenta completamente un efecto aleatorio de la irregularidad del muestreo. Por lo tanto, la expectativa de una FT de magnitud al cuadrado, es decir, PSD por su definición, se evalúa para tener un aunque más completo en el espectro de una señal muestreada irregularmente: & # 13

& # 13 El término de suma en (13) después de la sustitución (5) es & # 13

, tenemos debido a la presunta independencia: & # 13

& # 13 donde hemos asumido una función de probabilidad simétrica en (12) y por lo tanto una función característica de valor real.

, la expectativa (15) es 1 y se puede expandir formalmente como & # 13

en (14) y teniendo en cuenta (15), (16), la suma (14) se simplifica a & # 13

& # 13 El término de suma en (17) puede reconocerse como el FT de una sinusoide muestreada uniformemente multiplicada por la ventana triangular. Para facilitar la legibilidad de la fórmula, este FT se indicará como & # 13

& # 13 La fórmula final se convierte en & # 13

& # 13 A diferencia de (11), la expresión PSD (19), además de los picos relacionados con la sinusoide, incluye una porción de espectro suave relacionada con la aleatoriedad de los instantes de tiempo. Desde

, esta porción lisa exhibe un valle cerca de la frecuencia verdadera de la sinusoide compleja. La composición de PSD a partir de los términos presentes en (19) se aclara en la Figura 1.


4. Patrón de muestreo generado por el modelo IPFM

El modelo de modulación de frecuencia de pulso integral (IPFM) es ampliamente aceptado como el modelo de sincronización del nodo sinoauricular para modelar la variabilidad de la frecuencia cardíaca. Es un modelo de integración y disparo, también utilizado para describir una actividad neuronal. El modelo IPFM asume que se integra una señal de modulación y se genera un pulso de activación de latido cuando la señal integrada alcanza un umbral (Figura 2). Los tiempos de ocurrencia del latido (que representan los tiempos de muestreo)

generado por el modelo IPFM se puede definir como [12] & # 13

es el período cardíaco medio, y

es una señal de modulación adimensional interpretada como el cambio fraccionario de la frecuencia cardíaca instantánea en relación con


El patrón de muestreo se puede escribir como una secuencia de pulsos de Dirac: & # 13

& # 13 La transformada de Fourier de (21) se conoce como el espectro de cuentas (SPC) [12]: & # 13

& # 13 El espectro de una señal muestreada

está relacionado con el espectro real

por la convolución en el dominio de la frecuencia: & # 13

& # 13 El espectro de recuentos se ha analizado extensamente en varios trabajos [12-14]. Resumiremos aquí los hechos necesarios para evaluar el efecto del aliasing espectral. El SPC definido por (22) se puede reescribir en términos de la señal moduladora como & # 13

& # 13 Curiosamente, el espectro de la señal moduladora

se incluye directamente en (24). Además del pulso de Dirac ubicado en el origen de la frecuencia, un número infinito de términos de frecuencia modulada (FM) convolucionados con

están presentes. Las portadoras de estos términos modulados en FM que producen picos en (24) están ubicadas en múltiplos de la frecuencia de muestreo media, lo que explica el origen de los componentes con alias en el espectro (23). Para la señal de modulación sinusoidal: & # 13

es la frecuencia de modulación, el espectro (24) se puede ampliar utilizando las funciones de Bessel del primer tipo [15]. Las amplitudes espectrales de las portadoras siguen la función de Bessel de orden 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. The choice

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.

Referencias

  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, no. 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,” Astrophysical Journal, 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,” Avisos mensuales 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

Copyright

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.


Ver el vídeo: Practico5 parte2 - Biplot componentes principales y periodograma (Diciembre 2022).