Введение

ЭКГ регистрируют электрическую активность сердца человека в течение определенного периода времени. Врачи используют ЭКГ, чтобы визуально определить, является ли сердцебиение пациента нормальным или нерегулярным.

Мерцательная аритмия (AFib) - это тип нерегулярного сердцебиения, который возникает, когда верхние ёмкости сердца, предсердие, бьются вне координации с нижними ёмкостями, желудочками.

Этот пример использует данные ЭКГ из PhysioNet 2017 Challenge [1], [2], [3], который доступен в https://physionet.org/challenge/2017/. Данные состоят из набора сигналов ЭКГ, дискретизированных на частоте 300 Гц и разделенных группой экспертов на четыре различных класса: Normal (N), AFib (A), Other Rhythm (O) и Noisy Recording (~). В этом примере показано, как автоматизировать процесс классификации с помощью глубокого обучения. Процедура исследует двоичный классификатор, который может дифференцировать сигналы Normal ECG от сигналов, имеющих признаки AFib.

Этот пример использует сети долгой краткосрочной памяти (LSTM), тип рекуррентной нейронной сети (RNN), хорошо подходящий для изучения данных последовательности и timeseries. Сеть LSTM может изучать долгосрочные зависимости между временными шагами последовательности. LSTM- слоя (lstmLayer) может смотреть на временную последовательность в прямом направлении, в то время как двунаправленный слой LSTM (bilstmLayer) может посмотреть на временную последовательность как в прямом, так и в обратном направлениях. Этот пример использует двунаправленный слой LSTM.

Чтобы ускорить процесс обучения, запустите этот пример на машине с графическим процессором. Если ваша машина имеет графический процессор и Parallel Computing Toolbox™, то MATLAB ® автоматически использует GPU для обучения; в противном случае используется центральный процессор.

Загрузка и исследование данных

Запуск ReadPhysionetData скрипт, чтобы загрузить данные с сайта PhysioNet и сгенерировать MAT-файл (PhysionetData.mat), который содержит сигналы ЭКГ в соответствующем формате. Загрузка данных может занять несколько минут. Используйте условный оператор, который запускает скрипт только если PhysionetData.mat не существует в текущей папке.

if ~isfile('PhysionetData.mat')
    ReadPhysionetData         
end
load PhysionetData

Операция загрузки добавляет две переменные в рабочую область: Signals и Labels. Signals - массив ячеек, который содержит сигналы ECG. Labels является категориальным массивом, который содержит соответствующие метки основной истины сигналов.

Signals(1:5)

ans=5×1 cell array
    {1×9000  double}
    {1×9000  double}
    {1×18000 double}
    {1×9000  double}
    {1×18000 double}

Labels(1:5)

ans = 5×1 categorical
     N 
     N 
     N 
     A 
     A 

Используйте summary функция, чтобы увидеть, сколько сигналов AFib и сигналов Normal содержится в данных.

summary(Labels)

     A       738 
     N      5050 

Сгенерируйте гистограмму длин сигнала. Большинство сигналов имеют длину 9000 выборки.

L = cellfun(@length,Signals);
h = histogram(L);
xticks(0:3000:18000);
xticklabels(0:3000:18000);
title('Signal Lengths')
xlabel('Length')
ylabel('Count')

Визуализируйте сегмент одного сигнала от каждого класса. Сердцебиение AFib разнесено с нерегулярными интервалами, в то время как сердцебиение Normal происходит регулярно. AFib сигналов сердцебиения также часто не хватает P волны, которая импульсирует перед комплексом QRS в сигнале Normal heartbeat. График сигнала Normal показывает волну P и комплекс QRS.

normal = Signals{1};
aFib = Signals{4};

subplot(2,1,1)
plot(normal)
title('Normal Rhythm')
xlim([4000,5200])
ylabel('Amplitude (mV)')
text(4330,150,'P','HorizontalAlignment','center')
text(4370,850,'QRS','HorizontalAlignment','center')

subplot(2,1,2)
plot(aFib)
title('Atrial Fibrillation')
xlim([4000,5200])
xlabel('Samples')
ylabel('Amplitude (mV)')

Подготовка данных к обучению

Во время обучения trainNetwork функция разделяет данные на мини-пакеты. Затем функция заполняет или обрезает сигналы в том же мини-пакете, так что все они имеют одинаковую длину. Слишком большое количество заполнения или усечения может негативно эффект на эффективность сети, поскольку сеть может неправильно интерпретировать сигнал на основе добавленной или удаленной информации.

Чтобы избежать чрезмерного заполнения или усечения, примените segmentSignals функция к сигналам ЭКГ таким образом они все 9000 выборок длинные. Функция игнорирует сигналы с менее чем 9000 выборками. Если сигнал имеет более 9000 выборки, segmentSignals разбивает его на как можно больше 9000 сегментов выборки и игнорирует оставшиеся выборки. Для примера сигнал с выборками 18500 становится двумя 9000-выборочными сигналами, а оставшиеся 500 выборки игнорируются.

[Signals,Labels] = segmentSignals(Signals,Labels);

Просмотр первых пяти элементов Signals массив, чтобы убедиться, что каждая запись теперь имеет длину 9000 выборок.

Signals(1:5)

ans=5×1 cell array
    {1×9000 double}
    {1×9000 double}
    {1×9000 double}
    {1×9000 double}
    {1×9000 double}

Обучите классификатор, используя данные необработанного сигнала

Для разработки классификатора используйте необработанные сигналы, сгенерированные в предыдущем разделе. Разделите сигналы на набор обучающих данных для обучения классификатора и набор проверки для проверки точности классификатора на новых данных.

Используйте summary функция, чтобы показать, что отношение сигналов AFib к сигналам Normal составляет 718:4937 или приблизительно 1:7.

summary(Labels)

     A       718 
     N      4937 

Поскольку около 7/8 сигналов являются Нормальными, классификатор узнает, что он может достичь высокой точности просто путем классификации всех сигналов как Нормальных. Чтобы избежать этого смещения, увеличьте данные AFib, дублируя сигналы AFib в наборе данных так, чтобы было одинаковое количество сигналов Normal и AFib. Это дублирование, обычно называемое избыточной дискретизацией, является одной из форм увеличения данных, используемых в глубоком обучении.

Разделите сигналы согласно их классу.

afibX = Signals(Labels=='A');
afibY = Labels(Labels=='A');

normalX = Signals(Labels=='N');
normalY = Labels(Labels=='N');

Далее используйте dividerand разделить цели из каждого класса случайным образом на наборы для обучения и тестирования.

[trainIndA,~,testIndA] = dividerand(718,0.9,0.0,0.1);
[trainIndN,~,testIndN] = dividerand(4937,0.9,0.0,0.1);

XTrainA = afibX(trainIndA);
YTrainA = afibY(trainIndA);

XTrainN = normalX(trainIndN);
YTrainN = normalY(trainIndN);

XTestA = afibX(testIndA);
YTestA = afibY(testIndA);

XTestN = normalX(testIndN);
YTestN = normalY(testIndN);

Теперь для обучения существует 646 сигналов AFib и 4443 сигнала Normal. Чтобы достичь одинакового количества сигналов в каждом классе, используйте первые 4438 сигналов Normal, а затем используйте repmat повторить первые 634 AFib сигнала семь раз.

Для проверки существует 72 сигнала AFib и 494 сигнала Normal. Используйте первые сигналы 490 Normal, а затем используйте repmat чтобы повторить первые 70 сигналов AFib семь раз. По умолчанию нейронная сеть случайным образом тасует данные перед обучением, гарантируя, что смежные сигналы не все имеют одну и ту же метку.

XTrain = [repmat(XTrainA(1:634),7,1); XTrainN(1:4438)];
YTrain = [repmat(YTrainA(1:634),7,1); YTrainN(1:4438)];

XTest = [repmat(XTestA(1:70),7,1); XTestN(1:490)];
YTest = [repmat(YTestA(1:70),7,1); YTestN(1:490);];

Распределение между сигналами Normal и AFib теперь равномерно сбалансировано как в набор обучающих данных, так и в проверку наборе.

summary(YTrain)

     A      4438 
     N      4438 

summary(YTest)

     A      490 
     N      490 

Определите сетевую архитектуру LSTM

Сети LSTM могут изучать долгосрочные зависимости между временными шагами данных последовательности. Этот пример использует двунаправленный слой LSTM bilstmLayer, поскольку он рассматривает последовательность как в прямом, так и в обратном направлениях.

Поскольку входные сигналы имеют по одной размерности, задайте размер входа, чтобы быть последовательностями размера 1. Задайте двунаправленный слой LSTM с размером выходом 100 и выведите последний элемент последовательности. Эта команда предписывает двунаправленному слою LSTM сопоставить входные временные ряды с 100 функциями, а затем подготавливает выход для полностью подключенного уровня. Наконец, задайте два класса путем включения полносвязного слоя размера 2, за которым следуют слой softmax и слой классификации.

layers = [ ...
    sequenceInputLayer(1)
    bilstmLayer(100,'OutputMode','last')
    fullyConnectedLayer(2)
    softmaxLayer
    classificationLayer
    ]

layers = 
  5x1 Layer array with layers:

     1   ''   Sequence Input          Sequence input with 1 dimensions
     2   ''   BiLSTM                  BiLSTM with 100 hidden units
     3   ''   Fully Connected         2 fully connected layer
     4   ''   Softmax                 softmax
     5   ''   Classification Output   crossentropyex

Далее задайте опции обучения для классификатора. Установите 'MaxEpochs' 10, чтобы позволить сети сделать 10 проходов через обучающие данные. A 'MiniBatchSize' 150 направляет сеть на просмотр 150 обучающих сигналов за раз. Система координат 'InitialLearnRate' 0,01 помогает ускорить процесс обучения. Задайте 'SequenceLength' 1000, чтобы разбить сигнал на меньшие части, чтобы машина не исчерпала память, рассматривая слишком много данных в одно время. Установите 'GradientThreshold'to 1, чтобы стабилизировать процесс обучения путем предотвращения слишком больших градиентов. Задайте 'Plots' как 'training-progress' чтобы сгенерировать графики, которые показывают изображение процесса обучения с увеличениями количества итераций. Задайте 'Verbose' на false чтобы подавить выход таблицы, который соответствует данным, показанным на графике. Если вы хотите увидеть эту таблицу, задайте 'Verbose' на true.

Этот пример использует решатель адаптивной оценки момента (ADAM). ADAM работает лучше с RNN, такими как LSTM, чем стохастический градиентный спуск по умолчанию с импульсом (SGDM) решателя.

options = trainingOptions('adam', ...
    'MaxEpochs',10, ...
    'MiniBatchSize', 150, ...
    'InitialLearnRate', 0.01, ...
    'SequenceLength', 1000, ...
    'GradientThreshold', 1, ...
    'ExecutionEnvironment',"auto",...
    'plots','training-progress', ...
    'Verbose',false);

Обучите сеть LSTM

Обучите сеть LSTM с заданными опциями обучения и архитектурой слоя при помощи trainNetwork. Поскольку набор обучающих данных большая, процесс обучения может занять несколько минут.

net = trainNetwork(XTrain,YTrain,layers,options);

Верхний подграфик графика процесса обучения представляет точность обучения, которая является точностью классификации для каждого мини-пакета. Когда обучение успешно прогрессирует, это значение обычно увеличивается до 100%. В нижней подграфике отображаются потери обучения, которые являются потерями перекрестной энтропии на каждом мини-пакете. Когда обучение успешно прогрессирует, это значение обычно уменьшается к нулю.

Если обучение не сходится, графики могут колебаться между значениями, не стремясь в определенном направлении вверх или вниз. Это колебание означает, что точность обучения не улучшается и потеря обучения не уменьшается. Такая ситуация может возникнуть с начала обучения, или графики могут появиться после некоторого предварительного повышения точности обучения. Во многих случаях изменение опций обучения может помочь сети достичь сходимости. Уменьшение MiniBatchSize или уменьшение InitialLearnRate это может привести к увеличению времени обучения, но может помочь сети лучше учиться.

Точность обучения классификатора колеблется между примерно 50% и примерно 60%, а в конце 10 эпох на обучение ушло уже несколько минут.

Визуализация точности обучения и проверки

Вычислите точность обучения, которая представляет точность классификатора на сигналах, на которых он был обучен. Во-первых, классифицируйте обучающие данные.

trainPred = classify(net,XTrain,'SequenceLength',1000);

В задачах классификации матрицы неточностей используются, чтобы визуализировать эффективность классификатора на наборе данных, для которых известны истинные значения. Класс Target является меткой основная истина сигнала, а класс Output - меткой, назначенной сигналу сетью. Метки осей представляют метки классов, AFib (A) и Normal (N).

Используйте confusionchart команда для вычисления общей точности классификации для предсказаний тестовых данных. Задайте 'RowSummary' как 'row-normalized' отображение истинных положительных частот и ложных положительных частот в сводных данных строк. Кроме того, задайте 'ColumnSummary' как 'column-normalized' отображение положительных прогностических значений и частот ложного обнаружения в сводных данных столбцов.

LSTMAccuracy = sum(trainPred == YTrain)/numel(YTrain)*100

LSTMAccuracy = 61.7283

figure
confusionchart(YTrain,trainPred,'ColumnSummary','column-normalized',...
              'RowSummary','row-normalized','Title','Confusion Chart for LSTM');

Теперь классифицируйте тестовые данные с той же сетью.

testPred = classify(net,XTest,'SequenceLength',1000);

Вычислите точность проверки и визуализируйте эффективность классификации как матрицу неточностей.

LSTMAccuracy = sum(testPred == YTest)/numel(YTest)*100

LSTMAccuracy = 66.2245

figure
confusionchart(YTest,testPred,'ColumnSummary','column-normalized',...
              'RowSummary','row-normalized','Title','Confusion Chart for LSTM');

Улучшите эффективность с редукцией данных

Редукция данных из данных может помочь улучшить точность обучения и проверки классификатора. Чтобы решить, какие функции извлечь, этот пример адаптирует подход, который вычисляет частотно-временные изображения, такие как спектрограммы, и использует их для обучения сверточных нейронных сетей (CNNs) [4], [5].

Визуализируйте спектрограмму каждого типа сигнала.

fs = 300;

figure
subplot(2,1,1);
pspectrum(normal,fs,'spectrogram','TimeResolution',0.5)
title('Normal Signal')

subplot(2,1,2);
pspectrum(aFib,fs,'spectrogram','TimeResolution',0.5)
title('AFib Signal')

Поскольку этот пример использует LSTM вместо CNN, важно преобразовать подход, чтобы он применялся к одномерным сигналам. Частотно-временные (TF) моменты извлекают информацию из спектрограмм. Каждый момент может использоваться как одномерная функция для ввода в LSTM.

Исследуйте два момента TF во временном интервале:

The instfreq функция оценивает зависящую от времени частоту сигнала как первый момент спектрограммы степени. Функция вычисляет спектрограмму, используя кратковременные преобразования Фурье с временными окнами. В этом примере функция использует 255 временных окон. Временные выходы функции соответствуют центрам временных окон.

Визуализируйте мгновенную частоту для каждого типа сигнала.

[instFreqA,tA] = instfreq(aFib,fs);
[instFreqN,tN] = instfreq(normal,fs);

figure
subplot(2,1,1);
plot(tN,instFreqN)
title('Normal Signal')
xlabel('Time (s)')
ylabel('Instantaneous Frequency')

subplot(2,1,2);
plot(tA,instFreqA)
title('AFib Signal')
xlabel('Time (s)')
ylabel('Instantaneous Frequency')

Использование cellfun для применения instfreq функция для каждой камеры в наборах обучения и тестирования.

instfreqTrain = cellfun(@(x)instfreq(x,fs)',XTrain,'UniformOutput',false);
instfreqTest = cellfun(@(x)instfreq(x,fs)',XTest,'UniformOutput',false);

Спектральная энтропия измеряет, насколько колючим плоским является спектр сигнала. Сигнал с колючим спектром, как и сумма синусоидов, имеет низкую спектральную энтропию. Сигнал с плоским спектром, как и белый шум, имеет высокую спектральную энтропию. The pentropy функция оценивает спектральную энтропию на основе спектрограммы степени. Как и в случае мгновенной оценки частоты, pentropy использует 255 временных окон для вычисления спектрограммы. Временные выходы функции соответствуют центру временных окон.

Визуализируйте спектральную энтропию для каждого типа сигнала.

[pentropyA,tA2] = pentropy(aFib,fs);
[pentropyN,tN2] = pentropy(normal,fs);

figure

subplot(2,1,1)
plot(tN2,pentropyN)
title('Normal Signal')
ylabel('Spectral Entropy')

subplot(2,1,2)
plot(tA2,pentropyA)
title('AFib Signal')
xlabel('Time (s)')
ylabel('Spectral Entropy')

Использование cellfun для применения pentropy функция для каждой камеры в наборах обучения и тестирования.

pentropyTrain = cellfun(@(x)pentropy(x,fs)',XTrain,'UniformOutput',false);
pentropyTest = cellfun(@(x)pentropy(x,fs)',XTest,'UniformOutput',false);

Объедините функции так, чтобы каждая камера в новых наборах обучения и тестирования имела две размерности или два признака.

XTrain2 = cellfun(@(x,y)[x;y],instfreqTrain,pentropyTrain,'UniformOutput',false);
XTest2 = cellfun(@(x,y)[x;y],instfreqTest,pentropyTest,'UniformOutput',false);

Визуализируйте формат новых входов. Каждая камера больше не содержит один сигнал 9000 выборок -long; теперь он содержит две функции 255 выборок -long.

XTrain2(1:5)

ans=5×1 cell array
    {2×255 double}
    {2×255 double}
    {2×255 double}
    {2×255 double}
    {2×255 double}

Стандартизация данных

Мгновенная частота и спектральная энтропия имеют средства, которые различаются почти на один порядок величины. Кроме того, текущее среднее значение частоты может быть слишком высоким для того, чтобы LSTM эффективно обучался. Когда сеть помещается на данные с большим средним и большой областью значений значений, большие входы могут замедлить обучение и сходимость сети [6].

mean(instFreqN)

ans = 5.5615

mean(pentropyN)

ans = 0.6326

Используйте средний набор обучающих данных и стандартное отклонение для стандартизации наборов обучения и тестирования. Стандартизация, или z-оценка, является популярным способом улучшить эффективность сети во время обучения.

XV = [XTrain2{:}];
mu = mean(XV,2);
sg = std(XV,[],2);

XTrainSD = XTrain2;
XTrainSD = cellfun(@(x)(x-mu)./sg,XTrainSD,'UniformOutput',false);

XTestSD = XTest2;
XTestSD = cellfun(@(x)(x-mu)./sg,XTestSD,'UniformOutput',false);

Показать средства стандартизированной мгновенной частотной и спектральной энтропии.

instFreqNSD = XTrainSD{1}(1,:);
pentropyNSD = XTrainSD{1}(2,:);

mean(instFreqNSD)

ans = -0.3211

mean(pentropyNSD)

ans = -0.2416

Изменение сетевой архитектуры LSTM

Теперь, когда каждый из сигналов имеет две размерности, необходимо изменить архитектуру сети, задав размер входной последовательности как 2. Задайте двунаправленный слой LSTM с размером выходом 100 и выведите последний элемент последовательности. Задайте два класса путем включения полносвязного слоя размера 2, а затем слоя softmax и слоя классификации.

layers = [ ...
    sequenceInputLayer(2)
    bilstmLayer(100,'OutputMode','last')
    fullyConnectedLayer(2)
    softmaxLayer
    classificationLayer
    ]

layers = 
  5x1 Layer array with layers:

     1   ''   Sequence Input          Sequence input with 2 dimensions
     2   ''   BiLSTM                  BiLSTM with 100 hidden units
     3   ''   Fully Connected         2 fully connected layer
     4   ''   Softmax                 softmax
     5   ''   Classification Output   crossentropyex

Задайте опции обучения. Установите максимальное количество эпох равное 30, чтобы сеть могла сделать 30 проходов через обучающие данные.

options = trainingOptions('adam', ...
    'MaxEpochs',30, ...
    'MiniBatchSize', 150, ...
    'InitialLearnRate', 0.01, ...
    'GradientThreshold', 1, ...
    'ExecutionEnvironment',"auto",...
    'plots','training-progress', ...
    'Verbose',false);

Обучите сеть LSTM с частотно-временными функциями

Обучите сеть LSTM с заданными опциями обучения и архитектурой слоя при помощи trainNetwork.

net2 = trainNetwork(XTrainSD,YTrain,layers,options);

Есть большое улучшение точности обучения. Значение потери перекрестной энтропии трендов к 0. Кроме того, время, необходимое для обучения, уменьшается, потому что моменты TF короче, чем необработанные последовательности.

Визуализация точности обучения и проверки

Классификация обучающих данных с помощью обновленной сети LSTM. Визуализируйте эффективность классификации как матрицу неточностей.

trainPred2 = classify(net2,XTrainSD);
LSTMAccuracy = sum(trainPred2 == YTrain)/numel(YTrain)*100

LSTMAccuracy = 83.5962

figure
confusionchart(YTrain,trainPred2,'ColumnSummary','column-normalized',...
              'RowSummary','row-normalized','Title','Confusion Chart for LSTM');

Классификация тестовых данных по обновленной сети. Постройте матрицу неточностей, чтобы изучить точность проверки.

testPred2 = classify(net2,XTestSD);

LSTMAccuracy = sum(testPred2 == YTest)/numel(YTest)*100

LSTMAccuracy = 80.1020

figure
confusionchart(YTest,testPred2,'ColumnSummary','column-normalized',...
              'RowSummary','row-normalized','Title','Confusion Chart for LSTM');

Заключение

Этот пример показывает, как создать классификатор для обнаружения фибрилляции предсердий в сигналах ЭКГ с использованием сети LSTM. Процедура использует избыточную дискретизацию, чтобы избежать смещения классификации, которое происходит, когда человек пытается обнаружить ненормальные условия в населениях, состоящих в основном из здоровых пациентов. Обучение сети LSTM с использованием данных необработанного сигнала приводит к плохой точности классификации. Настройка сети с помощью двух функций времени-частотного момента для каждого сигнала значительно улучшает эффективность классификации, а также уменьшает время обучения.

Ссылки

[1] AF Classification from a Short Single Lead ECG Recording: the PhysioNet/Computing in Cardiology Challenge, 2017. https://physionet.org/challenge/2017/

[2] Клиффорд, Гари, Чэнъю Лю, Бенджамин Муди, Ли-вэй Х. Лехман, Икаро Силва, Цяо Ли, Алистер Джонсон и Роджер Г. Марк. "Классификация AF из короткой одной ведущей записи ЭКГ Вычисления в кардиологии (Rennes: IEEE). Том 44, 2017, с. 1-4.

[3] Гольдбергер, А. Л., Л. А. Н. Амарал, Л. Гласс, Ж. М. Хаусдорф, П. Ч. Иванов, Р. Г. Марк, Ж. Э. Миетус, Г. Б. Муди, К.-К. Пэн и Х. Э. Стэнли. PhysioBank, PhysioToolkit и PhysioNet: компоненты нового исследовательского ресурса комплексных физиологических сигналов. Циркуляция. Том 101, № 23, 13 июня 2000 года, стр. e215-e220. http://circ.ahajournals.org/content/101/23/e215.full

[4] Понс, Жорди, Томас Лиди и Ксавье Серра. Эксперименты с музыкально мотивированными сверточными нейронными сетями. 14-е Международное рабочее совещание по индексации мультимедиа на основе контента (CBMI). Июнь 2016 года.

[5] Wang, D. «Глубокое обучение reinvents the слуховой аппарат», IEEE Spectrum, Vol. 54, № 3, March 2017, pp. 32-37. doi: 10.1109/MSPEC.2017.7864754.

[6] Браунли, Джейсон. Как масштабировать данные для длинных краткосрочных сетей памяти на Python. 7 июля 2017 года. https://machinelearningmastery.com/how-to-scale-data-for-long-short-term-memory-networks-in-python/.