下述代码分别是对原始光谱数据进行一系列处理,最终得到了WDP数据,以及进行了原始数据的熔深数据提取(用空间密度趋势进行了深度曲线的提取,进行了可视化);另一个main和子函数则是对WDP数据依次进行了百分位滤波、LOF+百分位滤波、DBSCAN+百分位滤波,并进行了结果的可视化;请将2段代码进行合理拼接,因为中间的WDP数据是共享的呀,也就是第一段程序得到WDP后,第二段程序直接对WDP进行处理,最终呢呈现3个图:光谱数据处理过程图;A-scan曲线及包络提取图;空间密度趋势、百分位滤波、LOF+百分位滤波、DBSCAN+百分位滤波的熔深提取的对比图;
%% 探测器记录的是干涉光谱,通过物理关系:不同深度z的反射面会产生不同的相位延迟
%% 这个相位差会调制出不同频率的余弦振荡项在光谱上
%% 对光谱做傅里叶变换,就是把振荡频率转换成位置
%% 原始光谱-去直流项(去除S(k)背景)-k空间重采样-均匀间隔-加窗抑制旁瓣
%% ifft对齐0频-fft将k域信号转换为z域-abs取模得到包络-峰值对应反射层
%% 图2中有空间密度趋势,这个密度趋势是统计了B-scan 的趋势,由此进行统计,得到常见峰的位置;
%% 内置函数findpeaks用于提取显著反射层
%% ==================== 系统参数 ====================
lambda_center = 840e-9; % 中心波长 (m)
bandwidth_FWHM = 40e-9; % FWHM 带宽 (m)
sigma_lambda = bandwidth_FWHM / (2 * sqrt(2 * log(2))); % 高斯标准差计算
N_pixels = 2048; % 光谱仪像素数
lateral_range = 40e-3; % 扫描总长度 (m)
lateral_resolution = 20.5e-6; % 单步分辨率
lateral_step = 2 * lateral_resolution; % 采样间隔 ~41 μm
N_AScans = max(1, round(lateral_range / lateral_step));
% 构造波长向量(就把波长范围构造出来了)
lambda_min = lambda_center - 3 * sigma_lambda;
lambda_max = lambda_center + 3 * sigma_lambda;
%linsapce是生成一个包含N个元素的行向量,且在指定区间范围内均匀分布
lambda = linspace(lambda_min, lambda_max, N_pixels)';
% k空间基础(2pi/波长,进行重采样得基础)
k_vec = 2 * pi ./ lambda;
% 添加轻微非线性畸变(模拟 spectrometer 色散不均)
distortion = 1e4 * (k_vec - mean(k_vec)).^2 / max(k_vec)^2;
k_vec = k_vec + distortion;
% 深度轴(对称,ZPD 在中心)
c = 3e8;
n_eff = 1.0;
dz = pi / (n_eff * (max(k_vec) - min(k_vec)));
z_axis_full = (-N_pixels/2 : N_pixels/2 - 1)' * dz;
z_mm = z_axis_full * 1e3; % mm
x_mm = (0:N_AScans-1)' * lateral_step * 1e3;
% 噪声参数
read_noise_std = 0.01;
shot_noise_factor = 1.0; % 控制散粒噪声强度
source_instability = 0.03;
min_peak_height_ratio = 0.05;
min_peak_distance_base = round(N_pixels / 50);
fprintf('🔍 波长范围: %.2f ~ %.2f nm\n', min(lambda)*1e9, max(lambda)*1e9);
fprintf('📏 轴向分辨率: %.2f μm\n', (2*log(2)/pi) * lambda_center^2 / bandwidth_FWHM * 1e6);
fprintf('🔄 正在模拟 %d 条 A-scan...\n', N_AScans);
%% ==================== 数据源选择 ====================
use_real_data = false; % <<<<< 用户可改为 true 以启用真实数据
if use_real_data
try
load('user_spectrum.mat', 'lambda', 'raw_spectra_matrix');
if size(raw_spectra_matrix, 1) ~= N_pixels
error('光谱行数必须等于 N_pixels = %d', N_pixels);
end
if size(raw_spectra_matrix, 2) == 1
raw_spectra_matrix = repmat(raw_spectra_matrix, 1, N_AScans);
elseif size(raw_spectra_matrix, 2) < N_AScans
raw_spectra_matrix = raw_spectra_matrix(:, 1:N_AScans);
warning('原始光谱列数不足,已截取前%d列', N_AScans);
else
raw_spectra_matrix = raw_spectra_matrix(:, 1:N_AScans);
end
fprintf('✅ 成功加载用户光谱数据:%d pixels × %d A-scans\n', size(raw_spectra_matrix,1), size(raw_spectra_matrix,2));
catch ME
fprintf('❌ 加载用户光谱失败: %s\n', ME.message);
fprintf('🔁 切换至模拟光谱模式...\n');
use_real_data = false;
end
end
%% ==================== 初始化变量 ====================
bscan_intensity = zeros(N_pixels, N_AScans);
all_detected_points = cell(1, N_AScans);
example_spectrum_raw = [];
example_spectrum_proc = [];
example_ascan = [];
%% ==================== 主循环:模拟采集 ====================
for idx = 1:N_AScans
try
%% Step 1: 构建物理结构
surface_depth = -25e-6 + 50e-6 * rand(); % ±25μm 浮动
% 熔池散射点:密度随深度递减(指数分布)
%randi用于生成伪随机整数的函数,基于均匀分布
num_scatterers = randi([5, 12]);
%exprand生成均值为u的指数分布随机数
depths_norm = exprnd(1.8e-3, [num_scatterers, 1]);
depths_pool = surface_depth + 1.0e-3 + depths_norm;
%rand生成随机数,后面是数组,因此定义生成的是列向量,且每一个均在0-1之间
reflections_pool = 0.08 + 0.07 * rand(num_scatterers, 1);
bottom_depth = surface_depth + 2.8e-3 + 0.4e-3*rand();
depths_all = [surface_depth; depths_pool; bottom_depth];
reflections_all = [0.35; reflections_pool; 0.4]; % 表面强反射
% 5%概率添加虚假深层点
if rand() < 0.05
fake_deep = bottom_depth + 0.1e-3 + 0.3e-3*rand();
depths_all = [depths_all; fake_deep];
reflections_all = [reflections_all; 0.05];
end
%% Step 2: 生成或获取光谱
if use_real_data
raw_spectrum = raw_spectra_matrix(:, idx); % 使用实测数据(应为列向量)
else
raw_spectrum = generate_spectrum(lambda, k_vec, lambda_center, sigma_lambda, ...
depths_all, reflections_all, read_noise_std, shot_noise_factor, source_instability);
end
% 确保为列向量
raw_spectrum = raw_spectrum(:);
% 记录第一条有效的原始光谱用于可视化
if isempty(example_spectrum_raw)
example_spectrum_raw = raw_spectrum; % 已经是列向量
end
%% Step 3: 处理链
dc_removed = remove_dc(raw_spectrum);
[k_uniform, spec_resampled] = resample_to_uniform_k(k_vec, dc_removed);
if any(isnan(spec_resampled))
error('NaN detected in resampling');
end
spec_shaped = apply_window(spec_resampled, 'hanning');
ascan_envelope = compute_ascan(k_uniform, spec_shaped);
if isempty(example_spectrum_proc)
example_spectrum_proc = spec_resampled(:);
example_ascan = ascan_envelope(:);
end
bscan_intensity(:, idx) = ascan_envelope;
%% Step 4: 提取反射点(仅 z ≥ 0)
%这里是对z进行大于0的逻辑判断,因此返回的是一个真假的向量
valid_positive = (z_mm >= 0);
%zmm还是原来的,但这就提取了只有大于0的了
z_pos = z_mm(valid_positive);
%提取位于正深度区域的包络强度值
a_pos = ascan_envelope(valid_positive);
%设置2个相邻的峰之间的最小距离(就是峰不能太密)
min_dist = max(5, round(length(z_pos)/50));
%峰的阈值、间隔点、距离、阈值
[peaks, locs] = findpeaks(a_pos, 'MinPeakHeight', max(a_pos)*0.05, 'MinPeakDistance', min_dist);
%从数组索引转换为物理深度值
detected_z = z_pos(locs);
%去除靠近图像检测边缘的检测结果
detected_z = detected_z(detected_z <= max(z_mm)*0.95); % 截断边缘
%将本次扫描提取到的所有反射深度保存到单元数组中
all_detected_points{idx} = detected_z;
catch ME
fprintf('❌ A-scan %d 失败: %s\n', idx, ME.message);
bscan_intensity(:, idx) = eps;
all_detected_points{idx} = [];
end
end
fprintf('✅ 采集完成!\n');
%% ==================== 数据后处理:构建 pts_clean 并导出 CSV ====================
valid_mask = ~cellfun(@isempty, all_detected_points);
if any(valid_mask)
% 提取所有 z 值
pts_cat = vertcat(all_detected_points{valid_mask});
% 创建对应的 x_mm 向量(按每条 A-scan 复制)
x_repeated = arrayfun(@(i) repmat(x_mm(i), numel(all_detected_points{i}), 1), ...
find(valid_mask), 'UniformOutput', false);
x_pts = vertcat(x_repeated{:});
% 组合为 Nx2 [x_mm, z_mm]
pts_full = [x_pts, pts_cat];
% 过滤不合理深度
pts_clean = pts_full(pts_full(:,2) >= 0 & pts_full(:,2) <= max(z_mm)*0.95, :);
% 保存为 CSV
writematrix(pts_clean, 'WDP.csv');
fprintf('📁 已导出所有检测点至 bscan_detected_points.csv (共 %d 个点)\n', size(pts_clean,1));
else
pts_clean = [];
fprintf('⚠️ 无有效检测点,跳过CSV导出。\n');
end
%% ==================== 图1:五阶段处理链 ====================
figure('Name', '图1:五阶段处理链', 'Position', [100, 100, 1100, 700]);
subplot(3,2,1); plot(lambda*1e9, example_spectrum_raw); title('① 原始光谱'); xlabel('\lambda (nm)');
subplot(3,2,2); plot(k_vec, example_spectrum_raw); title('② k空间(有畸变,因为不均匀)'); xlabel('k (rad/m)');
dc_removed = remove_dc(example_spectrum_raw); subplot(3,2,3); plot(k_vec, dc_removed); title('③ 去DC');
[k_uni, spec_r] = resample_to_uniform_k(k_vec, dc_removed); subplot(3,2,4); plot(k_uni, spec_r); title('④ 均匀k(线性插值)');
w_spec = apply_window(spec_r, 'hanning'); subplot(3,2,5); plot(k_uni, w_spec); title('⑤ 加窗去旁瓣');
env = compute_ascan(k_uni, w_spec); subplot(3,2,6); plot(z_mm, env, 'c'); ax = gca; ax.XAxisLocation='origin'; grid on; title('⑥ A-scan包络');
try exportgraphics(gcf, 'figure1_chain.png', 'Resolution', 300); end
%% ==================== 图2:A-scan 分析 ====================
figure('Name', '图2:反射密度分析', 'Position', [200, 200, 900, 500]);
valid_region = (z_mm >= 0);
z_pos = z_mm(valid_region);
a_pos = example_ascan(valid_region);
min_peak_height = max(a_pos) * 0.1;
[~, pk_idx] = findpeaks(a_pos, 'MinPeakHeight', min_peak_height, 'MinPeakDistance', 5);
yyaxis left
plot(z_mm, example_ascan, 'b', 'LineWidth', 1.2);
hold on;
if ~isempty(pk_idx)
plot(z_pos(pk_idx), a_pos(pk_idx), 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 6);
end
ylabel('反射强度');
ylim auto;
bin_width = 0.05;
bins = 0:bin_width:max(z_mm);
[den, ~] = histcounts(pts_clean, bins);
den_smooth = smoothdata(den, 'gaussian', 5);
bin_centers = bins(1:end-1) + bin_width/2;
yyaxis right
plot(bin_centers, den_smooth, 'k--', 'LineWidth', 2);
ylabel('反射点密度(平滑)');
title('A-scan 包络 + 检测峰值 + 空间密度分布');
xlabel('深度 z (mm)');
legend('A-scan包络', '检测到的峰', '空间密度趋势', 'Location', 'best');
grid on;
hold off;
%% ==================== 图3:B-scan 与 熔深 ====================
figure('Name', '图3:B-scan 与 熔深提取', 'Position', [100, 100, 1200, 800]);
img_log = log10(bscan_intensity + 1e-10);
subplot(2,1,1); imagesc(x_mm, z_mm, img_log); axis image; colorbar; colormap(flipud(gray));
title('B-scan: 对数反射强度'); ylabel('深度 z (mm)'); hold on;
for idx = 1:N_AScans
x = x_mm(idx);
pts = all_detected_points{idx};
if ~isempty(pts)
plot(repmat(x, length(pts), 1), pts, 'y.', 'MarkerSize', 3);
end
end
% 提取熔深:取每个 A-scan 最深的有效峰
melt_depth_curve = NaN(size(x_mm));
for idx = 1:N_AScans
pts = all_detected_points{idx};
if isempty(pts), continue; end
[vals, ~] = sort(pts, 'descend');
for v = vals'
if v > 1.0 && v < max(z_mm)*0.95
melt_depth_curve(idx) = v;
break;
end
end
end
% 插值修复 NaN
valid_md = ~isnan(melt_depth_curve);
if sum(valid_md) >= 2
melt_depth_curve = interp1(x_mm(valid_md), melt_depth_curve(valid_md), x_mm, 'linear');
else
melt_depth_curve = fillmissing(melt_depth_curve, 'nearest');
end
plot(x_mm, melt_depth_curve, 'r-', 'LineWidth', 2.5);
legend('反射点', '估计熔深', 'Location', 'best');
hold off;
subplot(2,1,2);
plot(x_mm, melt_depth_curve, 'b-', 'LineWidth', 2);
xlabel('x (mm)'); ylabel('熔深 (mm)'); title('熔深曲线'); grid on;
% 安全设置 ylim
valid_vals = melt_depth_curve(~isnan(melt_depth_curve));
if isempty(valid_vals)
ylim([0, 1]); % 默认范围
else
ylim([0, max(valid_vals)*1.1]);
end
sgtitle('📌 B-scan 与 熔深提取:ZPD对齐表面,熔深为正');
%% ==================== 输出统计与保存 ====================
fprintf('\n📊 统计信息:\n');
fprintf(' 总A-scan数量: %d\n', N_AScans);
avg_pts = mean(cellfun(@numel, all_detected_points));
fprintf(' 平均每A-scan检测到 %.1f 个反射点\n', avg_pts);
if all(isnan(melt_depth_curve))
fprintf('⚠️ 未能提取任何有效熔深!\n');
else
range_md = [min(melt_depth_curve, [], 'omitnan'), max(melt_depth_curve, [], 'omitnan')];
fprintf(' 提取熔深范围: %.2f ~ %.2f mm\n', range_md(1), range_md(2));
fprintf(' 平均熔深: %.2f mm\n', mean(melt_depth_curve, 'omitnan'));
end
save('sd_oct_simulation_result.mat', 'bscan_intensity', 'melt_depth_curve', 'x_mm', 'z_mm', 'all_detected_points');
fprintf('💾 结果已保存至 sd_oct_simulation_result.mat\n');
disp('🎉 SD-OCT 高保真仿真完成!');
%% ========== 子函数1:生成光谱 ==========
function spectrum = generate_spectrum(lambda, k_vec, lambda_center, sigma_lambda, depths, reflections, noise_std, shot_factor, instability)
N = length(lambda);
I0 = exp(-(lambda - lambda_center).^2 / (2*sigma_lambda^2));
interference = zeros(N,1);
for i = 1:length(depths)
phase = 2 * k_vec * depths(i);
interference = interference + reflections(i) * cos(phase);
end
I_total = I0 + interference;
I_total = max(I_total, 0);
noise_readout = noise_std * randn(N,1);
noise_shot = shot_factor * sqrt(I_total) .* randn(N,1);
noise_drift = instability * I0 * randn();
spectrum = I_total + noise_readout + noise_shot + noise_drift;
end
%% ========== 子函数2:去直流 ==========
function y = remove_dc(x)
background = mean(x(1:2048));
y = x - background;
end
%% ========== 子函数3:k空间重采样 ==========
function [k_new, y_new] = resample_to_uniform_k(k, y)
k_new = linspace(min(k), max(k), length(k))';
y_new = interp1(k, y, k_new, 'pchip', 'extrap');
end
%% ========== 子函数4:加窗抑制旁瓣 ==========
function y = apply_window(x, type_str)
w = window(@hann, length(x));
y = x .* w(:);
end
%% ========== 子函数5:计算A-scan包络 ==========
function envelope = compute_ascan(k_uniform, spectrum)
%spectrum先确保为列向量
%ifft将数据“绕一圈”,使得0频成分移到首位,满足FFT的周期性假设
spec_shifted = ifftshift(spectrum(:));
%fft执行快速傅里叶变换,输出一个复数向量,模表示反射强度,相位可用于Doppler分析;
ascan_complex = fft(spec_shifted);
%取模得到包络信号,即看到的山峰,每一个山峰代表一个反射界面
envelope = abs(ascan_complex);
end
%%上边是第一个代码,接下来是第二个代码:
function main(defaultPath)
% weldingDataAnalysis_SingleWDP - 仅处理 WDP 文件的熔深滤波对比分析
%
% 功能:
% - 查找指定路径下唯一的 WDP 文件
% - 读取其熔深列
% - 执行三种滤波策略:
% 1. 百分位滤波(原始)
% 2. LOF 异常去除 + 插值 + 百分位滤波
% 3. DBSCAN 异常去除 + 插值 + 百分位滤波
% - 可视化三者结果对比
%% 参数初始化
if nargin < 1 || isempty(defaultPath) || ~ischar(defaultPath)
defaultPath = 'D:\OCT\20250919data'; % 默认路径
end
extensions = {'*.csv', '*.txt', '*.dat'};
keyword = 'WDP'; % 只关心 WDP 文件
fprintf('🔍 开始分析 WDP 文件...\n');
%% 步骤1:选择目录并查找 WDP 文件
folderPath = step0selectFolder(defaultPath);
targetFile = findWDPFile(folderPath, extensions, keyword);
if isempty(targetFile)
error('❌ 未找到包含 "%s" 的数据文件,请检查路径或文件命名!', keyword);
end
%% 步骤2:安全读取 WDP 数据
rawData = step2readDataSafely(targetFile);
if isempty(rawData)
error('❌ 无法加载 WDP 数据!');
end
fprintf('✅ 已加载 WDP 数据:%s [%d × %d]\n', targetFile, size(rawData,1), size(rawData,2));
% 提取时间/索引 和 熔深值
if size(rawData, 2) >= 2
t_full = rawData(:, 1);
y_full = rawData(:, 2);
else
t_full = (1:length(rawData))';
y_full = rawData(:);
end
% 截取中间段用于分析(避免首尾噪声)
subLength_analysis = 5000;
len = length(y_full);
startIdx = floor((len - subLength_analysis) / 2) + 1;
endIdx = startIdx + subLength_analysis - 1;
x_seg = t_full(startIdx:endIdx);
y_seg = y_full(startIdx:endIdx);
%% 步骤3:执行百分位滤波获取参数(弹窗一次)
[~, ~, filterParams] = step6performPercentileFiltering(y_seg, 'ShowPlot', false);
pct_val = filterParams.Percentile;
win_len = filterParams.WindowLength;
%% 步骤4:执行 LOF + 百分位滤波
try
[~, lofOutlierMask, ~, ~] = step7performLOFOnWDP(...
y_seg, ...
'PreExtracted', true, ...
'TimeIndex', x_seg, ...
'KFactor', 0.02, ...
'ShowPlot', false);
catch ME
warning('LOF 失败: %s。跳过该方法。', ME.message);
lofOutlierMask = false(size(y_seg));
end
% LOF 去噪插值(使用整数索引)
y_lof_clean = y_seg;
y_lof_clean(lofOutlierMask) = NaN;
idx_all = 1:length(y_seg); % 关键:使用索引而不是 x_seg
validIdx = ~isnan(y_lof_clean);
y_lof_interp = interp1(idx_all(validIdx), y_lof_clean(validIdx), idx_all, 'pchip');
% 补充边界 NaN
nanIdx = isnan(y_lof_interp);
if any(nanIdx)
y_lof_interp(nanIdx) = interp1(idx_all(~nanIdx), y_lof_interp(~nanIdx), ...
idx_all(nanIdx), 'nearest');
end
%% 步骤5:执行 DBSCAN + 百分位滤波
try
[dbscanOutlierMask_temp, ~, ~, ~] = step9performDBSCANOnWDP(...
y_seg, ...
'PreExtracted', true, ...
'TimeIndex', x_seg, ...
'MinPts', 10, ...
'Eps', 0.8, ...
'ShowPlot', false);
dbscanOutlierMask = logical(dbscanOutlierMask_temp);
catch ME
warning('DBSCAN 失败: %s。跳过该方法。', ME.message);
dbscanOutlierMask = false(size(y_seg));
end
% DBSCAN 去噪插值
y_dbscan_clean = y_seg;
y_dbscan_clean(dbscanOutlierMask) = NaN;
validIdx = ~isnan(y_dbscan_clean);
y_dbscan_interp = interp1(idx_all(validIdx), y_dbscan_clean(validIdx), idx_all, 'pchip');
% 补全
nanIdx = isnan(y_dbscan_interp);
if any(nanIdx)
y_dbscan_interp(nanIdx) = interp1(idx_all(~nanIdx), y_dbscan_interp(~nanIdx), ...
idx_all(nanIdx), 'nearest');
end
%% 步骤6:统一执行滑动百分位滤波
filtered_original = step8_1slidingPercentileFilter(y_seg, win_len, pct_val);
filtered_lof = step8_1slidingPercentileFilter(y_lof_interp, win_len, pct_val);
filtered_dbscan = step8_1slidingPercentileFilter(y_dbscan_interp, win_len, pct_val);
%% 步骤7:三重滤波结果可视化
figure('Position', [100, 100, 900, 600], 'Name', '【三重滤波对比】WDP 熔深处理结果', 'Color', 'white');
set(gcf, 'DefaultAxesFontName', 'Microsoft YaHei');
set(gcf, 'GraphicsSmoothing', 'on');
hold off;
subplot(3,1,1);
plot(x_seg, y_seg, 'Color', [0.7 0.7 0.7], 'LineWidth', 0.8); hold on;
plot(x_seg, filtered_original, 'k-', 'LineWidth', 1.8, 'DisplayName', sprintf('%g%% 滤波', pct_val));
title(sprintf('1. 原始信号 → %g%% 百分位滤波(窗口=%d)', pct_val, win_len));
xlabel('时间 / 索引'); ylabel('熔深 (mm)');
legend('show', 'Location', 'bestoutside'); grid on; box on;
subplot(3,1,2);
plot(x_seg, y_lof_interp, 'Color', [0.9 0.9 0.9], 'LineWidth', 0.8); hold on;
plot(x_seg, filtered_lof, 'b-', 'LineWidth', 1.8, 'DisplayName', ['LOF+' sprintf('%g%%滤波', pct_val)]);
scatter(x_seg(lofOutlierMask), y_seg(lofOutlierMask), 30, 'r', 'filled', 'MarkerFaceAlpha', 0.7, 'DisplayName', 'LOF异常点');
title('2. LOF去噪 → 百分位滤波');
xlabel('时间 / 索引'); ylabel('熔深 (mm)');
legend('show', 'Location', 'bestoutside'); grid on; box on;
subplot(3,1,3);
plot(x_seg, y_dbscan_interp, 'Color', [0.9 0.9 0.9], 'LineWidth', 0.8); hold on;
plot(x_seg, filtered_dbscan, 'g-', 'LineWidth', 1.8, 'DisplayName', ['DBSCAN+' sprintf('%g%%滤波', pct_val)]);
scatter(x_seg(dbscanOutlierMask), y_seg(dbscanOutlierMask), 30, 'm', 'diamond', 'filled', 'MarkerFaceAlpha', 0.7, 'DisplayName', 'DBSCAN异常点');
title('3. DBSCAN去噪 → 百分位滤波');
xlabel('时间 / 索引'); ylabel('熔深 (mm)');
legend('show', 'Location', 'bestoutside'); grid on; box on;
% 提取文件名(不含路径和扩展名)
[~, fileNameOnly, ~] = fileparts(targetFile);
sgtitle({'📈 WDP 熔深信号:三种滤波策略对比'; ...
sprintf('文件: %s', fileNameOnly)}, ...
'FontSize', 12, 'FontWeight', 'bold');
drawnow;
disp('🎉 所有处理完成!已生成三重滤波对比图。');
end
function filePath = findWDPFile(folderPath, extensions, keyword)
% 查找第一个匹配 *WDP* 的文件
filePath = '';
fprintf('🔍 正在搜索 WDP 文件...\n');
for e = 1:length(extensions)
ext = extensions{e};
list = dir(fullfile(folderPath, ext));
for k = 1:length(list)
fname = list(k).name;
if contains(fname, keyword, 'IgnoreCase', true)
filePath = fullfile(folderPath, fname);
fprintf('✅ 找到 WDP 文件: %s\n', fname);
return;
end
end
end
if isempty(filePath)
fprintf('❌ 未找到包含 "%s" 的文件。\n', keyword);
end
end
%% ========== 子函数 0: step0selectFolder ==========
function folderPath = step0selectFolder(defaultPath)
%% 参数处理与默认值
%判断输入参数个数,或者输入是否为空或者非字符串类型
if nargin < 1 || isempty(defaultPath) || ~ischar(defaultPath)
%pwd是内置函数,用于使用当前工作目录为默认路径,确保程序有可用的起点
defaultPath = pwd;
%isfolder是函数,用于检查指定路径是否存在,并且判断是不是一个文件夹
elseif ~isfolder(defaultPath)
warning('⚠️ 指定的默认路径不存在,将使用当前工作目录:\n%s', defaultPath);
defaultPath = pwd;
end
fprintf('🔍 正在打开文件夹选择器,默认路径: %s\n', defaultPath);
titleStr = '请选择焊接数据所在文件夹';
%uigetdir内置函数:打开一个图形化的选择文件夹对话框,允许鼠标点击浏览并选择目录
folderPath = uigetdir(defaultPath, titleStr);
if folderPath == 0
error('❌ 用户取消了文件夹选择,程序终止。');
end
% 规范化路径分隔符
%strrep就是将原路径下的符号更改为指定符号
folderPath = strrep(folderPath, '\', '/');
fprintf('✅ 已选择工作目录: %s\n', folderPath);
end
%% ========== 子函数 2: step2readDataSafely ==========
function data = step2readDataSafely(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
% 兼容旧版本写法
%detectImoportOptions内置函数:分析指定文本文件的内容结构,并生成一个导入选项对象
opts = detectImportOptions(filePath);
%readmatrix内置函数:用于读取纯数值矩阵数据的函数
data = readmatrix(filePath, opts); % 自动处理分隔符、标题行等
%isnumeric判断数据是否为数值类型
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';
%readtable函数:读取为表格结构/类型
tbl = readtable(filePath, opts);
% 检查所有变量是否可转为数值
%这里也是相当于引出一个变量作为逻辑值的载体,表示某种含义与否
allNumeric = true;
%width函数:输入的列数的读取
for i = 1:width(tbl)
%tb1{:,i};提取第i列的所有行,并尝试将其内容转换为一个数值向量或字符数组
var = tbl{:,i};
%any()函数:判断至少是否至少含有某一个指定量
if ~isnumeric(var) || any(isnan(var)) % 注意:字符转数字失败会引入 NaN
allNumeric = false;
break;
end
end
if allNumeric
%table2array函数:将表格类型变量转换为普通的数值数组
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
%fileparts函数:分解文件路径字符串的内置函数,分解为路径、名字、类型
[~, ~, ext] = fileparts(filePath);
%strcmpi函数:2个字符串是否相等,且不区分类型的大小写
if strcmpi(ext, '.csv')
%csvread函数:从指定的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);
%% ========== 子函数 6: step6performPercentileFiltering ==========
function [filteredY, x, finalParams] = step6performPercentileFiltering(rawData, varargin)
% performPercentileFiltering - 对 WDP 数据进行百分位滤波并可视化
%
% 输入:
% rawData - Nx1 或 Nx2 数值矩阵
% 名值对参数:
% 'Percentile' - 百分位数 (0~100),默认弹窗输入
% 'WindowLength' - 滑动窗口长度(奇数),默认弹窗输入
% 'ShowPlot' - 是否绘图,true/false,默认 true
%
% 输出:
% filteredY - 滤波后的 y 序列
% x - 对应的时间/索引序列
%% 默认参数解析
%创建一个inputParser对象,后续所有参数都添加到这个对象中进行管理
p = inputParser;
%addOptional给对象p添加可选位置参数
addOptional(p, 'Percentile', [], @isnumeric);
%这里对名称(名称对应的位置),名称的类型都作了规定,下面几条语句也是的
addOptional(p, 'WindowLength', [], @isnumeric);
%silogical表明传入的一定要是逻辑型
addOptional(p, 'ShowPlot', true, @islogical);
%将变长输入参数varargin拆开并传给parse函数,inputParser会自动按顺序匹配这些参数到前面定义的addOptional列表中
%执行验证,失败则抛出错误,成功后可通过p.Results获取结果
parse(p, varargin{:});
%将通过inputParser成功解析的输入参数,从p.Results结构体中取出并幅值给局部变量
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] = step6_1getFilterParamsFromDialog(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
%isscalar是用于判断是不是标量的,mod函数则可以用于判断输入量的第一位小数是不是0,以此判断是不是整数
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
%% 执行滑动窗口百分位滤波(镜像边界扩展)
%floor就是用于向下取整
halfWin = floor(windowLength / 2);
%下面这个是将整体拼接为一个向量,前向镜像填充+原始数据+后向镜像填充
%函数flipud就是进行列向量的上下翻转,fliplr就是进行行向量的上下翻转
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);
%prctile用于计算百分数,最终返回第输入百分位
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
% 返回最终使用的参数,便于后续复用
finalParams.Percentile = percentileValue;
finalParams.WindowLength = windowLength;
% 日志输出
fprintf('🔧 已完成 %.1f%% 百分位滤波处理,窗口大小=%d,数据长度=%d\n', ...
percentileValue, windowLength, n);
end
%% ========== 子函数 7: step7performLOFOnWDP==========
function [lofScores, isOutlier, x_sub, y_sub] = step7performLOFOnWDP(rawData, varargin)
% performLOFOnWDP_Optimized - 基于 y 和变化率的 LOF 异常检测(向量化加速 + 多阈值策略)
%
% 输入:
% rawData - Nx1 或 Nx2 数值矩阵,列2为熔深值;
% 或当 'PreExtracted'==true 时,可直接传入已提取的 y 向量(N×1)
%
% 名值对参数:
% 'TargetLength' - (default: 5000) 自动截取长度(中间段),仅在 PreExtracted=false 时有效
% 'UseMedianIQR' - (default: true) 是否使用中位数+IQR进行鲁棒阈值判定
% 'KFactor' - (default: 0.01) k = floor(n * KFactor),控制邻居数量
% 'ShowPlot' - (default: true) 是否绘制结果图
% 'PreExtracted' - (default: false) 是否已预提取信号段(rawData 即为 y_sub)
% 'TimeIndex' - (optional) 若 PreExtracted=true,则提供对应的时间轴(长度=N)
%
% 输出:
% lofScores - 每个样本的 LOF 得分(列向量)
% isOutlier - 逻辑数组,true 表示异常点(列向量)
% x_sub - 使用的时间索引或时间轴(列向量,便于绘图)
% y_sub - 实际用于分析的熔深值(列向量)
%% 参数解析
p = inputParser;
addOptional(p, 'TargetLength', 5000, @(x) isnumeric(x) && isscalar(x) && x > 0);
addOptional(p, 'UseMedianIQR', true, @islogical);
addOptional(p, 'KFactor', 0.01, @(x) isnumeric(x) && x > 0 && x <= 0.1);
addOptional(p, 'ShowPlot', true, @islogical);
addOptional(p, 'PreExtracted', false, @islogical);
addOptional(p, 'TimeIndex', [], @(x) isempty(x) || (isvector(x) && length(x) == numel(x)));
parse(p, varargin{:});
targetNum = p.Results.TargetLength;
useMedianIQR = p.Results.UseMedianIQR;
kFactor = p.Results.KFactor;
showPlot = p.Results.ShowPlot;
preExtracted = p.Results.PreExtracted;
timeIndex = p.Results.TimeIndex;
%% 输入验证
if isempty(rawData)
warning('⚠️ 输入数据为空!');
lofScores = []; isOutlier = []; x_sub = []; y_sub = [];
return;
end
if ~isnumeric(rawData) || ~ismatrix(rawData)
error('❌ rawData 必须是数值矩阵。');
end
% 提取熔深列:若非预提取模式且有多列,则取第2列
if size(rawData, 2) >= 2 && ~preExtracted
y_full = rawData(:, 2);
else
y_full = rawData(:); % 视为已提取的 y 向量
end
nTotal = length(y_full);
if nTotal < 3
error('❌ 数据点太少,无法进行 LOF 分析(至少需要3个点)。');
end
%% 【核心逻辑】根据是否预提取决定处理方式
if preExtracted
% 直接使用输入作为分析段
y_sub = y_full;
if isempty(timeIndex)
x_sub = (1:nTotal)'; % 默认使用索引
xlabelStr = '时间索引';
usePhysicalTime = false;
else
if length(timeIndex) ~= nTotal
error('❌ ''TimeIndex'' 长度必须与 rawData 一致 (%d vs %d)', length(timeIndex), nTotal);
end
x_sub = timeIndex(:);
xlabelStr = '时间 (s)';
usePhysicalTime = true;
end
fprintf('📊 执行 LOF 分析:使用用户提供的信号段,样本数=%d\n', nTotal);
else
% 原有逻辑:自动截取中间段
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_full(idxStart:idxEnd);
xlabelStr = '时间索引';
usePhysicalTime = false;
fprintf('📊 执行 LOF 分析:自动截取区间 [%d, %d],样本数=%d\n', idxStart, idxEnd, length(y_sub));
end
n = length(y_sub);
%% 特征工程:[熔深值, 一阶差分, 二阶差分]
dy = diff([y_sub; y_sub(end)]); % 补全最后一个 dy
ddy = diff([dy; dy(end)]); % 补全 ddy
features = [y_sub, dy, ddy]; % Nx3 特征矩阵
%% 归一化:Z-score 标准化
mu = mean(features, 1);
sigma = std(features, 0, 1);
sigma(sigma == 0) = 1; % 防止除零
data_norm = (features - mu) ./ sigma;
%% 自适应 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), rk_distance(j))
r_k_dist = D(sub2ind(size(D), repmat((1:n)', 1, k), kNN_idx));
reachDistMat = zeros(n, k);
for i = 1:n
for j_idx = 1:k
j = kNN_idx(i, j_idx);
reachDistMat(i, j_idx) = max(D(i, j), r_k_dist(j, j_idx));
end
end
% 平均可达距离
avgReachDist = mean(reachDistMat, 2);
%% 局部可达密度 LRD = 1 / avgReachDist
LRD = 1 ./ (avgReachDist + eps);
%% LOF 得分:LOF(i) = mean(LRD of kNN(i)) / LRD(i)
lofScores = zeros(n, 1);
for i = 1:n
neighbor_LRDs = LRD(kNN_idx(i, :));
lofScores(i) = mean(neighbor_LRDs) / LRD(i);
end
%% 自适应阈值判定
Q1 = prctile(lofScores, 25);
Q3 = prctile(lofScores, 75);
IQR = Q3 - Q1;
if useMedianIQR
threshold = median(lofScores) + 3 * IQR; % 更鲁棒
else
threshold = Q3 + 1.5 * 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', 'LineWidth', 2, ...
'DisplayName', sprintf('异常点(n=%d)', numOutliers));
xlabel(xlabelStr);
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), '%.5f'), ' ~ ', num2str(max(lofScores), '%.5f')]);
if numOutliers == 0
warning('🟡 未发现明显异常点,请检查信号是否过于平稳或需降低灵敏度');
end
end
%% ========== 子函数8_1: step8_1slidingPercentileFilter ==========
function result = step8_1slidingPercentileFilter(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 思路(若无工具箱则手动向量化)
%exist函数:检查类型文件是否存在
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
%% ========== 子函数 9: step9performDBSCANOnWDP ==========
function [isOutlier, x_sub, y_sub, eps_used] = step9performDBSCANOnWDP(rawData, varargin)
% performDBSCAN_OnWDP - 基于熔深和变化率的 DBSCAN 异常检测(支持预提取信号段)
%
% 输入:
% rawData - Nx1 或 Nx2 数值矩阵;
% 或当 'PreExtracted'==true 时,为已提取的 y 向量
%
% 名值对参数:
% 'TargetLength' - (default: 5000) 自动截取长度(仅 PreExtracted=false 时有效)
% 'Eps' - (optional) 邻域半径;若未指定,则自动估算
% 'MinPts' - (default: 5) 最小簇点数
% 'ShowPlot' - (default: true) 是否绘图
% 'PreExtracted' - (default: false) 是否已预提取信号段
% 'TimeIndex' - (optional) 若 PreExtracted=true,则提供对应时间轴
%
% 输出:
% isOutlier - 逻辑数组,true 表示异常点(噪声点)
% x_sub - 时间索引或时间轴(列向量,用于绘图)
% y_sub - 实际分析的熔深值(列向量)
% eps_used - 实际使用的 eps 值(便于记录)
%% 参数解析
p = inputParser;
addOptional(p, 'TargetLength', 5000, @(x) isnumeric(x) && isscalar(x) && x > 0);
addOptional(p, 'Eps', [], @(x) isempty(x) || (isnumeric(x) && x > 0));
addOptional(p, 'MinPts', 5, @(x) isnumeric(x) && isscalar(x) && x >= 1);
addOptional(p, 'ShowPlot', true, @islogical);
addOptional(p, 'PreExtracted', false, @islogical);
addOptional(p, 'TimeIndex', [], @(x) isempty(x) || (isvector(x) && length(x) == numel(x)));
parse(p, varargin{:});
targetNum = p.Results.TargetLength;
userEps = p.Results.Eps;
minPts = p.Results.MinPts;
showPlot = p.Results.ShowPlot;
preExtracted = p.Results.PreExtracted;
timeIndex = p.Results.TimeIndex;
%% 输入验证
if isempty(rawData)
warning('⚠️ 输入数据为空!');
isOutlier = []; x_sub = []; y_sub = []; eps_used = [];
return;
end
if ~isnumeric(rawData) || ~ismatrix(rawData)
error('❌ rawData 必须是数值矩阵。');
end
% 提取熔深列
if size(rawData, 2) >= 2 && ~preExtracted
y_full = rawData(:, 2);
else
y_full = rawData(:);
end
nTotal = length(y_full);
if nTotal < minPts + 1
error('❌ 数据点太少,无法运行 DBSCAN(至少需要 %d 个点)', minPts + 1);
end
%% 【核心逻辑】是否使用预提取信号段?
if preExtracted
y_sub = y_full(:); % 列向量
if isempty(timeIndex)
x_sub = (1:nTotal)';
xlabelStr = '时间索引';
usePhysicalTime = false;
else
if length(timeIndex) ~= nTotal
error('❌ ''TimeIndex'' 长度必须与输入 rawData 一致 (%d vs %d)', length(timeIndex), nTotal);
end
x_sub = timeIndex(:);
xlabelStr = '时间 (s)';
usePhysicalTime = true;
end
fprintf('📊 执行 DBSCAN 分析:使用用户提供的信号段,样本数=%d\n', nTotal);
else
% 自动截取中间段
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_full(idxStart:idxEnd);
y_sub = y_sub(:);
xlabelStr = '时间索引';
usePhysicalTime = false;
fprintf('📊 执行 DBSCAN 分析:自动截取区间 [%d, %d],样本数=%d\n', idxStart, idxEnd, length(y_sub));
end
n = length(y_sub);
if n < minPts + 1
error('❌ 截取后样本数不足,无法运行 DBSCAN(当前:%d,最少:%d)', n, minPts + 1);
end
%% 特征工程:[y, dy, ddy]
dy = diff([y_sub; y_sub(end)]);
ddy = diff([dy; dy(end)]);
features = [y_sub, dy, ddy];
%% 归一化:Z-score
mu = mean(features, 1);
sigma = std(features, 0, 1);
sigma(sigma == 0) = 1;
data_norm = (features - mu) ./ sigma;
%% 自动估算 eps(k-distance 图启发式)
if isempty(userEps)
% 计算每个点到其第 minPts 近邻的距离
D = pdist2(data_norm, data_norm);
[sortedD, ~] = sort(D, 2);
kDist = sortedD(:, minPts + 1); % 第 k+1 小的距离(排除自己)
kDistSorted = sort(kDist, 'descend');
% 简单启发:取下降最陡处附近(二阶差分最大)
diff1 = diff(kDistSorted);
diff2 = diff(diff1);
[~, idx] = max(abs(diff2)); % 找拐点
eps_est = kDistSorted(max(2, idx));
eps_used = round(eps_est * 1000) / 1000; % 保留三位小数
else
eps_used = userEps;
end
fprintf('🔍 使用 eps=%.4f, MinPts=%d 进行 DBSCAN 聚类...\n', eps_used, minPts);
%% 执行 DBSCAN
[idx_labels, coreIdx] = dbscan(data_norm, eps_used, minPts);
% DBSCAN 返回:
% idx_labels == -1 表示噪声点(异常)
% 其他为簇标签
isOutlier = (idx_labels == -1);
numOutliers = sum(isOutlier);
fprintf('✅ DBSCAN 完成:共发现 %d 个异常点(占比=%.1f%%)\n', numOutliers, 100*numOutliers/n);
%% 可视化
if showPlot
fig = figure('Position', [500, 300, 800, 500], ...
'Name', '【DBSCAN】异常检测结果', ...
'Color', 'white');
set(fig, 'DefaultAxesFontName', 'Microsoft YaHei');
set(fig, 'DefaultTextFontName', 'Microsoft YaHei');
set(fig, 'GraphicsSmoothing', 'on');
hold on; grid on; box on;
% 正常点(按簇着色)
uniqueLabels = unique(idx_labels(idx_labels ~= -1));
colors = lines(length(uniqueLabels));
for i = 1:length(uniqueLabels)
lbl = uniqueLabels(i);
clusterIdx = (idx_labels == lbl);
plot(x_sub(clusterIdx), y_sub(clusterIdx), '.', 'Color', colors(i,:), ...
'MarkerSize', 6, 'DisplayName', ['簇 #%d', num2str(lbl)]);
end
% 异常点(红色空心圈)
if any(isOutlier)
plot(x_sub(isOutlier), y_sub(isOutlier), 'ro', 'MarkerSize', 8, ...
'MarkerFaceColor', 'none', 'LineWidth', 2, ...
'DisplayName', sprintf('异常点(n=%d)', numOutliers));
end
xlabel(xlabelStr);
ylabel('熔深值 (mm)');
title(sprintf('DBSCAN 异常检测结果\neps=%.4f, MinPts=%d, 异常占比=%.1f%%', ...
eps_used, minPts, 100*numOutliers/n));
legend('show', 'Location', 'bestoutside');
hold off;
drawnow;
end
%% 输出统计
disp(['📈 DBSCAN 标签分布: 正常点=', num2str(n - numOutliers), ', 异常点=', num2str(numOutliers)]);
if numOutliers == 0
warning('🟡 未发现异常点,请检查 eps 是否过大或数据过于平稳');
end
end
最新发布