Документация

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

Этот пример использует данные о ECG, полученные из трех групп или классов, людей: люди с сердечной аритмией, люди с застойной сердечной недостаточностью и люди с нормальными ритмами пазухи. Пример использует 162 записи ECG от трех Физиосетевых баз данных: База данных Аритмии 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. Data 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);

Существует 113 записей в наборе trainData и 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