В этом примере показано, как классифицировать данные электрокардиограммы сердцебиения (ЭКГ) из вызова PhysioNet 2017 с использованием глубокого обучения и обработки сигналов. В частности, в примере используются сети долговременной кратковременной памяти и частотно-временной анализ.
ЭКГ регистрируют электрическую активность сердца человека в течение определенного периода времени. Врачи используют ЭКГ, чтобы визуально обнаружить, является ли сердцебиение пациента нормальным или нерегулярным.
Мерцательная аритмия (AFib) - это тип нерегулярного сердцебиения, которое возникает, когда верхние камеры сердца, предсердия, бьются из координации с нижними камерами, желудочками.
В этом примере используются данные ЭКГ из PhysioNet 2017 Challenge [1], [2], [3], которые доступны по адресу https://physionet.org/challenge/2017/. Данные состоят из набора ЭКГ-сигналов, дискретизированных на частоте 300 Гц и разделенных группой экспертов на четыре различных класса: Normal (N), AFib (A), Other Rhythm (O) и Noisy Recording (~). В этом примере показано, как автоматизировать процесс классификации с помощью глубокого обучения. Процедура исследует двоичный классификатор, который может отличать сигналы нормальной ЭКГ от сигналов, имеющих признаки AFib.
В этом примере используются сети с длительной кратковременной памятью (LSTM), тип рецидивирующей нейронной сети (RNN), хорошо подходящей для изучения последовательности и данных временных рядов. Сеть LSTM может изучать долгосрочные зависимости между временными шагами последовательности. Уровень LSTM (lstmLayer (Deep Learning Toolbox)) может смотреть на временную последовательность в прямом направлении, в то время как двунаправленный уровень LSTM (bilstmLayer (Deep Learning Toolbox)) может смотреть на временную последовательность как в прямом, так и в обратном направлениях. В этом примере используется двунаправленный уровень LSTM.
Чтобы ускорить процесс обучения, выполните этот пример на машине с графическим процессором. Если компьютер имеет графический процессор и Toolbox™ параллельных вычислений, MATLAB ® автоматически использует графический процессор для обучения; в противном случае он использует ЦП.
Запустить ReadPhysionetData сценарий для загрузки данных с веб-сайта PhysioNet и создания MAT-файла (PhysionetData.mat), который содержит сигналы ЭКГ в соответствующем формате. Загрузка данных может занять несколько минут. Использовать условный оператор, выполняющий сценарий, только если PhysionetData.mat не существует в текущей папке.
if ~isfile('PhysionetData.mat') ReadPhysionetData end load PhysionetData
Операция загрузки добавляет в рабочую область две переменные: Signals и Labels. Signals - матрица ячеек, которая содержит сигналы ЭКГ. 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 и нормальных сигналов, содержащихся в данных.
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 разнесены с нерегулярными интервалами, в то время как нормальные сердцебиения происходят регулярно. Сигналы пульса AFib также часто не имеют P-волны, которая подает импульсы перед комплексом QRS в сигнале пульса 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 к нормальным сигналам составляет 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 нормальных сигналов, а затем используйте 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 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'к 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 с указанными возможностями обучения и архитектурой уровней с помощью trainNetwork. Поскольку тренировочный набор большой, тренировочный процесс может занять несколько минут.
net = trainNetwork(XTrain,YTrain,layers,options);
Верхний подграф графика training-progress представляет точность обучения, которая является точностью классификации для каждой мини-партии. Когда обучение проходит успешно, это значение обычно увеличивается до 100%. Нижняя часть графика отображает потери при обучении, которые представляют собой потери при перекрестной энтропии для каждой мини-партии. Когда обучение проходит успешно, это значение обычно уменьшается до нуля.
Если тренировка не сходится, графики могут колебаться между значениями без тренда в определенном восходящем или нисходящем направлении. Это колебание означает, что точность тренировки не улучшается и потери тренировки не уменьшаются. Эта ситуация может возникнуть с начала обучения, или участки могут стать плато после некоторого предварительного улучшения точности обучения. Во многих случаях изменение вариантов обучения может помочь сети достичь сходимости. Уменьшение MiniBatchSize или уменьшение InitialLearnRate может привести к увеличению времени обучения, но это может помочь сети лучше учиться.
Точность обучения классификатора колеблется между примерно 50% и примерно 60%, а в конце 10 эпох на тренировку уже ушло несколько минут.
Вычислите точность обучения, которая представляет точность классификатора по сигналам, по которым он обучался. Во-первых, классифицируйте данные обучения.
trainPred = classify(net,XTrain,'SequenceLength',1000);В проблемах классификации матрицы путаницы используются для визуализации производительности классификатора на наборе данных, для которых известны истинные значения. Целевой класс является меткой «земля-истина» сигнала, а выходной класс - меткой, назначенной сигналу сетью. Метки осей представляют метки класса, 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');
Извлечение признаков из данных может помочь улучшить точность обучения и тестирования классификатора. Чтобы решить, какие функции извлекать, в этом примере адаптируется подход, который вычисляет частотно-временные изображения, такие как спектрограммы, и использует их для обучения сверточных нейронных сетей (CNN) [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')
Поскольку в этом примере вместо CNN используется LSTM, важно перевести подход таким образом, чтобы он применялся к одномерным сигналам. Моменты временной частоты (TF) извлекают информацию из спектрограмм. Каждый момент может использоваться в качестве одномерной функции для ввода в LSTM.
Изучите два момента TF во временной области:
Мгновенная частота (instfreq)
Спектральная энтропия (pentropy)
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);
Спектральная энтропия измеряет, насколько острым является спектр сигнала. Сигнал с острым спектром, как и сумма синусоид, имеет низкую спектральную энтропию. Сигнал с плоским спектром, как белый шум, имеет высокую спектральную энтропию. 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 выборок; теперь он содержит два элемента длиной 255 образцов.
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-scoring, является популярным способом повышения производительности сети во время обучения.
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
Теперь, когда каждый сигнал имеет два измерения, необходимо изменить архитектуру сети, указав размер входной последовательности как 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 с указанными возможностями обучения и архитектурой уровней с помощью 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 из короткой записи ECG с одним свинцом: физиоNet/Computing in Cardiology Challenge, 2017. https://physionet.org/challenge/2017/
[2] Клиффорд, Гари, Чэнъю Лю, Бенджамин Муди, Ли-вэй Х. Леман, Икаро Силва, Цяо Ли, Алистер Джонсон и Роджер Г. Марк. "AF-классификация из короткой записи ЭКГ с одиночным лидом: физиология Вычисления в кардиологии (Ренн: 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] Ван, Д. «Глубокое обучение заново изобретает слуховой аппарат», IEEE Spectrum, том 54, № 3, март 2017, стр. 32-37. дои: 10.1109/MSPEC.2017.7864754.
[6] Браунли, Джейсон. Как масштабировать данные для сетей долговременной памяти в Python. 7 июля 2017 года. https://machinelearningmastery.com/how-to-scale-data-for-long-short-term-memory-networks-in-python/.
instfreq | pentropy | bilstmLayer (инструментарий глубокого обучения) | lstmLayer (инструментарий для глубокого обучения) | trainingOptions (инструментарий для глубокого обучения) | trainNetwork (инструментарий для глубокого обучения)