ES2855979T3 - Descomposición de señales no estacionarias en componentes funcionales - Google Patents

Abstract

Un método implementado por ordenador (100) para descomponer una señal no estacionaria representativa de un fenómeno fisiológico en un conjunto de componentes funcionales autosimilares con el QGK (núcleo cuasi- Gaussiano) como elemento base, y para estimar los parámetros que caracterizan los respectivos componentes funcionales, el método que comprende las etapas de: aplicar filtros de media móvil de múltiples pasadas para estimar las tendencias dinámicas y de línea de base en el transcurso del tiempo de la señal (etapa 102), usar las tendencias para extraer una función de guía indicativa de los fragmentos de la señal sobre los que se desarrollan los componentes funcionales, y segmentar de forma adaptativa la función de guía mediante el uso de cruces por cero (ZC) como los puntos de segmentación; dividir la señal no estacionaria en las regiones de desarrollo de componentes, al segmentar la señal no estacionaria en puntos correspondientes a los puntos de segmentación de la función de guía (etapa 103); compensar la superposición de los QGK que coinciden con los componentes funcionales anteriores para formar un ECK (núcleo de componente empírico) (etapa 104), las subetapas recursivas de: transformar el ECK en un dominio de la frecuencia mediante el uso de un algoritmo de función de base similar (SBF) (etapa 105); estimar los parámetros κ (peso) y σ (dispersión) (etapa 106) y validar la fiabilidad del parámetro en términos de grados categóricos E, G, S, W y U, en donde: E - excelente, G - bueno, S - satisfactorio, W - regular y U - no satisfactorio al encontrar y evaluar la mejor coincidencia entre el espectro de amplitud normalizado del ECK y la plantilla del dominio de la frecuencia con un perfil Gaussiano; estimar el tiempo de inicio τ del QGK con los parámetros identificados κ y σ de la parte imaginaria del espectro complejo del ECK al encontrar puntos característicos que corresponden a los valores definidos del espectro de fase (etapa 107); expandir el modelo de la señal no estacionaria compuesta a partir de la suma de los QGK identificados en los ciclos anteriores de recursión al añadir el QGK recién identificado (etapa 108); y reorganizar los componentes identificados después de eliminar los QGK desfavorables mediante el uso de un parámetro categórico ξ para evaluar los resultados de la identificación (etapa 109); crear modelos parciales de la señal no estacionaria mediante la selección de conjuntos QGK cuyos parámetros pertenecen a clases predefinidas (etapa 110).

Description

DESCRIPCIÓN

Descomposición de señales no estacionarias en componentes funcionales

Campo de la invención

La presente invención se refiere, en general, al procesamiento de señales y, en particular, a la descomposición de señales biomédicas no estacionarias en componentes funcionales, y a la estimación de los parámetros que caracterizan cada uno de esos componentes funcionales.

Antecedentes

Las señales no estacionarias se producen por sistemas con atributos que varían con el tiempo. La necesidad de comprender y analizar dichas señales no estacionarias desde los puntos de vista teórico, de investigación y de diagnóstico surge en casi todos los campos de investigación biomédica. Sin embargo, las señales no estacionarias no están soportadas por la teoría general de señales, y los métodos matemáticos adecuados para tratar con señales no estacionarias son bastante limitados.

Los métodos de procesamiento de señales actuales se encaminan generalmente en la teoría de sistemas estacionarios lineales. En consecuencia, se ignora frecuentemente el carácter no estacionario de los datos reales o se presume que los efectos son insignificantes. Aunque las suposiciones correspondientes de la linealidad e invarianza en el tiempo proporcionan resultados aceptables para algunas aplicaciones, muchos fenómenos fisiológicos son altamente variables y no lineales, y por lo tanto no admiten una aproximación lineal razonablemente exacta. Mientras tanto, una gran cantidad de evidencias científicas y clínicas indican que se espera que un análisis más completo de las principales señales fisiológicas no estacionarias, tales como el electrocardiograma (ECG) y el electroencefalograma (EEG), conduzca a una comprensión más profunda de los mecanismos y estados de varios sistemas fisiológicos bajo diferentes condiciones normales y patológicas. En este contexto, los desarrollos de algoritmos avanzados del procesamiento de señales no estacionarias emergen como enfoques capaces para una serie de aplicaciones cada vez más desafiantes, tales como imágenes del cerebro y el cuerpo humanos, monitoreo de la salud implementado por ordenador e interfaces cerebro-computadora basadas en EEG. Al mismo tiempo, los continuos avances de los sistemas biomédicos computarizados, en términos de velocidad, tamaño, y costo, hicieron que el desarrollo de algoritmos sofisticados fuera práctico y efectivo.

El análisis de señales no estacionarias no está soportado por marcos teóricos y computacionales adecuados, tales como en el caso de la teoría de sistemas estacionarios lineales. Generalmente se acepta un enfoque heurístico, más que un enfoque matemáticamente riguroso, en el análisis de tipos particulares de señales fisiológicas. Por ejemplo, al analizar señales de EEG, dichas señales se pueden analizar mediante el uso de la noción de estacionariedad local (ver Barlow, J.S. "Methods of analysis of non-stationary EEGs, with emphasis on segmentation techniques: a comparative review". J. Clin. Neurophysiol. 1985, volumen 2, páginas 267-304). Esta noción supone que, aunque la señal es inherentemente no estacionaria, dentro de pequeños intervalos de tiempo el patrón de la señal se aparta solo ligeramente de la estacionariedad (ver Priestley, M.B.: "Non-Linear and Non-Stationary Time Series Analysis". Prensa Académica, Londres, 1988). Dichos elementos cuasi-estacionarios se asocian a diferentes componentes de la señal que se pueden producir por distintos sistemas no lineales.

En Melkonian y otros "Chaotic Dynamics of ECG Waveforms at Global and Microscopic Scales: Theory and Applications"; The Open Cybernetics & Systemics Joural, volumen 5, número 1, 24 de mayo 2011, páginas 32-44, XP055339387, ISSN: 1874-110X se divulga una teoría para describir el electrocardiograma (ECG) en términos de procesos celulares subyacentes de transporte de iones. La evolución del ECG a lo largo del tiempo se considera una secuencia de potenciales transitorios autosimilares parcialmente superpuestos, con el potencial de masa genérico (GMP) como el elemento básico. Mediante el uso de ecuaciones del proceso de nacimiento y muerte no homogéneo (BDP), se deduce un modelo de partículas de GMP en forma de BDP caótico. Sin embargo, este documento no dice nada sobre la validación de la fiabilidad de los parámetros en términos de grados categóricos.

Hasta ahora, el problema de detectar componentes de señales no estacionarias se ha abordado mediante el uso de metodologías de procesamiento de señales basadas en datos y basadas en modelos.

La ventaja de los métodos de procesamiento de señales basados en datos es su capacidad para describir la dinámica de las formas de onda empíricas complejas en detalle. Uno de dichos métodos basados en datos, la descomposición empírica en modos (EMD), se ha desarrollado por Huang y otros. (Huang y otros: "The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis". Proc. R. Soc. London, Ser. A, 1998, volumen 454, páginas 903-995) para descomponer sobre una base adaptativa una señal no estacionaria en la suma de funciones oscilatorias de "buen comportamiento", denominadas funciones de modo intrínseco. En los siguientes desarrollos el EMD se ha soportado por la transformada de Hilbert que proporcionó frecuencias instantáneas para cada función de modo intrínseco. Mediante el uso de estas medidas, la técnica describe el sistema en términos de cambios de frecuencia instantáneos. Sin embargo, la solución no se define únicamente porque esto depende de una serie de parámetros empíricos, cuya elección no está soportada por criterios estrictos. Sin una base matemática firme del EMD, este enfoque puramente empírico e intuitivo crea una serie de problemas metodológicos difíciles. Por ejemplo, la función de modo intrínseco se define solo numéricamente como la salida de un algoritmo iterativo, sin definición analítica. Bajo estas circunstancias, las opiniones sobre diferentes aspectos de la metodología provienen de comparaciones caso por caso realizadas empíricamente.

El enfoque basado en modelos para el procesamiento de la señal biomédica requiere modelos precisos de los sistemas involucrados. Un aspecto específico del problema es que las señales fisiológicas son las variables de escala global que se producen por múltiples fuentes de escala microscópica. La creación de un modelo físico requiere una identificación estricta de los elementos relevantes de escala microscópica y la síntesis, sobre estas bases, el modelo global. Sin embargo, la transición desde la escala global a la microscópica mediante el uso de estas nociones estrictamente deterministas no tiene una solución única porque es imposible delinear exactamente un patrón distintivo de los fenómenos celulares y moleculares involucrados. Ante tal incertidumbre, numerosos esfuerzos para crear modelos de la maquinaria celular que da lugar a las señales de escala global se soportan por una notable variedad de enfoques heurísticos que difieren ampliamente no solo en los detalles fisiológicos y anatómicos de los modelos, sino también en las herramientas matemáticas básicas. Se desconoce el grado en el que los modelos están en contradicción.

En consecuencia, todavía existe una necesidad, no satisfecha por las metodologías actuales, de identificar efectivamente los componentes de una señal de serie temporal no estacionaria, de manera que los modelos correspondientes se puedan usar eficazmente para clasificar los atributos de la señal.

Resumen

Es un objeto de la presente invención, cuyo alcance se define por las reivindicaciones adjuntas 1-7, superar sustancialmente, o al menos mejorar, una o más desventajas de los arreglos existentes.

De acuerdo con un primer aspecto de la presente invención, se proporciona un método para descomponer una señal no estacionaria en el conjunto de componentes funcionales autosimilares con el núcleo cuasi-Gaussiano (QGK) como elemento base, y para estimar parámetros que caracterizan respectivos componentes funcionales, el método que comprende las etapas de:

estimar tendencias dinámicas y de línea base en el transcurso de tiempo de la señal no estacionaria, al extraer la función de guía y su segmentación adaptativa;

dividir la señal no estacionaria en las regiones del desarrollo de componentes;

compensar recursivamente la superposición de las QGK emparejadas con los componentes funcionales anteriores para formar un núcleo de componente empírico (ECK) en cada ciclo de recursión, las subetapas de:

transformar el ECK en un dominio de frecuencia;

estimar los parámetros del QGK y su validación; y

la creación de modelos globales y parciales de la señal no estacionaria.

De acuerdo con otro aspecto de la presente invención se proporciona un aparato para implementar el método descrito anteriormente.

Breve Descripción de los Dibujos

Una o más de las modalidades de la presente invención se describirán ahora con referencia a los dibujos, en los cuales:

La Figura 1A ilustra una forma típica de un potencial eléctrico medido entre dos sensores que se colocan en la superficie de un cuero cabelludo;

La Figura 1B ejemplifica el carácter altamente no estacionario de una señal de EEG;

La Figura 2A ilustra una señal de típica ERP provocada en un paradigma auditivo excéntrico;

La Figura 2B ilustra los resultados típicos de sustituir cada segmento mostrado en la Figura 2A con un núcleo cuasi-Gaussiano que tiene parámetros que coinciden con ese segmento;

La Figura 2C ilustra el modelo obtenido mediante la suma de los núcleos cuasi-gaussianos componentes, así como también la señal original ERP mostrada en la Figura 2A para comparación;

La Figura 3 muestra un diagrama de flujo esquemático de un método de análisis de forma de onda que descompone de forma adaptativa una señal fisiológica en componentes funcionales, y después determina los parámetros que caracterizan cada componente funcional;

La Figura 4A muestra un diagrama de flujo esquemático de la estimación de tendencia dinámica y de línea de base aplicada en el método que se muestra en la Figura 3;

La Figura 4B muestra una señal fisiológica típica, una tendencia de línea de base resultante de aplicar el promedio de ventana móvil de múltiples pasadas a la señal fisiológica discreta, y una señal de línea de base sin tendencia. Las Figuras 5A y 5B ilustran puntos de segmentación que resultan de cruces por cero y puntos de mínimos globales y locales respectivamente;

La Figura 5C muestra un ejemplo de segmentación de señales de EEG con cruces por cero y mínimos ilustrados por líneas verticales continuas y de puntos, respectivamente.

La Figura 6A muestra una señal inicial y una función de guía que proporcionan puntos de segmentación con alta resolución temporal;

Las Figuras 6B y 6C replican la función de guía y la señal inicial de la Figura 6A con puntos de segmentación que se indican por líneas verticales;

La Figura 7A muestra la señal inicial de la Figura 6A y la función de guía que proporciona la segmentación de señal con baja resolución temporal;

Las Figuras 7B y 7C replican la función de guía y la función inicial de la Figura 7A con puntos de segmentación que se indican por líneas verticales;

La Figura 8A ilustra los componentes Gaussianos que se usan para la composición de un núcleo cuasi-Gaussiano; La Figura 8B ilustra un ejemplo de núcleo cuasi-Gaussiano;

La Figura 9A ejemplifica una descomposición de un interpolante en la suma de elementos finitos;

La Figura 9B muestra la descomposición de una función de sombrero en la suma de tres funciones de base triangular;

La Figura 9C muestra una función de base triangular usada en un algoritmo de función de base similar;

La Figura 9D muestra las partes real e imaginaria del espectro complejo de la función de base triangular;

La Figura 10A ilustra la plantilla de amplitud usada para estimar un parámetro a.

La Figura 10B ejemplifica el espectro de amplitud empírica normalizada y las diferentes posiciones T1, T2 y T3 de la plantilla de amplitud durante la estimación del parámetro a;

La Figura 11A ilustra el resultado de la coincidencia de la plantilla con el espectro de amplitud normalizado que pertenece a un grado E;

La Figura 11B ilustra el resultado de la coincidencia de la plantilla con el espectro de amplitud normalizado que pertenece a un grado G;

La Figura 12 ejemplifica varios resultados de la plantilla de amplitud que coinciden con los espectros de amplitud normalizados de diferentes segmentos de la señal de EEG mostrada en el panel central de la figura;

La Figura 13A muestra un ejemplo de segmentación de EEG mediante el uso de baja resolución temporal;

La Figura 13B muestra los espectros de amplitud de los 10 segmentos mostrados en la Figura 13A;

La Figura 13C muestra la superposición de los espectros de amplitud de la Figura 13B después de la normalización tanto de la amplitud como de la frecuencia;

Las Figuras 14A muestra los puntos a y b en la escala de frecuencia de abscisas de acuerdo con la cual se obtienen las estimaciones de un parámetro t;

La Figura 14B ejemplifica una situación en la que el punto b en la escala de abscisas se define por el mínimo de la parte imaginaria del espectro complejo;

La Figura 15A muestra un ejemplo de una forma de onda respiratoria de una persona al hablar;

La Figura 15B muestra las formas de onda de 54 QGK que se identifican en la señal de la Figura 15A;

La Figura 15C muestra la precisión de la reproducción de la señal de la Figura 15A mediante el modelo que se sintetiza como la suma de 54 QGK identificadas.

Descripción detallada

Cuando se hace referencia en uno o más de los dibujos acompañantes a las etapas y/o características, que tienen los mismos números de referencia, esas etapas y/o características tienen para los fines de esta descripción la(s) misma(s) función(es) u operación(es), a menos que aparezca la intención contraria.

Para describir una señal no estacionaria, una etapa crucial es crear un modelo. Los enfoques convencionales de este problema se soportan en modelos "basados en la física": el objetivo es describir la señal en términos de un modelo físico explícito de sus fuentes (verCandy J.V. "Model-based signal processing". John Wiley & Sons, Nueva Jersey, 2006). Sin embargo, las señales no estacionarias se producen generalmente por múltiples fuentes no lineales, cuya funcionalidad sigue sin estar clara.

Por el contrario de estos enfoques basados en la física, el método divulgado en la presente memoria se soporta con un enfoque de modelado basado en el caos (CBM). Esto tiene en cuenta de forma combinada tanto de los factores deterministas como de los probabilísticos que definen el comportamiento de la fuente. El marco general de este enfoque se introduce mediante referencia a los mecanismos generales que subyacen a la génesis de las señales biomédicas.

Un aspecto general del análisis de las señales biomédicas es que el proceso fisiológico es una variable de escala global producida por múltiples fuentes elementales que actúan en la escala macroscópica. En este contexto, los procesos fisiológicos reflejan el desempeño externo de sistemas monumentalmente complejos de los elementos moleculares y celulares, cuyo desempeño se rige tanto por factores deterministas como estocásticos. Una dificultad importante radica en el hecho de que los modelos de conjuntos celulares relevantes no son factibles, ya que el número de los elementos funcionales requeridos es demasiado alto y los detalles de la maquinaria celular son demasiado intrincados. Sin embargo, numerosos estudios teóricos y computacionales del efecto de masa producido por múltiples variables de escala microscópica arrojan algo de luz sobre aquellos aspectos generales de la maquinaria celular que son significativos en la escala global de aquellos que no lo son. Por lo tanto, un concepto físico general es que una activación sincrónica de un gran número de fuentes elementales en un volumen de escala microscópica local es una condición necesaria para producir un cambio detectable en la señal medida a escala global. Los aspectos principales de estos procesos se han estudiado ampliamente para tales señales biomédicas clínicamente importantes como el electroencefalograma (EEG), el electrocardiograma (ECG) y el electromiograma (EMG). El marco general que sugiere una forma de onda pico de escala global como marcador de un componente fisiológico distinto puede ilustrarse con referencia a los mecanismos subyacentes a la génesis del EEG.

Como se ilustra en la Figura 1A, el EEG se puede medir mediante el uso de sensores colocados en la superficie de un cuero cabelludo. El potencial eléctrico se produce por la actividad de las neuronas en las áreas corticales del cerebro más abajo del cuero cabelludo. La actividad de una sola neurona no es lo suficientemente alta para producir un potencial medible entre los sensores colocados en la superficie del cuero cabelludo. En cambio, cuando grandes conjuntos de neuronas corticales localizadas cerca trabajan en sincronía, su actividad total produce dipolos extracelulares a partir de los cuales resulta el potencial eléctrico medido.

Las fluctuaciones y la dinámica de estos potenciales de masa registrados en la superficie del cuero cabelludo siguen la actividad de los procesadores neurales corticales sin periodos de retraso. En este contexto las señales de EEG sirven como uno de los marcadores objetivos más importantes para el estudio en tiempo real de los mecanismos básicos del procesamiento de la información en el cerebro humano. Un ejemplo de una señal de EEG mostrada en la Figura 1B ilustra que la actividad eléctrica del cerebro es típicamente tanto altamente irregular como no estacionaria, y produce una variedad de patrones de actividad.

Una clave para la comprensión del contexto de procesamiento de la información de las señales de EEG se proporciona mediante la detección de cambios en las señales de EEG, y su relación en el tiempo con estímulos particulares durante una tarea cognitiva. Dichos fragmentos específicos de señales de EEG se denominan potenciales relacionados con eventos (ERP). La Figura 2A ilustra un ERP típico obtenido en un paradigma auditivo excéntrico. El paradigma auditivo excéntrico es una prueba durante la cual el carácter de la señal del EEG se afecta por la aplicación de dos tipos de tonos auditivos, específicamente los estímulos de fondo y objetivo. En la señal de EEG ilustrada en la Figura 2A, el estímulo objetivo se aplica ent= 0. A este evento le sigue una sucesión de deflexiones de forma de onda positivas y negativas. Sobre la base de evidencia experimental y clínica significativa, los picos de forma de onda que aparecen en intervalos de tiempo específicos después de la aplicación del estímulo se consideran entidades funcionales que se asocian con varias redes neuronales corticales. En el paradigma auditivo excéntrico, tales potenciales pico se conocen como latencia tardía ERP, denominados como N100, P200, N200 y P300 respectivamente. El número es representativo de un tiempo aproximado (en milisegundos) de la aparición de un pico negativo (N) o positivo (P) después de la aplicación del estímulo. Un aspecto específico de P300 es que al examinar los promedios de los ERP repetidos, este componente de latencia tardía aparece como un pico monolítico. Sin embargo, al examinar una sola forma de onda de ERP, es decir, la señal de EEG original, la latencia tardía ERP P300 aparece generalmente como un potencial complejo que consta de dos subcomponentes denominados P3a y P3b.

Las distintas fluctuaciones positivas-negativas etiquetadas en tal manera aparecen como marcadores objetivos de procesos activados en serie que están íntimamente relacionados con sistemas fisiológicos importantes. El concepto de una forma de onda pico como marcador de un componente producido por un sistema funcional con una organización de estructura-función específica se aplica convencionalmente a casi todos los procesos fisiológicos de importancia científica y clínica significativa. Dado este enfoque ampliamente aceptado, la variedad de diferentes técnicas de medición manuales o implementadas por ordenador (procedimientos de selección de pico) reducen la señal fisiológica a las amplitudes pico y latencias pico. Una limitación crítica de estas técnicas en gran parte empíricas es que la reducción de una señal a amplitudes y latencias pico en puntos de tiempo aislados no puede caracterizar la dinámica de las formas de onda biológicas que contienen información funcional y clínica importante. En un contexto computacional, el inconveniente de todas estas técnicas es su incapacidad para resolver la superposición temporal de los componentes. Sin correcciones de superposición, las mediciones de los parámetros de amplitud y latencia se pueden ver afectadas por errores no controlados significativos.

La Figura 3 muestra un diagrama de flujo esquemático de un método 100 de análisis de la forma de onda que considera que un componente no es solo un pico en la forma de onda, sino una desviación (positiva o negativa) con una forma específica. Más particularmente, el método 100 descompone de forma adaptativa una señal no estacionaria v(t) en componentes funcionales, y después aplica análisis de tiempo-frecuencia y procedimientos de estimación de parámetros para identificar cada componente funcional en términos de modelo dinámico y parámetros relevantes. La señal no estacionaria v(t) puede ser una señal de EEG, una señal de electrocardiograma (ECG), una señal de electromiograma (EMG), una señal derivada de la respiración, etc.

El método 100 se puede implementar mediante el uso de software, tal como uno o más programas de aplicación, ejecutados dentro de un sistema informático. El método 100 se puede implementar alternativamente en hardware dedicado, tal como uno o más circuitos integrados, que realizan las funciones o subfunciones del método 100. Dicho hardware dedicado puede incluir procesadores de señales digitales, o uno o más microprocesadores y memorias asociadas.

El método 100 comienza en la etapa 101 donde la señal no estacionaria v(t), que es una señal analógica, se condiciona para formar una secuencia de entrada discreta que consta de valores muestreados vm=v(tm). Más particularmente, la señal no estacionaria v(t) se amplifica dentro de un rango de voltajes estándar, y se muestrea a intervalos de muestreo regulares tm=mA, donde m es un entero y A es el intervalo de muestreo. Esto produce un conjunto de muestras v = {vo, ..., vm,..., vm}. La señal no estacionaria v(t) se puede considerar como una señal de banda limitada, y el intervalo de muestreo A se elige para exceder la velocidad de Nyquist, al ser la velocidad de muestreo más baja que es necesaria para la restauración precisa de la señal no estacionaria muestreada v(t).

Como se describe en detalle más abajo, el método 100 después opera en la secuencia de entrada discreta v = {v0, ..., vm, ..., vm}, al ser la representación discreta de la señal no estacionaria v(t), para dividir la secuencia de entrada v = {v⁰, ..., vm, ..., vm} en segmentos, para aplicar el análisis de tiempo-frecuencia a cada segmento e identificar para cada segmento el modelo dinámico con los parámetros adaptados al segmento. Se forma así una secuencia de parámetros, y esa secuencia de parámetros se usa para clasificar cada segmento de acuerdo con los criterios definidos por el usuario y también sintetizar la señal no estacionaria v(t) o sus segmentos definidos.

La señal no estacionaria v(t), y por lo tanto también su contraparte discreta v = {v⁰,...,vm,...,vM}, típicamente contiene ruido de alta frecuencia y actividad de fondo de baja frecuencia. Estas actividades irrelevantes generalmente se eliminan de análisis posteriores mediante el uso de filtros de paso bajo y paso alto. Sin embargo, debido a la naturaleza altamente irregular y no estacionaria tanto de los componentes funcionales significativos como de las actividades irrelevantes, la aplicación de filtros con características de frecuencia invariante en el tiempo puede distorsionar de manera significativa e impredecible la dinámica de los componentes funcionales, con la pérdida de información funcionalmente importante.

Una solución general a este problema es usar un filtro adaptativo. Dicho filtro adaptativo debe adaptar su rendimiento a las tendencias cambiantes en la dinámica de la señal. El filtro digital más comúnmente usado para la estimación de tendencias, el filtro de media móvil o boxcar, es un filtro digital con respuesta de impulso finita que reemplaza cada muestra con el promedio de la muestra en cuestión y un número determinado de puntos consecutivos anteriores y posteriores. Sin embargo, la respuesta de frecuencia del filtro boxcar está lejos de ser óptima y el promedio de movimiento simple puede causar distorsiones sustanciales en las trayectorias de tendencia. Para evitar tales distorsiones, la etapa 102 usa un filtro de media móvil de múltiples pasadas (MPMA) para eliminar el ruido de alta frecuencia y la actividad de fondo de baja frecuencia de la secuencia de entrada v. Más particularmente, se aplican dos filtros de MPMA a la secuencia de entrada v para estimar las tendencias dinámicas y de línea base en la secuencia de entrada v.

El MPMA es un proceso iterativo, del cual cada iteración aplica un promedio de ventana móvil a una secuencia de entrada, al ser la secuencia de entrada de la primera iteración la secuencia de entrada v, y la secuencia de entrada de las iteraciones subsiguientes es el resultado de la iteración anterior. En consecuencia, en la iteración k-th, las secuencias de entrada y salida son

y

respectivamente, donde el número de muestras es M+1, k de 1 a K es el número de la iteración, y v 0 = v.

La secuencia de salida vk+1 para la muestra mth viene dada por lo siguiente:

donde P define el ancho de la ventana, el ancho de la ventana es Nw= 2P+1.

Resulta claro de la Ecuación (1) que

se define por los puntos anteriores consecutivos P y los puntos posteriores consecutivos P. Sin embargo, de acuerdo con la definición de secuencia v, la primera muestra de los procesos en cuestión es Vq+1

Por consiguiente, los valores

con m+j<0 se excluyen de los cálculos. Esto se puede considerar como relleno de ceros que consiste en extender la secuencia v con valores cero vn=0 para n de -1 a -P.

Mientras que el promedio de ventana móvil convencional produce distorsiones de alta frecuencia en un espectro de la secuencia resultante, el MPMA que se aplica en la etapa 102 elimina gradualmente tales distorsiones del espectro.

^{Las residuales} á ^A v* ^jfc+i = v¡ ^l ¡ ⁺ ) ^l -Yj ^k g

con m desde P a M-P disminuyen después de cada iteración. Una mejora en la estimación de tendencias después de la iteraciónk-th se evalúa al medir el cambio de los valores muestreados entre las iteraciones. La precisión se expresa en términos del siguiente error cuadrático medio:

Las iteraciones terminan cuando se logra una mejora deseada. Cuando se aplica a señales electrofisiológicas típicas, como EEG, ERP, EMG y ECG, las iteraciones se terminan después de la iteración k-th si:

dondeMv es el valor medio de la señal de salida vk, y Tres el umbral adimensional. Dado 0,05 como umbral razonable Tren estas aplicaciones, las iteraciones se terminan típicamente después de 5 iteraciones.

La Figura 4A muestra un diagrama de flujo esquemático de la estimación de tendencia dinámica y de línea base que se aplica en la etapa 102. Se aplican dos MPMA separados en la secuencia de entrada v en las unidades MPMA 1 y MPMA 2 para formar secuencias filtradas: tendencia dinámica d= {do,...,dm,...,dM} y tendencia de referencia b = {bo,...,bm,...,bM} MPMA 1 y MPMA 2 tienen diferentes parámetros de ventana, NW1y NW2, que debe satisfacer la condición Nw2 s Nw1 +2. La secuencia de entradavy las secuencias filtradasdyb se combinan luego para producir dos secuencias de salida g y u.

La secuencia u = {u0,...,Um,...,UM} con um = vm-bm, representa una señal de la línea de base sin tendencia, es decir, la secuencia primaria v después de la sustracción de la tendencia de la línea de base b.

La Figura 4B ejemplifica el cálculo de la señal fisiológica sin tendencia de la línea base u. La línea continua es una porción de 201 puntos de una señal EEG típicavmuestreada en A = 0,004 s. La línea de puntos es la línea de basebresultante de aplicar el MPMA 2 con una ventana de 31 puntos y mediante el uso de 5 iteraciones. Después de eliminar la línea de basebde la señal de EEGv,la secuencia sin tendenciause muestra mediante la línea discontinua.

La secuencia g = {g0,...,gm,...,gM} con gm = dm- bm., se denomina función de guía. La función de guía se obtiene mediante el uso de la eliminación de componentes irrelevantes de alta y baja frecuencia del transcurso de tiempo de la señal. Estas operaciones de purificación soportan el papel de secuenciagcomo guía para la búsqueda de componentes funcionales significativos.

Haciendo referencia nuevamente a la Figura 3, después de la eliminación de la tendencia en la etapa 102, el método continúa a la etapa 103 donde se realiza la segmentación en la función de guía g. El objetivo de la segmentación realizada en la etapa 103 es estimar los puntos de segmentación que definen los segmentos de señal asociados con los componentes funcionales. Los puntos de segmentación se definen como cruces por cero (ZC) ilustrados en la Figura 5A y puntos de mínimos globales y locales del módulo de la función de guía ilustrados en la Figura 5B. La Figura 5C muestra un ejemplo de la función de guía gcon cruces por cero ilustrados en líneas verticales continuas y mínimos ilustrados en líneas verticales discontinuas.

Más particularmente, el punto tmse califica como un cruce por cero si:

El punto tm especifica un mínimo global o local de |g| = {|g0|,...,|gm|,...,|gM|} si:

Una secuencia de puntos de segmentación (xo,..., x ¡,...,xn} se forma así, con las longitudes -^i ,...,A,...,^n de los segmentos que se definen por los puntos de segmentación {xo,...,x ¡,..,Xn} siendo:

Generalmente, la función de guía g se puede considerar como una aproximación a la señal sin tendencia u. La precisión de la aproximación se puede modificar mediante el uso de diferentes parámetros del MPMA 1. En consecuencia, la resolución temporal del procedimiento de segmentación se rige por una elección adecuada de los parámetros del MPMA 1.

Las Figuras 6A y 7A ilustran dos funciones de guía (líneas discontinuas) calculadas para una y la misma señal de EEG sin tendenciat (líneas grises) mediante el uso de diferentes parámetros para MPMA 1. Más particularmente, las funciones de guía g1 y g2, que se muestran separadamente en la Figura 6By la Figura 7B, se han calculado mediante el uso de ventanas de 3 y 7 puntos, respectivamente. Ambos casos se han calculado mediante el uso de 2 iteraciones. Los resultados mostrados en la Figura 6A indican una discrepancia menor entre la secuencia sin tendencia u y la función de guía g1.. Por lo tanto la función de guía g1 reproduce la señal sin tendencia en detalle, y en consecuencia, se puede considerar como una reproducción prácticamente precisa de la secuencia sin tendencia u. Los puntos de segmentación indicados por líneas verticales incluyen 11 cruces por cero y 8 mínimos.

Por el contrario, los resultados mostrados en la Figura 7A, y también la comparación de la Figura 6B con la Figura 7B, muestran que el mismo procedimiento de MPMA con una ventana más grande elimina los mínimos de los puntos de segmentación y como resultado reemplaza las formas de onda pico complejas por formas de onda monolíticas. En algunos casos, la importancia del proceso es que elimina actividades funcionalmente irrelevantes, y enfatiza las tendencias dinámicas dominantes en el dominio del tiempo de una señal fisiológica de escala global como EEG, EMG, ECG, etc. Sin embargo, las fluctuaciones irregulares de la señal de alta frecuencia se pueden ser producir por fuentes físicamente significativas. Por ejemplo, en el caso de las actividades de EEG con frecuencias superiores a 30 Hz se conocen como ritmos gamma. Por tanto, en función del carácter de la aplicación, tanto la segmentación de alta resolución como la segmentación de baja resolución, así como también sus diferentes combinaciones proporcionan potentes herramientas computacionales para la investigación de la composición de componentes de señales complejas no estacionarias.

Aunque las transformaciones de señal aplicadas en la etapa 103 se realizan en el dominio del tiempo, éstas también se pueden asociar con nociones del dominio de la frecuencia. Tal como se desprende de las Figuras 6A y 7A, la ventana estrecha usada en MPMA 1 realiza un filtrado de paso bajo con un ancho de banda de frecuencia más alto, mientras que la ventana ancha usada en MPMA 2 realiza un filtrado de paso bajo con un ancho de banda de frecuencia relativamente baja. Sin embargo, dado que la señal original es un proceso no estacionario, el empleo de bandas de frecuencia convencionales con frecuencias de "corte" fijas es ambiguo.

Con referencia nuevamente a la Figura 3, después de la segmentación en la etapa 103, el método 100 realiza las etapas 104 a 108 como una recursión, cuyos ciclos consecutivos tratan la sucesión de segmentos de señal definidos. Para cada segmento se identifica un modelo de componente, cuyos parámetros se coinciden con el fragmento de señal correspondiente. Dado que el modelo de componente respectivo emparejado con la señal en ese segmento se extiende a los siguientes segmentos, el objetivo de las etapas de recursión desde 104 a 108 es también resolver la superposición de componentes temporal. La resolución de la superposición de componentes temporal es un problema urgente en muchos campos del procesamiento de señales que no se soporta por métodos comunes. Una gran dificultad es que cualquier enfoque de este problema debe involucrar la creación del modelo de componentes para describir estructuras predecibles y patrones esperados en el desarrollo del componente, y pronosticar los valores futuros del componente a partir de su trayectoria pasada.

El inconveniente de un enfoque "basado en la física" (verCandy J.V. "Model-based signal processing". John Wiley & Sons, Nueva Jersey, 2006) radica en el hecho de que generalmente no se dispone de suficiente información con relación a la señal en términos de modelos físicos explícitos de sus fuentes. Por el contrario de los enfoques basados en la física, el método divulgado en la presente memoria introduce una metodología única de modelado basado en el caos (CBM) que define los componentes principales de una señal no estacionaria en términos de caos determinista. Este enfoque de modelado basado en el caos se ha motivado por el hallazgo empírico del efecto de normalización (verMelkonian D, Blumenthal T, Gordon E. "Numerical Fourier transform spectroscopy of EMG half-waves: fragmentary-'decomposition-based approach to nonstationary signal analysis". Biol Cybem., 1999, volumen 81, páginas 457-467). El hallazgo indica que debido a una interferencia de factores deterministas y aleatorios la actividad de la fuente de la escala microscópica se convierte a escala global en medidas estadísticas que reflejan regularidades probabilísticas en lugar de una naturaleza física de las fuentes elementales. Por lo tanto, un componente fisiológico de escala global se desarrolla como un agregado estadístico acumulativo de múltiples fuentes elementales.

Para formular un modelo de componentes de escala global universal sobre esta base heurística, el método se refiere al teorema del límite central como una herramienta que define el comportamiento limitante de conjuntos de variables aleatorias (verGnedenko B, Kolmogorov A. "Limit distributions for sums of independent random variables". Addison-Wesley Publishing Company, Reading, Mass., 1954). Dada la suma de un cierto número de variables aleatorias f-if2,...,fn con medios n1,n2,.,fn y variaciones

_C "_Jj 1 _*C?2 1 _!*••!« , _{„ s}

el comportamiento limitante de sus transformadas exponenciales de Fourier F1,(/w),F2(iw),...Fn(/w) convergen en condiciones bastante generales a una función Gaussiana como n+ « (verPapoulis A. "The Fourier integral and its applications". McGraw-Hill Book Co, Nueva York, 1962)

donde F(io>) es una transformada exponencial de Fourier (función de variable compleja

),

n = ni + n ² +...+ n y

Por lo tanto, los efectos de fuentes únicas se transforman en parámetros de resumenoynque tienen significado como variables de escala global.

En términos físicos,

se puede considerar como el espectro de amplitud de un proceso temporal de escala global producido por la multiplicidad de fuentes de escala microscópica.

El soporte empírico de esta predicción teórica se ha proporcionado por el estudio de los potenciales relacionados con eventos auditivos (verMelkonian D, Gordon E, Rennie C, Bahramali H. Dynamic spectral analysis of event-related potentials. Electroencephalography and Clinical Neurophysiology, 1998, volumen 108, páginas 251-259). El análisis espectral dinámico reveló que una función Gaussiana proporciona una descripción universal de los perfiles de los espectros de amplitud de los principales componentes de latencia tardía de los potenciales relacionados con el evento auditivo.

La contraparte en el dominio del tiempo del espectro complejo de dominio de la frecuencia de la Ecuación (7) se define en la literatura relevante en un intervalo de escala de tiempo infinito. De manera similar a la Ecuación (8), la solución en el dominio del tiempo también tiende a una función Gaussiana de la siguiente forma:

Por el contrario, las rutinas de procesamiento de señales del método divulgado tratan con funciones de tiempo causales que se pueden considerar transitorias en sistemas físicamente realizables.

Para reubicar el análisis en el dominio del tiempo de una escala de tiempo infinita a semi-infinita, el método divulgado introduce una función causal w(t) que se define enf>0 como la suma de dos funciones Gaussianas

donde

)

en el que U(t) es una función de etapa unitaria, y e es una constante adimensional.

Al comparar la Ecuación (7) con las Ecuaciones (11) y (12) se muestra que el parámetro^(medio) de la Ecuación (7) se expresa en las Ecuaciones (11) y (12) como la fracción del parámetroo(dispersión) mediante el uso de la constantes.. El valor de la constantesdepende del tipo de señal no estacionaria. Por lo tanto, el análisis y las simulaciones por ordenador de EEG de humanos y animales mediante el uso del método divulgado sugieren s=2 como constante universal para estas aplicaciones. Esta estimación se considera como valor predeterminado de la constante e.

Las funciones gA(t) y - gB(t) con el parámetro o=1 y la constante s=2 se muestran en la Figura 8A mediante el guión y las líneas de puntos, respectivamente. Los extremos de estas funciones están en t = so y t = -so.

Como se muestra mediante la línea continua en la Figura 8A, la función g (t) es una función impar de la extensión infinita. La transformada de Fourier de la función g(/ es:

donde

m

Gc(m)= Jg(Ocos(fi*)ífr (15)

Ya que g(t) es una función impar, Ge(w) = 0 y

Gs(fti)= 2j g(*)sin(íHí)f/í=2j w(f)sin(¿af)ííf. (17)

o o

La última integral trigonométrica de la Ecuación (17) es la transformada sinusoidal de Fourier de w(t) para la cual la denotación

se usa.

Después de reemplazar w(t) por la Ecuación (10), esta integral se expresa analíticamente en la siguiente forma (verErdelyi, A. (Ed). "Tables of integral transforms", Volumen 1. McGraw-Hill, Nueva York, 1954, página 78):

M 2

Ws(o)=e 2 sín^ oT®) (19)

Esta expresión se puede considerar como la parte imaginaria del espectro complejo.

w (¡« > )= w í< » )r ', w (²⁰ )

donde

H 1

W(ffl)=e ^ (21)

es el espectro de amplitud y

es el espectro de fase.

Por lo tanto, el dúo de las funciones de frecuencia de las Ecuaciones (21) y (22) se puede considerar como las contrapartes en el dominio de la frecuencia que definen completamente la señal en el dominio del tiempo de w(t) en t>0. Esta función se muestra en la Figura 8B. Un valor aproximado en el que w(t) alcanza el máximo es t~2. El valor del máximo w(2)«0,399.

Físicamente, w(t) se puede considerar como un proceso transitorio inducido en t=0 mediante la activación de un conjunto de múltiples fuentes de escala microscópica. El comportamiento de fuentes separadas se rige por factores tanto probabilísticos como deterministas y se puede caracterizar como el caos determinista. Un rasgo característico principal de dicha forma de comportamiento es la no linealidad de los sistemas que producen el caos determinista. En este contexto, el proceso descrito por la Ecuación (11) se puede considerar como una respuesta de impulso de un sistema dinámico regido por el siguiente sistema de ecuaciones diferenciales no lineales:

donde g&A(t) y g&B(t) denotan las derivadas de gA(t) y gB(t).

El modelo de un proceso compuesto no estacionario se define en la forma:

donde k¡ y t son los parámetros libres denominados peso y tiempo de inicio.

Un modelo de componentes, denominado núcleo cuasi-Gaussiano (QGK), se define como:

que se caracteriza por un conjunto de tres parámetros, específicamente ( k ¡, Oí, r¡ ) .

En un contexto general, la estocasticidad de vn(Z) se deriva del hecho de que los parámetros (k , oí, t) son variables aleatorias. Un aspecto específico de la no estacionariedad es que la composición estructural del sistema que genera v(t) cambia en los instantes de tiempo t cuando las nuevas fuentes se conectan al sistema subyacente. Un aspecto crucial del caos determinista se expone por hecho de que aunque las dependencias Gaussianas de las Ecuaciones (10)-(12) son procesos deterministas, su evolución temporal se desarrolla como el comportamiento limitante de los sucesos aleatorios elementales subyacentes. Los aspectos probabilísticos de las transformaciones relevantes en la escala microscópica se han formulado en términos de procesos de Markov, específicamente procesos de nacimiento y muerte no homogéneos (verMelkonian D, Blumenthal T. Dynamics of the eyeblink EMG at global and microscopic scales. The Open Cybernetics and Systemics Journal, 2013, volumen 7, páginas.11-22).

El método 100 realiza las etapas de la 104 a la 108 para identificar los parámetros (k , oí, t í) de cada QGK Yn ( t) de la Ecuación (24). El algoritmo de identificación se organiza como una recursión cuyo principio rector proviene de la siguiente relación entre los términos "n" y "n-1 " de la Ecuación (24):

V»W = ¿ í= Ti i (0 = V ⁱM+Y*W (26)

Dado que la función vn(t) es conceptualmente un modelo, esta función vn(t) se puede considerar como una aproximación a la secuencia sin tendencia u(t), es decir, vn(t) = u(t) en t> 0. La contribución del QGKYn(t) al aproximado comienza desde el inicio de Yn(t) en t = m. Esto significa que

« ( í)a v H (í) para 0 < t iL

UW = ^ ( 0 ^ ( 0 para

Dado 71=0,

r t( í)= u (í) para

(28}

Para n>1,

T .W M O -V iM Para

(29)

La Ecuación (28) corresponde al primer ciclo de recursión en el que el impacto de los términos anteriores está ausente. La Ecuación (29) corresponde al ciclo nth de recursión donde el impacto de los términos anteriores se contabiliza por la función vn^-1(t). Esta función se define analíticamente por los parámetros principales (kí,o¡,t) (í=1,..,n-1) obtenido en los ciclos anteriores de recursión. El papel de vn-1(t) es la corrección de la superposición del componente temporal. En este contexto, la parte derecha de la Ecuación (29) representa una contraparte semi-empírica corregida por superposición del QGK analítico Yn(t).

Los parámetros de QGK Yn(t) se interrelacionan en el dominio del tiempo por complejas relaciones no lineales. La ventaja del dominio de la frecuencia es que las relaciones en el dominio de la frecuencia proporcionan los medios para separar estas relaciones y estimar la sucesión de Kn, On y Tn parámetros por procedimientos separados realizados en un orden consecutivo. Un aspecto importante de las correspondientes manipulaciones en el dominio de la frecuencia es que las estimaciones de Kn y Onse definen por los perfiles de los espectros de amplitud que son invariantes a los cambios de tiempo de la función de tiempo bajo el análisis. Esto proporciona los medios para usar el valor del punto de segmentación ^n-1como una estimación provisional (durante las etapas de la estimación Kn y On) del parámetroTn. En consecuencia, las Ecuaciones (28) y (29) se presentan a los fines de las transformaciones de tiempo a frecuencia en las siguientes formas:

Para respaldar las transformaciones del dominio del tiempo a frecuencia mediante procedimientos estándar, el elemento de forma de onda de la Ecuación (31) se cambia porXn^-1 y se describe como una sola función:

en el intervalo [0, Án], donde Án = Xn - Xn^-1. Esta función se denomina núcleo de componente empírico (ECK).

En la etapa 105 el ECKyn(t) se transforma del dominio del tiempo al dominio de la frecuencia mediante el uso de la estimación numérica de las transformaciones del seno y coseno finitos de Fourier. Más particularmente, se usa el algoritmo de función de base similar (SBF) divulgado en Melkonian (2010) (verMelkonian D. "Similar basis function algorithm for numerical estimation of Fourier integrals". Numer Algor., 2010, volumen 54, páginas 73-100).

El algoritmo SBF es una versión original de los métodos tipo Filon que proporcionan una precisión máxima en la estimación numérica de integrales trigonométricas mediante el uso de la interpolación polinomial. Por el contrario del análisis espectral digital convencional realizado por la Transformada Rápida de Fourier (FFT) y sus modificaciones, el algoritmo SBF está libre de la necesidad de lidiar con cuadrículas estrictas de puntos computacionales tanto en los dominios de tiempo como en los de frecuencia. Lo mismo es particularmente ventajoso en el presente caso donde las longitudes de los segmentos de señal sometidos a la transformación de tiempo a frecuencia en la etapa 105 son variables. En el dominio de la frecuencia, la principal ventaja de la elección flexible de los puntos computacionales es la posibilidad de expresar las dependencias del dominio de la frecuencia mediante el uso de una escala de frecuencia logarítmica.

Dado que el algoritmo SBF se realiza mediante procedimientos universales, el subíndice "n" (número del ciclo en la recursión) se omite en las siguientes descripciones del algoritmo.

El algoritmo SBF se dirige a la estimación de las siguientes transformaciones de seno y coseno finitos de Fourier:

usadas en estas ecuaciones suponen que y(t) = 0 en t > Á. Estas son las partes reales e imaginarias del espectro complejo de y(t)

donde

El espectro de amplitud correspondiente es:

La transcripción de y(t) a una forma digital especifica y(t) por sus valores muestreados y, = y (ti) en un número finito de puntos nodales ti (i = 0,1,..,N) con t0 = 0 y tN = Á. No es necesario espaciar los nodos igualmente. Por lo tanto, los cálculos se dirigen al conjunto de puntos de datos {(y0,t0),...,(y/,ti),...,(yw,tw)}. Mediante el uso de estos puntos de datos, y(t) se interpola sobre sus muestras mediante un polinomio lineal por partes h(t). El interpolante debe satisfacer la siguiente condición de interpolación:

donde h,= h(t,).

Suponiendo que y, = 0 para i > N, el algoritmo SBF crea para el intervalo [0, A] el interpolante como la suma de funciones de base similares:

donde a,son los coeficientes de ponderación, y 0¡(f) es una función de base similar (SBF). La SBF se define por la relación de similitud:

Esta escala de tiempo parametrizada simple produce la familia de SBF a partir de un elemento finito básico

fl- / si

1 ⁽ 0 ⁼ [ 0 de c

denominada "función de base triangular" (TBF). La TBF es un triángulo rectángulo unitario representado en la Figura 9C.

Mediante el uso de elementos finitos, el interpolante h(t) se presenta en la forma:

en donde 0,(t) es una función de sombrero. La función de sombrero $/(t) se define como:

en la malla t⁰ < ti <... < ti <... < tN. La Figura 9A ejemplifica una descomposición del interpolante h(t) como la suma del TBF adjunto a la muestra "0" seguida de la sucesión (muestras 1, 2, 3, 4) de las funciones de sombrero $;(t) de la manera expresada en la Ecuación (41).

Al comparar las Ecuaciones (37) y (41) se muestra que las funciones de sombrero de la Ecuación (41) se reemplazan en la Ecuación (37) por TBF. Los principios geométricos que resultan útiles en relación con esto se ilustran en la Figura 9B que muestra la descomposición de la función de sombrero (triángulo abc) en la suma de tres TBF (triángulos ohc, fod y goa). La aplicación de esta descomposición a todas las funciones de sombrero en la Ecuación (37) da lugar a la siguiente fórmula para la estimación de los coeficientes de interpolación:

donde

En el contexto de las soluciones numéricas, la principal ventaja de la construcción de la Ecuación (42) es que las integrales de seno y coseno de Fourier del TBF se expresan en una forma analítica como

Rs(®) = 3 sM4l} = Jr{f)siiifljtáf = — ™®. (46)

La Figura 9D muestra estas funciones denotadas como Re (Ecuación (45)) e Im (Ecuación (46)).

De acuerdo con la teoría de las transformadas de Fourier (teorema de la similitud), la compresión de la abscisa en el dominio del tiempo corresponde a la expansión de la abscisa más la contracción de la ordenada en el dominio de la frecuencia. Estas operaciones establecen las siguientes relaciones de similitud entre las integrales de Fourier de SBF y TBF:

Mediante el uso de estas ecuaciones, las integrales de coseno y seno del interpolante de la Ecuación (38) se expresan mediante las siguientes ecuaciones analíticas:

Dado que h(t) es un interpolante a y(t), estas transformadas de Fourier se consideran las estimaciones de las integrales finitas de Fourier de las Ecuaciones (33) y (34). Por consiguiente,

En este contexto, el reemplazo de las funciones Rc y Rs en las Ecuaciones (49) y (50) por sus expresiones analíticas (45) y (46) conduce a las siguientes ecuaciones computacionales:

Después de la estimación de los coeficientes de ponderación de la Ecuación (43), las Ecuaciones (52) y (53) se aplican directamente a los puntos de datos.

La ventaja del algoritmo SBF es que las características de frecuencia se pueden calcular en los conjuntos arbitrarios de las frecuencias angulares. Los siguientes procedimientos de identificación de parámetros emplean la escala de frecuencia logarítmica. En consecuencia, las muestras Yc(w,) y Ys(w,) de las características de frecuencia se estiman en los puntos

donde C>1 es una constante que define la velocidad de muestreo en una escala de abscisas logarítmica, ywAes la frecuencia inicial.

Los valores típicos de C que corresponden a 10, 20, 50 y 100 muestras por década son 1,2589, 1,122, 1,0471 y 1,0233, respectivamente.

El espectro de amplitud se calcula en los mismos puntos de la siguiente manera:

Dado el ciclo nth de recursión, el espectro de amplitud correspondiente Yn(w) sirve para la identificación de Kn (peso) y an (dispersión) parámetros. Para este propósito, el método 100 (Figura 3) después procede a la etapa 106.

Con la condición de que la frecuencia angular wA se selecciona de tal manera que Yn(wA) « Yn(0), el peso se define como

Con respecto a esto son decisivas las pruebas de adecuación del espectro de amplitud Gaussiana de la Ecuación (21) como modelo del espectro de amplitud empírico. Las pruebas se aplican al espectro de amplitud normalizado calculado a partir de la Ecuación (55), y usan la estimación del parámetro Kn de la Ecuación (56):

La comparación de este espectro con el modelo de la Ecuación (21) se realiza en una escala de frecuencia logarítmica mediante el uso del espectro Gaussiano como una plantilla. Este último se define analíticamente como

donde v es una variable adimensional real. Esta función se muestra en la Figura 10A. Para cálculos numéricos, la plantilla se define por sus valores muestreados en los puntos,

donde C>1 es la misma constante que en la Ecuación (54), y M+1 es el número de muestras de la plantilla.

Por lo tanto, la plantilla se extiende sobre el segmento [va, vb], donde v b = v aCm. Los puntos límite v a y v b se indican en la Figura 10A mediante las flechas A y B.

Las muestras del espectro de amplitud normalizado empírico con el que se compara la plantilla corresponden a los siguientes puntos en la escala de frecuencia angular logarítmica:

El intervalo de frecuencia (desde wAa wi = waC1) y la velocidad de muestreo se definen por los parámetros lyC. Una elección adecuada de estos parámetros proporciona los medios para una representación suficientemente precisa del espectro de amplitud en una forma digital. El intervalo de frecuencia del espectro de amplitud debe exceder el intervalo de frecuencia sobre el que se define la plantilla. En consecuencia, I>M.

La coincidencia de plantillas se realiza al mover la plantilla de amplitud Gaussiana a lo largo de la escala de abscisas. Este es un proceso de etapa a etapa durante el cual las muestras de la plantilla permanecen sin cambios mientras sus ubicaciones en la escala de abscisas se desplazan conjuntamente al multiplicar las frecuencias por la constante C. Por lo tanto, el primer desplazamiento se realiza mediante el cambio de la secuencia de puntos de la Ecuación (57) a la secuencia de puntos vm = vACm+1 (m = 0,1,...,M). La secuencia de los puntos de frecuencia después del desplazamiento kth es vm = vACm+k (m = 0,1,...,M). La Figura 10B ejemplifica el espectro de amplitud normalizado empírico ñn(w) (línea negra continua) y las plantillas de amplitud Gaussiana establecidas en diferentes desplazamientos como los gráficos T1, t 2 y T3.

La discrepancia entre la plantilla y el espectro de amplitud empírica se evalúa en cada posición de la plantilla mediante el uso del error cuadrático medio (MSE):

donde k es el número de la posición (k=1,...,I-M+1).

En términos de la frecuencia angular, la posición inicial de la plantilla cubre el intervalo de frecuencia [wa,wb] del espectro empírico, donde wb = waCm. Después k cambia el intervalo de frecuencia cubierto por la plantilla es

[o i,0 |J donde

y í l | = mñC k,

La mejor coincidencia de plantilla se indica mediante el valor mínimo del MSE. Dado K como el número del desplazamiento correspondiente, los parámetros I p M ]

definen en la escala de frecuencia angular logarítmica la posición de la plantilla en la que el MSE es mínimo. Para decidir si esta coincidencia fundamenta el espectro de amplitud de la Ecuación (21) como un modelo adecuado del espectro de amplitud empírico de la Ecuación (55), el residual entre estas funciones se considera una función de error adimensional £ (w) definido por sus muestras

donde j tiene valores desde 0 a M, ® j =Q AfCi ,

Aquí, ambos Yn(w) y T(v) son funciones adimensionales con valores iniciales iguales Yn(0) = T(0) = 1. En consecuencia, se presume que e(0)=0. El cambio de la función de error con el aumento de la frecuencia se evalúa mediante el uso del valor umbral £Tpara estimar el intervalo de frecuencia [0,wJ]sobre el cual

El valor final de Yn(w) en la frecuencia angular límite Wjse denota como Yj= Yn(wj). Esta estimación se usa para validar la calidad de la coincidencia de la plantilla mediante el uso del parámetro de categorización ^que toma valores enteros de 1 a 5 de acuerdo con el siguiente esquema:

4 “ 1 (grado E) si 0 £ Y ,£ 0 ,2 ,

\ ^{= 2 (grado G}} _si ^{0,2 < Y , S 0,4,}

£ = 3 (grado S) _si ^{0,4 <} % ^áQ,6, (64) £ = 4 (grado W) si 0,6 < Yj í 0,8,

^{£ = 5 (grado U) si} % ^{> 0,8.}

Dado el grado, es decir, el valor de ,^ el resultado del modelo se califica de la siguiente manera: E - excelente, G -bueno, S - satisfactorio, W - regular, U - no satisfactorio.

En función de la directiva, en el caso del grado W o U, se ignora el peso, es decir Kn=0.

Las Figuras 11A y 11B ilustran las funciones de error para dos ejemplos de coincidencias de plantilla mediante el uso del umbral sT= 0,05. Las frecuencias límite correspondientes se indican en la Figura 11A y la Figura 11B mediante líneas de puntos verticales. Estas líneas cruzan los espectros de amplitud en la Figura 11Ay la Figura 11B en los valores 0,065 y 0,301, respectivamente. En consecuencia, los modelos correspondientes pertenecen a las clases E y G, respectivamente.

La mejor coincidencia define alguna posición fija de la plantilla en la que las muestras de la plantilla se pueden asociar con entidades físicas: los valores de la frecuencia angular en las muestras correspondientes de Yn(oj) en la escala de abscisas logarítmica. En este contexto, las muestras de la plantilla se expresan como funciones de la frecuencia angular en la siguiente forma

donde el subíndice "F" indica la posición final de la plantilla, andenota el parámetro a de la Ecuación (21) en el ciclo de recursión nésima.

Dado la arbitraria wm, el parámetro requerido an se calcula como

Para una precisión estable de los cálculos, la frecuencia angular se selecciona del intervalo de frecuencias sobre el que Tp(wm) toma valores de 0,4 a 0,8.

La Figura 12 ejemplifica los resultados típicos de la coincidencia de la plantilla. El panel central muestra una pieza de EEG con los puntos de segmentación indicados por las líneas de puntos verticales. Los espectros de amplitud de 6 fragmentos de EEG consecutivos tomados entre los puntos de segmentación se muestran mediante líneas continuas en los paneles superior e inferior. Las plantillas correspondientes se representan por líneas grises. Cada ejemplo de coincidencia de plantilla se complementa por la siguiente información: el grado del modelo (E o B), la frecuencia Angular

y el valor

Para la misma señal redibujada en la Figura 13A, los espectros de amplitud de los fragmentos de señal entre los puntos de segmentación se muestran en la Figura 13B. La Figura 13C muestra la superposición de estas características espectrales después de la normalización tanto de la amplitud como de la frecuencia. Estos gráficos son típicos en el sentido de que la concordancia con la plantilla de amplitud Gaussiana (línea negra) es excelente para los intervalos de frecuencia baja e intermedia. Las simulaciones por ordenador han demostrado que las desviaciones de alta frecuencia se producen por actividades irrelevantes mezcladas con los componentes funcionales.

Al hacer referencia nuevamente al método 100 (Figura 3), al seguir la etapa 106, el método 100 procede a la etapa 107 donde se estima el tiempo de inicio Tn. Este procedimiento usa el espectro de amplitud normalizado Yn(w) y la parte imaginaria correspondiente del espectro complejo después de su normalización de una manera similar a la normalización del espectro de amplitud mediante el uso de la Ecuación (55):

donde el subíndice n en el número del ciclo de recursión.

Las Figuras 14A y 14B ejemplifican una forma típica de estas funciones denominadas AS (espectro de amplitud normalizado Yn(w)) e Im (parte imaginaria normalizada del espectro complejo Ysn(iw)).. La relación general entre estas funciones incluye el espectro de fase ^ (w) en términos de los cuales

De manera similar al espectro de fase de la Ecuación (22), ^n(w) se expresa como dependencia lineal en la frecuencia angular en la forma

_{V » = A < a} (69)

La diferencia entre los espectros de fase de las Ecuaciones (22) y (69) se explica por el coeficiente de desplazamiento

' . “ A " * » (70)

Físicamente, este cambio indica la diferencia entre el tiempo de inicio del componente Tn y el valor del punto de segmentación Xn-¹. En este contexto,

Los procedimientos del algoritmo SBF en los ciclos previos de recursión toman el punto de segmentación como un sustituto del tiempo de inicio desconocido del ECK bajo análisis. En consecuencia, el punto de segmentación sirve como el instante de tiempo a partir del cual comienza la integración numérica. Esta suposición tiene un efecto insignificante en los valores de los parámetros previamente definidos porque el espectro de amplitud, la fuente de los parámetros Knyon, es invariante a los cambios de tiempo de la función de tiempo de contraparte.

La estimación Tn se refiere a la Ecuación (71) de donde se deduce que

Para usar esta ecuación, el parámetro desconocido fin se debe estimar. La referencia a las ecuaciones (68) y (69) indica que el parámetro fin se define completamente por la siguiente relación entre las funciones Ysn(w) y Yn(w) como: Para encontrar el parámetro in a partir de esta ecuación, es suficiente tener estimaciones de Ysn(w) y Yn(w) para un solo punto de frecuencia angular. Sin embargo, para reducir los errores de cálculo, las dos estimaciones se calculan en diferentes frecuencias angulares w1n y w2n. Dado w2n > wm, la elección de estas frecuencias supone que en el intervalo desde w 1n hasta w2n Ysn(w) disminuye con el aumento de w.

Las líneas de puntos verticales en las Figuras 14A y 14B indicadas por la letra a en la escala de abscisas, definen la frecuencia angular w 1n como el punto en el que Ysn(wm) = Yn(wm). La referencia a la Ecuación (73) muestra que en este punto fin = n/2w1n. La inserción de este valor en la Ecuación (72) proporciona la primera estimación del tiempo de inicio:

De acuerdo con los modelos en el dominio de la frecuencia en la forma de las Ecuaciones (21) y (22), la forma esperada de Ysn(w) se caracteriza por un cruce por cero seguido de un mínimo local como se ilustra en la Figura 14A. Dada dicha situación, la frecuencia angular w2nse define como el punto donde Ysn(w) cruza la escala de abscisas, es decir, Ysn(W2n) = 0.

En la Figura 14A w2n se indica por la línea de puntos vertical que se origina desde b en la escala de abscisas. El valor del espectro de fase ^n(w2n) = n define i = n/w2n.

La Ecuación (72) da la segunda estimación del cambio de fase como

En algunos casos los elementos irrelevantes del ECK distorsionan el perfil de dominio de frecuencia esperado de Ysn(w). El gráfico Im en la Figura 14B es un ejemplo característico de Ysn(w) que alcanza el mínimo en lugar de cruzar la escala de abscisas. La frecuencia angular correspondiente se indica por la línea de puntos vertical que se origina desde b en la escala de abscisas. Esta frecuencia reemplaza W²n en la Ecuación (75) para obtener la segunda estimación del tiempo de inicio.

Dado Tn[1]y Tn[2], la elección final del tiempo de inicio es en base a en la evaluación de diferentes valores de Tn las discrepancias entre Ysn(w) y Ws(w) de la Ecuación (19) con o=on. La discrepancia se evalúa por el error cuadrático medio

donde las muestras Ysn(wi) se calculan mediante el uso de las estimaciones Tn[1] o Tn[2] de Tn, y N es un entero seleccionado de tal manera que wn = w 1n.

Dado que MSE[1]y MSE[2] son relevantes para Tn[1] y Tn[2], Tn =Tn[t] si MSE[1]< MSE[2]. De cualquier otra manera Tn = Tn[2].

Por tanto, las etapas 106 y 107 estiman los parámetros (Kn,On,Tn) para el segmento en consideración. La etapa 108 estima el modelo parcial

mediante el uso de la Ecuación (26) antes del método 100 después regresa a la etapa 104 desde donde las etapas 104 a 108 estiman los parámetros (Kn+i,on+i,Tn+i) del siguiente segmento.

La Figura 2B ilustra los resultados típicos de sustituir cada componente funcional mostrado en la Figura 2A (N100, p200, N200, P3a y P3b) con el QGK que tiene parámetros ( k,o ,t ) que coincida con el ECK correspondiente. Los segmentos se indican con los mismos símbolos N100, P200, N200, P3a y P3b que sus contrapartes empíricas en la Figura 2A. La Figura 2C ilustra la suma de los componentes QGK, es decir, un modelo global, así como también la señal ERP original que se muestra en la Figura 2A para comparación. Este es un grado de concordancia satisfactorio, ya que los parámetros de la descomposición de la señal y los componentes identificados se derivaron completamente de los procedimientos generales del método inventivo, sin ningún ajuste para adecuarlos a los fenómenos a los que se aplicaron subsecuentemente.

Se ha obtenido más evidencia a favor de esta conclusión mediante el uso de señales no estacionarias producidas por diferentes fuentes. Para ejemplificar este punto, la Figura 15A muestra un segmento de 60 segundos de la forma de onda respiratoria de una persona al hablar muestreada en A=0,056 s. El gráfico representa la señal sin tendencia del sensor de respiración mediante el uso de 5 pasadas para filtros de ventana de 3 y 47 puntos. Las unidades del eje y son las salidas del sensor de respiración sin procesar en bits.

La Figura 15B representa la superposición de 54 QGK identificadas. En la Figura 15C se muestra un modelo global de la forma de onda respiratoria sintetizada como la suma de todos los componentes del modelo mediante la línea continua. En el mismo gráfico, la línea de puntos es la forma de onda respiratoria original reproducida de la Figura 15A.

De manera similar al EEG, las formas de onda respiratorias son señales altamente no estacionarias. Sin embargo, por el contrario del EEG producido por las fuentes de electricidad (neuronas corticales), los movimientos respiratorios dependen de una serie de factores biomecánicos que rigen los movimientos de la pared torácica. Por lo tanto, los procesos físicos que subyacen a los movimientos respiratorios son bastante diferentes de la electrodinámica de las ondas cerebrales. Sin embargo, parece que el método descrito en la presente descripción proporciona herramientas de identificación de componentes igualmente efectivas para ambas aplicaciones. Los ajustes notablemente precisos de QCjK a formas de onda empíricas son típicos para una variedad de señales electrofisiológicas (EEG, Em G, ECG) y formas de onda respiratorias registradas de diferentes especies y bajo diferentes condiciones experimentales. Estos resultados apoyan fuertemente el enfoque de modelado basado en el caos (CBM) de la metodología inventiva.

Una vez procesados todos los segmentos y los parámetros £,k,o ,t se han estimado para cada uno de los segmentos formados en la etapa 103, el método 100 continúa hasta la etapa 109 en la que se examinan los resultados del modelado mediante el uso del parámetro £ como marcador. El valor del parámetro £ para cada segmento se compara con un valor de umbral definido £ t. Si £n s £ t, el modelo del componente correspondiente se elimina de las etapas de procesamiento posteriores. Los valores típicos de £ t son 4 y 5.

Después de eliminar los componentes inapropiados N i , los componentes restantes se vuelven a numerar en orden consecutivo desde n=1 hasta M =Ns-Ni, donde Ns-1 es el número de puntos de segmentación. El conjunto de los parámetros correspondientes se indica mediante:

$RM = {£*lcrfT} » (78)

donde k={k-i ,...,k/,...,km}, o={o i,...,ct,...,om}, t={ti,...,t/,...,tm}, son vectores, cada uno de los cuales contiene valores M del parámetro correspondiente.

Los elementos de los vectores k y o definen las formas de onda de los componentes del modelo mientras que los elementos del vector t definen un conjunto finito de puntos de tiempo a partir de los cuales se inician los desarrollos de los componentes definidos del modelo. En este contexto, el vector Kt = {k, t} proporciona el conjunto de puntos de datos {(k i,ti),...,(k/,7/),...,(km,tm)}, que representan la secuencia de los valores del parámetro k como la serie de tiempo. Igualmente, el vector aT = {o ,t} proporciona la serie de tiempo de los valores del parámetro o.

Las series de tiempo son la forma estándar de los datos numéricos a los que se adaptan numerosos métodos de procesamiento de señales. Por lo tanto, el método divulgado 100 que descompone una señal en conjuntos de parámetros estándar proporciona los medios para aplicar al procesamiento de señales no estacionarias métodos tradicionales del análisis de las series de tiempo tales como estandarización y transformación de datos, análisis de componentes principales, autocorrelación, modelos autorregresivos, agrupación de datos, pronóstico, etc.

Mediante el uso de los resultados de la descomposición de la señal no estacionaria en los conjuntos de los modelos de componentes (análisis de la señal), el método 100 ofrece una herramienta para resolver un problema opuesto, que es la síntesis de señal. Esto implica crear compuestos de componentes que se asignan a clases predefinidas.

La recopilación de los componentes para la síntesis comienza desde el procedimiento de clasificación de los parámetros realizado por la etapa 110. El tiempo de inicio se mantiene sin cambios, y la clasificación se aplica a los parámetros £,k y o para cada modelo de componente identificado. Más particularmente, un subconjunto de parámetros se determina a partir del conjunto completo de parámetros del modelo PM de acuerdo con las siguientes condiciones:

en donde ^i,ki,gi, y £f,kf,of son constantes que definen los valores límite (inicial y final) del parámetro correspondiente. De la Ecuación (78), el subconjunto de parámetros se denota

4 = { £ * , < * 4 (80)

donde un overhat (A) denota los elementos seleccionados, y L es el número de componentes. Obviamente, L<M. Finalmente, el método 100 termina en la etapa 110 donde el subconjunto de parámetros

de la etapa 110 se usa para la síntesis de un modelo parcial de la señal no estacionaria en la forma

La clasificación y la síntesis son procedimientos multipropósitos que pueden perseguir diferentes objetivos. Por ejemplo, la exclusión de los componentes con amplitudes enormemente pequeñas o grandes puede eliminar los elementos irrelevantes de la señal. La selección de los valores específicos de k y o puede enfatizar los componentes oscilatorios, etc.

El conjunto completo de los parámetros del modelo

y/o el conjunto de los parámetros seleccionados

se puede almacenar en un medio de almacenamiento del sistema informático, o comunicar a un sistema informático externo, antes del procesamiento adicional. El procesamiento adicional puede incluir diferentes métodos de procesamiento de señales basados en modelos, tales como la deconvolución óptima, los algoritmos de agrupamiento heurístico, la estimación Bayesiana recursiva, etc.

Claims (7)

1. REIVINDICACIONES

   i. Un método implementado por ordenador (100) para descomponer una señal no estacionaria representativa de un fenómeno fisiológico en un conjunto de componentes funcionales autosimilares con el QGK (núcleo cuasi-Gaussiano) como elemento base, y para estimar los parámetros que caracterizan los respectivos componentes funcionales, el método que comprende las etapas de:

   aplicar filtros de media móvil de múltiples pasadas para estimar las tendencias dinámicas y de línea de base en el transcurso del tiempo de la señal (etapa 102), usar las tendencias para extraer una función de guía indicativa de los fragmentos de la señal sobre los que se desarrollan los componentes funcionales, y segmentar de forma adaptativa la función de guía mediante el uso de cruces por cero (ZC) como los puntos de segmentación;

   dividir la señal no estacionaria en las regiones de desarrollo de componentes, al segmentar la señal no estacionaria en puntos correspondientes a los puntos de segmentación de la función de guía (etapa 103); compensar la superposición de los QGK que coinciden con los componentes funcionales anteriores para formar un ECK (núcleo de componente empírico) (etapa 104), las subetapas recursivas de:

   transformar el ECK en un dominio de la frecuencia mediante el uso de un algoritmo de función de base similar (SBF) (etapa 105);

   estimar los parámetros k (peso) y a (dispersión) (etapa 106) y validar la fiabilidad del parámetro en términos de grados categóricos E, G, S ,W y U, en donde: E - excelente, G - bueno, S - satisfactorio, W - regular y U - no satisfactorio al encontrar y evaluar la mejor coincidencia entre el espectro de amplitud normalizado del ECK y la plantilla del dominio de la frecuencia con un perfil Gaussiano; estimar el tiempo de inicio t del QGK con los parámetros identificados k y a de la parte imaginaria del espectro complejo del ECK al encontrar puntos característicos que corresponden a los valores definidos del espectro de fase (etapa 107);

   expandir el modelo de la señal no estacionaria compuesta a partir de la suma de los QGK identificados en los ciclos anteriores de recursión al añadir el QGK recién identificado (etapa 108); y reorganizar los componentes identificados después de eliminar los QGK desfavorables mediante el uso de un parámetro categórico 5 para evaluar los resultados de la identificación (etapa 109); crear modelos parciales de la señal no estacionaria mediante la selección de conjuntos QGK cuyos parámetros pertenecen a clases predefinidas (etapa 110).
2. 2. Un método (100) de acuerdo con la reivindicación 1, en donde el parámetro de categorización 5 toma un valor entero de 1 a 5 de acuerdo con lo siguiente:

   ^= 1 (grado E) si 0 < Yj < 0,2,

   5 = 2 (grado G) si 0,2 <Yj < 0,4,

   5 = 3 (grado S) si 0,4 <Yj < 0,6,

   5 = 4 (grado W) si 0.6 < V j <0,8,

   5 = 5 (grado U) si Yj > 0,8,

   en donde: E - excelente, G - bueno, S - satisfactorio, W - regular y U - no satisfactorio.
3. 3. Un método (100) de acuerdo con la reivindicación 1, en donde la tendencia dinámica se estima al aplicar un filtro de respuesta de impulso infinito de paso bajo a la señal no estacionaria.
4. 4. Un método de acuerdo con la reivindicación 1, en donde la segmentación adaptativa de la función de guía usa cruces por cero y mínimos globales y locales del módulo de la función de guía como puntos de segmentación.
5. 5. Un método (100) de acuerdo con la reivindicación 1, en donde el inicio de cada QGK emparejado con el respectivo ECK coincide con el punto de segmentación en un inicio del respectivo ECK.
6. 6. Un método (100) de acuerdo con la reivindicación 5, en donde la etapa de estimación comprende además la subetapa de estimar un inicio revisado de cada QGK que coincide con el ECK respectivo al encontrar una mejor coincidencia entre la fase del espectro complejo del ECK y la fase del espectro del complejo Gaussiano.
7. 7. Aparato para su uso en un método de acuerdo con las reivindicaciones de la 1 a la 6, el aparato que comprende:

   medios para estimar las tendencias dinámicas y de línea de base en el transcurso del tiempo de la señal, extraer y segmentar la función de guía;

   medios para dividir la señal no estacionaria en los segmentos sobre los que se desarrollan los componentes funcionales;

   medios para compensar la superposición del QGK que coincide con los componentes funcionales anteriores para formar un ECK, que comprende:

   medios para transformar el ECK en un dominio de la frecuencia;

   medios para estimar los parámetros de peso y dispersión del QGK, y validación de la fiabilidad del parámetro;

   medios para estimar el tiempo de inicio del QGK;

   medios para la expansión del modelo de la señal no estacionaria después de cada ciclo de recursión; medios para eliminar los QGK desfavorables y reorganizar los QGK restantes; y

   medios para crear modelos parciales de la señal no estacionaria con parámetros que pertenecen a clases predefinidas.