获取一个数字的上下百分之5的浮动值

public static void main(String[] args) {
        Random r = new Random();

        int s = 56;

        for (int i = 0; i < 10; i++) {
            // 第一种方法获取区间
            int num = (int) (Math.random() * 5 + 5);
            // 第二种方法获取区间
            int a = r.nextInt(5) + 10;

            double xiao = num * 0.1;
            double da = a * 0.1;

            double x = s * xiao;
            double d = s * da;

            System.out.println("最小值:" + (int)Math.rint(x) +"===值:" + (int)Math.rint(s) +"---最大值:"+ (int)Math.rint(d));

        }
    }

 

下述代码分别是对原始光谱据进行一系列处理,最终得到了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
最新发布
10-24
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值