В этом примере показано, как построить полный рабочий процесс оценки оставшегося срока полезного использования (RUL), включающий в себя этапы предварительной обработки, выбора модельных функций, построения индикатора работоспособности путем слияния датчиков, обучения оценщикам сходства RUL и проверки производительности прогностических функций. В примере используются обучающие данные из набора данных PHM2008 challenge из https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/ [1].
Поскольку набор данных мал, можно загрузить в память все данные ухудшения качества. Загрузите и распакуйте набор данных из https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/ в текущий каталог. Используйте helperLoadData вспомогательная функция для загрузки и преобразования обучающего текстового файла в массив ячеек расписаний. Обучающие данные содержат 218 имитаций пробега до отказа. Эта группа измерений называется «ансамблем».
degradationData = helperLoadData('train.txt');
degradationData(1:5)ans = 5×1 cell array
{223×26 table}
{164×26 table}
{150×26 table}
{159×26 table}
{357×26 table}
Каждый участник ансамбля представляет собой таблицу с 26 колонками. Столбцы содержат данные по идентификатору машины, отметке времени, 3 рабочим условиям и 21 измерению датчика.
head(degradationData{1})ans=8×26 table
id time op_setting_1 op_setting_2 op_setting_3 sensor_1 sensor_2 sensor_3 sensor_4 sensor_5 sensor_6 sensor_7 sensor_8 sensor_9 sensor_10 sensor_11 sensor_12 sensor_13 sensor_14 sensor_15 sensor_16 sensor_17 sensor_18 sensor_19 sensor_20 sensor_21
__ ____ ____________ ____________ ____________ ________ ________ ________ ________ ________ ________ ________ ________ ________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________
1 1 10.005 0.2501 20 489.05 604.13 1499.5 1310 10.52 15.49 394.88 2318.9 8770.2 1.26 45.4 372.15 2388.1 8120.8 8.6216 0.03 368 2319 100 28.58 17.174
1 2 0.0015 0.0003 100 518.67 642.13 1584.5 1404 14.62 21.61 553.67 2388 9045.8 1.3 47.29 521.81 2388.2 8132.9 8.3907 0.03 391 2388 100 38.99 23.362
1 3 34.999 0.8401 60 449.44 555.42 1368.2 1122.5 5.48 8 194.93 2222.9 8343.9 1.02 41.92 183.26 2387.9 8063.8 9.3557 0.02 334 2223 100 14.83 8.8555
1 4 20.003 0.7005 0 491.19 607.03 1488.4 1249.2 9.35 13.65 334.82 2323.8 8721.5 1.08 44.26 314.84 2388.1 8052.3 9.2231 0.02 364 2324 100 24.42 14.783
1 5 42.004 0.8405 40 445 549.52 1354.5 1124.3 3.91 5.71 138.24 2211.8 8314.6 1.02 41.79 130.44 2387.9 8083.7 9.2986 0.02 330 2212 100 10.99 6.4025
1 6 20.003 0.7017 0 491.19 607.37 1480.5 1258.9 9.35 13.65 334.51 2323.9 8711.4 1.08 44.4 315.36 2388.1 8053.2 9.2276 0.02 364 2324 100 24.44 14.702
1 7 42 0.84 40 445 549.57 1354.4 1131.4 3.91 5.71 139.11 2211.8 8316.9 1.02 42.09 130.16 2387.9 8082 9.3753 0.02 331 2212 100 10.53 6.4254
1 8 0.0011 0 100 518.67 642.08 1589.5 1407.6 14.62 21.61 553.48 2388.1 9050.4 1.3 47.5 521.74 2388 8133.3 8.4339 0.03 391 2388 100 38.98 23.234
Разбейте данные деградации на набор данных обучения и набор данных проверки для последующей оценки эффективности.
rng('default') % To make sure the results are repeatable numEnsemble = length(degradationData); numFold = 5; cv = cvpartition(numEnsemble, 'KFold', numFold); trainData = degradationData(training(cv, 1)); validationData = degradationData(test(cv, 1));
Укажите группы переменных, представляющих интерес.
varNames = string(degradationData{1}.Properties.VariableNames);
timeVariable = varNames{2};
conditionVariables = varNames(3:5);
dataVariables = varNames(6:26);Визуализация выборки данных ансамбля.
nsample = 10;
figure
helperPlotEnsemble(trainData, timeVariable, ...
[conditionVariables(1:2) dataVariables(1:2)], nsample)
Как показано в предыдущем разделе, отсутствует четкая тенденция, показывающая процесс деградации при каждом измерении от пробега до поломки. В этом и следующем разделах рабочие условия будут использоваться для извлечения более четких тенденций деградации из сигналов датчиков.
Обратите внимание, что каждый участник ансамбля содержит 3 условия работы: «op_setting_1,» «op_setting_2,» и «op_setting_3.». Сначала извлекем таблицы из каждой ячейки и объединим их в одну таблицу.
trainDataUnwrap = vertcat(trainData{:});
opConditionUnwrap = trainDataUnwrap(:, cellstr(conditionVariables));Визуализация всех рабочих точек на графике рассеяния 3D. Он ясно показывает 6 режимов, и точки в каждом режиме находятся в очень непосредственной близости.
figure helperPlotClusters(opConditionUnwrap)
Давайте используем методы кластеризации, чтобы найти 6 кластеров автоматически. Здесь используется алгоритм K-средства. K-means является одним из самых популярных алгоритмов кластеризации, но это может привести к локальной оптимизации. Рекомендуется повторять алгоритм кластеризации K-means несколько раз с различными исходными условиями и выбирать результаты с наименьшей стоимостью. В этом случае алгоритм работает 5 раз, и результаты идентичны.
opts = statset('Display', 'final'); [clusterIndex, centers] = kmeans(table2array(opConditionUnwrap), 6, ... 'Distance', 'sqeuclidean', 'Replicates', 5, 'Options', opts);
Replicate 1, 1 iterations, total sum of distances = 0.279547. Replicate 2, 1 iterations, total sum of distances = 0.279547. Replicate 3, 1 iterations, total sum of distances = 0.279547. Replicate 4, 1 iterations, total sum of distances = 0.279547. Replicate 5, 1 iterations, total sum of distances = 0.279547. Best total sum of distances = 0.279547
Визуализация результатов кластеризации и идентифицированных центроидов кластера.
figure helperPlotClusters(opConditionUnwrap, clusterIndex, centers)
Как показано на графике, алгоритм кластеризации успешно находит 6 рабочих режимов.
Давайте проведем нормализацию измерений, сгруппированных по различным рабочим режимам. Во-первых, вычислить среднее и стандартное отклонение каждого измерения датчика, сгруппированного по рабочим режимам, определенным в последнем разделе.
centerstats = struct('Mean', table(), 'SD', table()); for v = dataVariables centerstats.Mean.(char(v)) = splitapply(@mean, trainDataUnwrap.(char(v)), clusterIndex); centerstats.SD.(char(v)) = splitapply(@std, trainDataUnwrap.(char(v)), clusterIndex); end centerstats.Mean
ans=6×21 table
sensor_1 sensor_2 sensor_3 sensor_4 sensor_5 sensor_6 sensor_7 sensor_8 sensor_9 sensor_10 sensor_11 sensor_12 sensor_13 sensor_14 sensor_15 sensor_16 sensor_17 sensor_18 sensor_19 sensor_20 sensor_21
________ ________ ________ ________ ________ ________ ________ ________ ________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________
489.05 604.92 1502.1 1311.4 10.52 15.493 394.32 2319 8784.1 1.26 45.496 371.44 2388.2 8133.9 8.6653 0.03 369.74 2319 100 28.525 17.115
518.67 642.71 1590.7 1409.4 14.62 21.61 553.3 2388.1 9062.3 1.3 47.557 521.36 2388.1 8140.9 8.4442 0.03 393.27 2388 100 38.809 23.285
462.54 536.87 1262.8 1050.6 7.05 9.0275 175.4 1915.4 8014.9 0.93989 36.808 164.56 2028.3 7878 10.916 0.02 307.39 1915 84.93 14.262 8.5552
445 549.72 1354.7 1128.2 3.91 5.7158 138.62 2212 8327 1.0202 42.163 130.53 2388 8088.4 9.3775 0.02 331.12 2212 100 10.584 6.3515
491.19 607.59 1486 1253.6 9.35 13.657 334.46 2324 8729.1 1.0777 44.466 314.85 2388.2 8065 9.2347 0.022299 365.45 2324 100 24.447 14.67
449.44 555.82 1366.9 1131.9 5.48 8.0003 194.43 2223 8355.2 1.0203 41.995 183.01 2388.1 8071.1 9.3344 0.02 334.29 2223 100 14.827 8.8966
centerstats.SD
ans=6×21 table
sensor_1 sensor_2 sensor_3 sensor_4 sensor_5 sensor_6 sensor_7 sensor_8 sensor_9 sensor_10 sensor_11 sensor_12 sensor_13 sensor_14 sensor_15 sensor_16 sensor_17 sensor_18 sensor_19 sensor_20 sensor_21
__________ ________ ________ ________ __________ _________ ________ ________ ________ __________ _________ _________ _________ _________ _________ __________ _________ _________ __________ _________ _________
1.4553e-11 0.47617 5.8555 8.3464 1.1618e-12 0.0047272 0.65536 0.094487 18.057 1.51e-13 0.25089 0.52838 0.096183 16.012 0.037598 9.9929e-16 1.4723 0 0 0.14522 0.086753
4.059e-11 0.48566 5.9258 8.8223 3.7307e-14 0.0011406 0.86236 0.068654 20.061 1.2103e-13 0.25937 0.71239 0.069276 17.212 0.036504 9.9929e-16 1.5046 0 0 0.17565 0.10515
2.7685e-11 0.35468 5.2678 6.9664 2.043e-14 0.0043301 0.45074 0.2743 14.741 0.0010469 0.21539 0.33985 0.28932 13.624 0.044142 8.6744e-16 1.3083 0 5.9833e-12 0.11212 0.06722
0 0.44169 5.6853 7.5741 1.9763e-13 0.0049401 0.44198 0.30713 18.389 0.0014951 0.23584 0.3432 0.33144 16.917 0.037135 2.231e-15 1.4174 0 0 0.10778 0.063991
4.4456e-11 0.46992 5.7664 7.8679 8.9892e-13 0.0046633 0.59984 0.13032 17.983 0.0042388 0.2362 0.49055 0.13285 15.792 0.038156 0.0042084 1.4428 0 0 0.13352 0.079731
4.3489e-11 0.44341 5.7224 7.4842 3.7485e-13 0.0017642 0.4734 0.28889 17.608 0.0017978 0.23242 0.38243 0.31 15.91 0.038624 8.5009e-16 1.3974 0 0 0.11311 0.069348
Статистические данные в каждом режиме могут использоваться для нормализации данных обучения. Для каждого члена ансамбля извлеките рабочие точки каждого ряда, вычислите его расстояние до каждого центра кластера и найдите ближайший центр кластера. Затем для каждого измерения датчика вычитайте среднее значение и делите его на стандартное отклонение этого кластера. Если стандартное отклонение близко к 0, установите нормализованное измерение датчика на 0, потому что почти постоянное измерение датчика не полезно для оставшейся оценки полезного срока службы. Дополнительные сведения о функции нормализации см. в последнем разделе «Вспомогательные функции».
trainDataNormalized = cellfun(@(data) regimeNormalization(data, centers, centerstats), ... trainData, 'UniformOutput', false);
Визуализация данных, нормализованных рабочим режимом. Тренды деградации для некоторых измерений датчиков теперь обнаруживаются после нормализации.
figure helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(1:4), nsample)
Теперь выберите наиболее модные измерения датчиков для построения индикатора состояния для прогнозирования. Для каждого измерения датчика оценивается линейная модель деградации и ранжируются наклоны сигналов.
numSensors = length(dataVariables); signalSlope = zeros(numSensors, 1); warn = warning('off'); for ct = 1:numSensors tmp = cellfun(@(tbl) tbl(:, cellstr(dataVariables(ct))), trainDataNormalized, 'UniformOutput', false); mdl = linearDegradationModel(); % create model fit(mdl, tmp); % train mode signalSlope(ct) = mdl.Theta; end warning(warn);
Сортируйте наклоны сигнала и выберите 8 датчиков с наибольшими наклонами.
[~, idx] = sort(abs(signalSlope), 'descend');
sensorTrended = sort(idx(1:8))sensorTrended = 8×1
2
3
4
7
11
12
15
17
Визуализация выбранных модельных измерений датчика.
figure helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(sensorTrended(3:6)), nsample)
Обратите внимание, что некоторые из наиболее модных сигналов показывают положительные тенденции, в то время как другие показывают негативные тенденции.
Этот раздел посвящен слиянию измерений датчиков в единый индикатор состояния, с помощью которого обучается модель на основе подобия.
Предполагается, что все данные о сбое начинаются со здорового состояния. Состоянию работоспособности в начале присваивается значение 1, а состоянию работоспособности при сбое - значение 0. Предполагается, что состояние здоровья линейно ухудшается от 1 до 0 с течением времени. Эта линейная деградация используется для предохранения значений датчиков. Более сложные методы слияния датчиков описаны в литературе [2-5].
for j=1:numel(trainDataNormalized) data = trainDataNormalized{j}; rul = max(data.time)-data.time; data.health_condition = rul / max(rul); trainDataNormalized{j} = data; end
Визуализация состояния здоровья.
figure
helperPlotEnsemble(trainDataNormalized, timeVariable, "health_condition", nsample)
Состояние здоровья всех участников ансамбля изменяется от 1 до 0 с различной скоростью деструкции.
Теперь установите линейную регрессионную модель состояния здоровья с наиболее трендовыми измерениями датчиков в качестве регрессоров:
Состояние здоровья ~ 1 + Sensor2 + Sensor3 + Sensor4 + Sensor7 + Sensor11 + Sensor12 + Sensor15 + Sensor17
trainDataNormalizedUnwrap = vertcat(trainDataNormalized{:});
sensorToFuse = dataVariables(sensorTrended);
X = trainDataNormalizedUnwrap{:, cellstr(sensorToFuse)};
y = trainDataNormalizedUnwrap.health_condition;
regModel = fitlm(X,y);
bias = regModel.Coefficients.Estimate(1)bias = 0.5000
weights = regModel.Coefficients.Estimate(2:end)
weights = 8×1
-0.0308
-0.0308
-0.0536
0.0033
-0.0639
0.0051
-0.0408
-0.0382
Создайте один индикатор работоспособности путем умножения измерений датчика на соответствующие веса.
trainDataFused = cellfun(@(data) degradationSensorFusion(data, sensorToFuse, weights), trainDataNormalized, ... 'UniformOutput', false);
Визуализация слитого показателя здоровья для данных обучения.
figure helperPlotEnsemble(trainDataFused, [], 1, nsample) xlabel('Time') ylabel('Health Indicator') title('Training Data')
Данные от нескольких датчиков объединяются в единый индикатор состояния. Индикатор работоспособности сглаживается фильтром скользящего среднего. Посмотрите, что помощник функционирует «degradationSensorFusion» в последнем разделе для получения дополнительной информации.
Повторите процесс нормализации режима и слияния датчиков с набором данных проверки.
validationDataNormalized = cellfun(@(data) regimeNormalization(data, centers, centerstats), ... validationData, 'UniformOutput', false); validationDataFused = cellfun(@(data) degradationSensorFusion(data, sensorToFuse, weights), ... validationDataNormalized, 'UniformOutput', false);
Визуализация индикатора работоспособности для данных проверки.
figure helperPlotEnsemble(validationDataFused, [], 1, nsample) xlabel('Time') ylabel('Health Indicator') title('Validation Data')
Теперь создайте модель остаточного подобия RUL с использованием данных обучения. В этом параметре модель пытается подогнать все слитые данные полиномом 2-го порядка.
Расстояние между данными и данными вычисляется по 1-норме остатка
=||yj-yˆj,i||1
где - индикатор работоспособности машины , - оценочный индикатор работоспособности машины с использованием модели полинома 2-го порядка, идентифицированной в машине .
Оценка подобия вычисляется по следующей формуле:
i, j) 2)
Учитывая одного члена ансамбля в наборе данных проверки, модель найдет ближайшие 50 членов ансамбля в наборе данных обучения, подберет распределение вероятностей на основе 50 членов ансамбля и использует медиану распределения в качестве оценки RUL.
mdl = residualSimilarityModel(... 'Method', 'poly2',... 'Distance', 'absolute',... 'NumNearestNeighbors', 50,... 'Standardize', 1); fit(mdl, trainDataFused);
Чтобы оценить модель сходства RUL, используйте 50%, 70% и 90% данных проверки образца для прогнозирования его RUL.
breakpoint = [0.5, 0.7, 0.9];
validationDataTmp = validationDataFused{3}; % use one validation data for illustrationИспользуйте данные проверки перед первой точкой останова, которая составляет 50% срока службы.
bpidx = 1; validationDataTmp50 = validationDataTmp(1:ceil(end*breakpoint(bpidx)),:); trueRUL = length(validationDataTmp) - length(validationDataTmp50); [estRUL, ciRUL, pdfRUL] = predictRUL(mdl, validationDataTmp50);
Визуализируйте данные проверки, усеченные на 50%, и ближайшие соседи.
figure compare(mdl, validationDataTmp50);
Визуализируйте оцененный RUL по сравнению с истинным RUL и распределение вероятности оцененного RUL.
figure helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)
Существует относительно большая ошибка между оцененным значением RUL и истинным значением RUL, когда машина находится на промежуточной стадии работоспособности. В этом примере наиболее похожие 10 кривых близки в начале, но раздваиваются, когда они приближаются к состоянию отказа, что приводит примерно к двум режимам распределения RUL.
Используйте данные проверки перед второй точкой останова, которая составляет 70% срока службы.
bpidx = 2; validationDataTmp70 = validationDataTmp(1:ceil(end*breakpoint(bpidx)), :); trueRUL = length(validationDataTmp) - length(validationDataTmp70); [estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp70); figure compare(mdl, validationDataTmp70);
figure helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)
Когда наблюдается больше данных, оценка RUL улучшается.
Используйте данные проверки перед третьей точкой останова, которая составляет 90% срока службы.
bpidx = 3; validationDataTmp90 = validationDataTmp(1:ceil(end*breakpoint(bpidx)), :); trueRUL = length(validationDataTmp) - length(validationDataTmp90); [estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp90); figure compare(mdl, validationDataTmp90);
figure helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)
Когда машина близка к отказу, оценка RUL в этом примере еще более улучшается.
Теперь повторите ту же процедуру оценки для всего набора данных проверки и вычислите ошибку между расчетным значением RUL и истинным значением RUL для каждой точки останова.
numValidation = length(validationDataFused); numBreakpoint = length(breakpoint); error = zeros(numValidation, numBreakpoint); for dataIdx = 1:numValidation tmpData = validationDataFused{dataIdx}; for bpidx = 1:numBreakpoint tmpDataTest = tmpData(1:ceil(end*breakpoint(bpidx)), :); trueRUL = length(tmpData) - length(tmpDataTest); [estRUL, ~, ~] = predictRUL(mdl, tmpDataTest); error(dataIdx, bpidx) = estRUL - trueRUL; end end
Визуализируйте гистограмму ошибки для каждой точки останова вместе с ее распределением вероятности.
[pdf50, x50] = ksdensity(error(:, 1)); [pdf70, x70] = ksdensity(error(:, 2)); [pdf90, x90] = ksdensity(error(:, 3)); figure ax(1) = subplot(3,1,1); hold on histogram(error(:, 1), 'BinWidth', 5, 'Normalization', 'pdf') plot(x50, pdf50) hold off xlabel('Prediction Error') title('RUL Prediction Error using first 50% of each validation ensemble member') ax(2) = subplot(3,1,2); hold on histogram(error(:, 2), 'BinWidth', 5, 'Normalization', 'pdf') plot(x70, pdf70) hold off xlabel('Prediction Error') title('RUL Prediction Error using first 70% of each validation ensemble member') ax(3) = subplot(3,1,3); hold on histogram(error(:, 3), 'BinWidth', 5, 'Normalization', 'pdf') plot(x90, pdf90) hold off xlabel('Prediction Error') title('RUL Prediction Error using first 90% of each validation ensemble member') linkaxes(ax)
Постройте график ошибки прогнозирования в рамочном графике для визуализации медианы, 25-75 квантилей и отклонений.
figure boxplot(error, 'Labels', {'50%', '70%', '90%'}) ylabel('Prediction Error') title('Prediction error using different percentages of each validation ensemble member')
Вычислите и визуализируйте среднее и стандартное отклонение ошибки прогнозирования.
errorMean = mean(error)
errorMean = 1×3
-5.8944 3.1359 3.3555
errorMedian = median(error)
errorMedian = 1×3
-4.8538 5.3763 3.6580
errorSD = std(error)
errorSD = 1×3
26.4916 20.0720 18.0313
figure errorbar([50 70 90], errorMean, errorSD, '-o', 'MarkerEdgeColor','r') xlim([40, 100]) xlabel('Percentage of validation data used for RUL prediction') ylabel('Prediction Error') legend('Mean Prediction Error with 1 Standard Deviation Eror bar', 'Location', 'south')
Показано, что ошибка становится более концентрированной вокруг 0 (меньше отклонений), поскольку наблюдается больше данных.
[1] А. Саксена и К. Гебель (2008). «Набор данных проблемы PHM08», НАСА Репозиторий данных Предзнаменований Эймса (http://ti.arc.nasa.gov/project/prognostic-data-repository), Исследовательский центр Эймса, Моффетт-Филд, Калифорния
[2] Рёмер, Майкл Дж., Грегори Дж. Качпжиньский и Майкл Х. Шёллер. «Улучшенные диагностические и прогностические оценки с использованием слияния информации управления здоровьем». Разбирательство AUTOTESTCON, 2001. Конференция по технологиям готовности систем IEEE. IEEE, 2001.
[3] Гебель, Кай и Пьеро Бониссоне. «Прогностическое слияние информации для систем с постоянной нагрузкой». Информационный обмен, 2005 8-я Международная конференция по. Том 2. IEEE, 2005.
[4] Ван, Пэн и Дэвид У. Койт. «Прогнозирование надежности на основе моделирования деградации для систем с несколькими показателями деградации». Надежность и ремонтопригодность, ежегодный симпозиум 2004 года - РАМН. IEEE, 2004.
[5] Джардин, Эндрю КС, Даминг Лин и Драган Баневич. «Обзор диагностики оборудования и прогнозирования, реализующих обслуживание на основе состояния». Механические системы и обработка сигналов 20.7 (2006): 1483-1510.
function data = regimeNormalization(data, centers, centerstats) % Perform normalization for each observation (row) of the data % according to the cluster the observation belongs to. conditionIdx = 3:5; dataIdx = 6:26; % Perform row-wise operation data{:, dataIdx} = table2array(... rowfun(@(row) localNormalize(row, conditionIdx, dataIdx, centers, centerstats), ... data, 'SeparateInputs', false)); end function rowNormalized = localNormalize(row, conditionIdx, dataIdx, centers, centerstats) % Normalization operation for each row. % Get the operating points and sensor measurements ops = row(1, conditionIdx); sensor = row(1, dataIdx); % Find which cluster center is the closest dist = sum((centers - ops).^2, 2); [~, idx] = min(dist); % Normalize the sensor measurements by the mean and standard deviation of the cluster. % Reassign NaN and Inf to 0. rowNormalized = (sensor - centerstats.Mean{idx, :}) ./ centerstats.SD{idx, :}; rowNormalized(isnan(rowNormalized) | isinf(rowNormalized)) = 0; end function dataFused = degradationSensorFusion(data, sensorToFuse, weights) % Combine measurements from different sensors according % to the weights, smooth the fused data and offset the data % so that all the data start from 1 % Fuse the data according to weights dataToFuse = data{:, cellstr(sensorToFuse)}; dataFusedRaw = dataToFuse*weights; % Smooth the fused data with moving mean stepBackward = 10; stepForward = 10; dataFused = movmean(dataFusedRaw, [stepBackward stepForward]); % Offset the data to 1 dataFused = dataFused + 1 - dataFused(1); end