Этот пример показов, как классифицировать сигналы электрокардиограммы (ЭКГ) человека с помощью вейвлета рассеяния времени и классификатора машины опорных векторов (SVM). При вейвлет данные распространяются через ряд вейвлет, нелинейностей и усреднения, чтобы получить низкодисперсные представления временных рядов. Вейвлет рассеяние приводит к представлениям сигнала, нечувствительным к сдвигам в входном сигнале, не жертвуя различимостью классов. Для запуска этого примера необходимо иметь Wavelet Toolbox™ и Statistics and Machine Learning Toolbox™. Данные, используемые в этом примере, являются общедоступными от PhysioNet. Вы можете найти подход глубокого обучения к этой задаче классификации в этом примере Classify Time Series Using Wavelet Analysis and Deep Learning и подход машинного обучения в этом примере Signal Classification Using Wavelet-Based Features and машины опорных векторов.
Этот пример использует данные ЭКГ, полученные из трех групп или классов людей: людей с сердечной аритмией, людей с застойным сердечным отказом и людей с нормальными синусовыми ритмами. В примере используются 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
README.md
и ECGData.zip содержит
ECGData.mat
Modified_physionet_data.txt
License.txt.
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 только к обучающим данным (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 к преобразованиям рассеяния тестового набора задержки.
Анден, Дж., Маллат, С. 2014. Спектр глубокого рассеяния, Транзакции IEEE по обработке сигналов, 62, 16, стр. 4114-4128.
Baim DS, Colucci WS, Monrad ES, Smith HS, Wright RF, Lanoue A, Gauthier DF, Ransil BJ, Grossman W, Braunwald E. Выживаемость пациентов с тяжёлым застойным сердечным отказом, получавших J American College of Cardiology 1986 Mar; 7(3):661-670.
Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE. PhysioBank, PhysioToolkit и PhysioNet: компоненты нового исследовательского ресурса для сложных физиологических сигналов. Циркуляция. Том 101, № 23, 13 июня 2000 года, стр. e215-e220. http://circ.ahajournals.org/content/101/23/e215.full
Маллат, С., 2012. Групповое инвариантное рассеяние. Коммуникации в чистой и прикладной математике, 65, 10, стр. 1331-1398.
Moody GB, Mark RG. Влияние базы данных аритмии MIT-BIH. IEEE Eng in Med and Biol 20 (3): 45-50 (май-июнь 2001). (PMID: 11446209)
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