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

Данные, используемые в этом примере, были получены из открытого хранилища данных 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.

Текстовый файл trafflevibrationdata.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) и «Uncracked» (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). Если необходимо пропустить сетевое обучение, можно перейти непосредственно в раздел «Классификация и анализ». Там можно загрузить обученную сеть GRU, а также предварительно обработанные наборы обучения и тестирования.

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

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

Загрузите обученные стробированные рекуррентные модули (GRU) и 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». Сеть не перекрывала, чтобы классифицировать временные ряды как «Uncracked», несмотря на дисбаланс.

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

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» является плохим. Количество ложных срабатываний данных «Crack» довольно велико. Это именно то, чего вы ожидаете с несбалансированным набором данных, когда модель не научилась хорошо.

Протестируйте 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 с необработанными последовательностями работает лучше, чем сеть GRU, но не так хорошо, как сверточная сеть, обученная последовательностям вейвлет. Кроме того, уменьшение по временной размерности данных с помощью преобразования рассеяния привело к созданию сети, которая примерно в 8,5 раз меньше сверточной сети, обученной на необработанных данных.

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

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

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

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

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