山东省第五届省赛J题 Weighted Median(思维题 sort给结构体排序)

本文介绍了一种在O(n)最坏情况下计算加权中位数的方法。通过定义元素及其权重,使用结构体存储并根据数值大小进行排序,进而找到满足条件的加权中位数。

Description
For n elements x1, x2, …, xn with positive integer weights w1, w2, …, wn. The weighted median is the element xk satisfying
这里写图片描述 and 这里写图片描述 , S indicates 这里写图片描述
Can you compute the weighted median in O(n) worst-case?

Input
There are several test cases. For each case, the first line contains one integer n(1 ≤  n ≤ 10^7) — the number of elements in the sequence. The following line contains n integer numbers xi (0 ≤ xi ≤ 10^9). The last line contains n integer numbers wi (0 < wi < 10^9).

Output
One line for each case, print a single integer number— the weighted median of the sequence.

Sample Input
7
10 35 5 10 15 5 20
10 35 5 10 15 5 20
Sample Output
20

**【题意】:**给你n个数然后第一行是这些数xi,第二行是他们分别对应的权值wi。s是权值总和。然后让你找一个数xk满足公式,(x1xk)求和>s/2,(xkxn)求和小于等于s/2;可看题目中的公式因为两边求和都不存在xk的权值,则只要找到一半小于s/2则另一半一定小于等于s/2。

**【思路】:**将其存在一个结构体中,然后按x的大小进行排序(代码按降序排序),然后开一个变量ans=0遍历排完序的序列,当ans>=s/2时,记录下标k跳出循环即可,最后输出k所对应的x即可。

【注】:这里的数据量很大 所以结构体数组要开到1e7+10
并且 下面求和的两个变量都需要开long long 或double 不然会RE

【代码】:

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
using namespace std;
const int N=1e7+10;
struct F
{
  int x;
  int w;
}a[N];

bool cmp(struct F aa,struct F bb)
{
   return aa.x>bb.x;
}

int main()
{
    ios::sync_with_stdio(false);
    int n;
    while(cin>>n)
    {
       long long s=0;
       for(int i=1;i<=n;i++)
       {
           cin>>a[i].x;
       }
       for(int i=1;i<=n;i++)
       {
           cin>>a[i].w;
           s+=a[i].w;
       }
       s/=2;
       sort(a+1,a+n+1,cmp);
       long long ans=0;
       int k=0;
       for(int i=1;i<=n;i++)
       {
           ans+=a[i].w;
           if(ans>=s)//寻找xk
           {
               k=i;
               break;
           }
       }
       cout<<a[k].x<<endl;
    }
    return 0;
}
%% ==================== 完整融合版: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;
最新发布
10-24
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值