Описание данных

Этот пример использует данные ЭКГ, полученные из трех групп или классов людей: людей с сердечной аритмией, людей с застойным сердечным отказом и людей с нормальными синусовыми ритмами. В примере используются 162 записи ЭКГ из трех баз данных PhysioNet: База данных аритмии MIT-BIH [3] [5], База данных нормального синусового ритма MIT-BIH [3] и База данных застойной сердечной недостаточности BIDMC [2] [ Всего насчитывается 96 записей от лиц с аритмией, 30 записей от лиц с застойным сердечным отказом и 36 записей от лиц с нормальными синусовыми ритмами. Цель состоит в том, чтобы обучить классификатор для различения аритмии (ARR), застойного сердечного отказа (CHF) и нормального синусового ритма (NSR).

Загрузка данных

Первый шаг - загрузка данных из репозитория GitHub. Чтобы загрузить данные, нажмите Code и выберите Download ZIP. Сохраните файл physionet_ECG_data-main.zip в папке, в которой у вас есть разрешение на запись. Инструкции для этого примера предполагают, что вы загрузили файл во временную директорию, (tempdir в MATLAB). Измените последующие инструкции для распаковки и загрузки данных, если вы решите загрузить данные в папку, отличную от tempdir.

Файл physionet_ECG_data-main.zip содержит

и ECGData.zip содержит

ECGData.mat содержит данные, используемые в этом примере. Файл .txt, Modified_physionet_data.txt, требуется политикой копирования PhysioNet и предоставляет атрибуты источника для данных, а также описание шагов предварительной обработки, применяемых к каждой записи ЭКГ.

Загрузка файлов

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

unzip(fullfile(tempdir,'physionet_ECG_data-main.zip'),tempdir)
unzip(fullfile(tempdir,'physionet_ECG_data-main','ECGData.zip'),...
    fullfile(tempdir,'ECGData'))

После разархивации файла ECGData.zip загрузите данные в MATLAB.

load(fullfile(tempdir,'ECGData','ECGData.mat'))

ECGData массив структур с двумя полями: Data и Labels. Data представляет собой 162 на 65536 матрицу, где каждая строка является дискретизацией записи ЭКГ в 128 герц. Каждые временные ряды ЭКГ имеют общую длительность 512 секунд. Labels - массив ячеек 162 на 1 с диагностическими метками, по одному для каждой строки Data. Три диагностические категории: 'ARR' (аритмия), 'CHF' (застойный сердечный отказ) и 'NSR' (нормальный синусовый ритм).

Создайте обучающие и тестовые данные

Случайным образом разделите данные на два набора - обучающие и тестовые данные. Функция помощника helperRandomSplit выполняет случайное разделение. helperRandomSplit принимает требуемый процент разделения для обучающих данных и ECGData. The helperRandomSplit функция выводит два набора данных вместе с набором меток для каждого. Каждая строка trainData и testData является сигналом ЭКГ. Каждый элемент trainLabels и testLabels содержит метку класса для соответствующей строки матриц данных. В этом примере мы случайным образом присваиваем 70% процентов данных в каждом классе набору обучающих данных. Оставшиеся 30% удерживаются для проверки (предсказания) и присваиваются тестовому набору.

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

В trainData 113 записей установил и 49 записи в testData. По проекту обучающие данные содержат 69,75% (113/162) данных. Напомним, что класс ARR представляет 59,26% данных (96/162), класс CHF представляет 18,52% (30/162), и класс NSR представляет 22,22% (36/162). Исследуйте процент каждого класса в наборах обучения и тестов. Проценты в каждом соответствуют общим процентам классов в наборе данных.

Ctrain = countcats(categorical(trainLabels))./numel(trainLabels).*100

Ctrain = 3×1

   59.2920
   18.5841
   22.1239

Ctest = countcats(categorical(testLabels))./numel(testLabels).*100

Ctest = 3×1

   59.1837
   18.3673
   22.4490

График выборок

Постройте первые несколько тысяч выборок четырех случайным образом выбранных записей из ECGData. Функция помощника helperPlotRandomRecords делает это. helperPlotRandomRecords принимает ECGData и случайный seed как вход. Начальный seed устанавливается равным 14, так что строится по меньшей мере одна запись от каждого класса. Можно выполнять helperPlotRandomRecords с ECGData как единственный входной параметр столько раз, сколько вы хотите получить представление о разнообразии форм волны ЭКГ, сопоставленных с каждым классом. Исходный код для этого и всех вспомогательных функций можно найти в разделе Вспомогательные функции в конце этого примера.

helperPlotRandomRecords(ECGData,14)

Вейвлет вейвлет-времени

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

N = size(ECGData.Data,2);
sn = waveletScattering('SignalLength',N, 'InvarianceScale',150,'SamplingFrequency',128);

Можно визуализировать вейвлет в двух банках фильтров с помощью следующего.

[fb,f,filterparams] = filterbank(sn);
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');

Чтобы продемонстрировать шкалу инвариации, получите обратное преобразование Фурье функции масштабирования и центрируйте ее во времени 0 секунд. Две вертикальные черные линии отмечают контуры -75 и 75 секунд. Кроме того, постройте график действительной и мнимой частей самого грубого вейвлета (самой низкой частоты) от первой группы фильтров. Обратите внимание, что вейвлет самой грубой шкалы не превышает инвариантную шкалу, определяемую временной поддержкой функции масштабирования. Это важное свойство вейвлет времени.

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'})
hold off

После построения сети рассеяния получите коэффициенты рассеяния для обучающих данных в виде матрицы. Когда бегаешь featureMatrix с несколькими сигналами каждый столбец обрабатывается одним сигналом.

scat_features_train = featureMatrix(sn,trainData');

Область выхода featureMatrix при этом 409 на 16 на 113. Каждая страница тензора, scat_features_train, является преобразованием рассеяния одного сигнала. Вейвлет критически уменьшается по времени на основе полосы пропускания функции масштабирования. В этом случае это приводит к 16 временным окнам для каждого из 409 путей рассеяния.

В порядок получения матрицы, совместимой с классификатором SVM, измените форму многосигнального преобразования рассеяния на матрицу, где каждый столбец соответствует пути рассеяния, и каждая строка является временным окном рассеяния. В этом случае вы получаете 1808 строк, потому что для каждого из 113 сигналов в обучающих данных существует 16 временных окон.

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),[]);

Повторите процесс для тестовых данных. Первоначально, scat_features_test 409 на 16 на 49, потому что в наборе обучающих данных 49 сигналов ЭКГ. После изменения модели классификатора SVM матрица функции равна 784 на 416.

scat_features_test = featureMatrix(sn,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),[]);

Потому что для каждого сигнала мы получили 16 окон рассеяния, нам нужно создать метки, чтобы соответствовать количеству окон. Функция помощника createSequenceLabels делает это на основе количества окон.

[sequence_labels_train,sequence_labels_test] = createSequenceLabels(Nwin,trainLabels,testLabels);

Перекрестная валидация

Для классификации выполняются два анализа. Первый использует все данные рассеяния и подбирает мультикласс SVM с квадратичным ядром. Всего в набор данных в целом 2592 последовательности рассеяния, 16 для каждого из 162 сигналов. Вероятность ошибки, или потери, оценивается с помощью 5-кратной перекрестной валидации.

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

Вычислите потерю и матрицу неточностей. Отобразите точность.

predLabels = kfoldPredict(kfoldmodel);
loss = kfoldLoss(kfoldmodel)*100;
confmatCV = confusionmat(allLabels_scat,predLabels)

confmatCV = 3×3

        1535           0           1
           1         479           0
           0           0         576

fprintf('Accuracy is %2.2f percent.\n',100-loss);

Accuracy is 99.92 percent.

Точность составляет 99,88%, что довольно хорошо, но фактические результаты лучше, учитывая, что здесь каждое временное окно классифицируется отдельно. Для каждого сигнала существует 16 отдельных классификаций. Используйте простое большинство голосов, чтобы получить одно предсказание класса для каждого представления рассеяния.

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

Определите фактическую точность перекрестной проверки на основе режима предсказаний классов для каждого набора временных окон рассеяния. Если режим не уникален для заданного аппарата, helperMajorityVote возвращает ошибку классификации, обозначенную 'NoUniqueMode'. Это приводит к дополнительному столбцу в матрице неточностей, которая в этом случае является нулями, потому что для каждого набора предсказаний рассеяния существует уникальный режим.

CVaccuracy = sum(eq(ClassVotes,categorical([trainLabels; testLabels])))/162*100;
fprintf('True cross-validation accuracy is %2.2f percent.\n',CVaccuracy);

True cross-validation accuracy is 100.00 percent.

MVconfmatCV = confusionmat(categorical([trainLabels; testLabels]),ClassVotes);
MVconfmatCV

MVconfmatCV = 4×4

    96     0     0     0
     0    30     0     0
     0     0    36     0
     0     0     0     0

Рассеяние правильно классифицировало все сигналы в перекрестной проверенной модели. Если вы осмотрите ClassCountsВы видите, что 2 неправильно классифицированных временных окна в confmatCV относятся к 2 сигналам, где 15 из 16 окон рассеяния были классифицированы правильно.

Классификация SVM

Для следующего анализа мы подбираем мультикласс квадратичный SVM только к обучающим данным (70%), а затем используем эту модель, чтобы получить предсказания на 30% данных, сохраненных для проверки. В тестовом наборе 49 записей данных. Используйте большинство голосов для отдельных окон рассеяния.

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

The test accuracy is 97.96 percent. 

confusionchart(categorical(testLabels),TestVotes)

Точность классификации на тестовом наборе данных составляет приблизительно 98%. Матрица неточностей показывает, что одна запись CHF неправильно классифицирована как ARR. Все 48 других сигналов правильно классифицированы.

Сводные данные

Этот пример использовал рассеяние вейвлета времени и классификатор SVM, чтобы классифицировать формы волны ECG в один из трех диагностических классов. Вейвлет оказалось мощным экстрактором функций, который требовал только минимального набора заданных пользователем параметров, чтобы получить набор устойчивых функций для классификации. Сравните это с примером Signal Classification Using Wavelet-Based Функций and Машин опорных векторов, который требовал значительных знаний для функций ручной работы для использования в классификации. При вейвлет-временном рассеянии от вас требуется только задать шкалу временной инвариации, количество банков фильтров (или вейвлет-преобразования) и количество вейвлеты на октаву. Комбинация преобразования вейвлет и классификатора SVM дала 100% классификацию по перекрестно проверенной модели и 98% правильную классификацию при применении SVM к преобразованиям рассеяния тестового набора задержки.

Ссылки

Вспомогательные функции

helperPlotRandomRecords Строит четыре сигнала ЭКГ, случайным образом выбранных из ECGData.

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

helperMajorVote Находит режим в предсказанных метках классов для каждого набора временных окон рассеяния. Функция возвращает как режимы меток классов, так и количество предсказаний классов для каждого набора временных окон рассеяния. Если уникального режима нет, helperMajorityVote возвращает метку класса «error», чтобы указать, что набор окон рассеяния является ошибкой классификации.

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