AplicacindelmtododeHilbert-Huangdescripcinymetodologa

Vista previa del material en texto

1 
 
Aplicación del método de Hilbert-Huang a señales biológicas en el 
campo de la neurología: descripción y aspectos metodológicos. 
Mario Estévez Báez1, Calixto Machado Curbelo2, Eduardo Arrufat-Pié3, Aisel Santos Santos4 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
La Habana, Cuba Octubre de 2017 
 
 
Dirección para la correspondencia: Instituto de Neurología y Neurocirugía, 
Departmento de Neurofisiología Clínica, 29 y D, Vedado, La Habana 10400, Cuba. 
Tel./fax: +837-834 5578 
 
e-mail address: marioestevez@infomed.sld.cu (M. Estévez Báez) 
 
 
 
1 Doctor en Medicina, Especialista de Fisiología de Primer y Segundo Grados, Profesor e Investigador Titular, Doctor en Ciencias Médicas, 
Académico Titular de la IAA, Instituto de Neurología y Neurocirugía. 
 
2 Doctor en Medicina, Especialista de Fisiología Normal y Patológica de Primer Grado, Especialista de Neurología de Segundo Grado, Profesor e 
Investigador Titular, Doctor en Ciencias, Investigador de Mérito, Académico Titular de la ACC, Instituto de Neurología y Neurocirugía. 
 
3 Doctor en Medicina, Residente de Fisiología de Tercer Año, Instituto Superior de Ciencias Médicas de La Habana, “Victoria de Girón”. 
 
4 Doctor en Medicina, Especialista de Primer Grado en MGI y Neurología, Master en Ciencias, Profesor Auxiliar, Instituto de Neurología y 
Neurocirugía. 
 
Resumen: 
El análisis espectral constituye un método extensamente utilizado para el procesamiento 
digital de diferentes señales de interés como el EEG, el ECG y la variabilidad de la 
frecuencia cardiaca. Los procedimientos más utilizados están asociados a los métodos de 
Fourier, que estrictamente hablando, requieren que los procesos estudiados sean 
gaussianos, lineales y estacionarios, lo que realmente no se cumple para estas señales. El 
método de Hilbert-Huang (HH) constituye un novedoso abordaje que brinda una solución 
metodológica a las anteriores limitaciones y que en los últimos años ha comenzado a ser 
introducido en el campo de las neurociencias. En el trabajo se describen los métodos 
alternativos que más han sido empleados para afrontar la falta de linealidad y 
estacionariedad de las señales, señalando sus posibilidades y limitaciones, y luego se 
describe en detalles el método de HH, abordándose los aspectos teóricos y de carácter 
metodológico que se han de tomar en cuenta para el uso adecuado del mismo para la 
evaluación del contenido espectral de señales biológicas de tan alto significado clínico-
fisiológico para las neurociencias. 
Palabras clave: 
Empirical mode decomposition 
Hilbert-Huang transform 
Spectral analysis 
EEG 
Nonlinear methods 
mailto:marioestevez@infomed.sld.cu
2 
 
 Introducción 
El procesamiento digital de señales biológicas generadas 
por diferentes estructuras del sistema nervioso del hombre y los 
animales, ha hecho posible en gran medida el desarrollo alcanzado por 
las neurociencias. Múltiples escollos han tenido que ser debidamente 
evitados y las posibles limitaciones en el alcance de los resultados 
alcanzados, podrían en cierta medida estar vinculados a los métodos de 
análisis utilizados y especialmente a sus posibilidades y limitaciones. 
Probablemente, uno de los problemas que ha incidido con 
mayor relevancia, esté dado por la naturaleza de los procesos 
biológicos que son responsables de la información que portan muchas 
de estas señales biológicas. Múltiples evidencias existen que 
demuestran que la mayoría de los procesos biológicos generados por el 
sistema nervioso poseen gran complejidad. Por solo mencionar algunas 
de ellas: se ha reportado que el electroencefalograma (EEG) muestra 
en su estructura la existencia de componentes que reflejan procesos no 
lineales (Pritchard et al. , 2000, Li et al. , 2009, Pradhan et al. , 2012, 
Wang et al. , 2015), ha sido cuestionada la estructura gaussiana de la 
misma (McEwen et al. , 1975, Persson, 1977, Bender et al. , 1992, 
Kipinski et al. , 2011, Wang et al. , 2015), así como que esta señal 
contiene elementos no estacionarios (McEwen et al. , 1975, Persson, 
1977, Pritchard et al. , 2000, Guarín et al. , 2011, Kipinski et al. , 2011). 
De igual manera, son variadas las evidencias de no 
linealidad (Acharya et al. , 2004, Faes et al. , 2009, Chen et al. , 2010, 
Silva et al. , 2017) y no estacionariedad (Weber et al. , 1992, Faes et 
al. , 2009, Magagnin et al. , 2011) en la variabilidad de la frecuencia 
cardiaca (VFC), una señal que se deriva del electrocardiograma (ECG) 
y que está asociada con la actividad de la regulación cronotrópica por 
parte del sistema nervioso autónomo al sistema cardiovascular y 
específicamente al corazón. 
El estudio de las señales biológicas a las que antes nos 
hemos referido, ha logrado sus resultados más significativos en los 
estudios realizados en el dominio de la frecuencia, donde sin dudas, el 
método del análisis de Fourier ha resultado la piedra angular. No 
obstante, como han enfatizado Norden Huang y sus colaboradores 
(Huang et al. , 1998) este método “…posee cruciales restricciones: el 
sistema en estudio tiene que ser lineal y los datos que se someten al 
análisis tienen también que ser estrictamente periódicos o 
estacionarios; de no ser así, los espectros resultantes tendrían poco 
sentido físico”. 
Para evitar estas restricciones, se han tenido que crear 
diferentes convenciones, tales como la definición de la duración de 
tiempo de la señal que puede ser considerada al menos como casi 
estacionaria, la estricta observancia de determinadas condiciones 
estandarizadas para las sesiones de registro, la comprobación de la 
estacionariedad de las señales registradas mediante pruebas específicas 
antes de poder aplicar el método de Fourier (Weber et al. , 1992, Guarín 
et al. , 2011, Kipinski et al. , 2011), e incluso la utilización de métodos 
alternativos como el uso de espectros de estadígrafos de orden superior 
(en lengua inglesa “higher order statistics/spectra”) (Chua et al. , 2010, 
Shahid et al. , 2011, Pradhan et al. , 2012). 
1.2 Breve descripción de algunos métodos alternativos. 
1.2.1 El espectrograma 
Varios métodos han sido estudiados y aplicados para 
enfrentar la presencia de no estacionariedad de estas señales. Uno de 
los primeros intentos fue el método del espectrograma, o transformada 
de Fourier de corta duración (short-time Fourier transform), o también 
llamado método de Gabor. Del total de muestras de una señal que va a 
ser analizada, se selecciona una ventana con una duración obviamente 
inferior a la duración total del registro y se calcula para la misma su 
espectro usando el algoritmo que se desee ― generalmente se utiliza el 
periodograma clásico modificado por una ventana como la de 
Hamming. Posteriormente, se desplaza el inicio de la ventana por un 
corto espacio de tiempo de las muestras originales del registro y se 
calcula nuevamente su espectro. Sucesivamente se van obteniendo así 
espectros que están separados unos de otros por un mismo periodo de 
tiempo: el desplazamiento que se ha seleccionado emplear. Con ellos 
se crea un diagrama tridimensional donde se pueden observar las 
fluctuaciones espectrales (frecuencia y amplitud o potencia 
espectrales) contra el tiempo. 
 
 
Figura 1. Espectrograma de las fluctuaciones espectrales de la variabilidad de la frecuencia cardiaca evaluada en una paciente durante los treinta 
minutos ulteriores a haber sido declarada oficialmente en estado de muerte encefálica. La figura ilustra el empleo del método del espectrograma, o 
método de Gabor, para poner en evidencia los cambios de las frecuencias espectrales en función del tiempo. Este método fue el primero de los 
llamados métodos espectrales de tiempo-frecuencia. Este resultado fue objeto de publicación en Machado, C., Estevez, M., Perez-Nellar, J., & 
Schiavi, A. (2015). Residual vasomotor activity assessed by heartrate variability in a brain-dead case. BMJ Case Rep, 2015. 
3 
 
 
Figura 3. Ilustración del empleo del espectrograma para mostrar el efecto sobre el sistema nervioso autónomo del Zolpidem en siete pacientes en 
estado vegetativo persistente. Las flechas en cada diagrama indican el momento de administración del fármaco. En todos los casos puede apreciarse 
el efecto vagolítico del medicamento. Estos resultados fueron objeto de publicación en el trabajo Machado, C., Estévez, M., Rodríguez, R., Pérez 
Nellar, J., Chinchilla, M., Defina, P., et al. (2014). Zolpidem Arousing Effect in Persistent Vegetative State Patients: Autonomic, EEG and 
Behavioral Assessment. Current Pharmaceutical Design, 20(26), 4185-4202. 
 
En las Figuras 1, 2 y 3 se muestran algunos de estos 
diagramas que han sido parte de trabajos publicados por nuestro 
colectivo de autores (Estévez et al. , 2012, Machado et al. , 2014, 
Machado et al. , 2015). Los resultados que son capaces de ser puestos 
en evidencia mediante el uso del espectrograma dependen en gran 
medida de los valores seleccionados por el usuario para el proceso, 
quien debe decidir la duración de la ventana y con ello dejar fijadas la 
resolución espectral y las posibles bandas de frecuencia a estudiar, el 
tiempo o número de muestras a usar para desplazar la ventana, lo que 
determinará la posible presencia o no de eventos en el tiempo, el 
algoritmo para el cálculo del espectro (paramétrico, no paramétrico, u 
otros más modernos), y otros más. En fin, un usuario con poco 
conocimiento técnico del método puede dejar de advertir cambios que 
sí existieron durante el registro, o mostrar resultados que no tienen 
sentido para la clínica o la fisiología. 
1.2.2 Uso de distribuciones tiempo-frecuencia. 
Con posterioridad al anterior método se comenzaron a 
utilizar otros procedimientos de tiempo-frecuencia que están basados 
en otras distribuciones tales como la de Wigner-Ville, la de Kirkwood-
Rihaczek, la de Levin, la exponencial o de Choi-Williams, la Sinc o de 
Born-Jordan y otras. Muchas de ellas son llamadas distribuciones de 
Cohen. En la práctica, y en opinión nuestra, en el campo de señales en 
neurociencias no han aportado nada nuevo con respecto al 
espectrograma, por lo cual no profundizaremos en su descripción. Para 
los interesados hay una monografía al respecto que incluye amplia 
información (Boashash, 2003). 
1.2.3 Pequeñas ondulaciones (“wavelets”). 
Otro método de mayores potencialidades dentro del grupo 
de los llamados de tiempo-frecuencia es el de las pequeñas 
ondulaciones u ondículas (“wavelets” en lengua inglesa). En nuestra 
lengua es de uso generalizado llamar al procedimiento con su 
denominación original de wavelets y así lo haremos en este documento. 
El gran inconveniente del método de Fourier para el análisis 
de frecuencias está dado por el hecho de que al producirse la 
transformación de los datos originales (que tienen como dimensiones 
al tiempo y los valores de amplitud o intensidad de la señal) al dominio 
de la frecuencia, se pierde la dimensión tiempo, quedando entonces 
solo la dimensión amplitud-intensidad, que aparece como valores de 
potencia y las frecuencias espectrales discretas a que las mismas 
corresponden. Usando las wavelets es posible lograr una 
descomposición de la señal original sin perder la información de la 
dimensión tiempo. 
Una wavelet no es más que una onda de duración limitada 
y cuyo valor promedio de voltaje es cero. La transformada continua de 
las wavelets puede definirse como la suma para todo el tiempo que dura 
la señal de sus valores multiplicados por las versiones escaladas y 
4 
 
desplazadas de la wavelet seleccionada. Este proceso crea coeficientes 
que son función de escala y de posición. 
Hay que tener en cuenta que los resultados obtenidos van a 
depender entre otras muchas cosas de la wavelet escogida y es 
necesario decir que existen numerosísimas familias de wavelets que 
llevan el nombre de sus creadores (Daubechies, Coiflet, Meyer, etc.) o 
de la forma de las propias wavelets (por ejemplo, sombrero mejicano). 
Aunque a primera vista el método de las wavelets puede 
mostrarse sencillo e intuitivo, en el trasfondo subyace un complejo 
proceso matemático, que exige del investigador una preparación 
avanzada. Para algunas de las wavelets, como las de Morlet, su 
fundamento está directamente asociado al método de Fourier, por lo 
cual la interpretación física adecuada solo resulta posible si el 
fenómeno en estudio es estrictamente de naturaleza lineal. 
El software Matlab posee una serie de herramientas muy 
útiles que permiten con cierta rapidez al investigador introducirse en el 
uso de estos procedimientos, que debemos decir sin embargo, que no 
han sido muy usados por investigadores de las neurociencias en nuestro 
país. Nuestro grupo de trabajo solo las ha empleado ocasionalmente, 
aunque debemos reconocer que con excelentes resultados. 
De todos los métodos alternativos a los que nos referimos 
en este trabajo, es el que reúne las mejores condiciones para enfrentar 
la no estacionariedad y no linealidad de los procesos en las señales que 
más estudiamos (EEG, ECG, VFC), pero también es el más complejo 
de aplicar y de interpretar. 
1.2.3 Otros métodos. 
Han sido utilizados otros métodos tales como la estimación 
de tendencias por medio de mínimos cuadrados, el empleo de ventanas 
de promedios deslizantes y el uso de la diferenciación para tratar de 
obtener datos estacionarios. Estos métodos, no obstante, aunque útiles, 
resultan muy especializados y no han sido extensamente utilizados 
(Huang et al. , 1998). 
1.3 Método propuesto por Norden Huang y colaboradores. 
En 1998 fue publicada la descripción de un nuevo método 
concebido para poder ser utilizado en el análisis digital de señales en 
cuya generación biológica estuviesen implicados procesos no lineales 
y no estacionarios (Huang et al. , 1998). La publicación se realizó en 
una prestigiosa revista “Proceedings of the Royal Society A” de la 
Sociedad Real de Ciencias Matemáticas, Físicas y de Ingeniería de 
Gran Bretaña y los autores y coautores pertenecían a instituciones de 
alto prestigio, tales como el Centro de Vuelos Espaciales Goddard de 
la NASA, la Universidad John Hopkins, el Laboratorio de Procesos de 
la Hidrosfera de la NASA, el Laboratorio para Estudios Navales de los 
EEUU y el Centro de Guerra Naval de los EEUU, entre otros. Este 
trabajo describió un método, ya evidentemente desarrollado y aplicado 
en campos muy alejados al de la Biología, pero por su alcance generó 
“un antes y un después” (criterio personal nuestro) en el universo del 
procesamiento digital de señales. 
 Nuestro grupo de trabajo comenzó a familiarizarse con esta 
nueva manera de concebir el procesamiento de las señales con las que 
hemos trabajado en el transcurso de nuestra vida como investigadores 
de las neurociencias, a partir del año 2014 y presentó sus experiencias 
iniciales en el año 2015, en un evento internacional auspiciado por 
nuestro Instituto de Neurología y Neurocirugía: el VII Simposio 
Internacional sobre muerte cerebral y trastornos de conciencia, en su 
conferencia inaugural, la cual fue titulada “Incorporando nuevos 
métodos del procesamiento digital de señales a los métodos 
electrofisiológicos tradicionales en pacientes con trastornos de 
conciencia y muerte encefálica” (Estévez et al. , 2015). 
A partir de aquel momento se continuó profundizando en el 
desarrollo de diferentes materiales de software para el procesamiento 
de distintas señales específicas y acumulando resultados que ya han 
servido para la preparación y envío para publicación de diferentes 
trabajos a revistas internacionales (Estévez-Báez et al. , 2017, Machado 
et al. , 2017), así como para fundamentar el desarrollo de nuevos 
protocolos de investigación donde se puedan introducir los mismos 
para beneficiode nuestros pacientes y del desarrollo del conocimiento 
científico en el campo específico que aborda nuestra institución. 
2. Descomposición empírica modal 
2.1 Generalidades. 
La descomposición empírica modal es la piedra angular del 
método propuesto por Huang y sus colaboradores (Huang et al. , 1998). 
Para la aplicación del mismo se hace necesario que la señal que se 
analiza cumpla con tres supuestos: 
1. Que posea al menos dos extremos, uno máximo y otro mínimo. 
2. Que la escala de tiempo característica esté definida por el lapso 
de tiempo entre extremos. 
3. Que si los datos careciesen de extremos y contuviesen solamente 
puntos de inflexión, pudiesen ser sometidos al proceso de 
diferenciación una o más veces para poner en evidencia los 
extremos. 
Estos supuestos son perfectamente cumplidos por las 
señales biológicas que hemos explorado en nuestros estudios y que son 
el electroencefalograma, el electrocardiograma, el tacograma de las 
secuencias consecutivas de intervalos cardiacos R-R, el 
electronistagmograma y la fotopletismografía. 
La descomposición de la señal en componentes que reciben 
el nombre de funciones intrínsecas modales (en lengua inglesa intrinsic 
 
 
 
mode functions) se efectúa a través de un procedimiento que los 
creadores del método bautizaron como tamizado (sifting en lengua 
inglesa). 
Si denominamos como X(t) a nuestra señal, los pasos a 
seguir para la operación de tamizado implican los siguientes pasos: 
 Identificar los extremos en la señal X(t), tanto los superiores o 
máximos e inferiores o mínimos. 
 Obtener mediante interpolación de los extremos máximos y 
mínimos identificados (usando el método de las “cubic splines”) 
la envoltura superior e inferior de la señal, que podemos expresar 
como E(Mx) y E(Mn). 
 Calcular la media de los valores de la envoltura superior e inferior 
       2/MnEMxEtM 
 (2.1) 
 Extraer mediante sustracción el posible candidato a componente 
modal, D(t) 
     tMtXtD 
 (2.2) 
5 
 
 Si no se cumplen las condiciones para considerar al valor extraído 
como un componente modal, continuar iterando, aplicando los 
pasos anteriores, pero considerando en cada caso como señal, al 
valor de la señal original menos los valores sucesivamente 
extraídos. 
 Cuando se ha logrado extraer un componente modal que reúna las 
condiciones exigidas se efectúa la substracción del valor de ese 
componente de la señal original y se inicia de nuevo el proceso de 
tamizado para encontrar un próximo componente modal y así 
sucesivamente. 
 Al llegarse al final del proceso, se habrán extraído diferentes 
componentes modales que podríamos denominar como C1, C2,…, 
Cn y quedará un residuo, que podemos denominar como Rs. 
En la Figura 4 hemos ilustrado el proceso de tamizado para 
mejor comprensión. Hay aspectos que no han quedado debidamente 
definidos en el algoritmo del proceso de tamizado (sifting) que hemos 
anteriormente detallado y pasaremos a precisarlos en el próximo 
acápite. 
2.2 Componentes intrínsecos modales (IMFs). 
Un componente modal, o utilizando la terminología original 
de sus creadores, una función intrínseca modal, que llamaremos 
siguiendo las siglas inglesas como IMF, tiene que reunir dos 
condiciones para poder ser así considerado: 
1. Para todo el conjunto de datos que componen al componente 
modal la suma de los extremos (superiores e inferiores) tiene 
que ser igual al número de cruces por cero observados en la 
curva, o diferir cuando más en 1. 
2. En todos los puntos el valor medio de las envolventes 
calculadas a partir de los máximos y los mínimos tiene que 
ser cero. 
 
La primera condición resulta similar a los requerimientos 
tradicionales que debe poseer un proceso gaussiano estacionario de 
banda estrecha. La segunda condición fue una nueva idea de Huang et 
al. (Huang et al. , 1998) “… el mismo modifica el requerimiento clásico 
global a uno local; es necesario que la frecuencia instantánea no tenga 
fluctuaciones indeseadas inducidas por formas de onda asimétricas”. 
Como también expresaron literalmente estos autores, el 
nuevo método, a diferencia de casi todos los anteriores utilizados para 
analizar fenómenos no lineales y no estacionarios, resulta “intuitivo, 
directo, a posteriori y adaptable (adaptive en lengua inglesa) con la 
base de la descomposición basada en y derivada de los datos”. 
 
 
Figura 4. Ilustración del proceso de tamizado (sifting) que emplea el método de la descomposición empírica modal. En A se muestra la señal original 
correspondiente a 2048 muestras del tacograma de un sujeto sano, correspondientes a 5 minutos de registro del electrocardiograma. Al tacograma 
original se le ha eliminado la media, como parte del pre-procesamiento de esas señales, por lo cual se pueden observar valores positivos y negativos. 
En B se muestra la curva que contiene los valores de la envolvente o perfil (envelope en lengua inglesa) de la curva, calculada mediante el método 
de interpolación de los valores máximos de la señal (línea de color rojo), así como la envoltura o perfil calculado para los mínimos de la señal (línea 
de color azul). En C se muestra con color verde el resultado de promediar punto a punto los valores correspondientes a las envolventes de los 
valores máximos y de los mínimos en la envoltura o perfil de la señal. En D se muestra el resultado de la operación de restar los valores de la curva 
obtenida en C, de la señal original mostrada en A. 
 
El proceso de tamizado persigue dos grandes objetivos: la 
eliminación de ondas superpuestas (riding waves) y lograr que los 
perfiles de las ondas sean mucho más simétricos. Si utilizamos la 
nomenclatura que utilizaron Huang et al en su trabajo, la media de las 
envolventes (envelopes) superior e inferior se puede denominar m1 y la 
diferencia entre los datos originales y m1 sería el primer componente 
del tamizado que podría denominarse como h1, o sea, 
  111 hmtX  (2.3) 
6 
 
Al proseguir el tamizado se utiliza a h1 como los datos, 
obteniéndose que 
11111 hmh  (2.4) 
El proceso se repite k veces hasta que h1k cumpla los 
requisitos de una función intrínseca modal, o sea, 
  kkk hmh 1111  (2.5) 
Podemos entonces designar a este componente modal como 
c1 
khc 11  (2.6) 
Se puede separar entonces a c1 del resto de los datos 
  11 rctX  (2.7) 
Como este residuo r1 contiene aún información de otros 
componentes modales se le considera como los nuevos datos a procesar 
y se somete al tamizado. Al repetirse este proceso un número dado de 
veces el resultado será 
nnn rcrrcr  1221 ,..., (2.8) 
El proceso de tamizado se detiene cuando concurren dos 
posibles situaciones: el componente cn o el residuo rn son muy 
pequeños o que el residuo rn se convierte en una función monotónica 
de la cual no se puede ya extraer un nuevo componente modal. 
Tomando en cuenta las expresiones 2.7 y 2.8 se arriba a una importante 
expresión 
  


n
i
ni rctX
1 (2.9) 
La misma nos informa de la cualidad que tiene el proceso 
de reconstruir la señal original a partir de la sumatoria de los 
componentes modales extraídos con el residuo final. 
Otro aspecto que no habíamos dejado aclarado tampoco 
hasta ahora es el criterio del momento en el cual podemos detener el 
proceso de tamizado y considerar que se ha detectado un componente 
modal. Este criterio debe garantizar que los componentes modales 
conserven suficiente sentido físico tanto en las modulaciones 
observadas de amplitud y frecuencia. El criterio más usado y que fue 
originalmente descrito por Huang et al fue limitando el valor de la 
desviación estándar calculada entre dos tamizados sucesivos mediante 
la siguiente expresión, a valores entre 0.2 a 0.3. 
 

 








 

T
t k
kk
th
thth
DS
0
2
)1(1
2
1)1(1
)(
)()(( 2.10 
 
 
Figura 5. Ilustracióndel proceso de extracción de componentes modales, utilizando algoritmos propios, de una secuencia de intervalos R-R del 
electrocardiograma en un sujeto sano. Solo se muestran los primeros 5 componentes modales, llamados IMFs en la lámina, y no se incluyó tampoco 
por supuesto el residuo final del proceso. 
 
En un periodo inicial desarrollamos diferentes versiones 
algorítmicas en Matlab para efectuar el proceso de tamizado exigido 
para la descomposición empírica modal. Probamos diversos abordajes 
con mayor o menor éxito, pero logrando sobre todo con ello 
profundizar en las sutilezas de los problemas de implementación. 
Usamos incluso otros métodos para el proceso de interpolación con el 
fin de calcular las envolventes de las señales analizadas y entre ellos el 
uso del método cúbico de Hermite, hasta que pudimos comprender la 
conveniencia y diríamos mejor exigencia del empleo de cubic splines 
para el proceso de interpolación como habían ya señalado Huang et al. 
En las Figuras 5, 6 y 7 se muestran algunas ilustraciones 
acerca del proceso de extracción de componentes modales utilizados 
en una etapa inicial por nosotros utilizando algoritmos propios. En la 
Figura 6 se puede advertir la relación de la amplitud de los 
componentes extraídos respecto a los valores de amplitud de la señal 
original. En la Figura 7 se muestran en superposición todos los 5 
componentes modales extraídos de la señal original lo que 
complementa mejor lo que se deseaba mostrar en la figura anterior. 
2.3 Algoritmos de cálculo de la descomposición empírica modal. 
7 
 
Para el año 2004 ya se había demostrado que el proceso de 
descomposición modal se comportaba como un banco diádico de filtros 
(Flandrin et al. , 2004) y según nuestro criterio, un algoritmo en 
especial, desarrollado por D. Rilling y Holger Nahrstaedt en lenguaje 
Matlab, se convirtió en la base de casi todos los que luego se han 
utilizado con el fin de mejorar algunas características indeseables de 
los componentes modales extraídos. Estos autores desarrollaron una 
colección de software (toolbox) completa en Matlab, a la que 
denominaron “Empirical Mode Decomposition Toolbox, Version 1.3”, 
conteniendo múltiples herramientas de software y la hicieron 
disponible para su uso libre en un sitio web http://perso.ens-
lyon.fr/patrick.flandrin/emd.html. Nosotros hemos utilizado para 
nuestro trabajo actual muchos de los conceptos desarrollados por ese 
grupo de trabajo francés. 
Los problemas que se fueron detectando con la aplicación 
de diferentes algoritmos de la DEM podrían resumirse a dos: 
1. Presencia de oscilaciones de diferentes amplitudes en un 
mismo componente modal. 
2. Presencia de oscilaciones muy similares en diferentes 
componentes modales, a lo que en lengua inglesa se le dio en 
llamar mode mixing. 
Aunque no lo incluimos en la categoría de problemas 
metodológicos de la DEM, también deberíamos incluir a la eficiencia 
en tiempo de los diferentes algoritmos. Hay que tener en cuenta que los 
procesos de iteración para la extracción de un componente modal 
(IMF), pueden llegar a 3000 o 5000 y aunque las máquinas modernas 
son muy rápidas, la extracción de 10 o 12 componentes modales puede 
llevar en ocasiones varios minutos de ejecución. 
Por tanto, pronto comenzaron a surgir diferentes 
implementaciones para calcular la DEM y que han sido publicadas 
(Donghoh Kim et al. , 2009, Wu et al. , 2009, Torres et al. , 2011, 
Pustelnik et al. , 2012, Colominas et al. , 2014, Wang et al. , 2014, 
Colominas et al. , 2015, Colominas et al. , 2016). Algunas han sido 
ubicadas en la Internet para uso por quienes lo deseen. A continuación 
citaremos algunos de estos sitios: 
http://www-ljk.imag.fr/members/Thomas.Oberlin/EMDOS.tar.gz 
http://perso.ens-lyon.fr/nelly.pustelnik/Software/Prox-EMD_v1.0 
http://perso.ens-lyon.fr/patrick.flandrin/emd.html 
http://www.mit.edu/~gai/CODE/HRV/emd.m 
Uno de los procedimientos introducidos en el proceso de la 
DEM fue la adición de ruido blanco a la señal original en el proceso de 
tamizado con el marcado propósito de reducir el fenómeno de mode 
mixing, apareciendo algoritmos que han sido denominados en lengua 
inglesa “ensemble empirical mode decomposition” o EEMD por sus 
siglas. No obstante, tales algoritmos también han traído como 
consecuencia que la señal reconstruida por la sumatoria de los 
componentes modales y el residuo final, como describe la expresión 
2.9, contenga además ruido residual. 
 
 
 
 
 
Figura 6. Superposición sobre la señal original (color azul) de los componentes modales extraídos (color rojo) empleando algoritmos de cálculo 
desarrollados por nuestro grupo. 
 
http://perso.ens-lyon.fr/patrick.flandrin/emd.html
http://perso.ens-lyon.fr/patrick.flandrin/emd.html
http://www-ljk.imag.fr/members/Thomas.Oberlin/EMDOS.tar.gz
http://perso.ens-lyon.fr/nelly.pustelnik/Software/Prox-EMD_v1.0
http://perso.ens-lyon.fr/patrick.flandrin/emd.html
http://www.mit.edu/~gai/CODE/HRV/emd.m
8 
 
 
Figura 7. Superposición de los primeros 5 componentes modales extraídos de la señal original que se muestra en el diagrama superior de la figura 
empleando algoritmos de cálculo desarrollados por nuestro grupo. 
 
Como precisó M. Colominas (Colominas et al. , 2014) las 
diferentes realizaciones de la señal más ruido pueden producir un 
número diferente de componentes modales. Una solución orientada a 
solucionar este conflicto fue desarrollada por un algoritmo desarrollado 
por el grupo de la Universidad Nacional de Entre Ríos, de Argentina, 
a la que pertenece Colominas y que apareció publicada en el sumario 
de una conferencia científica en el año 2011 (Torres et al. , 2011). 
Llamaron a este algoritmo como “complete ensemble empirical mode 
decomposition with adaptive noise” y le dieron las siglas de 
CEEMDAN. 
Tres años más tarde se mejoró sensiblemente este 
procedimiento por el mismo grupo, ya que aún quedaba cierto grado de 
ruido residual que fue eliminado, haciendo desaparecer además 
algunos componentes modales espurios que se extraían en las etapas 
iniciales del tamizado. A este algoritmo optimizado sus autores le 
denominaron como “improved complete ensemble empirical mode 
decomposition with adaptive noise” con las siglas iCEEMDAN 
(Colominas et al. , 2014). 
Nuestro grupo de trabajo ha mantenido cordiales relaciones 
con este avanzado colectivo de investigadores dirigido por el profesor 
Gastón Shlotthauer, y muy en particular con el amigo Marcelo 
Colominas. Gran parte de los resultados que hemos ido acumulando 
han sido utilizando precisamente el iCEEMDAN como procedimiento 
para la extracción de componentes modales. 
Hay que decir, que el Dr. M. Colominas ha desarrollado un 
nuevo procedimiento para la extracción de componentes modales, que 
ha puesto a nuestra disposición, al que ha denominado “unconstrained 
optimization approach empirical mode decomposition” con las siglas 
UOA-EMD (Colominas et al. , 2015, Colominas et al. , 2016). En este 
trabajo no vamos a incluir los resultados preliminares que hemos ido 
obteniendo con ese software para la descomposición de señales, pero 
podemos afirmar que son muy prometedores y que el algoritmo resulta 
sensiblemente más rápido que cualquiera de las otros procedimientos 
de extracción de componentes modales excluyendo tal vez al método 
de D. Rilling y Holger Nahrstaedt ya antes mencionado y citado. 
Con esta información podemos ya pasar al segundo paso del 
método que propusieron originalmente Huang et al. (Huang et al. , 
1998) y que está asociado con el uso de la transformada de Hilbert a 
los componentes modales ya extraídos mediante el proceso de la DEM.
 
3 Utilización de la transformada de Hilbert 
3.1 Sobre la transformada de Hilbert. 
La transformada de Hilbert lleva el nombre de quien la creó, 
David Hilbert (1862-1943) y la misma fue introducida inicialmente 
para dar solución a un caso especialdel problema planteado por 
Riemann y el propio Hilbert para funciones holomórficas asociadas con 
funciones analíticas. La misma puede ser considerada como la 
convolución de una función de tiempo X(t) con la función 
 
t
th

1 . Como la función h(t) no resulta posible de integrar, las 
integrales que definen a la convolución no convergen. Se define por 
ello a la transformada de Hilbert utilizando el valor principal de Cauchy 
 
 
 
 
 
 
 
que convencionalmente denominaremos v.p. y entonces la 
transformada podría expresarse así: 
   
 
  





 



 d
t
X
pvdthtXpvtXH ..
1
..))((
 3.1 
Para el desarrollo del método (Huang et al. , 1998) resultaba 
importante la obtención de frecuencias instantáneas, un concepto 
controversial, pero que puede ser factible a partir del cálculo de las 
mismas a partir de una función analítica, las cuales pueden crearse 
precisamente utilizando la transformada de Hilbert. Para una serie de 
tiempo arbitraria X(t) se puede obtener su transformada de Hilbert, que 
podríamos llamar Y(t) y tendría la expresión 
9 
 
 
 




d
t
X
pvtY 




)(
..
1
 3.2 
X(t) y Y(t) conforman un par conjugado completo de una 
señal analítica Z(t) que podemos expresar de la siguiente manera 
      ;tiYtXtZ 
 3.3 
A partir de la señal analítica podemos obtener importantes 
valores instantáneos de la amplitud o potencia espectral, la frecuencia 
instantánea, así como la fase instantánea como se define en las 
expresiones siguientes: 
          
 
 
 
 
;;arctan;
22
td
td
t
tX
tY
ttYtXtP

 






 3.4 
donde P(t) serán las potencias instantáneas, θ(t) las frecuencias 
instantáneas y ω(t) las fases instantáneas. 
Como los componentes modales extraídos por la DEM 
poseen la propiedad de ser monotónicos o monocomponentes, se 
justifica el uso de la transformada de Hilbert para obtener valores 
instantáneos como los anteriormente mencionados. Podemos advertir 
entonces, la importancia que posee desde el punto de vista teórico la 
descomposición empírica modal. Es un paso básico para poder aplicar 
con rigurosidad la transformada de Hilbert a los componentes extraídos 
y así poder obtener los valores instantáneos necesarios para analizar las 
señales mediante herramientas tales como el espectro de Hilbert y el 
espectro marginal de Hilbert que veremos más adelante. 
El análisis espectral de Hilbert (AEH), adicionado al 
método de la descomposición empírica modal ha permitido la 
obtención de variaciones temporales de tres importantes elementos del 
dominio de la frecuencia en una señal: su amplitud o potencia, su 
frecuencia y su fase. Como recientemente han planteado Huang et al. 
(Huang et al. , 2016) el AEH introdujo dos ideas novedosas. En primer 
lugar, la expansión aditiva del proceso está fundamentalmente basada 
en una función modal intrínseca adaptable (recordemos que el término 
“adaptable” se debe entender en el sentido de que los valores surgen 
basados en y derivados de los datos como tal) obtenida a través de la 
DEM y se puede representar de la siguiente manera 
         
 
  
  


N
j
N
j
N
j
di
jjjj
t
j
etattatctX
1 1 1
Recos


 3.5 
donde la frecuencia se define como la derivada en el tiempo 
de la función de fase calculada “adaptablemente” θj(t). El resultado ya 
no resulta ser un valor medio del espacio de tiempo determinado por 
integración, sino que ahora posee valores instantáneos en diferentes 
momentos de tiempo. Recordemos como referencia la expresión del 
análisis de Fourier utilizando esta misma convención 
   


N
j
tfi
j
jetatX
1
2
Re

 3.6 
En segundo lugar, esta expansión adaptable ha permitido 
también la representación tanto de la frecuencia como de la amplitud 
aj(t) como funciones del tiempo. De tal manera, el método de análisis 
espectral desarrollado puede revelar variaciones en dependencia del 
tiempo, lo que lo convierte en un método estándar para procesos no 
estacionarios. 
 
 
Figura 8. Espectro de amplitud espectral de Hilbert correspondiente a cinco componentes modales extraídos del tacograma de intervalos R-R del 
electrocardiograma de un sujeto control. Se puede advertir que resultan más visibles los valores para los componentes 2, 3, 4 y 5 que corresponden 
a frecuencias espectrales más lentas, en tanto que las que corresponden al componente 1 que son las más rápidas y además consecuentemente menos 
intensas, resultan difíciles de detectar, siendo solo visibles algunos puntos entre los segundos 200 y 250 en las frecuencias entre 0.4 y 0.45 Hz. 
 
 
 
 
 
 
 
 
 
 
10 
 
 
Figura 9. Espectro de potencias espectrales de Hilbert del tacograma de intervalos R-R en un paciente comatoso donde pueden advertirse los 
elementos que se señalaron en la Figura 8, pero en este caso resultan de interés los valores observados del primer componente modal extraído que 
pueden ser identificados para frecuencias muy rápidas (por encima de 0.4 Hz) y que se corresponden con la llamada banda espectral de muy alta 
frecuencia de la variabilidad de la frecuencia cardiaca, en estudio por nuestro grupo y que forma parte de un trabajo en preparación para publicación. 
Los colores seleccionados son diferentes al de la Figura 8 y se han usado para buscar mejor contraste. 
 
 
 
Figura 10. Otra manera de representar el espectro de amplitudes de Hilbert, usando un fondo blanco en este caso para facilitar la identificación de 
componentes con baja intensidad de sus valores. Se trata igualmente de otro caso de paciente comatoso y el estudio se corresponde también al 
análisis del tacograma de intervalors R-R del electrocardiograma, aunque utilizando una duración de tiempo mayor de 420 segundos en este caso. 
 
Aunque estas generalizaciones son válidas, Huang et al. han 
enfatizado (Huang et al. , 2016) que las variaciones de la frecuencia 
están determinadas por portadores (carriers) de ondas rápidamente 
cambiantes ω(t). Con este abordaje las variaciones temporales de la 
amplitud a(t) y la frecuencia ω(t) podrían cubrir los procesos no 
estacionarios, pero las variaciones de la amplitud (o la energía 
proporcional al cuadrado de la amplitud) podrían aún quedar en 
términos de función del tiempo. Esto hace que el AEH asociado a la 
DEM sea extremadamente efectivo y poderoso para señales 
intramodales de frecuencia modulada. Sin embargo, no se ha tomado 
en cuenta el posible papel de las modulaciones de frecuencia y 
amplitud inter-componentes modales (AM y FM). Por todo ello, en esta 
mencionada publicación (Huang et al. , 2016) los autores han propuesto 
una modificación al método que han denominado como análisis 
espectral de Holo-Hilbert. El trabajo es muy reciente, no tiene 
importantes implicaciones para las señales con las cuales trabajamos y 
por tanto solo lo mencionamos sin profundizar en el mismo. Aún, en la 
revisión de la literatura en la base de datos Pubmed NLM no ha habido 
repercusiones con referencia a este perfeccionamiento del método. 
11 
 
Hay que decir que con mucha frecuencia se le denomina en 
la literatura al método que hemos venido analizando como método de 
Hilbert-Huang e incluso transformada de Hilbert-Huang. 
3.2 Espectro de Hilbert. 
 Una vez que tenemos los valores de amplitud o potencia 
espectral instantáneos calculados para cada uno de los componentes 
modales de interés podemos crear una matriz de datos y representarla 
en un gráfico tridimensional donde se puedan apreciar las fluctuaciones 
de la intensidad de la amplitud o potencia espectrales para cada 
frecuencia espectral instantánea en función de tiempo. 
Generalmente, el espectro de amplitudes o de potencias de 
Hilbert se representa de modo que la intensidad de los valores sea 
puesta en evidencia por una escala de colores, sobre un fondo de colorfijo como mostramos en las Figuras 8 y 9, aunque para que resulte más 
adecuada su inspección se prefiera en ocasiones que el fondo sea de 
color blanco como se muestra en la Figura 10. 
 
 
3.3 Espectro marginal de Hilbert. 
La definición formal del llamado espectro marginal de 
Hilbert se puede encontrar en la siguiente expresión 

T
dttHh
0
),()( 
 3.7 
El mismo representa a la amplitud o potencia acumulada 
durante todo el periodo de tiempo de la señal analizada. En ese sentido, 
aparentemente resulta un retroceso con respecto a las potencialidades 
del espectro de Hilbert, pero trataremos de explicar sus características 
y su superioridad al espectro de Fourier. 
Cuando se calcula el contenido espectral de una señal 
utilizando el método de Fourier la resolución espectral del proceso que 
se aplica queda “sellada” cuando se fija la duración de la secuencia que 
será objeto de análisis. La resolución espectral del proceso se puede 
calcular como el valor recíproco de la duración de la señal. Si la 
duración es digamos, de 250 segundos, la resolución más baja que 
puede detectarse será de 1/250 s = 0.004 Hz. 
Cuando se obtiene el espectro usando Fourier el primer 
valor se conoce como frecuencia cero, tiene un valor real y se considera 
que contiene el componente DC de la señal analizada, muy asociado a 
la media aritmética de los valores de la señal original. A partir de esta 
frecuencia cero, se obtiene el valor de la primera frecuencia espectral 
discreta, que en este caso será la de 0.004 Hz. La próxima frecuencia 
discreta será 0.004+0.004 = 0.008 Hz, la posterior será 0.008+0.004 = 
0.012 Hz y así sucesivamente hasta alcanzar el máximo valor posible 
de acuerdo con los parámetros que se estén utilizando. 
Para cada valor de estas frecuencias discretas, que podemos 
advertir que son valores múltiplos enteros de la resolución espectral, se 
obtienen los valores de potencia que les correspondan. Debemos notar, 
sin embargo, que si el fenómeno que nos interesa detectar posee una 
frecuencia de 0.006 Hz, entonces el proceso que estamos empleando 
no lo detectará de modo efectivo. Solamente se podrá advertir el efecto 
que el mismo pueda ejercer en las frecuencias cercanas de 0.004 y 
0.008 Hz. (u otras más por la extensión del efecto, que ha sido 
relacionado con los llamados fenómenos de “aliasing” y de fuga o 
“leakage”), pero con distorsión en sus características de magnitud de 
la potencia espectral. 
Si analizamos esa misma señal con el método de la DEM y 
la ulterior aplicación de la transformada de Hilbert, es decir, la 
aplicación del método de Hilbert-Huang, como resultado tendríamos 
un número dado de componentes modales y para cada uno de ellos los 
valores de amplitud-potencia, valores de frecuencia y de fase 
instantáneos para cada unidad de tiempo. Si por ejemplo, la señal de 
250 segundos ha sido muestreada a 200 Hz, resulta que tendremos para 
cada componente modal 50,000 de tales valores instantáneos. Las 
frecuencias instantáneas obtenidas no están sujetas a ninguna exigencia 
metodológica. Aparecen las que contiene el componente modal, pero 
tengamos en cuenta que a lo largo del tiempo que dura la señal de este 
componente modal específico, la frecuencia mostrará oscilaciones, e 
igualmente ocurrirá con los otros componentes instantáneos obtenidos 
mediante el proceso. Se obtiene así una distribución tiempo-
frecuencia-amplitud que nos permite detectar cambios con facilidad del 
proceso descrito por la señal en estudio. Este comportamiento queda 
claramente evidenciado en las Figuras 8, 9 y 10. 
Ahora bien, si deseáramos conocer cuáles fueron las 
frecuencias observadas durante todo el periodo de tiempo analizado de 
la señal y su intensidad respecto a todo ese periodo, sin importarnos el 
momento en el tiempo en que las mismas aparecieron, podríamos 
desarrollar un proceso que consistiría en obtener una especie de 
histograma, donde en la escala de las abscisas tendríamos los valores 
de las frecuencias, desde cero hasta un valor dado fijado por razones 
en las que no vamos a profundizar, pero digamos que hasta aquel valor 
de interés por parte del investigador, y en las ordenadas los 
correspondientes valores de amplitud o potencia espectral. Ahora bien, 
esta es la misma forma de representación gráfica del espectro de 
señales usando el tradicional método de Fourier. Bueno, pues es a esta 
representación, pero utilizando para construirla los valores calculados 
de frecuencias y potencias instantáneas obtenidas con el método de 
Hilbert-Huang, a lo que se le denomina espectro marginal de Hilbert. 
Siguiendo con nuestra hipotética señal muestreada a 200 Hz 
por 250 segundos, si así lo deseáramos podríamos tener el espectro 
marginal de Hilbert correspondiente a un componente modal 
específico. Para ello, necesitaríamos ordenar nuestros valores de 
frecuencia y potencia instantáneos de menor frecuencia hasta la mayor 
frecuencia, pero conservando los valores de potencia correspondientes 
a cada valor de frecuencia originalmente detectados. Nos podrá ocurrir, 
obviamente, que la misma frecuencia aparezca muchas veces, lo cual 
nos permitiría calcular la sumatoria (integración) de las potencias de 
todas ellas. También nos va a ocurrir que las frecuencias que surgieron 
del proceso no tienen que tener entre sí similares intervalos, es decir, 
una será 0.012345, la próxima podrá ser 0.012533 y la siguiente 
0.013145. 
Si tratamos de graficar nuestro espectro marginal utilizando 
todos los valores de frecuencias ascendentemente ordenadas en el eje 
de las abscisas, con los valores de potencias correspondientes 
obtendríamos un auténtico “churro”, y además, su valor sería nulo si 
deseásemos comparar este espectro marginal con el espectro de Fourier 
al que antes nos habíamos referido, y que como vimos estaba obligado 
a emplear como frecuencias discretas a valores estrictamente múltiplos 
de la resolución espectral, que en este caso era 0.004 Hz. 
Si tomamos la decisión de que nuestro espectro marginal 
posea frecuencias igualmente espaciadas y si deseamos incluso que 
estén espaciadas seleccionando un valor prefijado que nos permita 
efectuar comparaciones con espectros tradicionalmente obtenidos por 
los métodos basados en Fourier, solamente tendríamos que elaborar un 
procedimiento que nos permitiese calcular por ejemplo, mediante 
sumatoria, los valores de potencia de las frecuencias mayores a 0.004 
y menores o igual que 0.008 y ese valor adjudicarlo a la frecuencia 
0.008 Hz. Hacer lo mismo para las subsecuentes, y al final poseer para 
cada valor de frecuencia discreta de intervalo uniforme, sus 
12 
 
correspondientes valores de potencia espectral y luego representar esos 
valores en un gráfico que sería el espectro marginal de Hilbert. 
Si en vez de tomar los valores de un solo componente 
modal, tomamos todos aquellos que fueron extraídos y que muestran 
frecuencias discretas en el rango de interés de la señal que estamos 
estudiando, la disponibilidad de valores por supuesto se hace mucho 
mayor (en este caso hipotético 50,000 x número de componentes 
modales) y se abarcará obviamente todo el rango de frecuencia que se 
desea analizar. 
Si con esta representación se pierde por completo la 
dimensión tiempo, podría suponerse que los resultados no superarían a 
los que se obtienen utilizando el procedimiento tradicional de los 
métodos basados en Fourier, pero debemos tener en cuenta, que ahora 
sí habrán sido tomadas en cuenta y sin distorsiones, todas las posibles 
frecuencias discretas que puso en evidencia la transformada de Hilbert, 
y no solo aquellas que obligatoriamente está obligado a emplear el 
método tradicional. 
 
 
Figura 11. Representación gráfica de la densidad espectral de potencia de una secuencia de 420 segundos del tacograma de intervalos R-R del ECG 
en un sujetos sano, utilizando el método delperiodograma de Welch, basado en el método de Fourier (A), y del espectro marginal de Hilbert sin 
suavizado (B) y con suavizado (C). Las siglas utilizadas para identificar las bandas de frecuencia son: VLFB: banda de muy baja frecuencia; LFB: 
banda de baja frecuencia; MFB: banda de media frecuencia; HFB: banda de alta frecuencia; VHFB: banda de muy alta frecuencia. Los cálculos 
incluyeron el rango de frecuencias entre 0.02 y 0.6 Hz. 
 
En la Figura 11 mostramos en un collage las diferencias que 
pueden observarse del espectro de Fourier, así como del espectro 
marginal de Hilbert, calculado con la misma resolución espectral que 
el método de Fourier, sin suavizado (smoothing) ninguno y aplicando 
una ventana de suavizado del espectro, que hemos incluido en un 
software que hemos desarrollado en el Grupo de Trabajo, cuyo código 
desarrollado en Matlab se encuentra disponible en la URL que a 
continuación se señala: 
https://wwwresearchgatenet/publication/320161261_Matlab_Code_of
_the_Hilbert_Marginal_Spectrum_HMS. 
No caben dudas de que el espectro marginal de Hilbert, 
aunque se parece mucho al obtenido mediante el uso del periodograma 
clásico, resulta mucho más eficiente y real con respecto al contenido 
espectral de la señal estudiada. 
. 
 
Referencias 
 
Acharya UR, Kannathal N, Sing OW, Ping LY, Chua T. Heart rate 
analysis in normal subjects of various age groups. Biomed Eng Online. 
2004;3:24. 
Bender R, Schultz B, Schultz A, Pichlmayr I. Testing the Gaussianity of 
the human EEG during anesthesia. Methods of information in medicine. 
1992;31:56-9. 
Boashash B. Time Frequency Analysis: A comprehensive reference. 
First ed. Amsterdam-Boston-Heidelberg-London-New York-Oxford-
Paris-SanDiego-San Francisco-Singapore-Sydney-Tokyo: Elsevier Ltd; 
2003. 
Chen Z, Brown E, Barbieri R. Characterizing Nonlinear Heartbeat 
Dynamics within a Point Process Framework. IEEE transactions on bio-
medical engineering. 2010;10.1109/TBME.2010.2041002. 
Chua KC, Chandran V, Acharya UR, Lim CM. Application of higher 
order statistics/spectra in biomedical signals--a review. Med Eng Phys. 
2010;32:679-89. 
Colominas MA, Humeau-Heurtier A, Schlotthauer G. Orientation-
Independent Empirical Mode Decomposition for Images Based on 
Unconstrained Optimization. IEEE transactions on image processing : a 
publication of the IEEE Signal Processing Society. 2016;25:2288-97. 
Colominas MA, Schlotthauer G, Torres ME. Improved complete 
ensemble EMD: A suitable tool for biomedical signal processing. 
Biomedical Signal Processing and Control. 2014;14:19-29. 
Colominas MA, Schlotthauer G, Torres ME. An unconstrained 
optimization approach to empirical mode decomposition. Digital Signal 
Processing. 2015;DOI: 10.1016/j.dsp.2015.02.013. 
Donghoh Kim, Hee-Seok O. EMD: A Package for Empirical Mode 
Decomposition and Hilbert Spectrum. The R Journal 2009;1:40-6. 
Estévez-Báez M, Machado C, Montes-Brown J, Jas-García J, Leisman 
G, Schiavi A, et al. Very high frequency oscillations of heart rate 
variability in healthy humans and in patients with cardiovascular 
https://wwwresearchgatenet/publication/320161261_Matlab_Code_of_the_Hilbert_Marginal_Spectrum_HMS
https://wwwresearchgatenet/publication/320161261_Matlab_Code_of_the_Hilbert_Marginal_Spectrum_HMS
13 
 
autonomic neuropathy. BMC Biomedical Engineering (submitted to 
publication). 2017. 
Estévez M, Machado C. Incorporating novel methods of digital signal 
processing to traditional electrophysiological methods in patients with 
disorders of consciousness and brain death. VII International 
Symposium on Brain Death and Disorders of Conciousness Havana, 
Cuba DOI: 1013140/RG2113638484. Havana, Cuba2015. 
Estévez M, Montes Brown J, Hernández-Cruz A, Machado C. Heart 
Rate Variability Spectral Analysis and Orthostatic Tests: A Time-
Frequency Approach. Functional Neurology Rehabilitation and 
Ergonomics. 2012;2:51-2. 
Faes L, Zhao H, Chon KH, Nollo G. Time-varying surrogate data to 
assess nonlinearity in nonstationary time series: application to heart rate 
variability. IEEE transactions on bio-medical engineering. 2009;56:685-
95. 
Flandrin P, Rilling G, Gonçalvés P. Empirical Mode Decomposition as 
a Filter Bank. IEEE Signal Processing Letters 2004;11:112-4. 
Guarín D, Delgado E, Orozco A. Testing for nonlinearity in non-
stationary physiological time series. 2011. p. 2671-4. 
Huang NE, Hu K, Yang ACC, Chang HC, Jia D, Liang WK, et al. On 
Holo-Hilbert spectral analysis: a full informational spectral 
representation for nonlinear and non-stationary data. Phil Trans R Soc 
A. 2016;374: 20150206. 
Huang NE, Shen Z, Long SR, Wu MC, Shih HH, Zheng Q, et al. The 
empirical mode decomposition and the Hilbert spectrum for nonlinear 
and non-stationary time series analysis. Proc R Soc Lond. 
1998;454:903-95. 
Kipinski L, Konig R, Sieluzycki C, Kordecki W. Application of modern 
tests for stationarity to single-trial MEG data: transferring powerful 
statistical tools from econometrics to neuroscience. Biological 
cybernetics. 2011;105:183-95. 
Li C, Ding GH, Wu GQ, Poon CS. Fractal, entropic and chaotic 
approaches to complex physiological time series analysis: a critical 
appraisal. Conference proceedings : Annual International Conference 
of the IEEE Engineering in Medicine and Biology Society IEEE 
Engineering in Medicine and Biology Society Conference. 
2009;2009:3429-32. 
Machado C, Defina P, Estévez M, Leisman G. A reason for care in the 
clinical evaluation of function on the spectrum of consciousness. BMJ 
Case Reports (submitted to publication). 2017. 
Machado C, Estévez M, Pérez-Nellar J, Schiavi A. Residual vasomotor 
activity assessed by heart rate variability in a brain-dead case. BMJ Case 
Reports 04/2015; 2015(apr01 1), DOI:101136/bcr-2014-205677. 2015. 
Machado C, Estevez M, Rodriguez R, Perez-Nellar J, Chinchilla M, 
DeFina P, et al. Zolpidem arousing effect in persistent vegetative state 
patients: autonomic, EEG and behavioral assessment. Current 
pharmaceutical design. 2014;20:4185-202. 
Magagnin V, Bassani T, Bari V, Turiel M, Maestri R, Pinna GD, et al. 
Non-stationarities significantly distort short-term spectral, symbolic and 
entropy heart rate variability indices. Physiological Measurement. 
2011;32:1775-86. 
McEwen JA, Anderson GB. Modeling the stationarity and Gaussianity 
of spontaneous electroencephalographic activity. IEEE transactions on 
bio-medical engineering. 1975;22:361-9. 
Persson J. Comments on "Modeling the stationarity and Gaussianity of 
spontaneous electroencephalographic activity". IEEE transactions on 
bio-medical engineering. 1977;24:302. 
Pradhan C, Jena SK, Nadar SR, Pradhan N. Higher-order spectrum in 
understanding nonlinearity in EEG rhythms. Computational and 
mathematical methods in medicine. 2012;2012:206857. 
Pritchard WS, Stam CJ. Nonlinearity in human resting, eyes-closed 
EEG: an in-depth case study. Acta neurobiologiae experimentalis. 
2000;60:109-21. 
Pustelnik N, Borgnat P, Flandrin P. A multicomponent proximal 
algorithm for empirical mode decomposition. In: EUPSICO, editor. 
20th European Signal Processing Conference Proceedings pp 1880–
1884: EUPSICO; 2012. p. 1880–4. 
Shahid S, Prasad G. Bispectrum-based feature extraction technique for 
devising a practical brain-computer interface. Journal of neural 
engineering. 2011;8:025014. 
Silva LEV, Lataro RM, Castania JA, Aguiar Silva CA, Salgado HC, 
Fazan R, Jr., et al. Nonlinearities of heart rate variability in animal 
models of impaired cardiac control: contribution of different time scales. 
Journal of applied physiology. 2017:jap 00059 2017. 
Torres ME, Colominas MA, Schlotthauer G, Flandrin P. A complete 
ensemble empirical mode decomposition with adaptive noise. 36th 
IEEE Int Conf on Acoust, Speech and Signal Process. Prague, Czech 
Republic: ICASSP; 2011. p. 4144–7. 
Wang R, Wang J, Li S, Yu H, Deng B, Wei X. Multiple featureextraction and classification of electroencephalograph signal for 
Alzheimers' with spectrum and bispectrum. Chaos. 2015;25:013110. 
Wang YH, Yeh CH, Young HWV, Hu K, Lo MT. On the computational 
complexity of the empirical mode decomposition algorithm. Physica A: 
Statistical Mechanics and its Applications. 2014;400:159-67. 
Weber EJ, Molenaar PC, van der Molen MW. A nonstationarity test for 
the spectral analysis of physiological time series with an application to 
respiratory sinus arrhythmia. Psychophysiology. 1992;29:55-65. 
Wu Z, Huang NE. Ensemble empirical mode decomposition: a noise-
assisted data analysis method. Adv Adapt Data Anal. 2009;1:1-41.