Данные - описание и обязательные атрибуты

Данные, использованные в этом примере, были получены из репозитория открытых данных Mendalay Data [2]. Данные распространяются под лицензией Creative Commons (CC) BY 4.0. Загрузите данные и модели, использованные в этом примере, во временный каталог, указанный в MATLAB ®tempdir команда. Если вы решили загрузить данные в папку, отличную от tempdirизмените имя каталога в последующих инструкциях. Распакуйте данные в папку, указанную как TransverseCrackData.

dataURL = 'https://ssd.mathworks.com/supportfiles/wavelet/crackDetection/transverse_crack.zip';
saveFolder = fullfile(tempdir,'TransverseCrackData'); 
zipFile = fullfile(tempdir,'transverse_crack.zip');
websave(zipFile,dataURL);
unzip(zipFile,saveFolder)

После распаковки данных, TransverseCrackData папка содержит подпапку с именем transverse_crack_latest. Все последующие команды должны выполняться в этой папке, иначе эту папку можно поместить в папку matlabpath.

Текстовый файл, vilevibrationdata.rights, включенный в zip-файл, содержит текст лицензии CC BY 4.0. Данные были переупакованы из исходного формата Excel в MAT-файлы.

Сбор данных подробно описан в [1]. Было использовано двенадцать четырехметровых участков асфальта, содержащих центральную поперечную трещину, и двенадцать равноудаленных участков. Данные получены с трех отдельных дорог. Поперечные трещины имели ширину от 2 - 13 мм с промежутками трещин от 7 - 35 мм. Секции приводились в движение с тремя различными скоростями: 30 км/ч, 40 км/ч и 50 км/ч. Измерения вертикального ускорения от шарнира подвески переднего пассажира производятся с частотой дискретизации 1,28 кГц. Скорости 30, 40 и 50 км/ч соответствуют измерениям датчиков каждые 6,5 мм при 30 км/ч, 8,68 мм при 40 км/ч и 10,85 мм при 50 км/ч. Подробный вейвлет-анализ этих данных, отличных от анализа в этом примере, см. в [1].

В соответствии с условиями лицензии CC BY 4.0, мы отмечаем, что информация о скорости исходных данных не сохраняется в данных, используемых в этом примере. Информация о скорости и дороге сохраняется в road1.mat, road2.mat, road3.mat файлы данных, включенные в папку данных для полноты.

Данные - загрузка и импорт

Загрузите данные акселерометра и соответствующие метки. Имеется 327 записей акселерометров.

load(fullfile(saveFolder,"transverse_crack_latest","allroadData.mat"));
load(fullfile(saveFolder,"transverse_crack_latest","allroadLabel.mat"));

Данные - предварительная обработка

Получение длины всех временных рядов. Отображение гистограммы количества временных рядов на длину.

tslen = cellfun(@length,allroadData);
uLen = unique(tslen);
Ng = histcounts(tslen);
Ng = Ng(Ng > 0);
bar(uLen,Ng,0.5)
grid on
AX = gca;
AX.YLim = [0 max(Ng)+15];
text(uLen(1),Ng(1)+10,num2str(Ng(1)))
text(uLen(2),Ng(2)+10,num2str(Ng(2)))
text(uLen(3),Ng(3)+10,num2str(Ng(3)))
xlabel('Length in Samples')
ylabel('Number of Series')
title('Time Series Length')

Существует три уникальных длины в наборе данных: 369, 461 и 615 выборок. Транспортное средство движется с тремя различными скоростями, но пройденное расстояние и частота дискретизации постоянны, что приводит к различным длинам записи данных. Определите, сколько записей находится в «Cracked» (CR) и «Без трещин» (UNCR) классы.

countlabels(allroadLabel)

ans=2×3 table
    Label    Count    Percent
    _____    _____    _______

    CR        109     33.333 
    UNCR      218     66.667 

Этот набор данных значительно несбалансирован. Существует в два раза больше временных рядов без трещины (UNCR) в виде серии, содержащей трещину (CR). Это означает, что классификатор, который предсказывает «Uncracked» на каждой записи, достигнет точности 67% без какого-либо обучения.

Временные ряды также имеют разную длину. Для использования сети свертки необходима общая форма ввода. В повторяющихся сетях в качестве входных данных можно использовать временные ряды с неравной длиной, но все временные ряды в мини-пакете дополняются или усекаются на основе вариантов обучения. Это требует осторожности при создании мини-партий как для обучения, так и для тестирования для обеспечения правильного распределения дополненных последовательностей. Кроме того, это требует, чтобы вы не тасовали данные во время обучения. С помощью этого небольшого набора данных желательно перетасовать обучающие данные для каждой эпохи. Соответственно, используется общая длина временного ряда.

Наиболее распространенная длина - 461 образец. Кроме того, трещина, если она имеется, расположена в центре записи. Соответственно, мы можем симметрично расширить серию с 369 выборок до длины 461, отражая начальную и конечную 46 выборок. В записях с 615 образцами удалите исходные 77 и последние 77 образцов.

Обучение - извлечение функций и обучение работе с сетью

Следующие разделы генерируют обучающие и тестовые наборы, создают последовательности вейвлет-рассеяния и обучают как стробированные повторяющиеся блоки (GRU), так и 1-D сверточные сети. Во-первых, расширить или усечь временной ряд в обоих наборах, чтобы получить общую длину 461 выборки.

allroadData = equalLenTS(allroadData);
all(cellfun(@numel,allroadData)== 461)

ans = logical
   1

Теперь каждый временной ряд как в треснувших, так и без трещины наборах данных имеет 461 выборку. Разделите данные в обучающем наборе, состоящем из 80% временных рядов в каждом классе, и удерживайте оставшиеся 20% каждого класса для тестирования. Убедитесь, что несбалансированные пропорции сохранены в каждом наборе.

splt8020 = splitlabels(allroadLabel,0.80);
countlabels(allroadLabel(splt8020{1}))

ans=2×3 table
    Label    Count    Percent
    _____    _____    _______

    CR         87     33.333 
    UNCR      174     66.667 

countlabels(allroadLabel(splt8020{2}))

ans=2×3 table
    Label    Count    Percent
    _____    _____    _______

    CR        22      33.333 
    UNCR      44      66.667 

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

TrainData = allroadData(splt8020{1});
TrainLabels = allroadLabel(splt8020{1});
TestData = allroadData(splt8020{2});
TestLabels = allroadLabel(splt8020{2});

Перетасуйте данные и этикетки один раз перед тренировкой.

idxS = randperm(length(TrainData));
TrainData = TrainData(idxS);
TrainLabels = TrainLabels(idxS);
idxS = randperm(length(TrainData));
TrainData = TrainData(idxS);
TrainLabels = TrainLabels(idxS);

Вычислите последовательности рассеяния для каждого обучающего ряда. Последовательности рассеяния хранятся в матрице ячеек, чтобы быть совместимыми с сетью GRU. Эти последовательности впоследствии преобразуются в 4-D входной массив для использования с 1-D сверточной сетью.

XTrainSCAT = cell(size(TrainData));
for kk = 1:numel(TrainData)
    XTrainSCAT{kk} = helperscat(TrainData{kk});
end
npaths = cellfun(@(x)size(x,1),XTrainSCAT);
inputSize = npaths(1);

Обучение - Сеть ГРУ

Создайте сетевые уровни GRU. Используйте два слоя GRU с 30 скрытыми блоками каждый, а также два уровня отсева. Входной размер - это количество путей рассеяния. Чтобы обучить эту сеть необработанному временному ряду, измените inputSize 1 и транспонирование каждого временного ряда в вектор строки (1 на 461). Если вы хотите пропустить обучение сети, вы можете перейти непосредственно в раздел Классификация и анализ. Там можно загрузить обученную сеть ГРУ, а также предварительно обработанные обучающие и тестовые наборы.

numHiddenUnits1 = 30;
numHiddenUnits2 = 30;
numClasses = 2;

GRUlayers = [ ...
    sequenceInputLayer(inputSize,'Name','InputLayer',...
        'Normalization','zerocenter')
    gruLayer(numHiddenUnits1,'Name','GRU1','OutputMode','sequence')
    dropoutLayer(0.35,'Name','Dropout1')
    gruLayer(numHiddenUnits2,'Name','GRU2','OutputMode','last')
    dropoutLayer(0.2,'Name','Dropout2')
    fullyConnectedLayer(numClasses,'Name','FullyConnected')
    softmaxLayer('Name','smax');
    classificationLayer('Name','ClassificationLayer')];

Обучение сети ГРУ. Используйте мини-партию размером 15 с 150 эпохами.

maxEpochs = 150;
miniBatchSize = 15;

options = trainingOptions('adam', ...
    'ExecutionEnvironment','gpu', ...
    'L2Regularization',1e-3,...
    'MaxEpochs',maxEpochs, ...
    'MiniBatchSize',miniBatchSize, ...
    'Shuffle','every-epoch',...
    'Verbose',0,...
    'Plots','none');
iterPerEpoch = floor(length(XTrainSCAT)/miniBatchSize);
[scatGRUnet,infoGRU] = trainNetwork(XTrainSCAT,TrainLabels,GRUlayers,options);

Постройте график сглаженной точности обучения и потерь на итерацию.

figure
subplot(2,1,1)
smoothedACC = filter(1/iterPerEpoch*ones(iterPerEpoch,1),1,...
    infoGRU.TrainingAccuracy);
smoothedLoss = filter(1/iterPerEpoch*ones(iterPerEpoch,1),1,...
    infoGRU.TrainingLoss);
plot(smoothedACC)
title(['Training Accuracy (Smoothed) '...
    num2str(iterPerEpoch) ' iterations per epoch'])
ylabel('Accuracy (%)')
ylim([0 100.1])
grid on
xlim([1 length(smoothedACC)])
subplot(2,1,2)
plot(smoothedLoss)
ylim([-0.01 1])
grid on
xlim([1 length(smoothedLoss)])
ylabel('Loss')
xlabel('Iteration')

Получить преобразования вейвлет-рассеяния удержанных тестовых данных для классификации.

XTestSCAT = cell(size(TestData));
for kk = 1:numel(TestData)
    XTestSCAT{kk} = helperscat(TestData{kk});
end

Обучение - 1-D сверточная сеть

Обучить 1-D сверточную сеть последовательностями вейвлет-рассеяния. Последовательности рассеяния - 28 на 58, где 58 - количество временных шагов, а 28 - количество путей рассеяния. Чтобы использовать это в 1-D сверточной сети, транспонируйте и изменяйте форму последовательностей рассеяния, которые должны быть 58-by-1-by-28-by-Nsig, где Nsig - количество обучающих примеров. Если вы хотите пропустить обучение сети, вы можете перейти непосредственно к разделу Классификация и анализ. Там можно загрузить обученную сверточную сеть, а также предварительно обработанные обучающие и тестовые наборы.

scatCONVTrain = zeros(58,1,28,length(XTrainSCAT),'single');
for kk = 1:length(XTrainSCAT)
    scatCONVTrain(:,:,:,kk) = reshape(XTrainSCAT{kk}',[58,1,28]);
end

Постройте и обучите 1-D сверточную сеть.

conv1dLayers = [
    imageInputLayer([58 1 28]);
    convolution2dLayer([3 1],64,'stride',1);
    batchNormalizationLayer;
    reluLayer;
    maxPooling2dLayer([2 1]);
    convolution2dLayer([3 1],32);
    reluLayer;
    maxPooling2dLayer([2 1]);
    fullyConnectedLayer(500);
    fullyConnectedLayer(2);
    softmaxLayer;
    classificationLayer;
    ];

convoptions = trainingOptions('sgdm', ...
    'InitialLearnRate',0.01, ...
    'LearnRateSchedule','piecewise', ...
    'LearnRateDropFactor',0.5, ...
    'LearnRateDropPeriod',5, ...
    'Plots','none',...
    'MaxEpochs',50,...
    'Verbose',0,...
    'Plots','none',...
    'MiniBatchSize',15);
[scatCONV1Dnet,infoCONV] = ...
    trainNetwork(scatCONVTrain,TrainLabels,conv1dLayers,convoptions);

Постройте график сглаженной точности обучения и потерь на итерацию.

iterPerEpoch = floor(size(scatCONVTrain,4)/15);
figure
subplot(2,1,1)
smoothedACC = filter(1/iterPerEpoch*ones(iterPerEpoch,1),1,...
    infoCONV.TrainingAccuracy);
smoothedLoss = filter(1/iterPerEpoch*ones(iterPerEpoch,1),1,...
    infoCONV.TrainingLoss);
plot(smoothedACC)
title(['Training Accuracy (Smoothed) '...
    num2str(iterPerEpoch) ' iterations per epoch'])
ylabel('Accuracy (%)')
ylim([0 100.1])
grid on
xlim([1 length(smoothedACC)])
subplot(2,1,2)
plot(smoothedLoss)
ylim([-0.01 1])
grid on
xlim([1 length(smoothedLoss)])
ylabel('Loss')
xlabel('Iteration')

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

scatCONVTest = zeros(58,1,28,length(XTestSCAT));
for kk = 1:length(XTestSCAT)
    scatCONVTest(:,:,:,kk) = reshape(XTestSCAT{kk}',58,1,28);
end

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

Загрузить обученный стробируемый рекуррентный блок (ГРУ) и 1-D сверточные сети вместе с тестовыми данными и последовательностями рассеяния. Все данные, функции и сети были созданы в разделе Обучение - Извлечение функций и сетевое обучение.

load(fullfile(saveFolder,"transverse_crack_latest","TestData.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","TestLabels.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","XTestSCAT.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","scatGRUnet"))
load(fullfile(saveFolder,"transverse_crack_latest","scatCONV1Dnet.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","scatCONVTest.mat"))

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

load(fullfile(saveFolder,"transverse_crack_latest","TrainData.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","TrainLabels.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","XTrainSCAT.mat"))
load(fullfile(saveFolder,"transverse_crack_latest","scatCONVTrain.mat"))

Проверьте количество временных рядов на класс в тестовом наборе. Следует отметить, что набор тестов значительно несбалансирован, как описано в разделе «Данные - предварительная обработка».

countlabels(TestLabels)

ans=2×3 table
    Label    Count    Percent
    _____    _____    _______

    CR        22      33.333 
    UNCR      44      66.667 

XTestSCAT содержит последовательности вейвлет-рассеяния, вычисленные для необработанного временного ряда в TestData.

Отображение производительности модели GRU на тестовых данных, не используемых при обучении модели.

miniBatchSize = 15;
ypredSCAT = classify(scatGRUnet,XTestSCAT, ...
    'MiniBatchSize',miniBatchSize);
figure
confusionchart(TestLabels,ypredSCAT,'RowSummary','row-normalized',...
    'ColumnSummary','column-normalized');
title({'GRU network -- Wavelet Scattering Sequences'; ...
    'Confusion chart with precision and recall'});

Несмотря на большой дисбаланс в классах и малом наборе данных, точность и значения повторного вызова указывают на то, что сеть хорошо работает с тестовыми данными. В частности, значение точности для данных «Cracked» составляет около 98%, а коэффициент отзыва составляет приблизительно 96%. Эти значения особенно хороши, учитывая, что полные 67% записей в обучающем наборе были «Без трещины». Сеть не перекрыла себя, чтобы классифицировать временные ряды как «Uncracked», несмотря на дисбаланс.

Если установить inputSize = 1 и перенести временной ряд в данные обучения, можно переучить сеть ГРУ на необработанные данные временного ряда. Это было сделано на тех же данных в обучающем наборе. Можно загрузить эту сеть и проверить производительность на тестовом наборе.

load(fullfile(saveFolder,"transverse_crack_latest","tsGRUnet.mat"));
rawTest = cellfun(@transpose,TestData,'UniformOutput',false);
miniBatchSize = 15;
YPredraw = classify(tsGRUnet,rawTest, ...
    'MiniBatchSize',miniBatchSize);
confusionchart(TestLabels,YPredraw,'RowSummary','row-normalized',...
    'ColumnSummary','column-normalized');
title({'GRU network -- Raw Time Series'; ...
    'Confusion chart with precision and recall'});

Для этой сети производительность плохая. В частности, отзыв для данных «Cracked» является плохим. Количество ложных негативов для данных «Крэк» достаточно велико. Это именно то, что вы ожидаете с несбалансированным набором данных, когда модель плохо научилась.

Проверьте 1-D сверточную сеть, обученную последовательностям вейвлет-рассеяния.

miniBatchSize = 15;
YPredSCAT = classify(scatCONV1Dnet,scatCONVTest,...
    'MiniBatchSize',miniBatchSize);
figure
confusionchart(TestLabels,YPredSCAT,'RowSummary','row-normalized',...
    'ColumnSummary','column-normalized');
title({'1-D Convolutional Network-- Wavelet Scattering Sequences'; ...
    'Confusion chart with precision and recall'});

Производительность сверточной сети с последовательностями рассеяния превосходна и превосходит производительность сети ГРУ. Точность и отзыв о меньшинстве демонстрируют надежное обучение.

Обучать 1-D сверточную сеть на сыром наборе последовательностей inputSize = [461 1 1] в imageInputLayer. Вы можете загрузить и протестировать эту сеть.

load(fullfile(saveFolder,"transverse_crack_latest","tsCONV1Dnet.mat"));
XTestData = cell2mat(TestData');
XTestData = reshape(XTestData,[461 1 1 66]);
miniBatchSize = 15;
YPredRAW = classify(tsCONV1Dnet,XTestData,...
    'MiniBatchSize',miniBatchSize);
confusionchart(TestLabels,YPredRAW,'RowSummary','row-normalized',...
    'ColumnSummary','column-normalized');
title({'1-D Convolutional Network-- Raw Sequences'; ...
    'Confusion chart with precision and recall'});

Сеть свертки 1-D с необработанными последовательностями работает лучше, чем сеть ГРУ, но не так хорошо, как сеть свертки, обученная последовательностям вейвлет-рассеяния. Кроме того, уменьшение во временном измерении данных с преобразованием рассеяния привело к сети, которая приблизительно в 8,5 раз меньше, чем сверточная сеть, обученная на необработанных данных.

Вывод и анализ вейвлета

В этом разделе показано, как классифицировать один временной ряд с помощью вейвлет-анализа с предварительно подготовленной моделью. В качестве модели используется сверточная сеть 1-D, обученная последовательностям вейвлет-рассеяния. Загрузите обученную сеть и некоторые тестовые данные, если они еще не были загружены в предыдущем разделе.

load(fullfile(saveFolder,"transverse_crack_latest","scatCONV1Dnet.mat"));
load(fullfile(saveFolder,"transverse_crack_latest","TestData.mat"));

Создайте сеть вейвлет-рассеяния для преобразования данных. Выберите временной ряд из данных теста и классифицируйте данные. Если модель классифицирует временной ряд как «Треснувший», исследуйте ряд на предмет положения трещины в форме сигнала.

sf = waveletScattering('SignalLength',461,...
    'OversamplingFactor',1,'qualityfactors',[8 1],...
    'InvarianceScale',0.05,'Boundary','reflect','SamplingFrequency',1280);
idx = 1;
data = TestData{idx};
[smat,x] = featureVectors(data,sf);
PredictedClass = classify(scatCONV1Dnet,smat);
if isequal(PredictedClass,'CR') 
    fprintf('Crack detected. Computing wavelet transform modulus maxima.\n')
    wtmm(data,'Scaling','local');
end

Crack detected. Computing wavelet transform modulus maxima.

Метод максимумов модуля вейвлет-преобразования (WTMM) показывает линию максимумов, сходящуюся к самой точной шкале в выборке 205. Линии максимумов, сходящиеся к мелким масштабам, являются хорошей оценкой того, где сингулярности находятся во временном ряду. Это делает образец 205 хорошей оценкой местоположения трещины.

plot(x,data)
axis tight
hold on
plot([x(205) x(205)],[min(data) max(data)],'k')
grid on
title(['Crack located at ' num2str(x(205)) 'meters'])
xlabel('Distance (m)')
ylabel('Amplitude')

Вы можете повысить свою уверенность в этом местоположении, используя методы анализа множественных решений (MRA) и выявляя изменения в наклоне в серии длинных вейвлет MRA. Введение в методы MRA см. в разделе Практическое введение в анализ множественных решений (Wavelet Toolbox). В [1] разница в энергии между сериями «Cracked» и «Uncracked» произошла в низкочастотных диапазонах, в частности в интервале [10,20] Гц. Соответственно, следующая MRA сфокусирована на компонентах сигнала в частотных диапазонах от [10,80] Гц. В этих диапазонах определите линейные изменения данных. Постройте график точек изменения вместе с компонентами MRA.

[mra,chngpts] = helperMRA(data,x);

Анализ точек изменений на основе MRA помог подтвердить анализ WTMM при идентификации области около 1,8 метра в качестве вероятного местоположения трещины.

Резюме

Этот пример показал, как использовать последовательности вейвлет-рассеяния с повторяющимися и сверточными сетями для классификации временных рядов. Пример дополнительно продемонстрировал, как вейвлет-методы могут помочь локализовать признаки в том же пространственном (временном) масштабе, что и исходные данные.

Ссылки

[1] Ян, Цюнь и Шиши Чжоу. «Определение поперечного растрескивания асфальтового покрытия на основе анализа сигнала вибрации транспортного средства», Проект дорожных материалов и дорожного покрытия, 2020, 1-19. https://doi.org/10.1080/14680629.2020.1714699.

[2] Чжоу, шиши. «Данные о вибрации транспортного средства» https://data.mendeley.com/datasets/3dvpjy4m22/1. Данные используются в CC BY 4.0. Данные переупаковываются из исходного формата данных Excel в MAT-файлы. Метка скорости удалена и сохраняется только метка «crack» или «nocrack».

Приложение

Вспомогательные функции, используемые в этом примере.

function smat = helperscat(datain)
% This helper function is only in support of Wavelet Toolbox examples.
% It may change or be removed in a future release.
datain = single(datain);

sn = waveletScattering('SignalLength',length(datain),...
    'OversamplingFactor',1,'qualityfactors',[8 1],...
    'InvarianceScale',0.05,'Boundary','reflect','SamplingFrequency',1280);
smat = sn.featureMatrix(datain);

end

%-----------------------------------------------------------------------
function dataUL = equalLenTS(data)
% This function in only in support of Wavelet Toolbox examples.
% It may change or be removed in a future release.
N = length(data);
dataUL = cell(N,1);
for kk = 1:N
    L = length(data{kk});
    switch L
        case 461
            dataUL{kk} = data{kk};
        case 369
            Ndiff = 461-369;
            pad = Ndiff/2;
            dataUL{kk} = [flip(data{kk}(1:pad)); data{kk} ; ...
                flip(data{kk}(L-pad+1:L))];     
        otherwise
            Ndiff = L-461;
            zrs = Ndiff/2;
            dataUL{kk} = data{kk}(zrs:end-zrs-1);
    end
end       

end

%--------------------------------------------------------------------------
function [fmat,x] = featureVectors(data,sf)
% This function is only in support of Wavelet Toolbox examples.
% It may change or be removed in a future release.
data = single(data);
N = length(data);
dt = 1/1280;
if N < 461
    Ndiff = 461-N;
    pad = Ndiff/2;
    dataUL = [flip(data(1:pad)); data ; ...
                flip(data(N-pad+1:N))];   
     rate = 5e4/3600;
     dx = rate*dt;
     x = 0:dx:(N*dx)-dx;     
    
elseif N > 461
    Ndiff = N-461;
    zrs = Ndiff/2;
    dataUL = data(zrs:end-zrs-1);
    rate = 3e4/3600;
    dx = rate*dt;
    x = 0:dx:(N*dx)-dx;
else
    dataUL = data;
    rate = 4e4/3600;
    dx = rate*dt;
    x = 0:dx:(N*dx)-dx;
end
fmat = sf.featureMatrix(dataUL);
fmat = reshape(fmat',[58 1 28]);
end

%------------------------------------------------------------------------------
function [mra,chngpts] = helperMRA(data,x)
% This function is only in support of Wavelet Toolbox examples.
% It may change or be removed in a future release.
mra = modwtmra(modwt(data,'sym3'),'sym3');
mraLev = mra(4:6,:);
Ns = size(mraLev,1);
thresh = [2, 4, 8];
chngpts = false(size(mraLev));
% Determine changepoints. We want different thresholds for different
% resolution levels.
for ii = 1:Ns
    chngpts(ii,:) = ischange(mraLev(ii,:),"linear",2,"Threshold",thresh(ii));
end

for kk = 1:Ns
    idx = double(chngpts(kk,:));
    idx(idx == 0) = NaN;    
    subplot(Ns,1,kk)
    plot(x,mraLev(kk,:))
    if kk == 1
        title('MRA Components')
    end
    yyaxis right
    hs = stem(x,idx);
    hs.ShowBaseLine = 'off';
    hs.Marker = '^';
    hs.MarkerFaceColor = [1 0 0];
end
grid on
axis tight
xlabel('Distance (m)')
end