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

В этом примере используются данные ЭКГ, полученные из трех групп или классов людей: людей с сердечной аритмией, людей с застойной сердечной недостаточностью и людей с нормальным синусовым ритмом. В примере используются 162 записи ЭКГ из трех баз данных PhysioNet: база данных аритмии MIT-BIH [3] [5], база данных нормального синусового ритма MIT-BIH [3] и база данных застойной сердечной недостаточности BIDMC [2] [3]. Всего насчитывается 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. helperRandomSplit функция выводит два набора данных вместе с набором меток для каждого. Каждая строка trainData и testData - сигнал ЭКГ. Каждый элемент trainLabels и testLabels содержит метку класса для соответствующей строки матриц данных. В этом примере мы случайным образом присваиваем обучающему набору 70% данных в каждом классе. Оставшиеся 30% удерживаются для тестирования (прогнозирования) и назначаются тестовому набору.

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

Есть 113 записей в trainData набор и 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 и случайное начальное число в качестве входных данных. Начальное начальное число устанавливается равным 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 строк, поскольку в данных обучения имеется 16 временных окон для каждого из 113 сигналов.

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 сигналов. Частота ошибок или потерь оценивается с использованием пятикратной перекрестной проверки.

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-классификатор для классификации форм сигналов ЭКГ в один из трех диагностических классов. Вейвлет-рассеяние оказалось мощным экстрактором признаков, который требовал только минимального набора пользовательских параметров для получения набора надежных признаков для классификации. Сравните это с примером Классификация сигналов с использованием вейвлет-основанных функций и поддержки векторных машин, которые требовали значительного количества опыта для рукоделия для использования в классификации. При вейвлет-рассеянии требуется только указать масштаб инвариантности времени, количество блоков фильтров (или вейвлет-преобразований) и количество вейвлетов на октаву. Комбинация преобразования вейвлет-рассеяния и классификатора 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

helPerPotionVote Находит режим в прогнозируемых метках классов для каждого набора временных окон рассеяния. Функция возвращает как режимы метки класса, так и количество предсказаний класса для каждого набора временных окон рассеяния. При отсутствии уникального режима 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