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 wavelet para la clasificación de señal ECG

Este ejemplo muestra cómo clasificar las señales de electrocardiograma humano (ECG) mediante la dispersión de tiempo de wavelet y un clasificador de máquina de vectores de soporte (SVM). En la dispersión de Wavelet, los datos se propagan a través de una serie de transformaciones de Wavelet, no linealidades y promediación para producir representaciones de baja varianza de series de tiempo. La dispersión de tiempo wavelet produce representaciones de señal insensibles a los cambios en la señal de entrada sin sacrificar la discriminabilidad de clase. Para ejecutar este ejemplo, debe tener el cuadro de herramientas wavelet™ y el cuadro de herramientas estadísticas y aprendizaje automático™. Los datos utilizados en este ejemplo están disponibles públicamente.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 sinusales normales. El ejemplo utiliza 162 grabaciones de ECG de tres bases de datos de PhysioNet: [3] [5], [3] y [2] [3].MIT-BIH base de datos de arritmiasMIT-BIH base de datos de ritmo sinusal normalLa 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 sinusales normales. El objetivo es entrenar un clasificador para distinguir entre arritmias (ARR), insuficiencia cardíaca congestiva (CHF) y ritmo sinusal normal (NSR).

Descargar datos

El primer paso es descargar los datos de la.Repositorio de GitHub Para descargar los datos, haga clic y seleccione.Clone or downloadDownload ZIP Guarde el archivo en una carpeta en la que tenga permiso de escritura.physionet_ECG_data-master.zip Las instrucciones de este ejemplo suponen que ha descargado el archivo en el directorio temporal (en MATLAB).tempdir Modifique las instrucciones subsiguientes para descomprimir y cargar los datos si decide descargar los datos de la carpeta diferente.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 usando.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 pre-procesamiento aplicados a cada grabación ECG.

Cargue archivos

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

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 162-by-65536 en la que cada fila es una grabación de ECG muestreada a 128 Hertz.ECGDataDataLabelsData Cada serie de tiempo de ECG tiene una duración total de 512 segundos. es una matriz de celdas de diagnóstico de 162-by-1, 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 prueba

Divida aleatoriamente los datos en dos conjuntos: conjuntos de datos de entrenamiento y pruebas. 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 de 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 un 70% de los datos de cada clase al conjunto de entrenamiento. El 30% restante se mantiene para las pruebas (predicción) y se asignan al conjunto de pruebas.

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

Hay 113 registros en el conjunto y 49 registros en.trainDatatestData Por diseño, los datos de entrenamiento 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 generales 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  

Las muestras de parcela

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

helperPlotRandomRecords(ECGData,14) 

Dispersión de tiempo wavelet

Los parámetros clave que se especifican en una descomposición de dispersión de tiempo de wavelet son la escala de la invariante del tiempo, el número de transformaciones de wavelet y el número de longitudes de onda 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 longitudes de onda por octava en el primer banco de filtros y 1 wavelet por octava en el segundo banco de filtros. La escala de invarianza 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 filtro 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 invarianza, obtener la transformada inversa de Fourier de la función de escalado y centrar en 0 segundos en el tiempo. Las dos líneas negras verticales marcan los límites-75 y 75 segundos. Además, trace las partes reales e imaginarias de la wavelet 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 excede la escala invariable determinada por el tiempo de soporte de la función de escalado. Esta es una propiedad importante de la dispersión de tiempo de wavelet.

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 de una sola señal.featureMatrix

scat_features_train = featureMatrix(sf,trainData'); 

La salida en este caso es 416-by-16-by-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 ondulación de wavelet se muestrea de forma crítica 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 rutas de dispersión 416.

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 señales 113 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 de ECG en el conjunto de entrenamiento.scat_features_test Después de remodelar para el clasificador SVM, la matriz de características es 784-by-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 hemos obtenido 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. La primera 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 DataSet, 16 para cada una de las señales 162. 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); 

Calcule 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 por mayoría simple para obtener una predicción de clase única para cada representación de dispersión.

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

Determine la precisión de validación cruzada real basada en el 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 de validación cruzada. Si usted examina, usted ve que las 2 ventanas de tiempo mal clasificadas en son atribuibles a 2 señales donde 15 de las 16 ventanas de dispersión fueron clasificadas correctamente.ClassCountsconfmatCV

Clasificación de SVM

Para el siguiente análisis, ajustamos una SVM cuadrática de varias clases a los datos de entrenamiento solamente (70%) y, a continuación, utilice ese modelo para realizar predicciones sobre el 30% de los datos mantenidos 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 DataSet de prueba es aproximadamente del 98%. La matriz de confusión muestra que un registro CHF se clasifica erróneamente como ARR. Todas las 48 otras señales se clasifican correctamente.

Resumen

Este ejemplo utilizaba la dispersión de tiempo de wavelet y un clasificador de SVM para clasificar las formas de onda de ECG en una de tres clases de diagnóstico. La dispersión de wavelet demostró ser un potente extractor de características, que requería solo un conjunto mínimo de parámetros especificados por el usuario para producir un conjunto de características sólidas para la clasificación. Compare esto con el ejemplo, que requería una cantidad significativa de experiencia para que las características de artesanía se utilizen en la clasificación.Signal Classification Using Wavelet-Based Features and Support Vector Machines (Wavelet Toolbox) Con la dispersión de tiempo de Wavelet, solo se requiere especificar la escala de la invarianza de tiempo, el número de bancos de filtros (o transformaciones de wavelet) y el número de longitudes de onda por octava. La combinación de una transformación de dispersión de wavelet y un clasificador de SVM produjo una clasificación del 100% en un modelo validado de forma cruzada 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 Colegio Americano de Cardiología 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. Fisiobank, 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 la dispersión invariable. Comunicaciones en matemática pura y aplicada, 65, 10, PP. 1331-1398.

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

Funciones de soporte

Traza cuatro señales de ECG elegidas aleatoriamente.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 previstas 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