Данные — описание и необходимые приписывания

Данные, используемые в этом примере, были получены из Данных Mendeley открытый репозиторий данных [2]. Данные распределяются в соответствии с лицензией BY 4.0 Creative Commons (CC). Загрузите данные и модели, используемые в этом примере в вашей временной директории, заданной 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.

Текстовый файл, vehiclevibrationdata.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 выборок. Транспортное средство перемещается на трех различных скоростях, но пересеченное расстояние и частота дискретизации постоянно получившийся в различных длинах записи данных. Определите, сколько записей находится во "Взломанном" (CR) и "Невзломанный" (UNCRКлассы.

countlabels(allroadLabel)

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

    CR        109     33.333 
    UNCR      218     66.667 

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

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

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

Обучение — извлечение признаков и сетевое обучение

Следующие разделы генерируют наборы обучающих данных и наборы тестов, создают последовательности рассеивания вейвлета и обучают и закрытый текущий модуль (GRU) и 1D сверточные сети. Во-первых, расширьте или обрежьте временные ряды в обоих наборах, чтобы получить общую продолжительность 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);

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

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);

Обучение — сеть ГРУ

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

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

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

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

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

Создайте и обучите 1D сверточную сеть.

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) и 1D сверточные сети наряду с тестовыми данными и рассеивающимися последовательностями. Все данные, функции и сети были созданы в разделе Training — Feature Extraction и Network Training.

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"))

Исследуйте количество временных рядов в классе в наборе тестов. Обратите внимание, что набор тестов является значительно неустойчивым, как обсуждено в разделе Data — Preprocessing.

countlabels(TestLabels)

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

    CR        22      33.333 
    UNCR      44      66.667 

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

Покажите производительность модели ГРУ на тестовых данных, не используемых в обучении модели.

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'});

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

Если вы устанавливаете 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'});

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

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

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'});

Эффективность сверточной сети с рассеивающимися последовательностями превосходна и превышает эффективность сети ГРУ. Точность и отзыв на классе меньшинства демонстрируют устойчивое изучение.

Чтобы обучить 1D сверточную сеть на необработанных последовательностях устанавливает 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'});

1D сеть свертки с необработанными последовательностями выполняет лучше, чем сеть ГРУ, но не, а также сверточная сеть, обученная с последовательностями рассеивания вейвлета. Далее, сокращение вдоль измерения времени данных с рассеивающимся преобразованием привело к сети, которая приблизительно в 8.5 раз меньше, чем сверточная сеть, обученная на необработанных данных.

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

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

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. В [1] различие в энергии между серией "Cracked" и "Uncracked" произошло в низкочастотных полосах, в частности в интервале [10,20] Гц. Соответственно, следующий MRA фокусируется на компонентах сигнала в диапазонах частот от [10,80] Гц. В этих полосах идентифицируйте линейные изменения в данных. Постройте точки перехода наряду с компонентами MRA.

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

Основанный на MRA changepoint анализ помог подтвердить анализ 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-файлов. Удаленная метка Speed и только "раскалывается" или сохраненная метка "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