Введение

Электрическое действие в человеческом сердце может быть измерено как последовательность амплитуд далеко от базового сигнала. Для одного нормального цикла heartbeat сигнал ECG может быть разделен на следующую морфологию удара [1]:

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

Этот пример использует сигналы ECG от общедоступной Базы данных QT 3[] 4[]. Данные состоят примерно из 15 минут записей ECG, с частотой дискретизации 250 Гц, измеренных от в общей сложности 105 пациентов. Чтобы получить каждую запись, ревизоры поместили два электрода в другие места на груди пациента, приводящей к двухканальному сигналу. База данных обеспечивает метки области сигнала, сгенерированные автоматизированной экспертной системой [2]. Этот пример стремится использовать решение для глубокого обучения обеспечить метку для каждой выборки сигнала ECG согласно области, где выборка расположена. Этот процесс маркировки необходимых областей через сигнал часто упоминается как сегментация формы волны.

Чтобы обучить глубокую нейронную сеть классифицировать области сигнала, можно использовать сеть Long Short-Term Memory (LSTM). В этом примере показано, как методы предварительной обработки сигнала и частотно-временной анализ могут использоваться, чтобы улучшать производительность сегментации LSTM. В частности, пример использует synchrosqueezed преобразование Фурье, чтобы представлять неустановившееся поведение сигнала ECG.

Загрузите и подготовьте данные

Каждый канал 105 двухканальных сигналов ECG был помечен независимо автоматизированной экспертной системой и обработан независимо для в общей сложности 210 сигналов ECG, которые хранились вместе с метками области в 210 MAT-файлах. Файлы доступны в следующем местоположении: https://www.mathworks.com/supportfiles/SPT/data/QTDatabaseECGData.zip.

Загрузите файлы данных в свою временную директорию, местоположение которой задано tempdir MATLAB® команда. Если вы хотите поместить файлы данных в папку, отличающуюся от tempdir, измените имя каталога в последующих инструкциях.

% Download the data
dataURL = 'https://www.mathworks.com/supportfiles/SPT/data/QTDatabaseECGData.zip';
datasetFolder = fullfile(tempdir,'QTDataset');
zipFile = fullfile(tempdir,'QTDatabaseECGData.zip');
if ~exist(datasetFolder,'dir')
     websave(zipFile,dataURL);
end

% Unzip the data
unzip(zipFile,tempdir)

unzip операция создает QTDatabaseECGData папка в вашей временной директории с 210 MAT-файлами в нем. Каждый файл содержит сигнал ECG в переменной ecgSignal и таблица области помечает в переменной signalRegionLabels. Каждый файл также содержит частоту дискретизации сигнала в переменной Fs. В этом примере все сигналы имеют частоту дискретизации 250 Гц.

Создайте datastore сигнала, чтобы получить доступ к данным в файлах. Этот пример принимает, что набор данных хранился в вашей временной директории под QTDatabaseECGData папка. Если дело обстоит не так, изменение путь к данным в коде ниже. Задайте имена переменной сигнала, которые вы хотите считать из каждого файла с помощью SignalVariableNames параметр.

sds = signalDatastore(datasetFolder,'SignalVariableNames',["ecgSignal","signalRegionLabels"])

sds = 
  signalDatastore with properties:

                       Files:{
                             '/tmp/QTDataset/ecg1.mat';
                             '/tmp/QTDataset/ecg10.mat';
                             '/tmp/QTDataset/ecg100.mat'
                              ... and 207 more
                             }
    AlternateFileSystemRoots: [0×0 string]
                    ReadSize: 1
         SignalVariableNames: ["ecgSignal"    "signalRegionLabels"]

Datastore возвращает двухэлементный массив ячеек с сигналом ECG, и таблица области помечает каждый раз, когда вы вызываете read функция. Используйте preview функция datastore, чтобы видеть, что содержимое первого файла является 225,000 выборки длинный сигнал ECG и таблица, содержащая 3 385 меток области.

data = preview(sds)

data=2×1 cell array
    {225000×1 double}
    {  3385×2 table }

Посмотрите на первые несколько строк таблицы меток области и заметьте, что каждая строка содержит предельные индексы области и значение класса области (P, T, или QRS).

head(data{2})

ans=8×2 table
    ROILimits     Value
    __________    _____

     83    117     P   
    130    153     QRS 
    201    246     T   
    285    319     P   
    332    357     QRS 
    412    457     T   
    477    507     P   
    524    547     QRS 

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

displayWaveformLabels(data,true,1000)

Обычная процедура классификации машинного обучения следующая:

Сеть обучена с 70% данных и протестирована с остающимися 30%.

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

rng default
[trainIdx,~,testIdx] = dividerand(numel(sds.Files),0.7,0,0.3);
trainDs = subset(sds,trainIdx);
testDs = subset(sds,testIdx);

В этой проблеме сегментации вход к сети LSTM является сигналом ECG, и выход является последовательностью или маской меток с той же длиной как входной сигнал. Сетевая задача состоит в том, чтобы пометить каждую выборку сигнала именем области, которой это принадлежит. Поэтому необходимо преобразовать метки области на наборе данных к последовательностям, содержащим одну метку на выборку сигнала. Используйте преобразованный datastore и roi2mask функция помощника, чтобы преобразовать метки. Предварительно просмотрите преобразованный datastore, чтобы заметить, что он возвращает сигнальный вектор и вектор метки равных длин. Постройте первые 1 000 элементов категориального вектора маски. roi2mask функция добавляет категорию меток, n/a, пометить выборки, которые не принадлежат никакой необходимой области.

trainDs = transform(trainDs, @roi2mask);
testDs = transform(testDs, @roi2mask);

transformedData = preview(trainDs)

transformedData=1×2 cell array
    {224993×1 double}    {224993×1 categorical}

plot(transformedData{2}(1:1000))

Передача очень длинных входных сигналов в сеть LSTM может привести к ухудшению производительности оценки и чрезмерному использованию памяти. Чтобы избежать этих эффектов, повредите сигналы ECG и их соответствующие маски метки с помощью преобразованного datastore и resizeData функция помощника. Функция помощника создает как можно больше сегментов с 5000 выборками и отбрасывает остающиеся выборки. Предварительный просмотр выхода преобразованного datastore показывает, что первый сигнал ECG и его маска метки повреждаются в сегменты с 5000 выборками. Обратите внимание на то, что предварительный просмотр преобразованного datastore только показывает первые 8 элементов в противном случае floor(224993/5000) = 44 массива ячеек элемента, которые закончились бы, если бы мы вызвали datastore read функция.

trainDs = transform(trainDs,@resizeData);
testDs = transform(testDs,@resizeData);
preview(trainDs)

ans=8×2 cell array
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}

Введите необработанные сигналы ECG непосредственно в сеть LSTM

Во-первых, обучите сеть LSTM с помощью необработанных сигналов ECG от обучающего набора данных.

Задайте сетевую архитектуру перед обучением. Задайте sequenceInputLayer из размера 1, чтобы принять одномерные временные ряды. Задайте слой LSTM с 'sequence' режим вывода, чтобы обеспечить классификацию для каждой выборки в сигнале. Используйте 200 скрытых узлов в оптимальной производительности. Задайте fullyConnectedLayer с выходным размером 4, один для каждого из классов формы волны. Добавьте softmaxLayer и classificationLayer выводить предполагаемые метки.

layers = [ ...
    sequenceInputLayer(1)
    lstmLayer(200,'OutputMode','sequence')
    fullyConnectedLayer(4)
    softmaxLayer
    classificationLayer];

Выберите опции для учебного процесса, которые гарантируют хорошую производительность сети. Обратитесь к trainingOptions документация для описания каждого параметра.

options = trainingOptions('adam', ...
    'MaxEpochs',10, ...
    'MiniBatchSize',50, ...
    'InitialLearnRate',0.01, ...
    'LearnRateDropPeriod',3, ...
    'LearnRateSchedule','piecewise', ...
    'GradientThreshold',1, ...
    'Plots','training-progress',...
    'shuffle','every-epoch',...
    'Verbose',0,...
    'DispatchInBackground',true);

Поскольку целый обучающий набор данных умещается в памяти, возможно использовать tall функция datastore, чтобы преобразовать данные параллельно, если Parallel Computing Toolbox™ доступен, и затем собирает его в рабочую область. Обучение нейронной сети является итеративным. В каждой итерации datastore считывает данные из файлов и преобразовывает данные прежде, чем обновить сетевые коэффициенты. Если совпадения данных в память о вашем компьютере, импортируя данные в рабочую область включают более быстрое обучение, потому что данные считаны и преобразованы только однажды. Обратите внимание на то, что, если данные не умещаются в памяти, вы должны, чтобы передать datastore в учебную функцию, и преобразования выполняются в каждую учебную эпоху.

Создайте длинные массивы и для наборов обучающих данных и для наборов тестов. В зависимости от вашей системы количество рабочих в параллельном пуле, который создает MATLAB, может отличаться.

tallTrainSet = tall(trainDs);

Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 8).

tallTestSet = tall(testDs);

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

 trainData = gather(tallTrainSet);

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete

- Pass 1 of 1: Completed in 28 sec
Evaluation completed in 28 sec

 trainData(1,:)

ans=1×2 cell array
    {1×5000 double}    {1×5000 categorical}

 testData = gather(tallTestSet);

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 9.8 sec
Evaluation completed in 10 sec

Обучение сети

Используйте trainNetwork команда, чтобы обучить сеть LSTM. Из-за большого размера набора данных, этот процесс может занять несколько минут. Если ваша машина имеет графический процессор и Parallel Computing Toolbox™, то MATLAB автоматически использует графический процессор в обучении. В противном случае это использует центральный процессор.

Учебные подграфики точности и потери в фигуре отслеживают процесс обучения через все итерации. Используя необработанные данные сигнала, сеть правильно классифицирует приблизительно 77% выборок как принадлежащий волне P, комплексу QRS, волне T или N/A.

net = trainNetwork(trainData(:,1),trainData(:,2),layers,options);

Классифицируйте данные о тестировании

Классифицируйте данные о тестировании с помощью обученной сети LSTM и classify команда. Задайте мини-пакетный размер 50, чтобы совпадать с опциями обучения.

predTest = classify(net,testData(:,1),'MiniBatchSize',50);

Матрица беспорядка обеспечивает интуитивное, и информативное означает визуализировать производительность классификации. Используйте confusionchart команда, чтобы вычислить полную точность классификации для предсказаний данных о тестировании. Для каждого входа преобразуйте массив ячеек категориальных меток к вектору-строке. Задайте нормированное столбцом отображение, чтобы просмотреть результаты как проценты выборок для каждого класса.

confusionchart([predTest{:}],[testData{:,2}],'Normalization','column-normalized');

Используя необработанный сигнал ECG, как введено к сети, только приблизительно 60% выборок T-волны, 40% выборок P-волны и 60% QRS-комплексных выборок были правильны. Чтобы улучшать производительность, примените некоторое знание характеристик сигнала ECG до входа к нейронной сети для глубокого обучения, например, базовая линия, блуждающая вызванный дыхательным движением пациента.

Примените методы фильтрации, чтобы удалить базовое блуждание и высокочастотный шум

Три морфологии удара занимает различные диапазоны частот. Спектр комплекса QRS обычно имеет центральную частоту приблизительно 10-25 Гц, и ее компоненты лежат ниже 40 Гц. P и T-волны происходят на еще более низких частотах: P-волновые-компоненты ниже 20 Гц, и T-волновые-компоненты ниже 10 Гц [5].

Базовое блуждание является низкочастотным (<0,5 Гц) колебание, вызванное движением дыхания пациента. Это колебание независимо от морфологии удара и не предоставляет значимую информацию [6].

Спроектируйте полосовой фильтр с частотным диапазоном полосы пропускания [0.5, 40] Гц, чтобы удалить блуждание и любой высокочастотный шум. Удаление этих компонентов улучшает обучение LSTM, потому что сеть не изучает несоответствующие функции. Используйте cellfun на высоких массивах ячейки данных, чтобы отфильтровать набор данных параллельно.

% Bandpass filter design
hFilt = designfilt('bandpassiir', 'StopbandFrequency1',0.4215,'PassbandFrequency1', 0.5, ...
    'PassbandFrequency2',40,'StopbandFrequency2',53.345,...
    'StopbandAttenuation1',60,'PassbandRipple',0.1,'StopbandAttenuation2',60,...
    'SampleRate',250,'DesignMethod','ellip');

% Create tall arrays from the transformed datastores and filter the signals
tallTrainSet = tall(trainDs);
tallTestSet = tall(testDs);

filteredTrainSignals = gather(cellfun(@(x)filter(hFilt,x),tallTrainSet(:,1),'UniformOutput',false));

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete

- Pass 1 of 1: Completed in 19 sec
Evaluation completed in 19 sec

trainLabels = gather(tallTrainSet(:,2));

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 19 sec
Evaluation completed in 19 sec

filteredTestSignals = gather(cellfun(@(x)filter(hFilt,x),tallTestSet(:,1),'UniformOutput',false));

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 8.5 sec
Evaluation completed in 8.6 sec

testLabels = gather(tallTestSet(:,2));

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 8.3 sec
Evaluation completed in 8.5 sec

Постройте сырые данные и отфильтрованные сигналы для типичного случая.

trainData = gather(tallTrainSet);

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 19 sec
Evaluation completed in 19 sec

figure
subplot(2,1,1)
plot(trainData{95,1}(2001:3000))
title('Raw')
grid
subplot(2,1,2)
plot(filteredTrainSignals{95}(2001:3000))
title('Filtered')
grid

Обучите сеть с фильтрованными сигналами ECG

Обучите сеть LSTM на отфильтрованных сигналах ECG с помощью той же сетевой архитектуры как прежде.

filteredNet = trainNetwork(filteredTrainSignals,trainLabels,layers,options);

Предварительная обработка сигналов улучшает учебную точность до лучше, чем 80%.

Классифицируйте фильтрованные сигналы ECG

Классифицируйте предварительно обработанные тестовые данные с обновленной сетью LSTM.

predFilteredTest = classify(filteredNet,filteredTestSignals,'MiniBatchSize',50);

Визуализируйте производительность классификации как матрицу беспорядка.

figure
confusionchart([predFilteredTest{:}],[testLabels{:}],'Normalization','column-normalized');

Простая предварительная обработка улучшает классификацию T-волн приблизительно на 15% и QRS-комплексную классификацию на 10%. Точность классификации P-волн не изменилась.

Представление частоты времени сигналов ECG

Общий подход для успешной классификации данных timeseries должен извлечь функции частоты времени и накормить ими сеть вместо исходных данных. Сеть затем изучает шаблоны через время и частоту одновременно [7].

Synchrosqueezed преобразование Фурье (FSST) вычисляет спектр частоты для каждой выборки сигнала, таким образом, это идеально для проблемы сегментации под рукой, где мы должны обеспечить то же разрешение времени как исходные сигналы. Используйте fsst функция, чтобы смотреть преобразование одного из учебных сигналов. Задайте окно Кайзера длины 128, чтобы обеспечить соответствующее разрешение частоты.

data =  preview(trainDs);
figure
fsst(data{1,1},250,kaiser(128),'yaxis')

Вычислите FSST каждого сигнала в обучающем наборе данных по частотному диапазону интереса, [0.5, 40] Гц. Обработайте действительные и мнимые части FSST как отдельные функции и подайте оба компонента в сеть. Кроме того, стандартизируйте учебные функции путем вычитания среднего значения и деления на стандартное отклонение. Используйте преобразованный datastore, extractFSSTFeatures функция помощника и tall функция, чтобы обработать данные параллельно.

fsstTrainDs = transform(trainDs,@(x)extractFSSTFeatures(x,250));
fsstTallTrainSet = tall(fsstTrainDs);
fsstTrainData = gather(fsstTallTrainSet);

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete

- Pass 1 of 1: Completed in 2 min 35 sec
Evaluation completed in 2 min 35 sec

Повторите эту процедуру для данных о тестировании.

fsstTTestDs = transform(testDs,@(x)extractFSSTFeatures(x,250));
fsstTallTestSet = tall(fsstTTestDs);
fsstTestData = gather(fsstTallTestSet);

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 1 min 23 sec
Evaluation completed in 1 min 23 sec

Настройте сетевую архитектуру

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

size(fsstTrainData{1,1})

ans = 1×2

          40        5000

Задайте sequenceInputLayer из 40 входных функций. Сохраните остальную часть сетевых параметров неизменной.

layers = [ ...
    sequenceInputLayer(40)
    lstmLayer(200,'OutputMode','sequence')
    fullyConnectedLayer(4)
    softmaxLayer
    classificationLayer];

Обучите сеть с FSST сигналов ECG

Обучите обновленную сеть LSTM с преобразованным набором данных.

fsstNet = trainNetwork(fsstTrainData(:,1),fsstTrainData(:,2),layers,options);

Использование функций частоты времени улучшает учебную точность, которая теперь превышает 90%.

Классифицируйте тестовые данные с FSST

Используя обновленную сеть LSTM и извлеченные функции FSST, классифицируйте данные о тестировании.

predFsstTest = classify(fsstNet,fsstTestData(:,1),'MiniBatchSize',50);

Визуализируйте производительность классификации как матрицу беспорядка.

confusionchart([predFsstTest{:}],[fsstTestData{:,2}],'Normalization','column-normalized');

Используя частоту времени представление улучшает классификацию T-волн приблизительно на 20%, классификацию P-волн приблизительно на 40% и QRS-комплексную классификацию на 30%, когда по сравнению с необработанными данными заканчивается.

Используйте displayWaveformLabels функция помощника, чтобы сравнить сетевое предсказание с основной истиной помечает для одного сигнала ECG.

testData = gather(tall(testDs));

Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 8.3 sec
Evaluation completed in 8.5 sec

figure
subplot(2,1,1)
displayWaveformLabels({testData{1,1}(3000:4000),testData{1,2}(3000:4000)},false)
title('Ground Truth')

subplot(2,1,2)
displayWaveformLabels({testData{1,1}(3000:4000),predFsstTest{1}(3000:4000)},false)
title('Predicted')

Заключение

Этот пример показал, как предварительная обработка сигнала и частотно-временной анализ могут улучшать производительность сегментации формы волны LSTM. Фильтрация полосы пропускания и основанный на Фурье synchrosqueezing приводят к среднему улучшению через все выходные классы от 55% приблизительно до 85%.

Ссылки

[1] Макшарри, Патрик Э., и др. "Динамическая модель для генерации синтетических сигналов электрокардиограммы". IEEE® Transactions на Биоинженерии. Издание 50, № 3, 2003, стр 289–294.

[2] Laguna, Пабло, Рэймон Джейне и Пере Каминаль. "Автоматическое обнаружение контуров волны в мультиведущих сигналах ECG: Валидация с базой данных CSE". Компьютеры и Биомедицинское Исследование. Издание 27, № 1, 1994, стр 45–60.

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

[4] Laguna, Пабло, Роджер Г. Марк, Ари Л. Голдбергер и Джордж Б. Муди. "База данных для Оценки Алгоритмов для Измерения QT и Других Интервалов Формы волны в ECG". Компьютеры в Кардиологии. Vol.24, 1997, стр 673–676.

[5] Sörnmo, Леиф и Пабло Лагуна. "Электрокардиограмма (ECG) обработка сигналов". Вайли Энкиклопедия Биоинженерии, 2006.

[6] Колер, B-U., Карстен Хенниг и Райнхольд Орглмайстер. "Принципы обнаружения программного обеспечения QRS". Разработка IEEE в Журнале Медицины и Биологии. Издание 21, № 1, 2002, стр 42–57.

[7] Salamon, Джастин и Хуан Пабло Белло. "Глубокие сверточные нейронные сети и увеличение данных для экологической звуковой классификации". Буквы Обработки сигналов IEEE. Издание 24, № 3, 2017, стр 279–283.