%% ==================== 完整融合版：OCT光谱 → WDP → 多策略滤波对比 ====================
function fused_WDP_analysis_pipeline()
% fused_WDP_analysis_pipeline
% 功能：端到端仿真 OCT 数据采集 → 提取 WDP → 多方法熔深滤波对比
%
% 输出三张图：
% Figure 1: 光谱处理五步链
% Figure 2: A-scan 分析 + 空间密度趋势
% Figure 3: 四种熔深提取与滤波结果对比图
%% ==================== 系统参数 ====================
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;
lambda = linspace(lambda_min, lambda_max, N_pixels)';
% k空间基础（2pi/波长）
k_vec = 2 * pi ./ lambda;
% 添加轻微非线性畸变
distortion = 1e4 * (k_vec - mean(k_vec)).^2 / max(k_vec)^2;
k_vec = k_vec + distortion;
% 深度轴（ZPD 在中心）
dz = pi / (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;
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);
%% ==================== 初始化变量 ====================
bscan_intensity = zeros(N_pixels, N_AScans);
all_detected_points = cell(1, N_AScans);
example_spectrum_raw = [];
example_spectrum_proc = [];
example_ascan = [];
%% ==================== 主循环：模拟采集 ====================
% 初始化错误日志
error_log = struct('frame_idx', {}, 'message', {}, 'stack', {});
for idx = 1:N_AScans
try
%% Step 1: 构建物理结构
surface_depth = -25e-6 + 50e-6 * rand(); % ±25μm 浮动
num_scatterers = randi([5, 12]);
depths_norm = exprnd(1.8e-3, [num_scatterers, 1]);
depths_pool = surface_depth + 1.0e-3 + depths_norm;
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];
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
raw_spectrum = generate_spectrum(lambda, k_vec, lambda_center, sigma_lambda, ...
depths_all, reflections_all, read_noise_std, shot_noise_factor, source_instability);
raw_spectrum = raw_spectrum(:);
if isempty(example_spectrum_raw)
example_spectrum_raw = raw_spectrum;
end
dc_removed = remove_dc(raw_spectrum);
[k_uniform, spec_resampled] = resample_to_uniform_k(k_vec, dc_removed);
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;
valid_positive = (z_mm >= 0);
z_pos = z_mm(valid_positive);
a_pos = ascan_envelope(valid_positive);
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
new_error = struct(...
'frame_idx', idx, ...
'message', ME.message, ...
'identifier', ME.identifier, ...
'line', get(ME.stack, 'line'), ...
'file', fileparts(get(ME.stack, 'file')));
error_log(end+1) = new_error;
% 回退填充
bscan_intensity(:, idx) = eps;
all_detected_points{idx} = [];
end
end
% 结束后打印摘要
if ~isempty(error_log)
fprintf('⚠️ 共有 %d 条 A-scan 处理失败：\n');
for i = 1:min(5, length(error_log)) % 只显示前5个
e = error_log(i);
fprintf(' 帧 %d: %s [第%d行]\n', e.frame_idx, e.message, e.line);
end
if length(error_log) > 5
fprintf(' （还有 %d 个错误未列出）\n', length(error_log)-5);
end
else
fprintf('✅ 所有 %d 条 A-scan 处理成功！\n', N_AScans);
end
%% ==================== 后处理：构建 pts_clean 并提取初始熔深曲线 ====================
valid_mask = ~cellfun(@isempty, all_detected_points);
if any(valid_mask)
pts_cat = vertcat(all_detected_points{valid_mask});
x_repeated = arrayfun(@(i) repmat(x_mm(i), numel(all_detected_points{i}), 1), ...
find(valid_mask), 'UniformOutput', false);
x_pts = vertcat(x_repeated{:});
pts_full = [x_pts, pts_cat];
pts_clean = pts_full(pts_full(:,2) >= 0 & pts_full(:,2) <= max(z_mm)*0.95, :);
else
pts_clean = [];
fprintf('⚠️ 无有效检测点。\n');
end
% 提取初始熔深曲线（基于最深峰）
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_interp = interp1(x_mm(valid_md), melt_depth_curve(valid_md), x_mm, 'linear');
else
melt_depth_curve_interp = fillmissing(melt_depth_curve, 'nearest');
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包络');
%% ==================== 图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：四合一熔深提取与滤波对比 ====================
% 使用 melt_depth_curve_interp 作为原始信号
y_full = melt_depth_curve_interp;
x_full = x_mm;
% ✅ 自适应截取中间段（最多5000）
analysis_length = min(5000, length(y_full));
if analysis_length == length(y_full)
x_seg = x_full;
y_seg = y_full;
else
center = floor(length(y_full) / 2);
half = floor(analysis_length / 2);
startIdx = max(1, center - half + 1);
endIdx = min(length(y_full), center + half);
x_seg = x_full(startIdx:endIdx);
y_seg = y_full(startIdx:endIdx);
end
fprintf('📊 开始滤波分析：使用 %d 个数据点\n', length(y_seg));
% ✅ 【关键】不再依赖 GUI 弹窗，改为自动设置滤波参数
pct_val = 95; % 常用 90% 百分位保留高值
win_len_base = floor(length(y_seg) / 20); % 窗口约为总长度的 5%
win_len = max(11, win_len_base - mod(win_len_base-1, 2)); % 取奇数，最小11
fprintf('🔧 使用滑动百分位滤波：窗口=%d, 百分位=%.1f%%\n', win_len, pct_val);
% --- 方法1：仅百分位滤波 ---
filtered_original = slidingPercentileFilter(y_seg, win_len, pct_val);
% --- 方法2：LOF + 插值 + 百分位滤波（增强版）---
try
% ✅ 预处理：轻度平滑
y_seg_preproc = smoothdata(y_seg, 'movmedian', 5);
y_seg_preproc = smoothdata(y_seg_preproc, 'gaussian', 11);
% ✅ 执行 LOF
[lofOutlierMask_raw, lofScores, ~, ~] = performLOFOnWDP(...
y_seg_preproc, ...
'PreExtracted', true, ...
'TimeIndex', x_seg, ...
'KFactor', 0.03, ... % k ≈ 30
'ShowPlot', false);
% ✅ 类型与长度校正
lofOutlierMask_raw = logical(lofOutlierMask_raw(:));
if length(lofOutlierMask_raw) < length(y_seg)
pad = false(length(y_seg) - length(lofOutlierMask_raw), 1);
lofOutlierMask_raw = [lofOutlierMask_raw; pad];
elseif length(lofOutlierMask_raw) > length(y_seg)
lofOutlierMask_raw = lofOutlierMask_raw(1:length(y_seg));
end
% ✅ 动态阈值（IQR）
Q1 = prctile(lofScores, 25);
Q3 = prctile(lofScores, 75);
IQR = Q3 - Q1;
adaptive_threshold = Q3 + 1.5 * IQR;
lofOutlierMask = lofScores > adaptive_threshold;
% ✅ 限制最大异常比例
max_outlier_ratio = 0.1;
if sum(lofOutlierMask) > max_outlier_ratio * length(lofOutlierMask)
[~, idx_desc] = sort(lofScores, 'descend');
topN = floor(max_outlier_ratio * length(lofScores));
lofOutlierMask(:) = false;
lofOutlierMask(idx_desc(1:topN)) = true;
fprintf('💡 LOF 异常点过多，已限制为前 %d 个\n', topN);
end
% ✅ 插值修复
y_lof_clean = double(y_seg);
y_lof_clean(lofOutlierMask) = NaN;
idx_all = 1:length(y_seg);
validIdx = ~isnan(y_lof_clean);
if sum(validIdx) == 0
warning('🟡 LOF：无有效点，使用中位数填充');
y_lof_interp = median(y_seg, 'omitnan') * ones(size(y_seg));
elseif sum(validIdx) == 1
singleIdx = find(validIdx, 1);
y_lof_interp = interp1(singleIdx, y_lof_clean(singleIdx), idx_all, 'nearest');
else
y_lof_interp = interp1(idx_all(validIdx), y_lof_clean(validIdx), idx_all, 'pchip');
nanIdx = isnan(y_lof_interp);
if any(nanIdx) && sum(~nanIdx) >= 1
y_lof_interp(nanIdx) = interp1(idx_all(~nanIdx), y_lof_interp(~nanIdx), idx_all(nanIdx), 'nearest');
end
end
% ✅ 应用百分位滤波
filtered_lof = slidingPercentileFilter(y_lof_interp, win_len, pct_val);
fprintf('✅ LOF 完成：发现 %d 个受控异常点（自适应阈值=%.3f）\n', sum(lofOutlierMask), adaptive_threshold);
catch ME
fprintf('❌ LOF 失败: %s\n', ME.message);
warning('🔴 回退到原始信号');
filtered_lof = filtered_original;
y_lof_interp = y_seg;
end
% --- 方法3：DBSCAN + 插值 + 百分位滤波 ---
try
[dbscanOutlierMask_temp, ~, ~, ~] = performDBSCANOnWDP(...
y_seg, ...
'PreExtracted', true, ...
'TimeIndex', x_seg, ...
'MinPts', 5, ...
'Eps', 0.4, ...
'ShowPlot', false);
dbscanOutlierMask = logical(dbscanOutlierMask_temp);
y_dbscan_clean = y_seg;
y_dbscan_clean(dbscanOutlierMask) = NaN;
idx_all = 1:length(y_seg);
validIdx = ~isnan(y_dbscan_clean);
if sum(validIdx) < 2
warning('DBSCAN后有效点太少，跳过插值');
y_dbscan_interp = y_seg;
else
y_dbscan_interp = interp1(idx_all(validIdx), y_dbscan_clean(validIdx), idx_all, 'pchip');
nanIdx = isnan(y_dbscan_interp);
if any(nanIdx) && sum(~nanIdx) >= 2
y_dbscan_interp(nanIdx) = interp1(idx_all(~nanIdx), y_dbscan_interp(~nanIdx), idx_all(nanIdx), 'nearest');
end
end
filtered_dbscan = slidingPercentileFilter(y_dbscan_interp, win_len, pct_val);
catch ME
warning('🟡 DBSCAN 处理失败: %s', ME.message);
filtered_dbscan = filtered_original;
dbscanOutlierMask = false(size(y_seg));
y_dbscan_interp = y_seg;
end
% --- 绘制四合一结果 ---
figure('Position', [100, 100, 1000, 800], 'Name', '【四合一】熔深提取与滤波对比', 'Color', 'white');
set(gcf, 'DefaultAxesFontName', 'Microsoft YaHei');
set(gcf, 'GraphicsSmoothing', 'on');
% Subplot 1: 空间密度趋势法原始熔深（已包含线性插值）
subplot(4,1,1);
plot(x_mm, melt_depth_curve_interp, 'k-', 'LineWidth', 1.8);
xlabel('x (mm)'); ylabel('熔深 (mm)');
title('1. 空间密度趋势法提取的熔深曲线（线性插值补全）');
grid on; box on;
% Subplot 2: 百分位滤波
subplot(4,1,2);
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);
title(sprintf('2. %g%% 百分位滤波（窗口=%d）', pct_val, win_len));
xlabel('时间 / 位置 (mm)'); ylabel('熔深 (mm)');
legend('原始点', sprintf('%g%%滤波', pct_val), 'Location', 'bestoutside');
grid on; box on; hold off;
% Subplot 3: LOF + 百分位
subplot(4,1,3);
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);
scatter(x_seg(lofOutlierMask), y_seg(lofOutlierMask), 30, 'r', 'filled', 'MarkerFaceAlpha', 0.7);
title('3. LOF去噪 → 百分位滤波');
xlabel('时间 / 位置 (mm)'); ylabel('熔深 (mm)');
legend('插值后信号', 'LOF+滤波', 'LOF异常点', 'Location', 'bestoutside');
grid on; box on; hold off;
% Subplot 4: DBSCAN + 百分位
subplot(4,1,4);
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);
scatter(x_seg(dbscanOutlierMask), y_seg(dbscanOutlierMask), 30, 'm', 'diamond', 'filled', 'MarkerFaceAlpha', 0.7);
title('4. DBSCAN去噪 → 百分位滤波');
xlabel('时间 / 位置 (mm)'); ylabel('熔深 (mm)');
legend('插值后信号', 'DBSCAN+滤波', 'DBSCAN异常点', 'Location', 'bestoutside');
grid on; box on; hold off;
sgtitle({'📊 四种熔深提取与滤波策略对比'; ...
'蓝色=LOF, 绿色=DBSCAN, 黑色=百分位, 灰色=原始'}, ...
'FontSize', 12, 'FontWeight', 'bold');
drawnow;
disp('🎉 全流程完成！已生成三张图：');
disp(' ➤ 图1：光谱处理五步链');
disp(' ➤ 图2：A-scan与密度趋势');
disp(' ➤ 图3：四合一熔深对比图');
end % 主函数结束
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
function y = remove_dc(x)
background = mean(x(1:2048));
y = x - background;
end
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
function y = apply_window(x, type_str)
w = window(@hann, length(x));
y = x .* w(:);
end
function envelope = compute_ascan(k_uniform, spectrum)
spec_shifted = ifftshift(spectrum(:));
ascan_complex = fft(spec_shifted);
envelope = abs(ascan_complex);
end
function result = slidingPercentileFilter(signal, windowLen, pct)
n = length(signal);
result = zeros(n, 1);
halfWin = floor(windowLen / 2);
% 边界保持原样
leftPad = min(halfWin, n-1);
rightPad = min(halfWin, n-1);
result(1:leftPad) = signal(1:leftPad);
result(end-rightPad+1:end) = signal(end-rightPad+1:end);
% 中心区域
startCenter = halfWin + 1;
endCenter = n - halfWin;
if startCenter <= endCenter
idx_center = startCenter : endCenter;
windows = arrayfun(@(i) signal(i-halfWin:i+halfWin), idx_center, 'UniformOutput', false);
result(startCenter:endCenter) = cellfun(@(w) prctile(w, pct), windows);
end
end
function [filteredY, x, finalParams] = performPercentileFiltering(rawData, varargin)
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 size(rawData, 2) >= 2
x = rawData(:, 1);
y = rawData(:, 2);
else
x = (1:length(rawData))';
y = rawData(:);
end
n = length(y);
if isempty(percentileValue) || isempty(windowLength)
[p_val, w_len] = getFilterParamsFromDialog(percentileValue, windowLength);
if isempty(p_val) || isempty(w_len)
error('用户取消');
end
percentileValue = p_val;
windowLength = w_len;
end
if mod(windowLength, 2) == 0
warning('⚠️ 窗口长度应为奇数，已自动调整');
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
finalParams.Percentile = percentileValue;
finalParams.WindowLength = windowLength;
fprintf('🔧 已完成 %.1f%% 百分位滤波处理，窗口大小=%d\n', percentileValue, windowLength);
end
function [lofScores, isOutlier, x_sub, y_sub] = performLOFOnWDP(rawData, varargin)
p = inputParser;
addOptional(p, 'TargetLength', 5000, @(x) isnumeric(x) && isscalar(x) && x > 0);
addOptional(p, 'UseMedianIQR', true, @islogical);
addOptional(p, 'KFactor', 0.02, @(x) isnumeric(x) && x > 0 && x <= 0.1);
addOptional(p, 'ShowPlot', false, @islogical);
addOptional(p, 'PreExtracted', false, @islogical);
addOptional(p, 'TimeIndex', [], @(x) isempty(x) || (isvector(x) && length(x) == numel(x)));
parse(p, varargin{:});
preExtracted = p.Results.PreExtracted;
timeIndex = p.Results.TimeIndex;
if ~preExtracted
error('只支持 PreExtracted=true 模式');
end
y_sub = rawData(:);
if isempty(timeIndex)
x_sub = (1:length(y_sub))';
else
x_sub = timeIndex(:);
end
dy = diff([y_sub; y_sub(end)]);
ddy = diff([dy; dy(end)]);
features = [y_sub, dy, ddy];
mu = mean(features, 1); sigma = std(features, 0, 1); sigma(sigma==0)=1;
data_norm = (features - mu) ./ sigma;
k = min(20, max(3, floor(p.Results.KFactor * length(y_sub))));
D = pdist2(data_norm, data_norm);
[~, sortedIdx] = sort(D, 2);
kNN_idx = sortedIdx(:, 2:k+1);
r_k_dist = D(sub2ind(size(D), repmat((1:size(D,1))', 1, k), kNN_idx));
reachDistMat = zeros(size(kNN_idx));
for i = 1:size(reachDistMat,1)
for j = 1:k
reachDistMat(i,j) = max(D(i,kNN_idx(i,j)), r_k_dist(kNN_idx(i,j),j));
end
end
avgReachDist = mean(reachDistMat, 2);
LRD = 1 ./ (avgReachDist + eps);
lofScores = zeros(size(LRD));
for i = 1:length(LRD)
lofScores(i) = mean(LRD(kNN_idx(i,:))) / LRD(i);
end
Q1 = prctile(lofScores, 25); Q3 = prctile(lofScores, 75); IQR = Q3 - Q1;
threshold = median(lofScores) + 3 * IQR;
isOutlier = lofScores > threshold;
fprintf('✅ LOF 分析完成：共发现 %d 个异常点（阈值=%.3f）\n', sum(isOutlier), threshold);
end
function [isOutlier, x_sub, y_sub, eps_used] = performDBSCANOnWDP(rawData, varargin)
p = inputParser;
addOptional(p, 'MinPts', 5, @(x) isnumeric(x) && isscalar(x) && x >= 1);
addOptional(p, 'Eps', [], @(x) isempty(x) || (isnumeric(x) && x > 0));
addOptional(p, 'ShowPlot', false, @islogical);
addOptional(p, 'PreExtracted', true, @islogical);
addOptional(p, 'TimeIndex', [], @(x) isempty(x) || (isvector(x) && length(x)==numel(rawData)));
parse(p, varargin{:});
y_sub = rawData(:);
if isempty(p.Results.TimeIndex)
x_sub = (1:length(y_sub))';
else
x_sub = p.Results.TimeIndex(:);
end
dy = diff([y_sub; y_sub(end)]); ddy = diff([dy; dy(end)]);
features = [y_sub, dy, ddy];
mu = mean(features,1); sigma = std(features,0,1); sigma(sigma==0)=1;
data_norm = (features - mu)./sigma;
if isempty(p.Results.Eps)
D = pdist2(data_norm, data_norm);
[~, sortedD] = sort(D,2); kDist = sortedD(:, p.Results.MinPts+1);
kDistSorted = sort(kDist, 'descend');
diff2 = diff(diff(kDistSorted));
[~, idx] = max(abs(diff2)); eps_est = kDistSorted(max(2,idx));
eps_used = round(eps_est*1000)/1000;
else
eps_used = p.Results.Eps;
end
[labels,~] = dbscan(data_norm, eps_used, p.Results.MinPts);
isOutlier = (labels == -1);
fprintf('✅ DBSCAN 完成：共发现 %d 个异常点（eps=%.4f）\n', sum(isOutlier), eps_used);
end
1：对代码作出正确的理解；
2：对代码目前存在的问题/冗余进行详细的说明，并进行改进，改进后呈现完整的代码给我；
3：我的matlab版本是2024b；
最新发布
本文详细介绍了如何在Apache服务器中安装并配置mod_bandwidth模块，以实现对不同用户和域名的下载速度限制。通过设置BandWidth、LargeFileLimit等参数，管理员可以在不影响网站性能的情况下，灵活地控制网络流量。
扫一扫
2357
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
