%% ========== 子函数 1: selectFolder ==========
function folderPath = selectFolder(defaultPath)
% selectFolder - 交互式选择文件夹路径
% 输入:
% defaultPath: (可选) 默认打开的路径,若为空或无效则使用当前目录
% 输出:
% folderPath: 用户选择的文件夹路径(字符串)
%
% 示例:
% path = selectFolder('D:\OCT\20250919data');
%% 参数处理与默认值
if nargin < 1 || isempty(defaultPath) || ~ischar(defaultPath)
defaultPath = pwd; % 使用当前工作目录作为 fallback
elseif ~isfolder(defaultPath)
warning('⚠️ 指定的默认路径不存在,将使用当前工作目录。\n原路径: %s', defaultPath);
defaultPath = pwd;
end
%% 弹出文件夹选择器
fprintf('🔍 正在打开文件夹选择器,默认路径: %s\n', defaultPath);
titleStr = '请选择焊接数据所在文件夹'; % 自定义对话框标题
folderPath = uigetdir(defaultPath, titleStr);
%% 处理用户取消操作
if folderPath == 0
error('❌ 用户取消了文件夹选择,程序终止。');
end
% 确保路径为绝对路径且格式统一(使用正斜杠)
folderPath = fullfile(folderPath);
% 输出确认信息(可用于调试)
fprintf('✅ 已选择工作目录: %s\n', folderPath);
end
%% ========== 子函数 2: findTargetFiles ==========
function targetFiles = findTargetFiles(folderPath, extensions, keywords)
% 查找每个关键词对应的首个匹配文件(按扩展名遍历)
%
% 输入:
% folderPath - 目标文件夹路径(字符串)
% extensions - 扩展名通配符元胞数组,如 {'*.csv', '*.txt', '*.dat'}
% keywords - 关键词元胞数组,如 {'WD','WDP','WDS','WDSQ'}
% 输出:
% targetFiles - 元胞数组,每个元素对应找到的第一个匹配文件的完整路径;
% 未找到则为空字符串 ''
% cell就是用于创建一个空的元胞数据,这就相当于把关键词元胞数组的尺寸给复制下来
targetFiles = cell(size(keywords));
fprintf('🔍 正在搜索目录:%s\n', folderPath);
for i = 1:length(keywords)
% found标记作用:表示是否已找到匹配文件
found = false;
for e = 1:length(extensions)
% fullfile是MATLAB的函数,自动用操作系统兼容的分隔符进行路径拼接
searchPattern = fullfile(folderPath, extensions{e});
% dir(pattern)返回结构体数组,包含符合通配符的所有文件信息
files = dir(searchPattern);
for k = 1:length(files)
filename = files(k).name;
% contains用于判断字符串str是否包含字串pattern,IgnoreCase, true合起来是忽略大小写
if contains(filename, keywords{i}, 'IgnoreCase', true)
targetFiles{i} = fullfile(folderPath, filename);
fprintf('✅ 匹配第%d个文件: %s\n', i, filename);
% found标记找到
found = true;
break; % 跳出当前文件循环
end
end
if found
break; % 跳出当前扩展名循环(优先级高的先匹配到就停止)
end
end
% ~是逻辑非操作:如果仍未找到,发出警告
if ~found
warning('⚠️ 未找到包含 "%s" 的文件!', keywords{i});
targetFiles{i} = '';
end
end
end
%% ========== 子函数 3: readDataSafely ==========
function data = readDataSafely(filePath)
% readDataSafely - 安全读取数值型表格数据(支持 CSV/TXT/DAT/XLSX)
% 输入:
% filePath - 文件路径(字符串)
% 输出:
% data - 数值型矩阵(double),MxN
%
% 特性:
% - 自动检测分隔符和编码
% - 支持多种格式(包括 Excel)
% - 优先使用现代函数(readmatrix)
% - 兜底策略清晰
% - 提供详细错误信息
%% 输入验证
if nargin < 1 || isempty(filePath) || ~ischar(filePath)
error('❌ 输入参数 filePath 必须是非空字符串。');
end
if ~isfile(filePath)
error('❌ 指定的文件不存在:%s', filePath);
end
fprintf('📊 正在读取数据文件: %s\n', filePath);
%% 第一尝试:使用 readmatrix(推荐方式,自动处理分隔符)
try
data = readmatrix(filePath, 'DetectImportOptions', true);
if isnumeric(data) && ~isempty(data)
fprintf('✅ 成功使用 readmatrix 读取数据,尺寸: %d×%d\n', size(data));
return;
else
warning('⚠️ readmatrix 返回了空或非数值数据,尝试备用方法...');
end
catch ME1
fprintf('💡 readmatrix 失败: %s\n', ME1.message);
end
%% 第二尝试:readtable + 转换(兼容复杂表结构)
try
opts = detectImportOptions(filePath);
opts.VariableNamingRule = 'preserve';
tbl = readtable(filePath, opts);
% 检查所有变量是否可转为数值
allNumeric = true;
for i = 1:width(tbl)
var = tbl{:,i};
if ~isnumeric(var) || any(isnan(var)) % 注意:字符转数字失败会引入 NaN
allNumeric = false;
break;
end
end
if allNumeric
data = table2array(tbl);
fprintf('✅ 成功用 readtable 读取并转换为数值矩阵,尺寸: %d×%d\n', size(data));
return;
else
warning('⚠️ 表格中包含非数值列,跳过 readtable 方法...');
end
catch ME2
fprintf('💡 readtable 失败: %s\n', ME2.message);
end
%% 第三尝试:旧式 csvread(仅限 .csv 文件且纯数字)
try
[~, ~, ext] = fileparts(filePath);
if strcmpi(ext, '.csv')
data = csvread(filePath);
if isnumeric(data) && ~isempty(data)
fprintf('✅ 成功用 csvread 读取数据(兼容模式),尺寸: %d×%d\n', size(data));
return;
end
end
catch ME3
fprintf('💡 csvread 失败: %s\n', ME3.message);
end
%% 所有方法均失败
error('❌ 所有读取方法均失败,无法加载文件:\n%s', filePath);
%% ========== 子函数 4: subsampleData ==========
function result = subsampleData(data, step)
% subsampleData - 通用降采样函数(支持多种模式)
% 输入:
% data - 二维数值矩阵(MxN),每行为一个样本
% step - 正整数,采样间隔;即每 step 行取一行
% 输出:
% result - 降采样后的矩阵(KxN),K ≈ ceil(M/step)
%% 输入验证
if nargin < 2
error('❌ 缺少必要输入参数:需要 data 和 step。');
end
if isempty(data)
fprintf('🔍 输入数据为空,返回空矩阵。\n');
result = [];
return;
end
if ~isnumeric(data) || ndims(data) ~= 2
error('❌ data 必须是非空二维数值矩阵。');
end
if ~isscalar(step) || ~isnumeric(step) || step <= 0 || floor(step) ~= step
error('❌ step 必须是正整数。');
end
nRows = size(data, 1);
if nRows == 1
fprintf('📌 数据仅有一行,直接返回原数据。\n');
result = data;
return;
elseif step == 1
fprintf('💡 step=1,无需降采样,返回原始数据。\n');
result = data;
return;
elseif step >= nRows
fprintf('⚠️ step (%d) ≥ 数据行数 (%d),仅保留第一行。\n', step, nRows);
result = data(1, :);
return;
else
% 标准降采样:从第1行开始,每隔 step-1 行取一次
indices = 1:step:nRows;
result = data(indices, :);
fprintf('✅ 成功降采样:原 %d 行 → 采样后 %d 行(step=%d)\n', nRows, length(indices), step);
end
end
%% ========== 子函数 5: visualizeMultiChannel ==========
function fig = visualizeMultiChannel(allData, titles, ylabels, validIndices)
% visualizeMultiChannel - 多通道数据可视化(支持空数据+自动布局)
%
% 输入:
% allData - 元胞数组,每个元素是一个 Nx1 或 Nx2 数值矩阵
% titles - 元胞数组,每通道标题(字符串)
% ylabels - 元胞数组,每通道 Y 轴标签
% validIndices - 正整数向量,表示哪些通道有有效数据(如 [1,2,4])
%
% 输出:
% fig - 图形对象句柄
%% 输入验证
narginchk(4, 4);
if ~iscell(allData) || ~iscell(titles) || ~iscell(ylabels)
error('❌ allData, titles, ylabels 必须是元胞数组。');
end
numPlots = length(allData);
if length(titles) ~= numPlots || length(ylabels) ~= numPlots
error('❌ titles 和 ylabels 的长度必须与 allData 一致。');
end
if isempty(validIndices)
warning('⚠️ validIndices 为空,将尝试绘制所有通道...');
validIndices = 1:numPlots;
end
%% 创建图形窗口
fig = figure('Position', [100, 100, 800, min(1000, 220 * numPlots)], ...
'Name', '多通道数据可视化', ...
'Color', 'white');
% ✅ 正确方式:设置默认字体和字号(影响所有后续子对象)
set(fig, 'DefaultAxesFontName', 'Microsoft YaHei'); % 中文友好字体
set(fig, 'DefaultAxesFontSize', 10);
set(fig, 'DefaultTextFontName', 'Microsoft YaHei');
set(fig, 'DefaultTextFontSize', 10);
% 启用抗锯齿提升线条质量
set(fig, 'GraphicsSmoothing', 'on');
%% 循环绘制每个子图
for i = 1:numPlots
subplot(numPlots, 1, i);
% 判断该通道是否有有效数据
hasValidData = any(validIndices == i) && ~isempty(allData{i});
if ~hasValidData
text(0.5, 0.5, '无数据', 'HorizontalAlignment', 'center', ...
'VerticalAlignment', 'middle', ...
'FontSize', 14, 'Color', [0.8 0.2 0.2], 'FontWeight', 'bold');
title(titles{i}, 'Interpreter', 'none');
xlabel('时间/索引');
ylabel(ylabels{i});
box on;
grid on;
continue;
end
data = allData{i};
if ~isnumeric(data) || ndims(data) > 2
warning('⚠️ 通道 %d 数据格式异常,跳过...', i);
text(0.5, 0.5, '数据错误', 'HorizontalAlignment','center',...
'VerticalAlignment','middle','FontSize',12,'Color','r');
title(titles{i});
box on; grid on;
continue;
end
% 提取 x 和 y
if size(data, 2) >= 2
x = data(:, 1);
y = data(:, 2);
xLabel = '时间';
else
x = 1:size(data, 1);
y = data(:);
xLabel = '索引';
end
% 根据通道选择绘图类型
if i == 2 % WDP 特殊处理:散点图
scatter(x, y, 6, 'b', 'filled', 'MarkerFaceAlpha', 0.5);
xlabel(xLabel);
ylabel(ylabels{i});
title(titles{i}, 'Interpreter', 'none');
grid on; box on;
else % 其他通道:折线图
plot(x, y, 'Color', [0 0.4 1], 'LineWidth', 1.3);
xlabel(xLabel);
ylabel(ylabels{i});
title(titles{i}, 'Interpreter', 'none');
grid on; box on;
end
end
%% 添加总标题
sgtitle('【焊接过程多源数据可视化】', 'FontSize', 14, 'FontWeight', 'bold', 'Color', [0.1 0.3 0.6]);
% 立即刷新显示
drawnow;
fprintf('📊 已成功绘制 %d 个通道的数据可视化图表。\n', numPlots);
%% ========== 子函数 6: visualizeComparison ==========
function fig = visualizeComparison(d1, d2, s1, s2)
% visualizeComparison - 熔深预测曲线 vs 实测点 对比图
%
% 输入:
% d1 - Nx1 或 Nx2 数值矩阵,拟合结果(x可选)
% d2 - Mx1 或 Mx2 数值矩阵,实测数据(x可选)
% s1 - 正数,d1 的降采样步长(用于图例标注)
% s2 - 正数,d2 的降采样步长
%
% 输出:
% fig - 图形窗口句柄
%% 输入验证
narginchk(4, 4);
% 验证数据
if ~isempty(d1) && (~isnumeric(d1) || ~ismatrix(d1))
error('❌ d1 必须是数值矩阵或空数组。');
end
if ~isempty(d2) && (~isnumeric(d2) || ~ismatrix(d2))
error('❌ d2 必须是数值矩阵或空数组。');
end
% 验证步长
if ~isscalar(s1) || ~isnumeric(s1) || s1 <= 0
error('❌ s1 必须是正数标量。');
end
if ~isscalar(s2) || ~isnumeric(s2) || s2 <= 0
error('❌ s2 必须是正数标量。');
end
%% 创建图形窗口
fig = figure('Position', [200, 200, 800, 500], ...
'Name', '熔深预测 vs 实测点对比', ...
'Color', 'white');
% 设置默认字体(防中文乱码)
set(fig, 'DefaultAxesFontName', 'Microsoft YaHei');
set(fig, 'DefaultTextFontName', 'Microsoft YaHei');
set(fig, 'DefaultAxesFontSize', 10);
set(fig, 'GraphicsSmoothing', 'on'); % 启用抗锯齿
hold on; grid on;
%% 绘制拟合曲线(红色实线)
if ~isempty(d1)
if size(d1, 2) >= 2
x1 = d1(:, 1);
y1 = d1(:, 2);
plot(x1, y1, 'Color', [1, 0, 0], 'LineWidth', 2, ...
'DisplayName', sprintf('拟合曲线(间隔=%g)', s1));
else
y1 = d1(:);
plot(y1, 'Color', [1, 0, 0], 'LineWidth', 2, ...
'DisplayName', sprintf('拟合序列(间隔=%g)', s1));
end
else
text(0.5, 0.5, '无拟合数据', 'HorizontalAlignment', 'center', ...
'VerticalAlignment', 'middle', 'Color', [0.6 0.2 0.2], ...
'FontSize', 12, 'FontWeight', 'bold');
end
%% 绘制实测点(蓝色填充散点)
if ~isempty(d2)
if size(d2, 2) >= 2
x2 = d2(:, 1);
y2 = d2(:, 2);
scatter(x2, y2, 10, 'filled', 'MarkerFaceColor', [0, 0, 0.8], ...
'MarkerEdgeColor', [0, 0, 0.6], ...
'DisplayName', sprintf('实测点(间隔=%g)', s2), ...
'MarkerFaceAlpha', 0.7);
else
y2 = d2(:);
scatter(1:length(y2), y2, 10, 'filled', 'MarkerFaceColor', [0, 0, 0.8], ...
'MarkerEdgeColor', [0, 0, 0.6], ...
'DisplayName', sprintf('实测序列(间隔=%g)', s2), ...
'MarkerFaceAlpha', 0.7);
end
else
text(0.5, 0.5, '无实测数据', 'HorizontalAlignment', 'center', ...
'VerticalAlignment', 'middle', 'Color', [0.2 0.2 0.6], ...
'FontSize', 12, 'FontWeight', 'bold');
end
%% 标签与标题
xlabel('时间索引 / 时间戳');
ylabel('熔深值 (mm)');
title({'熔深:拟合曲线 vs 实测点'; ...
sprintf('降采样策略:预测间隔=%g,实测间隔=%g', s1, s2)}, ...
'FontWeight', 'bold');
% 图例
legend('show', 'Location', 'best', 'Box', 'off');
box on;
drawnow;
%% 日志输出(✅ 使用 if-else 替代 ?:)
if isempty(d1)
str_d1 = '空';
else
str_d1 = [num2str(size(d1,1)) '行'];
end
if isempty(d2)
str_d2 = '空';
else
str_d2 = [num2str(size(d2,1)) '行'];
end
fprintf('📈 已完成熔深对比图绘制:拟合数据(%s) + 实测数据(%s),采样步长=(%g, %g)\n', ...
str_d1, str_d2, s1, s2);
%% ========== 子函数 7: performPercentileFiltering ==========
function [filteredY, x] = performPercentileFiltering(rawData, varargin)
% performPercentileFiltering - 对 WDP 数据进行百分位滤波并可视化
%
% 输入:
% rawData - Nx1 或 Nx2 数值矩阵
% 名值对参数:
% 'Percentile' - 百分位数 (0~100),默认弹窗输入
% 'WindowLength' - 滑动窗口长度(奇数),默认弹窗输入
% 'ShowPlot' - 是否绘图,true/false,默认 true
%
% 输出:
% filteredY - 滤波后的 y 序列
% x - 对应的时间/索引序列
%% 默认参数解析
p = inputParser;
addOptional(p, 'Percentile', [], @isnumeric);
addOptional(p, 'WindowLength', [], @isnumeric);
addOptional(p, 'ShowPlot', true, @islogical);
parse(p, varargin{:});
showPlot = p.Results.ShowPlot;
percentileValue = p.Results.Percentile;
windowLength = p.Results.WindowLength;
%% 输入验证
if isempty(rawData)
warning('⚠️ 无法进行滤波:输入数据为空!');
if showPlot
figure('Name', '滤波提示', 'Color', 'w');
text(0.5, 0.5, '无数据可处理', 'HorizontalAlignment', 'center', ...
'FontSize', 14, 'Color', 'r', 'FontWeight', 'bold');
title('百分位滤波失败');
end
filteredY = [];
x = [];
return;
end
if ~isnumeric(rawData) || ~ismatrix(rawData)
error('❌ rawData 必须是数值矩阵。');
end
% 提取 x 和 y
if size(rawData, 2) >= 2
x = rawData(:, 1);
y = rawData(:, 2);
else
x = (1:length(rawData))';
y = rawData(:);
end
n = length(y);
if n < 3
error('❌ 数据点太少,无法应用滑动窗口滤波(至少需要3个点)。');
end
%% 参数获取方式:优先名值对,否则弹窗
if isempty(percentileValue) || isempty(windowLength)
[p_val, w_len] = getFilterParamsFromDialog(percentileValue, windowLength);
if nargin == 1 && (isempty(p_val) || isempty(w_len))
warning('跳过滤波:用户取消输入。');
filteredY = y;
return;
end
percentileValue = p_val;
windowLength = w_len;
end
% 参数校验
if percentileValue < 0 || percentileValue > 100
error('❌ 百分位必须在 0~100 范围内!');
end
if windowLength < 3 || ~isscalar(windowLength) || mod(windowLength,1) ~= 0
error('❌ 窗口长度必须是大于等于3的整数!');
end
if mod(windowLength, 2) == 0
warning('⚠️ 窗口长度应为奇数,已自动调整为 %d', windowLength + 1);
windowLength = windowLength + 1;
end
%% 执行滑动窗口百分位滤波(镜像边界扩展)
halfWin = floor(windowLength / 2);
y_ext = [flipud(y(1:halfWin)); y; flipud(y(end-halfWin+1:end))];
filteredY = zeros(n, 1);
for i = halfWin + 1 : n + halfWin
windowData = y_ext(i - halfWin : i + halfWin);
filteredY(i - halfWin) = prctile(windowData, percentileValue);
end
%% 可视化
if showPlot
fig = figure('Position', [300, 300, 800, 500], ...
'Name', sprintf('%g%% 百分位滤波结果', percentileValue), ...
'Color', 'white');
set(fig, 'DefaultAxesFontName', 'Microsoft YaHei');
set(fig, 'DefaultTextFontName', 'Microsoft YaHei');
set(fig, 'GraphicsSmoothing', 'on');
hold on; grid on;
scatter(x, y, 6, [0.5 0.5 1], 'filled', ...
'MarkerFaceAlpha', 0.5, 'DisplayName', '原始实测点');
plot(x, filteredY, 'Color', [1, 0, 0], 'LineWidth', 2, ...
'DisplayName', sprintf('%g%% 百分位滤波曲线', percentileValue));
xlabel('时间 / 索引');
ylabel('熔深值 (mm)');
title(sprintf('WDP 数据滤波结果(窗口=%d, 百分位=%.1f%%)', windowLength, percentileValue));
legend('show', 'Location', 'best', 'Box', 'off');
box on;
drawnow;
end
% 日志输出
fprintf('🔧 已完成 %.1f%% 百分位滤波处理,窗口大小=%d,数据长度=%d\n', ...
percentileValue, windowLength, n);
end
%% ========== 子函数 7-1: getFilterParamsFromDialog ==========
% 局部函数:getFilterParamsFromDialog
% 必须放在文件末尾!
% ===========================================================================
function [p_val, w_len] = getFilterParamsFromDialog(p_default, w_default)
prompt = {'请输入百分位数 (建议 0-10):', '请输入滑动窗口长度 (推荐奇数):'};
dlgTitle = '📊 百分位滤波参数设置';
defaults = {p_default, w_default};
if isempty(defaults{1}), defaults{1} = '5'; end
if isempty(defaults{2}), defaults{2} = '101'; end
answer = inputdlg(prompt, dlgTitle, 1, defaults);
if isempty(answer)
p_val = [];
w_len = [];
else
p_val = str2double(answer{1});
w_len = str2double(answer{2});
if isnan(p_val)
error('❌ 百分位输入无效,请输入合法数字!');
end
if isnan(w_len)
error('❌ 窗口长度输入无效,请输入整数!');
end
end
end
%% ========== 子函数 8: performPercentileFiltering ==========
function [lofScores, isOutlier, x_sub, y_sub] = performLOFOnWDP_Optimized(rawData, varargin)
% performLOFOnWDP_Optimized - 基于 y 和变化率的 LOF 异常检测(向量化加速 + 多阈值策略)
%
% 输入:
% rawData - Nx1 或 Nx2 数值矩阵,列2为熔深值
% 名值对参数:
% 'TargetLength' - 截取长度,默认 5000
% 'UseMedianIQR' - 是否使用中位数 IQR 抗噪,默认 true
% 'KFactor' - k = floor(n * KFactor),默认 0.01
% 'ShowPlot' - 是否绘图,默认 true
%
% 输出:
% lofScores - 每个样本的 LOF 得分
% isOutlier - 逻辑数组,标记是否为异常点
% x_sub - 时间索引(截取后)
% y_sub - 熔深值(截取后)
%% 参数解析
p = inputParser;
addOptional(p, 'TargetLength', 5000, @isnumeric);
addOptional(p, 'UseMedianIQR', true, @islogical);
addOptional(p, 'KFactor', 0.01, @(x) isnumeric(x) && x > 0 && x <= 0.1);
addOptional(p, 'ShowPlot', true, @islogical);
parse(p, varargin{:});
targetNum = p.Results.TargetLength;
useMedianIQR = p.Results.UseMedianIQR;
kFactor = p.Results.KFactor;
showPlot = p.Results.ShowPlot;
%% 输入验证
if isempty(rawData)
warning('⚠️ 无法执行 LOF:输入数据为空!');
lofScores = []; isOutlier = []; x_sub = []; y_sub = [];
return;
end
if ~isnumeric(rawData) || ~ismatrix(rawData)
error('❌ rawData 必须是数值矩阵。');
end
% 提取 y
if size(rawData, 2) >= 2
y = rawData(:, 2);
else
y = rawData(:);
end
nTotal = length(y);
if nTotal < 3
error('❌ 数据点太少,无法进行 LOF 分析(至少需要3个点)。');
end
%% 截取中间段数据(避免首尾噪声影响)
idxStart = 1; idxEnd = nTotal;
if nTotal > targetNum
mid = floor(nTotal / 2);
half = floor(targetNum / 2);
idxStart = max(1, mid - half + 1);
idxEnd = min(nTotal, mid + half);
end
x_sub = (idxStart:idxEnd)';
y_sub = y(idxStart:idxEnd);
% 补充:若长度仍不足3,报错
if length(y_sub) < 3
error('❌ 截取后数据过短,无法进行 LOF 计算。');
end
%% 特征工程:[熔深值, 一阶差分, 二阶差分]
dy = diff([y_sub; y_sub(end)]); % 向下填充最后一个值
ddy = diff([dy; dy(end)]); % 二阶差分
data_raw = [y_sub, dy, ddy]; % 三维特征空间
%% 归一化:Z-score 标准化(每列独立)
mu = mean(data_raw, 1);
sigma = std(data_raw, 0, 1);
sigma(sigma == 0) = 1;
data_norm = (data_raw - mu) ./ sigma;
n = size(data_norm, 1);
%% 自适应 k 值选择
k = min(20, max(3, floor(kFactor * n)));
fprintf('📊 使用 k = %d 进行 LOF 计算(样本数=%d)...\n', k, n);
%% 向量化计算距离矩阵(欧氏距离)
D = pdist2(data_norm, data_norm); % n×n 距离矩阵
%% 获取 k-近邻索引(排除自身)
[~, sortedIdx] = sort(D, 2);
kNN_idx = sortedIdx(:, 2:k+1); % 每行前 k 个最近邻(跳过自己)
%% 向量化计算可达距离 reachDist(i,j) = max(d(i,j), d(j, kNN_j_last))
% 先获取每个点的第k个邻居的距离 → 即 rk-distance
r_k_dist = D(sub2ind(size(D), repelem((1:n)', 1), kNN_idx(:, end))); % shape: n×1
% 构造完整的可达距离矩阵(广播)
reachDistMat = zeros(n, n);
for i = 1:n
j_list = kNN_idx(i, :);
for j_idx = 1:k
j = j_list(j_idx);
reachDistMat(i, j) = max(D(i, j), r_k_dist(j)); % reachDist(i,j)
end
end
% 提取每个点对其 k-近邻的平均可达距离
avgReachDist = zeros(n, 1);
for i = 1:n
avgReachDist(i) = mean(reachDistMat(i, kNN_idx(i, :)));
end
%% 计算局部可达密度 LRD = 1 / avgReachDist
LRD = 1 ./ (avgReachDist + eps);
%% 计算 LOF 得分:LOF(i) = mean(LRD(nn)) / LRD(i)
lofScores = zeros(n, 1);
for i = 1:n
neighbor_LRDs = LRD(kNN_idx(i, :));
lofScores(i) = mean(neighbor_LRDs) / LRD(i);
end
%% 自适应阈值判定(支持 Median-IQR 提升鲁棒性)
Q1 = prctile(lofScores, 25);
Q3 = prctile(lofScores, 75);
IQR = Q3 - Q1;
if useMedianIQR
threshold = median(lofScores) + 3 * IQR; % 更稳健
else
threshold = Q3 + 1.5 * IQR; % 经典 IQR
end
isOutlier = lofScores > threshold;
numOutliers = sum(isOutlier);
fprintf('🔍 LOF 分析完成:共发现 %d 个异常点(阈值=%.3f)\n', numOutliers, threshold);
%% 可视化
if showPlot
fig = figure('Position', [400, 400, 800, 500], ...
'Name', '【优化版】LOF 异常检测结果', ...
'Color', 'white');
set(fig, 'DefaultAxesFontName', 'Microsoft YaHei');
set(fig, 'DefaultTextFontName', 'Microsoft YaHei');
set(fig, 'GraphicsSmoothing', 'on');
hold on; grid on; box on;
% 正常点:蓝色小点
plot(x_sub(~isOutlier), y_sub(~isOutlier), 'b.', 'MarkerSize', 4, ...
'DisplayName', '正常点');
% 异常点:红色大圆圈(空心边框)
plot(x_sub(isOutlier), y_sub(isOutlier), 'ro', 'MarkerSize', 8, ...
'MarkerFaceColor', 'none', 'MarkerEdgeColor', 'r', 'LineWidth', 2, ...
'DisplayName', sprintf('异常点(n=%d)', numOutliers));
xlabel('时间索引');
ylabel('熔深值 (mm)');
title(sprintf('LOF 异常检测结果\nk=%d, 阈值=%.3f, 异常占比=%.1f%%', ...
k, threshold, 100*numOutliers/n));
legend('show', 'Location', 'bestoutside');
hold off;
drawnow;
end
%% 输出统计信息
disp(['📈 LOF得分范围: ', num2str(min(lofScores)), ' ~ ', num2str(max(lofScores))]);
if numOutliers == 0
warning('🟡 未发现明显异常点,请检查信号是否过于平稳或需降低灵敏度');
end
end
%% ========== 子函数9: comparePercentileFilteringWithLOF ==========
function [filtered_original, filtered_after_lof, lofOutlierMask, x_seg, y_seg] = ...
comparePercentileFilteringWithLOF(rawData, varargin)
% comparePercentileFilteringWithLOF_Optimized - 对比两种滤波策略:
%
% 方法1: 原始信号 → 百分位滤波
% 方法2: LOF去异常 + 插值 → 百分位滤波
%
% 输入:
% rawData - Nx1 或 Nx2 数值矩阵
% 名值对参数:
% 'TargetLength' - 截取长度,默认 5000
% 'Percentile' - 百分位数,默认 5
% 'WindowLength' - 滑动窗口长度,默认 101
% 'ShowPlot' - 是否绘图,默认 true
%
% 输出:
% filtered_original - 方法一的滤波结果
% filtered_after_lof - 方法二的滤波结果
% lofOutlierMask - 异常点逻辑数组
% x_seg - 时间索引(截取段)
% y_seg - 原始熔深值(截取段)
%% 参数解析
p = inputParser;
addOptional(p, 'TargetLength', 5000, @isnumeric);
addOptional(p, 'Percentile', 5, @(x) isnumeric(x) && x >= 0 && x <= 100);
addOptional(p, 'WindowLength', 101, @(x) isnumeric(x) && mod(x,2)==1 && x >= 3);
addOptional(p, 'ShowPlot', true, @islogical);
parse(p, varargin{:});
targetNum = p.Results.TargetLength;
pct_value = p.Results.Percentile;
windowLen = p.Results.WindowLength;
showPlot = p.Results.ShowPlot;
%% 输入验证
if isempty(rawData)
warning('⚠️ 无法执行对比分析:输入数据为空!');
filtered_original = []; filtered_after_lof = []; lofOutlierMask = [];
x_seg = []; y_seg = [];
return;
end
if ~isnumeric(rawData) || ~ismatrix(rawData)
error('❌ rawData 必须是数值矩阵。');
end
% 提取 y
if size(rawData, 2) >= 2
y_full = rawData(:, 2);
else
y_full = rawData(:);
end
x_full = (1:length(y_full))';
nTotal = length(y_full);
if nTotal < 3
error('❌ 数据点太少,无法进行分析(至少需要3个点)。');
end
%% 截取中间段
idxStart = 1; idxEnd = nTotal;
if nTotal > targetNum
mid = floor(nTotal / 2);
half = floor(targetNum / 2);
idxStart = max(1, mid - half + 1);
idxEnd = min(nTotal, mid + half);
end
x_seg = x_full(idxStart:idxEnd)';
y_seg = y_full(idxStart:idxEnd)';
n = length(y_seg);
fprintf('🔍 开始对比分析:处理区间 [%d, %d],样本数=%d\n', idxStart, idxEnd, n);
%% 步骤1:使用优化版 LOF 检测异常点(复用已有函数)
try
[~, isOutlier, ~, ~] = performLOFOnWDP_Optimized([x_seg, y_seg], ...
'TargetLength', targetNum, ...
'KFactor', 0.01, ...
'UseMedianIQR', true, ...
'ShowPlot', false);
catch ME
warningId = 'WDPAnalysis:LOF:ExecutionFailed';
warning(warningId, 'LOF 异常检测执行失败:%s。将跳过去噪步骤并使用原始信号进行滤波。', ME.message);
isOutlier = false(n, 1); % 默认无异常点
end
numOutliers = sum(isOutlier);
lofOutlierMask = isOutlier;
fprintf('📊 LOF检测完成,发现 %d 个异常点\n', numOutliers);
%% 步骤2:构建两组输入信号
% --- 方法一:原始完整信号 ---
y_original = y_seg;
% --- 方法二:去除异常点后插值补全 ---
y_clean = y_seg;
y_clean(isOutlier) = NaN;
% 使用 pchip 插值(保形,减少振荡)
y_interpolated = interp1(x_seg(~isOutlier), y_clean(~isOutlier), x_seg, 'pchip');
% 边界处理:若首尾为NaN,则用最近邻填充
if any(isnan(y_interpolated))
missingIdx = isnan(y_interpolated);
y_interpolated(missingIdx) = interp1(x_seg(~missingIdx), y_interpolated(~missingIdx), ...
x_seg(missingIdx), 'nearest');
end
%% 步骤3:分别执行滑动百分位滤波(调用优化函数)
filtered_original = slidingPercentileFilter(y_original, windowLen, pct_value);
filtered_after_lof = slidingPercentileFilter(y_interpolated, windowLen, pct_value);
%% 步骤4:绘制双子图对比
if showPlot
fig = figure('Position', [100, 100, 800, 600], ...
'Name', '【对比】百分位滤波 vs LOF+百分位滤波', ...
'Color', 'white');
set(fig, 'DefaultAxesFontName', 'Microsoft YaHei');
set(fig, 'DefaultTextFontName', 'Microsoft YaHei');
set(fig, 'GraphicsSmoothing', 'on');
% 子图1:原始 → 百分位滤波
subplot(2,1,1);
hold on;
plot(x_seg, y_original, 'Color', [0.7 0.7 0.7], 'LineWidth', 0.8, 'DisplayName', '原始实测点');
plot(x_seg, filtered_original, 'k-', 'LineWidth', 1.8, 'DisplayName', sprintf('%g%% 滤波', pct_value));
title(sprintf('方法一:直接对原始数据进行 %g%% 百分位滤波', pct_value));
legend('show', 'Location', 'bestoutside');
grid on; box on; hold off;
ylabel('熔深值 (mm)');
% 子图2:LOF去噪后 → 百分位滤波
subplot(2,1,2);
hold on;
plot(x_seg, y_interpolated, 'Color', [0.9 0.9 0.9], 'LineWidth', 0.8, 'DisplayName', 'LOF去噪+插值');
plot(x_seg, filtered_after_lof, 'b-', 'LineWidth', 1.8, 'DisplayName', ['LOF+' sprintf('%g%%滤波', pct_value)]);
scatter(x_seg(isOutlier), y_seg(isOutlier), 30, 'r', 'filled', ...
'DisplayName', sprintf('LOF异常点(n=%d)', numOutliers), 'MarkerFaceAlpha', 0.8);
title(sprintf('方法二:先LOF去噪再进行 %g%% 百分位滤波', pct_value));
legend('show', 'Location', 'bestoutside');
grid on; box on; hold off;
xlabel('时间索引');
ylabel('熔深值 (mm)');
sgtitle({'💡 LOF预处理提升滤波稳定性'; ...
sprintf('窗口=%d, 百分位=%g%%, 异常点=%d/%d', windowLen, pct_value, numOutliers, n)}, ...
'FontSize', 12, 'FontWeight', 'bold');
drawnow;
end
disp('✅ 双子图对比已完成:显示了是否使用LOF预处理的影响。');
end
%% ========== 子函数9—1: slidingPercentileFilter ==========
function result = slidingPercentileFilter(signal, windowLen, pct)
% slidingPercentileFilter_Vec - 向量化滑动窗口百分位滤波
%
% 输入:
% signal - 一维信号
% windowLen - 奇数窗口长度
% pct - 百分位数 (0~100)
%
% 输出:
% result - 滤波后信号
n = length(signal);
result = zeros(n, 1);
halfWin = floor(windowLen / 2);
% 边界保持原样
result(1:halfWin) = signal(1:halfWin);
result(end-halfWin+1:end) = signal(end-halfWin+1:end);
% 中心区域:使用 movprctile 思路(若无工具箱则手动向量化)
if exist('movprctile', 'file') % Statistics and Machine Learning Toolbox
result(halfWin+1:n-halfWin) = movprctile(signal, pct, windowLen);
else
% 手动构造滑动窗口(利用 arrayfun 向量化)
idx_center = (halfWin + 1):(n - halfWin);
windows = arrayfun(@(i) signal(i-halfWin : i+halfWin), idx_center, 'UniformOutput', false);
result(halfWin+1:n-halfWin) = cellfun(@(w) prctile(w, pct), windows);
end
end
最新发布