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

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

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

Первый шаг должен загрузить данные из репозитория GitHub. Чтобы загрузить данные, нажмите Clone or download и выберите Download ZIP. Сохраните файл physionet_ECG_data-master.zip в папке, где у вас есть разрешение записи. Инструкции для этого примера принимают, что вы загрузили файл на свою временную директорию, (tempdir в MATLAB). Измените последующие инструкции для того, чтобы разархивировать и загрузить данные, если вы принимаете решение загрузить данные в папке, отличающейся от tempdir. Если вы знакомы с Git, можно загрузить последнюю версию инструментов (Git) и получить данные из системной командной строки с помощью git clone https://github.com/mathworks/physionet_ECG_data/ .

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

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

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

Загрузите файлы

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

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

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

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

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

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

Случайным образом разделите данные в два набора - наборы тестовых данных и обучение. Функция помощника helperRandomSplit выполняет случайное разделение. helperRandomSplit принимает желаемый процент разделения для обучающих данных и ECGData. helperRandomSplit функционируйте выходные параметры два набора данных наряду с набором меток для каждого. Каждая строка trainData и testData сигнал ECG. Каждый элемент 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), класс швейцарского франка представляет 18,52% (30/162), и класс NSR представляет 22,22% (36/162). Исследуйте процент каждого класса в наборах обучающих данных и наборах тестов. Проценты в каждом сопоставимы с полными процентами класса в наборе данных.

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

Ctrain =

   59.2920
   18.5841
   22.1239


Ctest =

   59.1837
   18.3673
   22.4490

Постройте выборки

Постройте первую несколько тысяч выборок четырех случайным образом выбранных записей от ECGData. Функция помощника helperPlotRandomRecords делает это. helperPlotRandomRecords принимает ECGData и случайный seed, как введено. Начальный seed установлен в 14 так, чтобы по крайней мере одна запись от каждого класса была построена. Можно выполнить helperPlotRandomRecords с ECGData как единственный входной параметр так много раз, как вы хотите получить смысл множества форм волны ECG, сопоставленных с каждым классом. Можно найти исходный код для этого и всех функций помощника в разделе Supporting Functions в конце этого примера.

helperPlotRandomRecords(ECGData,14)

Время вейвлета, рассеиваясь

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

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

Можно визуализировать вейвлет, просачивается эти два набора фильтров со следующим.

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

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

scat_features_train = featureMatrix(sf,trainData');

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

Для получения матрицу, совместимую с классификатором SVM, изменитесь, рассеивание мультисигнала преобразовывают в матрицу, где каждый столбец соответствует рассеивающемуся пути, и каждая строка является рассеивающимся окном времени. В этом случае вы получаете 1 808 строк, потому что существует 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 416 16 49, потому что существует 49 форм волны ECG в наборе обучающих данных. После изменения для классификатора SVM матрица функции 784 416.

scat_features_test = featureMatrix(sf,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 с квадратичным ядром. Всего там 2 592 рассеивающихся последовательности в наборе данных в целом, 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)
fprintf('Accuracy is %2.2f percent.\n',100-loss);

confmatCV =

        1535           0           1
           1         479           0
           0           0         576

Accuracy is 99.92 percent.

Точность составляет 99,92%, который довольно хорош, но фактические результаты лучше полагают, что здесь каждый раз окно классифицируется отдельно. Существует 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);
MVconfmatCV = confusionmat(ClassVotes,categorical([trainLabels; testLabels]));
MVconfmatCV

True cross-validation accuracy is 100.00 percent.

MVconfmatCV =

    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);
testconfmat = confusionmat(TestVotes,categorical(testLabels));
testconfmat

The test accuracy is 97.96 percent. 

testconfmat =

    29     1     0     0
     0     8     0     0
     0     0    11     0
     0     0     0     0

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

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

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

Ссылки

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

Графики helperPlotRandomRecords четыре сигнала ECG, случайным образом выбранные из 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

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

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