В этом примере показано, как создать полный рабочий процесс оценки Остающегося срока полезного использования (RUL) включая шаги для предварительной обработки, выбора trendable функций, построения медицинского индикатора сплавом датчика, учебное подобие средства оценки RUL и проверка эффективности предзнаменований. Пример использует обучающие данные от набора данных проблемы PHM2008 от 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 столбцами. Столбцы содержат данные для ID машины, метки времени, 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-средних значений являются одним из самых популярных алгоритмов кластеризации, но он может привести к локальным оптимумам. Это - хорошая практика, чтобы повторить K-средних значений, кластеризирующих алгоритм несколько раз с различными начальными условиями и выбрать результаты с самой низкой ценой. В этом случае алгоритм запускается 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, потому что почти постоянное измерение датчика не полезно для остающейся оценки срока полезного использования. Обратитесь к последнему разделу, "Функции Помощника", для получения дополнительной информации о regimeNormalization функционируют.
trainDataNormalized = cellfun(@(data) regimeNormalization(data, centers, centerstats), ... trainData, 'UniformOutput', false);
Визуализируйте данные, нормированные на рабочий режим. Тренды ухудшения для некоторых измерений датчика теперь показаны после нормализации.
figure helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(1:4), nsample)
Теперь выберите самые trendable измерения датчика, чтобы создать медицинский индикатор для предсказания. Для каждого измерения датчика оценивается линейная модель ухудшения, и наклоны сигналов оцениваются.
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
Визуализируйте выбранные trendable измерения датчика.
figure helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(sensorTrended(3:6)), nsample)
Заметьте, что некоторые самые trendable сигналы показывают положительные тренды, в то время как другие показывают отрицательные тренды.
Этот раздел фокусируется на плавлении измерений датчика в один медицинский индикатор, с которым обучена основанная на подобии модель.
Все данные запуска к отказу приняты, чтобы запуститься со здорового условия. Состояние здоровья вначале присвоено, значение 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 нормой невязки
где медицинский индикатор машины , предполагаемый медицинский индикатор машины использование 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), Исследовательский центр Эймса, Поле Moffett, CA
[2] Roemer, Майкл Дж., Грегори Дж. Кацпрзынский и Майкл Х. Шоеллер. "Улучшенные диагностические и предвещающие оценки с помощью медицинского сплава информации управления". Продолжения AUTOTESTCON, 2001. Системная Технологическая Конференция по Готовности IEEE. IEEE, 2001.
[3] Goebel, Kai и Пьеро Бониссоне. "Предвещающий информационный сплав для постоянных систем загрузки". Информационный Fusion, 2 005 8-х Международных конференций по вопросам. Издание 2. IEEE, 2005.
[4] Ван, Пенг и Дэвид В. Койт. "Предсказание надежности на основе моделирования ухудшения для систем с несколькими мерами по ухудшению". Надежность и Поддерживаемость, 2 004 Ежегодных ПОРШНЯ Симпозиума. 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