Классификация сигнала, использующая основанные на вейвлете функции и машины опорных векторов

В этом примере показано, как классифицировать человеческую электрокардиограмму (ECG) сигналы с помощью основанного на вейвлете извлечения признаков и классификатора машины опорных векторов (SVM). Проблема классификации сигнала упрощена путем преобразования необработанных сигналов ECG в намного меньший набор функций, которые служат в агрегате, чтобы дифференцировать различные классы. У вас должны быть Wavelet Toolbox™, Signal Processing Toolbox™ и Statistics and Machine Learning Toolbox™, чтобы запустить этот пример. Данные, используемые в этом примере, общедоступны от PhysioNet.

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

Этот пример использует данные о ECG, полученные из трех групп или классов, людей: люди с сердечной аритмией, люди с застойной сердечной недостаточностью и люди с нормальными ритмами пазухи. Пример использует 162 записи ECG от трех Физиосетевых баз данных: База данных Аритмии MIT-BIH [3] [7], MIT-BIH Нормальная База данных Ритма Пазухи [3] и База данных Застойной сердечной недостаточности BIDMC [1] [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

  • README.md

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

  • ECGData.mat

  • Modified_physionet_data.txt

  • License.txt.

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 герц. 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 блоках каждого сигнала приблизительно одна минута в длительности (8 192 выборки):

  • Авторегрессивная модель (AR) коэффициенты порядка 4 [8].

  • Значения шенноновской энтропии (SE) для максимального перекрытия дискретный пакет вейвлета преобразовывают (MODPWT) на уровне 4 [5].

  • Мультифрактальные оценки лидера вейвлета второго cumulant масштабирующихся экспонент и области значений экспонент Держателя или спектра [4] сингулярности.

Кроме того, многошкальные оценки отклонения вейвлета извлечены для каждого сигнала по целой длине данных [6]. Объективная оценка отклонения вейвлета используется. Это требует, чтобы только уровни по крайней мере с одним коэффициентом вейвлета, незатронутым граничными условиями, использовались в оценках отклонения. Для длины сигнала 2^16 (65,536) и 'db2' вейвлет это приводит к 14 уровням.

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

Коэффициенты AR для каждого окна оцениваются с помощью метода Города, arburg. В [8], авторы использовали методы выбора порядка модели решить, что модель AR (4) обеспечила лучшее пригодное для форм волны ECG в подобной проблеме классификации. В [5], информация теоретическая мера, шенноновская энтропия, вычислялась на терминальных узлах пакетного дерева вейвлета и использовалась со случайным лесным классификатором. Здесь мы используем неподкошенный пакет вейвлета, преобразовывают, modwpt, вниз к уровню 4.

Определение шенноновской энтропии для неподкошенного пакета вейвлета преобразовывает следующий [5], дают: где количество соответствующих коэффициентов в j-th узле и является нормированными квадратами пакетных коэффициентов вейвлета в j-th терминальном узле.

Две фрактальных меры, оцененные методами вейвлета, используются в качестве функций. После [4], мы используем ширину спектра сингулярности, полученного из dwtleader как мера мультифрактальной природы сигнала ECG. Мы также используем второй cumulant масштабирующихся экспонент. Масштабирующиеся экспоненты являются основанными на шкале экспонентами, описывающими поведение закона степени в сигнале в различных разрешениях. Второй cumulant широко представляет отъезд масштабирующихся экспонент от линейности.

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

Всего существует 190 функций: 32 функции AR (4 коэффициента на блок), 128 значений шенноновской энтропии (16 значений на блок), 16 фрактальных оценок (2 на блок) и 14 оценок отклонения вейвлета.

helperExtractFeatures функция вычисляет эти функции и конкатенирует их в характеристический вектор для каждого сигнала. Можно найти исходный код для этой функции помощника в разделе Supporting Functions в конце этого примера.

timeWindow = 8192;
ARorder = 4;
MODWPTlevel = 4;
[trainFeatures,testFeatures,featureindices] = ...
    helperExtractFeatures(trainData,testData,timeWindow,ARorder,MODWPTlevel);

trainFeatures и testFeatures 113 190 и 49 190 матрицы, соответственно. Каждая строка этих матриц является характеристическим вектором для соответствующих данных о ECG в trainData и testData, соответственно. В создании характеристических векторов данные уменьшаются с 65 536 выборок до 190 векторов элемента. Это - значительное сокращение данных, но целью не является только сокращение данных. Цель состоит в том, чтобы уменьшать данные до намного меньшего набора функций, который получает различие между классами так, чтобы классификатор мог точно разделить сигналы. Индексы для функций, которые составляют оба trainFeatures и testFeatures содержатся в массиве структур, featureindices. Можно использовать эти индексы, чтобы исследовать функции группой. Как пример, исследуйте область значений экспонент Держателя в спектрах сингулярности впервые окно. Отобразите данные на графике для целого набора данных.

allFeatures = [trainFeatures;testFeatures];
allLabels = [trainLabels;testLabels];
figure
boxplot(allFeatures(:,featureindices.HRfeatures(1)),allLabels,'notch','on')
ylabel('Holder Exponent Range')
title('Range of Singularity Spectrum by Group (First Time Window)')
grid on

Можно выполнить односторонний дисперсионный анализ на этой функции и подтвердить то, что появляется в коробчатой диаграмме, а именно, что у ARR и групп NSR есть значительно большая область значений, чем группа швейцарского франка.

[p,anovatab,st] = anova1(allFeatures(:,featureindices.HRfeatures(1)),...
    allLabels);
c = multcompare(st,'display','off')
c =

    1.0000    2.0000    0.0176    0.1144    0.2112    0.0155
    1.0000    3.0000   -0.1591   -0.0687    0.0218    0.1764
    2.0000    3.0000   -0.2975   -0.1831   -0.0687    0.0005

Как дополнительный пример, считайте различие в отклонении во второй самой низкой частоте (вторая по величине шкала) поддиапазоном вейвлета для этих трех групп.

boxplot(allFeatures(:,featureindices.WVARfeatures(end-1)),allLabels,'notch','on')
ylabel('Wavelet Variance')
title('Wavelet Variance by Group')
grid on

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

Классификация сигнала

Теперь, когда данные уменьшались до характеристического вектора для каждого сигнала, следующий шаг должен использовать эти характеристические векторы в классификации сигналов ECG. Можно использовать приложение Classification Learner, чтобы быстро оценить большое количество классификаторов. В этом примере мультикласс используется SVM с квадратичным ядром. Выполняются два исследования. Сначала мы используем набор данных в целом (обучение и тестирующие наборы) и оцениваем misclassification уровень и матрицу беспорядка использование 5-кратной перекрестной проверки.

features = [trainFeatures; testFeatures];
rng(1)
template = templateSVM(...
    'KernelFunction','polynomial',...
    'PolynomialOrder',2,...
    'KernelScale','auto',...
    'BoxConstraint',1,...
    'Standardize',true);
model = fitcecoc(...
    features,...
    [trainLabels;testLabels],...
    'Learners',template,...
    'Coding','onevsone',...
    'ClassNames',{'ARR','CHF','NSR'});
kfoldmodel = crossval(model,'KFold',5);
classLabels = kfoldPredict(kfoldmodel);
loss = kfoldLoss(kfoldmodel)*100
[confmatCV,grouporder] = confusionmat([trainLabels;testLabels],classLabels);
loss =

    8.0247

5-кратная ошибка классификации составляет 8,02% (правильные 91,98%). Матрица беспорядка, confmatCV, показывает, какие записи были неправильно классифицированы. grouporder дает упорядоченное расположение групп. Две из группы ARR были неправильно классифицированы как швейцарский франк, восемь из группы швейцарского франка были неправильно классифицированы как ARR и один, как NSR, и два от группы NSR были неправильно классифицированы как ARR.

Точность, вспомните, и счет F1

В задаче классификации точность для класса является количеством правильных положительных результатов, разделенных на количество положительных результатов. Другими словами, всех записей, что классификатор присваивает данную метку, какая пропорция на самом деле принадлежит классу. Отзыв задан как количество правильных меток, разделенных на количество меток для данного класса. А именно, всех записей, принадлежащих классу, какая пропорция сделала нашу марку классификатора как тот класс. В оценке точности ваша система машинного обучения вы идеально хотите преуспеть и на точности и на отзыве. Например, предположите, что у нас был классификатор, который пометил каждую запись как ARR. Затем наш отзыв для класса ARR был бы 1 (100%). Все записи, принадлежащие классу ARR, были бы помечены ARR. Однако точность была бы низкой. Поскольку наш классификатор пометил все записи как ARR, будет 66 ложных положительных сторон в этом случае для точности 96/162, или 0.5926. Счет F1 является средним гармоническим точности и отзыва и поэтому обеспечивает одну метрику, которая обобщает производительность классификатора и в терминах отзыва и в терминах точности. Следующая функция помощника вычисляет точность, вспомните, и музыка F1 к этим трем классам. Вы видите как helperPrecisionRecall вычисляет точность, вспомните, и счет F1 на основе матрицы беспорядка путем исследования кода в разделе Supporting Functions.

CVTable = helperPrecisionRecall(confmatCV);

Можно отобразить таблицу, возвращенную helperPrecisionRecall со следующей командой.

disp(CVTable)
           Precision    Recall    F1_Score
           _________    ______    ________

    ARR     90.385      97.917         94 
    CHF     91.304          70     79.245 
    NSR     97.143      94.444     95.775 

И точность и отзыв хороши для ARR и классов NSR, в то время как отзыв значительно ниже для класса швейцарского франка.

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

model = fitcecoc(...
     trainFeatures,...
     trainLabels,...
     'Learners',template,...
     'Coding','onevsone',...
     'ClassNames',{'ARR','CHF','NSR'});
predLabels = predict(model,testFeatures);

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

correctPredictions = strcmp(predLabels,testLabels);
testAccuracy = sum(correctPredictions)/length(testLabels)*100
[confmatTest,grouporder] = confusionmat(testLabels,predLabels);
testAccuracy =

   97.9592

Точность классификации на тестовом наборе данных составляет приблизительно 98%, и матрица беспорядка показывает, что запись за один швейцарский франк была неправильно классифицирована как NSR.

Подобно тому, что было сделано в анализе перекрестной проверки, получите точность, вспомните, и музыка F1 к набору тестов.

testTable = helperPrecisionRecall(confmatTest);
disp(testTable)
           Precision    Recall    F1_Score
           _________    ______    ________

    ARR        100         100        100 
    CHF        100      88.889     94.118 
    NSR     91.667         100     95.652 

Классификация на необработанных данных и кластеризации

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

rawData = [trainData;testData];
Labels = [trainLabels;testLabels];
rng(1)
template = templateSVM(...
    'KernelFunction','polynomial', ...
    'PolynomialOrder',2, ...
    'KernelScale','auto', ...
    'BoxConstraint',1, ...
    'Standardize',true);
model = fitcecoc(...
    rawData,...
    [trainLabels;testLabels],...
    'Learners',template,...
    'Coding','onevsone',...
    'ClassNames',{'ARR','CHF','NSR'});
kfoldmodel = crossval(model,'KFold',5);
classLabels = kfoldPredict(kfoldmodel);
loss = kfoldLoss(kfoldmodel)*100
[confmatCVraw,grouporder] = confusionmat([trainLabels;testLabels],classLabels);
rawTable = helperPrecisionRecall(confmatCVraw);
disp(rawTable)
loss =

   33.3333

           Precision    Recall    F1_Score
           _________    ______    ________

    ARR        64          100     78.049 
    CHF       100       13.333     23.529 
    NSR       100       22.222     36.364 

misclassification уровень для необработанных данных временных рядов составляет 33,3%. При повторении точности вспомните, и анализ счета F1 показывает очень плохую музыку F1 и к швейцарскому франку (23.52) и к группам NSR (36.36). Получите коэффициенты дискретного преобразования Фурье (DFT) величины для каждого сигнала выполнить анализ в частотном диапазоне. Поскольку данные с действительным знаком, мы можем достигнуть некоторого снижения объема данных с помощью ДПФ путем использования того, что величины Фурье являются четной функцией.

rawDataDFT = abs(fft(rawData,[],2));
rawDataDFT = rawDataDFT(:,1:2^16/2+1);
rng(1)
template = templateSVM(...
    'KernelFunction','polynomial',...
    'PolynomialOrder',2,...
    'KernelScale','auto',...
    'BoxConstraint',1,...
    'Standardize',true);
model = fitcecoc(...
    rawDataDFT,...
    [trainLabels;testLabels],...
    'Learners',template,...
    'Coding','onevsone',...
    'ClassNames',{'ARR','CHF','NSR'});
kfoldmodel = crossval(model,'KFold',5);
classLabels = kfoldPredict(kfoldmodel);
loss = kfoldLoss(kfoldmodel)*100
[confmatCVDFT,grouporder] = confusionmat([trainLabels;testLabels],classLabels);
dftTable = helperPrecisionRecall(confmatCVDFT);
disp(dftTable)
loss =

   19.1358

           Precision    Recall    F1_Score
           _________    ______    ________

    ARR     76.423      97.917     85.845 
    CHF        100      26.667     42.105 
    NSR     93.548      80.556     86.567 

Используя ДПФ величины уменьшают misclassification уровень до 19,13%, но это - все еще более двух раз коэффициент ошибок, полученный с нашими 190 функциями. Эти исследования демонстрируют, что классификатор извлек выгоду из тщательного выбора функций.

Чтобы ответить на вопрос относительно роли классификатора, попытайтесь кластеризировать данные с помощью только характеристические векторы. Используйте k-средних значений, кластеризирующихся наряду со статистической величиной разрыва, чтобы определить и оптимальное количество кластеров и кластерное присвоение. Допускайте возможность 1 - 6 кластеров для данных.

rng default
eva = evalclusters(features,'kmeans','gap','KList',[1:6]);
eva
eva = 

  GapEvaluation with properties:

    NumObservations: 162
         InspectedK: [1 2 3 4 5 6]
    CriterionValues: [1.2777 1.3539 1.3644 1.3570 1.3591 1.3752]
           OptimalK: 3

Статистическая величина разрыва указывает, что оптимальное количество кластеров равняется трем. Однако, если вы смотрите на количество записей в каждом из этих трех кластеров, вы видите, что k-средних значений, кластеризирующиеся на основе характеристических векторов, сделали плохое задание разделения трех диагностических категорий.

countcats(categorical(eva.OptimalY))
ans =

    61
    74
    27

Вспомните, что существует 96 человек в классе ARR, 30 в классе швейцарского франка, и 36 в классе NSR.

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

Этот пример использовал обработку сигналов, чтобы извлечь функции вейвлета из сигналов ECG и использовал те функции, чтобы классифицировать сигналы ECG в три класса. Мало того, что извлечение признаков приводило к существенному количеству снижения объема данных, оно также получило различия между ARR, швейцарским франком и классами NSR, как продемонстрировано результатами перекрестной проверки и производительностью классификатора SVM на наборе тестов. Пример далее продемонстрировал, что применение классификатора SVM к необработанным данным привело к низкой производительности также, как и кластеризация характеристических векторов, не используя классификатор. Ни классификатор, ни одни только функции не были достаточны, чтобы разделить классы. Однако, когда извлечение признаков использовалось в качестве шага снижения объема данных до использования классификатора, эти три класса были хорошо разделены.

Ссылки

  1. Baim DS, WS Colucci, ES Monrad, Смит ХС, Wright RF, Lanoue A, Готье ДФ, Ransil BJ, Гроссман В, Браунвальд E. Выживание пациентов с тяжелой застойной сердечной недостаточностью отнеслось с устным milrinone. J американский Колледж Кардиологии 1 986 мартов; 7 (3):661-670.

  2. Engin, M., 2004. ECG разбил классификацию с помощью нейронечеткой сети. Буквы Распознавания образов, 25 (15), pp.1715-1722.

  3. Голдбергер ЭЛ, LAN Amaral, Стекло L, Гаусдорф ДЖМ, Иванов PCh, Марк РГ, Mietus JE, Капризный Гбайт, Пенг C-K, Стэнли ХЭ. PhysioBank, PhysioToolkit и PhysioNet: Компоненты Нового Ресурса Исследования для Комплексных Физиологических Сигналов. Циркуляция. Издание 101, № 23, 13 июня 2000, стр e215-e220. http://circ.ahajournals.org/content/101/23/e215.full

  4. Leonarduzzi, R.F., Schlotthauer, G. и Торрес. M.E. 2010. Вейвлет основанный на лидере мультифрактальный анализ изменчивости сердечного ритма во время миокардиальной ишемии. Разработка в Обществе Медицины и Биологии (EMBC), 2 010 Ежегодных Международных конференциях IEEE.

  5. Литий, T. и Чжоу, M., 2016. Классификация ECG с помощью пакета вейвлета энтропийные и случайные леса. Энтропия, 18 (8), p.285.

  6. Махарадж, Э.Э. и Алонсо, Утра 2014. Дискриминантный анализ многомерных временных рядов: Приложение к диагнозу на основе сигналов ECG. Вычислительная Статистика и Анализ данных, 70, стр 67-87.

  7. Капризный Гбайт, Марк РГ. Удар Базы данных Аритмии MIT-BIH. Инженер IEEE в Медиане и Biol 20 (3):45-50 (мочь-июнь 2001). (PMID: 11446209)

  8. Чжао, Q. и Чжан, L., 2005. Извлечение признаков ECG и классификация с помощью вейвлета преобразовывают и машины опорных векторов. Международная конференция IEEE по вопросам Нейронных сетей и Мозга, 2, стр 1089-1092.

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

Графики 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

Извлечения helperExtractFeatures функции вейвлета и коэффициенты AR для блоков данных заданного размера. Функции конкатенированы в характеристические векторы.

function [trainFeatures, testFeatures,featureindices] = helperExtractFeatures(trainData,testData,T,AR_order,level)
% This function is only in support of XpwWaveletMLExample. It may change or
% be removed in a future release.
trainFeatures = [];
testFeatures = [];

for idx =1:size(trainData,1)
    x = trainData(idx,:);
    x = detrend(x,0);
    arcoefs = blockAR(x,AR_order,T);
    se = shannonEntropy(x,T,level);
    [cp,rh] = leaders(x,T);
    wvar = modwtvar(modwt(x,'db2'),'db2');
    trainFeatures = [trainFeatures; arcoefs se cp rh wvar']; %#ok<AGROW>

end

for idx =1:size(testData,1)
    x1 = testData(idx,:);
    x1 = detrend(x1,0);
    arcoefs = blockAR(x1,AR_order,T);
    se = shannonEntropy(x1,T,level);
    [cp,rh] = leaders(x1,T);
    wvar = modwtvar(modwt(x1,'db2'),'db2');
    testFeatures = [testFeatures;arcoefs se cp rh wvar']; %#ok<AGROW>

end

featureindices = struct();
% 4*8
featureindices.ARfeatures = 1:32;
startidx = 33;
endidx = 33+(16*8)-1;
featureindices.SEfeatures = startidx:endidx;
startidx = endidx+1;
endidx = startidx+7;
featureindices.CP2features = startidx:endidx;
startidx = endidx+1;
endidx = startidx+7;
featureindices.HRfeatures = startidx:endidx;
startidx = endidx+1;
endidx = startidx+13;
featureindices.WVARfeatures = startidx:endidx;
end


function se = shannonEntropy(x,numbuffer,level)
numwindows = numel(x)/numbuffer;
y = buffer(x,numbuffer);
se = zeros(2^level,size(y,2));
for kk = 1:size(y,2)
    wpt = modwpt(y(:,kk),level);
    % Sum across time
    E = sum(wpt.^2,2);
    Pij = wpt.^2./E;
    % The following is eps(1)
    se(:,kk) = -sum(Pij.*log(Pij+eps),2);
end
se = reshape(se,2^level*numwindows,1);
se = se';
end


function arcfs = blockAR(x,order,numbuffer)
numwindows = numel(x)/numbuffer;
y = buffer(x,numbuffer);
arcfs = zeros(order,size(y,2));
for kk = 1:size(y,2)
    artmp =  arburg(y(:,kk),order);
    arcfs(:,kk) = artmp(2:end);
end
arcfs = reshape(arcfs,order*numwindows,1);
arcfs = arcfs';
end


function [cp,rh] = leaders(x,numbuffer)
y = buffer(x,numbuffer);
cp = zeros(1,size(y,2));
rh = zeros(1,size(y,2));
for kk = 1:size(y,2)
    [~,h,cptmp] = dwtleader(y(:,kk));
    cp(kk) = cptmp(2);
    rh(kk) = range(h);
end
end

helperPrecisionRecall возвращает точность, вспомните, и баллы F1 на основе матрицы беспорядка. Выводит результаты как таблицу MATLAB.

function PRTable = helperPrecisionRecall(confmat)
% This function is only in support of XpwWaveletMLExample. It may change or
% be removed in a future release.
precisionARR = confmat(1,1)/sum(confmat(:,1))*100;
precisionCHF = confmat(2,2)/sum(confmat(:,2))*100 ;
precisionNSR = confmat(3,3)/sum(confmat(:,3))*100 ;
recallARR = confmat(1,1)/sum(confmat(1,:))*100;
recallCHF = confmat(2,2)/sum(confmat(2,:))*100;
recallNSR = confmat(3,3)/sum(confmat(3,:))*100;
F1ARR = 2*precisionARR*recallARR/(precisionARR+recallARR);
F1CHF = 2*precisionCHF*recallCHF/(precisionCHF+recallCHF);
F1NSR = 2*precisionNSR*recallNSR/(precisionNSR+recallNSR);
% Construct a MATLAB Table to display the results.
PRTable = array2table([precisionARR recallARR F1ARR;...
    precisionCHF recallCHF F1CHF; precisionNSR recallNSR...
    F1NSR],'VariableNames',{'Precision','Recall','F1_Score'},'RowNames',...
    {'ARR','CHF','NSR'});

end