Введение

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

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

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

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

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

Загрузите и исследуйте данные

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

ReadPhysionetData
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 функция, чтобы видеть, что существует 738 сигналов AFib и 5 050 Нормальных сигналов.

summary(Labels)

     A       738 
     N      5050 

Сгенерируйте гистограмму длин сигнала. Заметьте, что большинство сигналов является 9 000 выборок долго.

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

Визуализируйте сегмент одного сигнала от каждого класса. Heartbeat AFib растянуты в неправильных интервалах, в то время как Нормальные heartbeat регулярно происходят. Сигналы heartbeat AFib также часто испытывают недостаток в волне P, какие импульсы, прежде чем QRS объединяют в Нормальном сигнале heartbeat. График Нормального сигнала показывает волну 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 функционируйте к сигналам ECG, таким образом, они - все 9 000 выборок долго. Функция игнорирует сигналы меньше чем с 9 000 выборок. Если сигнал имеет больше чем 9 000 выборок, segmentSignals пропуски это в как можно больше сегментов с 9000 выборками и игнорирует остающиеся выборки. Например, сигнал с 18 500 выборками становится двумя сигналами с 9000 выборками, и остающиеся 500 выборок проигнорированы.

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

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

Signals(1:5)

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

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

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

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

summary(Labels)

     A       718 
     N      4937 

Поскольку 87,3% сигналов Нормален, классификатор узнал бы, что это может достигнуть высокой точности просто путем классификации всех сигналов как Нормальные. Чтобы избежать этого смещения, увеличьте данные AFib путем дублирования сигналов AFib в наборе данных так, чтобы было то же количество сигналов 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 и 4 443 Нормальных сигнала для обучения. Чтобы достигнуть того же количества сигналов в каждом классе, используйте первые 4 438 Нормальных сигналов, и затем используйте repmat повторить первые 634 AFib сигнализирует семь раз.

Для тестирования существует 72 сигнала AFib и 494 Нормальных сигнала. Используйте первые 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);];

Распределение между Нормальным и сигналами 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 проходит через обучающие данные. 'MiniBatchSize' из 150 направляет сеть, чтобы посмотреть на 150 учебных сигналов за один раз. 'InitialLearnRate' из 0,01 помогает ускорить учебный процесс. Задайте 'SequenceLength' из 1 000, чтобы повредить сигнал в мелкие кусочки так, чтобы машина не исчерпывала память путем рассмотрения слишком большого количества данных одновременно. Установите 'GradientThreshold'к 1, чтобы стабилизировать учебный процесс, препятствуя тому, чтобы градиенты стали слишком большим. Задайте 'Plots' как 'training-progress' сгенерировать графики, которые показывают диаграмму процесса обучения как количество увеличений итераций. Установите 'Verbose' к false подавить таблицу выход, который соответствует данным, показанным в графике. Если вы хотите видеть эту таблицу, установите 'Verbose' к true.

Этот пример использует адаптивную оценку момента (ADAM) решатель. ADAM выполняет лучше с RNNs как LSTMs, чем стохастический градиентный спуск по умолчанию с импульсом (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);

В проблемах классификации матрицы беспорядка используются, чтобы визуализировать производительность классификатора на наборе данных, которыми известны истинные значения. Целевой Класс является меткой основной истины сигнала, и Выходной класс является меткой, присвоенной сигналу сетью. Метки осей представляют метки класса, AFib (A) и Нормальный (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');

Матрица беспорядка показывает, что 81,7% основной истины, сигналы AFib правильно классифицируются как AFib, в то время как 31,1% основной истины Нормальные сигналы правильно классифицируются как Нормальные. Кроме того, 54,2% сигналов, классифицированных как AFib, является на самом деле AFib, и 63,0% сигналов, классифицированных, когда Нормальный на самом деле Нормальны. Полная учебная точность составляет 56,4%.

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

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');

Эта матрица беспорядка похожа на учебную матрицу беспорядка. Полная точность тестирования между 50% к 60%.

Улучшайте производительность с извлечением признаков

Извлечение признаков из данных может помочь улучшить точность обучения и тестирования классификатора. Решить, какие функции извлечь, этот пример адаптирует подход, который вычисляет изображения частоты времени, такие как спектрограммы, и использует их, чтобы обучить сверточные нейронные сети (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 во временном интервале:

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);

Визуализируйте формат новых входных параметров. Каждая ячейка больше не содержит одни 9 000 выборок, долго сигнализируют; теперь это содержит два, 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-выигрыш, является популярным способом улучшать производительность сети во время обучения.

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);

Существует большое улучшение учебной точности, которая теперь больше 90%. Тренды потери перекрестной энтропии к 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');

Заключение

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

Ссылки

[1] Классификация AF от Короткой Одной Ведущей Записи ECG: Физиосетевое/вычислительное в проблеме Кардиологии, 2017. https://physionet.org/challenge/2017/

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

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

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

[5] Ван, D. "Глубокое обучение заново изобретает слуховой аппарат", Спектр IEEE, Издание 54, № 3, март 2017, стр 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/.