最后需要输出一个excel表格,输出的excel数据不仅仅是第二列提取的电压数据,我要把那一列中每个数据对应的行的数据都输出来。修改下面的代码,其他的都不变
%% 燃料电池数据自适应斜率采样分析
clc; clear; close all;
% 1. 读取Excel数据(带时间戳去重)
filename = 'FC1_Ageing_part.xlsx';
try
data = readtable(filename);
time = data{:,1}; % 第一列时间
voltage = data{:,2}; % 第二列电压
% 处理重复时间点
[time, uniqueIdx] = unique(time, 'stable');
voltage = voltage(uniqueIdx);
data = data(uniqueIdx, :);
catch ME
error('数据读取失败: %s\n请检查: 1)文件路径 2)Excel格式 3)列标题', ME.message);
end
% 2. 数据预处理
numPoints = length(time);
if numPoints < 10
error('数据点不足(最少需要10个点),当前仅%d点', numPoints);
end
% 3. 自适应斜率采样法提取关键点(5%数据)
targetPoints = max(3, round(0.05 * numPoints)); % 修改此处
selectedIndices = adaptiveSlopeSampling(time, voltage, targetPoints);
% 4. 确保采样点时间唯一性
sampled_time = time(selectedIndices);
sampled_voltage = voltage(selectedIndices);
% 5. 提取并输出完整行数据
sampledData = data(selectedIndices, :);
try
writetable(sampledData, 'sampled_data.xlsx');
disp('采样数据已保存: sampled_data.xlsx');
catch
warning('Excel写入失败,改为CSV格式输出');
writetable(sampledData, 'sampled_data.csv');
end
% 6. 使用稳健插值方法(PCHIP)
predicted_voltage = interp1(sampled_time, sampled_voltage, time, 'pchip', 'extrap');
% 7. 计算性能指标
[metrics, error] = calculateMetrics(voltage, predicted_voltage);
fprintf('性能指标:\nMSE=%.4e\nRMSE=%.4e\nMAE=%.4e\nR²=%.4f\n',...
metrics.mse, metrics.rmse, metrics.mae, metrics.r2);
%% 修改后的绘图部分代码
figure('Position', [100, 100, 1200, 400]);
% ===== 使用 tight_subplot 创建紧凑布局 =====
gap = [0.02, 0.06]; % 垂直间隙5%,水平间隙3%
marg_h = [0.12, 0.05]; % 下边距10%,上边距%
marg_w = [0.05, 0.015]; % 左右边距各5%
ha = tight_subplot(1, 2, gap, marg_h, marg_w);
% ---- 图1: 主对比图 ----
axes(ha(1)); % 激活第一个子图
hold on;
h1_raw = plot(time, voltage, 'k-', 'LineWidth', 1.5);
h1_pred = plot(time, predicted_voltage, 'r-', 'LineWidth', 1.5);
hold off;
xlabel('Time(h)');
ylabel('Voltage(V)');
title('Raw data Vs sampled data', 'FontSize', 12, 'FontWeight', 'bold');
grid on;
% 设置第一个子图的边框和坐标轴为黑色且加粗
set(gca, 'FontSize', 12, ...
'LineWidth', 1.5, ... % 边框线宽
'XColor', [0 0 0], ... % X轴颜色(黑色)
'YColor', [0 0 0], ... % Y轴颜色(黑色)
'Box', 'on', ... % 显示完整边框
'XTickLabelMode', 'auto', ... % 确保刻度标签显示
'YTickLabelMode', 'auto'); % 确保刻度标签显示
% 第一个子图图例位置控制(相对于子图的归一化坐标)
legend1_rel_pos = [0.05, 0.05, 0.3, 0.15]; % [左, 下, 宽, 高]
ax1_pos = get(gca, 'Position');
legend1_abs_pos = [
ax1_pos(1) + legend1_rel_pos(1)*ax1_pos(3),
ax1_pos(2) + legend1_rel_pos(2)*ax1_pos(4),
legend1_rel_pos(3)*ax1_pos(3),
legend1_rel_pos(4)*ax1_pos(4)
];
legend([h1_raw, h1_pred], {'Raw data', 'Sample data'}, 'Position', legend1_abs_pos);
% ---- 图2: 误差分布图 ----
axes(ha(2)); % 激活第二个子图
hold on;
h2_error = plot(time, error, 'Color', [0, 0.447, 0.741], 'LineWidth', 1.5);
h2_zero = plot([0, 1200], [0, 0], 'r--', 'LineWidth', 1.5, ...
'DisplayName', 'Zero-error line');
hold off;
xlabel('Time(h)');
ylabel('Error(V)');
title('Error distribution', 'FontSize', 12, 'FontWeight', 'bold');
grid on;
% 设置第二个子图的边框和坐标轴为黑色且加粗
set(gca, 'FontSize', 12, ...
'LineWidth', 1.5, ... % 边框线宽
'XColor', [0 0 0], ... % X轴颜色(黑色)
'YColor', [0 0 0], ... % Y轴颜色(黑色)
'Box', 'on', ... % 显示完整边框
'XTickLabelMode', 'auto', ... % 确保刻度标签显示
'YTickLabelMode', 'auto'); % 确保刻度标签显示
% 第二个子图图例位置控制(相对于子图的归一化坐标)
legend2_rel_pos = [0.05, 0.05, 0.3, 0.15]; % [左, 下, 宽, 高]
ax2_pos = get(gca, 'Position');
legend2_abs_pos = [
ax2_pos(1) + legend2_rel_pos(1)*ax2_pos(3),
ax2_pos(2) + legend2_rel_pos(2)*ax2_pos(4),
legend2_rel_pos(3)*ax2_pos(3),
legend2_rel_pos(4)*ax2_pos(4)
];
legend([h2_error, h2_zero], {'Error', 'Zero-error line'}, 'Position', legend2_abs_pos);
% ===== 文本框位置计算(基于第二个子图)=====
axPos = get(ha(2), 'Position'); % 获取归一化位置 [left, bottom, width, height]
% 定义文本框相对位置参数
box_rel_left = 0.40; % 距离左边5%
box_rel_bottom = 0.05; % 距离底部65%
box_rel_width = 0.25; % 宽度25%
box_rel_height = 0.25; % 高度25%
% 计算文本框绝对位置
box_left = axPos(1) + box_rel_left * axPos(3);
box_bottom = axPos(2) + box_rel_bottom * axPos(4);
box_width = box_rel_width * axPos(3);
box_height = box_rel_height * axPos(4);
% 定义性能指标文本内容
metrics_text = sprintf([
'RMSE = %.6f\n'......
'MSE = %.6f\n'...
'MAE = %.6f\n'...
'R² = %.4f'],...
metrics.mse, metrics.rmse, metrics.mae, metrics.r2);
% ===== 创建半透明背景的文本框(固定在误差分布图内)=====
annotation('textbox', [box_left, box_bottom, box_width, box_height], ...
'String', metrics_text, ...
'FitBoxToText', 'on', ...
'BackgroundColor', [1, 1, 1, 0.85], ... % 半透明白色背景
'EdgeColor', [0.5, 0.5, 0.5], ... % 灰色边框
'LineWidth', 1.5, ...
'FontSize', 10, ...
'FontWeight', 'bold', ...
'Margin', 8, ...
'VerticalAlignment', 'top', ...
'HorizontalAlignment', 'left');
% 9. 保存采样信息到Excel <<< 修改输出内容
samplingInfo = table(length(sampled_time), numPoints, 100*length(sampled_time)/numPoints, ...
'VariableNames', {'SampledPoints', 'TotalPoints', 'SamplingRatePercent'});
writetable(samplingInfo, 'sampling_info.xlsx');
disp('采样信息已保存: sampling_info.xlsx');
%% 自适应斜率采样函数
function indices = adaptiveSlopeSampling(time, data, targetPoints)
n = length(data);
indices = [1, n]; % 始终包含首尾点
% 计算初始斜率
slopes = diff(data) ./ diff(time);
% 计算斜率变化率 (二阶差分)
slopeChanges = abs(diff(slopes));
% 如果目标点数过少,使用固定阈值
if targetPoints <= 3
indices = sort([1, round(n/2), n]);
return;
end
% 动态调整阈值算法
maxIter = 100;
tol = 0.02 * targetPoints; % 2%容忍度
lowThresh = min(slopeChanges);
highThresh = max(slopeChanges);
iter = 0;
while iter < maxIter
currentThresh = (lowThresh + highThresh) / 2;
candidatePoints = find(slopeChanges > currentThresh) + 1; % +1补偿索引
% 包含首尾点并确保唯一性
candidatePoints = unique([1; candidatePoints(:); n]);
numSelected = length(candidatePoints);
if abs(numSelected - targetPoints) <= tol
indices = candidatePoints;
break;
elseif numSelected > targetPoints
lowThresh = currentThresh; % 提高阈值减少点数
else
highThresh = currentThresh; % 降低阈值增加点数
end
iter = iter + 1;
end
% 确保精确点数(避免重复)
if length(indices) > targetPoints
% 优先保留变化最大的点
[~, sortedIdx] = sort(slopeChanges, 'descend');
topPoints = sortedIdx(1:min(targetPoints-2, length(sortedIdx))) + 1;
indices = sort([1; topPoints; n]);
end
% 最终确保索引唯一性
indices = unique(indices);
end
%% 性能指标计算函数
function [metrics, error] = calculateMetrics(actual, predicted)
error = actual - predicted;
metrics.mse = mean(error.^2);
metrics.rmse = sqrt(metrics.mse);
metrics.mae = mean(abs(error));
ss_res = sum(error.^2);
ss_tot = sum((actual - mean(actual)).^2);
if ss_tot == 0
metrics.r2 = 1; % 处理常数值情况
else
metrics.r2 = 1 - (ss_res / ss_tot);
end
end