Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Dispersión de tiempo de onda para la clasificación de señales ECG

Este ejemplo muestra cómo clasificar las señales de electrocardiograma humano (ECG) mediante la dispersión de tiempo de onda y un clasificador de máquina vectorial de soporte (SVM). En la dispersión de ondas, los datos se propagan a través de una serie de transformaciones de wavelet, no linealidades y promedio para producir representaciones de baja varianza de series temporales. La dispersión del tiempo de onda produce representaciones de señal insensibles a los cambios en la señal de entrada sin sacrificar la discriminabilidad de la clase. Debe tener Wavelet Toolbox™ y Statistics and Machine Learning Toolbox™ para ejecutar este ejemplo. Los datos utilizados en este ejemplo están disponibles públicamente desde .PhysioNet Puede encontrar un enfoque de aprendizaje profundo para este problema de clasificación en este ejemplo y un enfoque de aprendizaje automático en este ejemplo.Classify Time Series Using Wavelet Analysis and Deep Learning (Wavelet Toolbox)Signal Classification Using Wavelet-Based Features and Support Vector Machines (Wavelet Toolbox)

Descripción de los datos

Este ejemplo utiliza datos de ECG obtenidos de tres grupos, o clases, de personas: personas con arritmia cardíaca, personas con insuficiencia cardíaca congestiva y personas con ritmos sinuales normales. El ejemplo utiliza 162 grabaciones ECG de tres bases de datos PhysioNet: [3][5], [3] y [2][3].Base de datos de arritmias del MIT-BIHBase de datos de ritmo sinusal normal MIT-BIHLa base de datos de insuficiencia cardíaca congestiva BIDMC En total, hay 96 grabaciones de personas con arritmia, 30 grabaciones de personas con insuficiencia cardíaca congestiva y 36 grabaciones de personas con ritmos sinuales normales. El objetivo es entrenar a un clasificador para distinguir entre arritmia (ARR), insuficiencia cardíaca congestiva (CHF) y ritmo sinusal normal (NSR).

Descargar datos

El primer paso es descargar los datos del archivo .Repositorio GitHub Para descargar los datos, haga clic y seleccione .Clone or downloadDownload ZIP Guarde el archivo en una carpeta donde tenga permiso de escritura.physionet_ECG_data-master.zip En las instrucciones de este ejemplo se supone que ha descargado el archivo en el directorio temporal, (en MATLAB).tempdir Modifique las instrucciones posteriores para descomprimir y cargar los datos si decide descargar los datos en una carpeta diferente de .tempdir Si está familiarizado con Git, puede descargar la última versión de las herramientas ( ) y obtener los datos de un símbolo del sistema mediante .Gitgit clone https://github.com/mathworks/physionet_ECG_data/

El archivo contienephysionet_ECG_data-master.zip

  • ECGData.zip

  • README.md

y ECGData.zip contiene

  • ECGData.mat

  • Modified_physionet_data.txt

  • License.txt.

ECGData.mat contiene los datos utilizados en este ejemplo. El archivo .txt, Modified_physionet_data.txt, es requerido por la política de copia de PhysioNet y proporciona las atribuciones de origen para los datos, así como una descripción de los pasos de preprocesamiento aplicados a cada registro ECG.

Cargar archivos

Si ha seguido las instrucciones de descarga de la sección anterior, escriba los siguientes comandos para descomprimir los dos archivos de archivado.

unzip(fullfile(tempdir,'physionet_ECG_data-master.zip'),tempdir) unzip(fullfile(tempdir,'physionet_ECG_data-master','ECGData.zip'),...     fullfile(tempdir,'ECGData')) 

Después de descomprimir el archivo ECGData.zip, cargue los datos en MATLAB.

load(fullfile(tempdir,'ECGData','ECGData.mat')) 

es una matriz de estructura con dos campos: y . es una matriz de 162 por 65536 donde cada fila es una grabación ECG muestreada a 128 hercios.ECGDataDataLabelsData Cada serie temporal de ECG tiene una duración total de 512 segundos. es una matriz de celdas de 162 por 1 de etiquetas de diagnóstico, una para cada fila de .LabelsData Las tres categorías de diagnóstico son: 'ARR' (arritmia), 'CHF' (insuficiencia cardíaca congestiva) y 'NSR' (ritmo sinusal normal).

Crear datos de entrenamiento y pruebas

Divida aleatoriamente los datos en dos conjuntos: conjuntos de datos de entrenamiento y prueba. La función auxiliar realiza la división aleatoria. acepta el porcentaje de división deseado para los datos de entrenamiento y .helperRandomSplithelperRandomSplitECGData La función genera dos conjuntos de datos junto con un conjunto de etiquetas para cada uno.helperRandomSplit Cada fila de y es una señal ECG.trainDatatestData Cada elemento de y contiene la etiqueta de clase para la fila correspondiente de las matrices de datos.trainLabelstestLabels En este ejemplo, asignamos aleatoriamente el 70% de los datos de cada clase al conjunto de entrenamiento. El 30% restante se mantiene para las pruebas (predicción) y se asigna al conjunto de pruebas.

percent_train = 70; [trainData,testData,trainLabels,testLabels] = ...     helperRandomSplit(percent_train,ECGData); 

Hay 113 registros en el set y 49 registros en .trainDatatestData Por diseño, los datos de formación contienen el 69,75% (113/162) de los datos. Recuerde que la clase ARR representa el 59,26% de los datos (96/162), la clase CHF representa el 18,52% (30/162) y la clase NSR representa el 22,22% (36/162). Examine el porcentaje de cada clase en los conjuntos de entrenamiento y prueba. Los porcentajes de cada uno son coherentes con los porcentajes de clase globales del conjunto de datos.

Ctrain = countcats(categorical(trainLabels))./numel(trainLabels).*100 Ctest = countcats(categorical(testLabels))./numel(testLabels).*100 
 Ctrain =     59.2920    18.5841    22.1239   Ctest =     59.1837    18.3673    22.4490  

Muestras de parcela

Trazar las primeras miles de muestras de cuatro registros seleccionados al azar de .ECGData La función auxiliar hace esto. acepta y una semilla aleatoria como entrada.helperPlotRandomRecordshelperPlotRandomRecordsECGData El valor inicial se establece en 14 para que se trace al menos un registro de cada clase. Puede ejecutar con como único argumento de entrada tantas veces como desee para obtener una idea de la variedad de formas de onda ECG asociadas a cada clase.helperPlotRandomRecordsECGData Puede encontrar el código fuente para esto y todas las funciones auxiliares en la sección Funciones auxiliares al final de este ejemplo.

helperPlotRandomRecords(ECGData,14) 

Dispersión del tiempo de Wavelet

Los parámetros clave que se pueden especificar en una descomposición de dispersión de tiempo de onda son la escala de la invariante de tiempo, el número de transformaciones de wavelet y el número de ondas por octava en cada uno de los bancos de filtros de wavelet. En muchas aplicaciones, la cascada de dos bancos de filtros es suficiente para lograr un buen rendimiento. En este ejemplo, construimos una descomposición de dispersión de tiempo de wavelet con los bancos de filtros predeterminados: 8 ondas por octava en el primer banco de filtroy 1 onda por octava en el segundo banco de filtros. La escala de invariancia se establece en 150 segundos.

N = size(ECGData.Data,2); sf = waveletScattering('SignalLength',N, 'InvarianceScale',150,'SamplingFrequency',128); 

Puede visualizar los filtros de wavelet en los dos bancos de filtros con lo siguiente.

[fb,f,filterparams] = filterbank(sf); subplot(211) plot(f,fb{2}.psift) xlim([0 128]) grid on title('1st Filter Bank Wavelet Filters'); subplot(212) plot(f,fb{3}.psift) xlim([0 128]) grid on title('2nd Filter Bank Wavelet Filters'); xlabel('Hz'); 

Para demostrar la escala de invariancia, obtenga la transformación inversa de Fourier de la función de escalado y celénela en 0 segundos en el tiempo. Las dos líneas negras verticales marcan los límites de -75 y 75 segundos. Además, trace las partes reales e imaginarias de la onda de escala más gruesa (frecuencia más baja) del primer banco de filtros. Tenga en cuenta que la wavelet de escala más gruesa no supera la escala invariable determinada por la compatibilidad de tiempo de la función de escalado. Esta es una propiedad importante de la dispersión del tiempo de las olas.

figure; phi = ifftshift(ifft(fb{1}.phift)); psiL1 = ifftshift(ifft(fb{2}.psift(:,end))); t = (-2^15:2^15-1).*1/128; scalplt = plot(t,phi); hold on grid on ylim([-1.5e-4 1.6e-4]); plot([-75 -75],[-1.5e-4 1.6e-4],'k--'); plot([75 75],[-1.5e-4 1.6e-4],'k--'); xlabel('Seconds'); ylabel('Amplitude'); wavplt = plot(t,[real(psiL1) imag(psiL1)]); legend([scalplt wavplt(1) wavplt(2)],{'Scaling Function','Wavelet-Real Part','Wavelet-Imaginary Part'}); title({'Scaling Function';'Coarsest-Scale Wavelet First Filter Bank'}) 

Después de construir el marco de descomposición de dispersión, obtenga los coeficientes de dispersión para los datos de entrenamiento como una matriz. Cuando se ejecuta con varias señales, cada columna se trata una sola señal.featureMatrix

scat_features_train = featureMatrix(sf,trainData'); 

La salida de en este caso es 416-por-16-por-113.featureMatrix Cada página del tensor, , es la transformación de dispersión de una señal.scat_features_train La transformación de dispersión de wavelet se reduce críticamente en el tiempo en función del ancho de banda de la función de escalado. En este caso, esto da como resultado 16 ventanas de tiempo para cada una de las 416 rutas de dispersión.

Para obtener una matriz compatible con el clasificador SVM, remodele la transformación de dispersión multiseñal en una matriz donde cada columna corresponde a una ruta de dispersión y cada fila es una ventana de tiempo de dispersión. En este caso, obtendrá 1808 filas porque hay 16 ventanas de tiempo para cada una de las 113 señales en los datos de entrenamiento.

Nwin = size(scat_features_train,2); scat_features_train = permute(scat_features_train,[2 3 1]); scat_features_train = reshape(scat_features_train,...     size(scat_features_train,1)*size(scat_features_train,2),[]); 

Repita el proceso para los datos de prueba. Inicialmente, es 416 por 16 por 49 porque hay 49 formas de onda ECG en el conjunto de entrenamiento.scat_features_test Después de remodelar para el clasificador SVM, la matriz de características es 784-por-416.

scat_features_test = featureMatrix(sf,testData'); scat_features_test = permute(scat_features_test,[2 3 1]); scat_features_test = reshape(scat_features_test,...     size(scat_features_test,1)*size(scat_features_test,2),[]); 

Porque para cada señal obtuvimos 16 ventanas de dispersión, necesitamos crear etiquetas para que coincidan con el número de ventanas. La función auxiliar hace esto en función del número de ventanas.createSequenceLabels

[sequence_labels_train,sequence_labels_test] = createSequenceLabels(Nwin,trainLabels,testLabels); 

Validación cruzada

Para la clasificación, se realizan dos análisis. El primero utiliza todos los datos de dispersión y se adapta a una SVM de varias clases con un kernel cuadrático. En total hay 2592 secuencias de dispersión en todo el conjunto de datos, 16 para cada una de las 162 señales. La tasa de error, o pérdida, se estima mediante la validación cruzada de 5 veces.

scat_features = [scat_features_train; scat_features_test]; allLabels_scat = [sequence_labels_train; sequence_labels_test]; rng(1); template = templateSVM(...     'KernelFunction', 'polynomial', ...     'PolynomialOrder', 2, ...     'KernelScale', 'auto', ...     'BoxConstraint', 1, ...     'Standardize', true); classificationSVM = fitcecoc(...     scat_features, ...     allLabels_scat, ...     'Learners', template, ...     'Coding', 'onevsone', ...     'ClassNames', {'ARR';'CHF';'NSR'}); kfoldmodel = crossval(classificationSVM, 'KFold', 5); 

Calcular la pérdida y la matriz de confusión. Visualice la precisión.

predLabels = kfoldPredict(kfoldmodel); loss = kfoldLoss(kfoldmodel)*100; confmatCV = confusionmat(allLabels_scat,predLabels) fprintf('Accuracy is %2.2f percent.\n',100-loss); 
 confmatCV =          1535           0           1            1         479           0            0           0         576  Accuracy is 99.92 percent. 

La precisión es 99.92%, que es bastante bueno, pero los resultados reales son mejor teniendo en cuenta que aquí cada ventana de tiempo se clasifica por separado. Hay 16 clasificaciones separadas para cada señal. Utilice un voto de mayoría simple para obtener una predicción de una sola clase para cada representación de dispersión.

classes = categorical({'ARR','CHF','NSR'}); [ClassVotes,ClassCounts] = helperMajorityVote(predLabels,[trainLabels; testLabels],classes); 

Determine la precisión real de la validación cruzada en función del modo de las predicciones de clase para cada conjunto de ventanas de tiempo de dispersión. Si el modo no es único para un conjunto determinado, devuelve un error de clasificación indicado por .helperMajorityVote'NoUniqueMode' Esto da como resultado una columna adicional en la matriz de confusión, que en este caso es todos ceros porque existe un modo único para cada conjunto de predicciones de dispersión.

CVaccuracy = sum(eq(ClassVotes,categorical([trainLabels; testLabels])))/162*100; fprintf('True cross-validation accuracy is %2.2f percent.\n',CVaccuracy); MVconfmatCV = confusionmat(ClassVotes,categorical([trainLabels; testLabels])); MVconfmatCV 
True cross-validation accuracy is 100.00 percent.  MVconfmatCV =      96     0     0     0      0    30     0     0      0     0    36     0      0     0     0     0  

La dispersión ha clasificado correctamente todas las señales en el modelo validado cruzado. Si examina , verá que las 2 ventanas de tiempo mal clasificadas en son atribuibles a 2 señales donde 15 de las 16 ventanas de dispersión se clasificaron correctamente.ClassCountsconfmatCV

Clasificación SVM

Para el siguiente análisis, ajustamos una SVM cuadrática multiclase solo a los datos de entrenamiento (70%) y luego usa ese modelo para hacer predicciones sobre el 30% de los datos retenidos para las pruebas. Hay 49 registros de datos en el conjunto de pruebas. Utilice un voto mayoritario en las ventanas de dispersión individuales.

model = fitcecoc(...      scat_features_train, ...      sequence_labels_train, ...      'Learners', template, ...      'Coding', 'onevsone', ...      'ClassNames', {'ARR','CHF','NSR'}); predLabels = predict(model,scat_features_test); [TestVotes,TestCounts] = helperMajorityVote(predLabels,testLabels,classes); testaccuracy = sum(eq(TestVotes,categorical(testLabels)))/numel(testLabels)*100; fprintf('The test accuracy is %2.2f percent. \n',testaccuracy); testconfmat = confusionmat(TestVotes,categorical(testLabels)); testconfmat 
The test accuracy is 97.96 percent.   testconfmat =      29     1     0     0      0     8     0     0      0     0    11     0      0     0     0     0  

La precisión de clasificación en el conjunto de datos de prueba es de aproximadamente el 98%. La matriz de confusión muestra que un registro CHF está clasificado erróneamente como ARR. Las otras 48 señales están correctamente clasificadas.

Resumen

Este ejemplo utiliza la dispersión del tiempo de las ondas y un clasificador SVM para clasificar las formas de onda ECG en una de las tres clases de diagnóstico. La dispersión de wavelet demostró ser un potente extractor de características, que solo requería un conjunto mínimo de parámetros especificados por el usuario para producir un conjunto de características robustas para la clasificación. Compare esto con el ejemplo, que requería una cantidad significativa de experiencia en las características artesanales para usar en la clasificación.Signal Classification Using Wavelet-Based Features and Support Vector Machines (Wavelet Toolbox) Con la dispersión de tiempo de las ondas, solo es necesario especificar la escala de la invariancia de tiempo, el número de bancos de filtros (u transformaciones de wavelet) y el número de wavelets por octava. La combinación de una transformación de dispersión de wavelet y un clasificador SVM produjo una clasificación del 100% en un modelo validado cruzado y una clasificación correcta del 98% al aplicar una SVM a las transformaciones de dispersión de un conjunto de pruebas de retención.

Referencias

  1. Anden, J., Mallat, S. 2014. Espectro de dispersión profunda, Transacciones IEEE en procesamiento de señales, 62, 16, pp. 4114-4128.

  2. Baim DS, Colucci WS, Monrad ES, Smith HS, Wright RF, Lanoue A, Gauthier DF, Ransil BJ, Grossman W, Braunwald E. Supervivencia de pacientes con insuficiencia cardíaca congestiva grave tratados con milrinona oral. J American College of Cardiology 1986 Mar; 7(3):661-670.

  3. Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE. PhysioBank, PhysioToolkit y PhysioNet: Componentes de un nuevo recurso de investigación para señales fisiológicas complejas. .Circulation Vol. 101, No 23, 13 de junio de 2000, pp. e215-e220.http://circ.ahajournals.org/content/101/23/e215.full

  4. Mallat, S., 2012. Agrupar dispersión invariable. Comunicaciones en Matemáticas Puras y Aplicadas, 65, 10, págs. 1331-1398.

  5. Moody GB, Mark RG. El impacto de la base de datos de arritmias del MIT-BIH. IEEE Eng en Med y Biol 20(3):45-50 (mayo-junio 2001). (PMID: 11446209)

Funciones de apoyo

Traza cuatro señales ECG elegidas aleatoriamente de .helperPlotRandomRecordsECGData

function helperPlotRandomRecords(ECGData,randomSeed) % This function is only intended to support the XpwWaveletMLExample. It may % change or be removed in a future release.  if nargin==2     rng(randomSeed) end  M = size(ECGData.Data,1); idxsel = randperm(M,4); for numplot = 1:4     subplot(2,2,numplot)     plot(ECGData.Data(idxsel(numplot),1:3000))     ylabel('Volts')     if numplot > 2         xlabel('Samples')     end     title(ECGData.Labels{idxsel(numplot)}) end  end 

Busca el modo en las etiquetas de clase predichas para cada conjunto de ventanas de tiempo de dispersión.helperMajorityVote La función devuelve los modos de etiqueta de clase y el número de predicciones de clase para cada conjunto de ventanas de tiempo de dispersión. Si no hay ningún modo único, devuelve una etiqueta de clase de "error" para indicar que el conjunto de ventanas de dispersión es un error de clasificación.helperMajorityVote

function [ClassVotes,ClassCounts] = helperMajorityVote(predLabels,origLabels,classes) % This function is in support of ECGWaveletTimeScatteringExample. It may % change or be removed in a future release.  % Make categorical arrays if the labels are not already categorical predLabels = categorical(predLabels); origLabels = categorical(origLabels); % Expects both predLabels and origLabels to be categorical vectors Npred = numel(predLabels); Norig = numel(origLabels); Nwin = Npred/Norig; predLabels = reshape(predLabels,Nwin,Norig); ClassCounts = countcats(predLabels); [mxcount,idx] = max(ClassCounts); ClassVotes = classes(idx); % Check for any ties in the maximum values and ensure they are marked as % error if the mode occurs more than once modecnt = modecount(ClassCounts,mxcount); ClassVotes(modecnt>1) = categorical({'error'}); ClassVotes = ClassVotes(:);  %------------------------------------------------------------------------- function modecnt = modecount(ClassCounts,mxcount) modecnt = Inf(size(ClassCounts,2),1); for nc = 1:size(ClassCounts,2)     modecnt(nc) = histc(ClassCounts(:,nc),mxcount(nc)); end  end  % EOF end