function [all_sequences, all_pri, remaining_toa] = optimized_CDIF(filename, max_level, a, tolerance, min_pulse_sequence, min_total_pulses)
% 读取数据并初始化
data = load(filename);
toa_original = data(:,4);
toa_original = sort(toa_original);
remaining_toa = toa_original;
% 计算初始参数
T_total = max(toa_original) - min(toa_original);
absolute_max_pri = 1000000; % 1000ms
PRI_max = min(T_total/2, absolute_max_pri);
% 存储提取的序列
all_sequences = {};
all_pri = [];
sequence_counter = 1;
% 存储所有序列的CDIF数据（按层级分组）
all_seq_cdif_data = cell(1, max_level); % 每个层级的序列CDIF数据
% 主循环：处理所有序列
while length(remaining_toa) >= min_total_pulses
[sequence, pri, level, remaining_toa, level_hists] = extract_next_sequence(...
remaining_toa, PRI_max, a, tolerance, ...
min_pulse_sequence, min_total_pulses, max_level);
if ~isempty(sequence)
all_sequences{end+1} = sequence;
all_pri(end+1) = pri;
sequence_counter = sequence_counter + 1;
fprintf('提取序列%d: PRI=%.1fus, 脉冲数=%d (分析层级:%d)\n', ...
sequence_counter-1, pri, length(sequence), level);
for lvl = 1:length(level_hists)
level_data = level_hists{lvl};
if lvl <= max_level
if isempty(all_seq_cdif_data{lvl})
all_seq_cdif_data{lvl} = struct(...
'seq_id', sequence_counter-1, ...
'hist', level_data.hist, ...
'threshold', level_data.threshold, ...
'pri', pri);
else
all_seq_cdif_data{lvl}(end+1) = struct(...
'seq_id', sequence_counter-1, ...
'hist', level_data.hist, ...
'threshold', level_data.threshold, ...
'pri', pri);
end
end
end
else
break;
end
end
% 绘制组合CDIF图
for lvl = 1:max_level
if ~isempty(all_seq_cdif_data{lvl})
plot_combined_cdif(all_seq_cdif_data{lvl}, lvl);
end
end
% 显示最终结果
fprintf('\n===== 分析完成 =====\n');
fprintf('提取序列总数: %d\n', length(all_sequences));
fprintf('剩余脉冲数: %d\n', length(remaining_toa));
end
function [sequence, pri, level, remaining_toa, level_hists] = extract_next_sequence(...
toa, PRI_max, a, tolerance, min_pulse_sequence, min_total_pulses, max_level)
level_hists = {};
cumulative_hist = []; % 累积直方图
T_total = max(toa) - min(toa);
% CDIF分析循环 (从1级到max_level)
for current_level = 1:max_level
% 计算当前级别的直方图
[level_hist, thresholds] = compute_level_hist(toa, current_level, PRI_max, T_total, a);
% 累积到前一级的直方图
if isempty(cumulative_hist)
cumulative_hist = level_hist;
else
% 确保直方图长度一致
max_len = max(length(cumulative_hist), length(level_hist));
cumulative_hist(end+1:max_len) = 0;
level_hist(end+1:max_len) = 0;
cumulative_hist = cumulative_hist + level_hist;
end
% 保存当前层级直方图
level_data.hist = cumulative_hist;
level_data.threshold = thresholds;
level_hists{end+1} = level_data;
% 寻找候选PRI（高于最优门限）
candidate_indices = find(cumulative_hist > thresholds);
% 过滤无效PRI候选值
candidate_indices(candidate_indices < 10) = []; % 过滤过小的PRI
candidate_indices(candidate_indices > PRI_max) = []; % 过滤过大的PRI
if ~isempty(candidate_indices)
% 按显著性排序（直方图值降序）
[sorted_values, sort_idx] = sort(cumulative_hist(candidate_indices), 'descend');
candidate_pri = candidate_indices(sort_idx);
% 尝试提取序列
for i = 1:length(candidate_pri)
tolerance_val = candidate_pri(i) * tolerance;
[sequence, new_toa] = extract_sequence(toa, candidate_pri(i), tolerance_val, min_pulse_sequence);
if ~isempty(sequence)
pri = candidate_pri(i);
level = current_level;
remaining_toa = new_toa;
return;
end
end
end
end
% 未找到序列
pri = [];
level = [];
sequence = [];
remaining_toa = toa;
level_hists = {};
end
function [level_hist, thresholds] = compute_level_hist(toa, level, PRI_max, T_total, a)
% 初始化直方图
level_hist = zeros(1, ceil(PRI_max));
if length(toa) > level
% 计算当前级别的dTOA
dtoa = toa(level+1:end) - toa(1:end-level);
% 过滤有效值
valid_indices = (dtoa >= 1) & (dtoa <= PRI_max);
dtoa = dtoa(valid_indices);
% 统计直方图
for i = 1:length(dtoa)
idx = min(ceil(dtoa(i)), length(level_hist));
if idx > 0
level_hist(idx) = level_hist(idx) + 1;
end
end
end
% 计算最优门限函数 D_{cdif} = a * T / t
t = 1:length(level_hist);
thresholds = a * T_total ./ t;
% 处理除零和无效值
thresholds(isinf(thresholds) | isnan(thresholds)) = max(level_hist);
thresholds(t == 0) = max(level_hist);
% 确保门限不低于最小值
min_threshold = 3; % 最小门限值
thresholds(thresholds < min_threshold) = min_threshold;
end
function plot_combined_cdif(cdif_data, level)
fig_name = sprintf('层级%d - 所有序列CDIF直方图', level);
figure('Name', fig_name, 'NumberTitle', 'off');
% 获取最大PRI范围
max_pri_len = 0;
for i = 1:length(cdif_data)
max_pri_len = max(max_pri_len, length(cdif_data(i).hist));
end
% 创建颜色映射
colors = lines(length(cdif_data));
hold on;
% 绘制所有序列的直方图
for i = 1:length(cdif_data)
data = cdif_data(i);
x = 1:length(data.hist);
% 绘制直方图
plot(x, data.hist, 'Color', colors(i, :), 'LineWidth', 1.5, ...
'DisplayName', sprintf('序列%d (PRI=%.1fµs)', data.seq_id, data.pri));
% 找到最接近实际PRI的索引位置
[~, idx] = min(abs(x - data.pri));
peak_value = data.hist(idx);
% 标记检测到的PRI - 使用实际峰值位置
plot(x(idx), peak_value, 'o', 'MarkerSize', 8, ...
'MarkerFaceColor', colors(i, :), 'MarkerEdgeColor', 'k', ...
'HandleVisibility', 'off');
end
% 绘制统一的门限曲线
if ~isempty(cdif_data)
x_thresh = 1:length(cdif_data(1).threshold);
plot(x_thresh, cdif_data(1).threshold, 'k--', 'LineWidth', 1.5, ...
'DisplayName', '最优门限');
end
hold off;
% 添加图例和标签
title(fig_name);
xlabel('PRI (\mu s)');
ylabel('累积计数');
legend('show', 'Location', 'best');
grid on;
% 设置坐标轴范围
xlim([1, min(max_pri_len, 100000)]);
% 添加最优门限说明
annotation('textbox', [0.15, 0.85, 0.2, 0.05], 'String', ...
sprintf('最优门限: D_{cdif} = \\alpha T_{total}/t'), ...
'FitBoxToText', 'on', 'BackgroundColor', 'white', 'EdgeColor', 'none');
end
function [sequence, remaining] = extract_sequence(toa, pri, tolerance_val, min_length)
sequence = [];
remaining = toa;
if length(toa) < min_length
return;
end
best_sequence = [];
best_length = 0;
% 尝试不同的起点
for start_idx = 1:min(20, length(toa))
current_sequence = toa(start_idx);
current_time = toa(start_idx);
last_index = start_idx;
lost_count = 0;
total_missed = 0;
% 向前搜索
for j = (start_idx+1):length(toa)
expected_time = current_time + pri;
% 检查当前脉冲是否在期望时间附近
if abs(toa(j) - expected_time) <= tolerance_val
current_sequence(end+1) = toa(j);
current_time = toa(j);
last_index = j;
lost_count = 0; % 重置丢失计数
% 检查是否可能丢失脉冲
elseif toa(j) > expected_time + tolerance_val
% 计算可能丢失的脉冲数
missed_count = floor((toa(j) - expected_time) / pri);
% 检查是否在容差范围内
next_expected = expected_time + missed_count * pri;
if abs(toa(j) - next_expected) <= tolerance_val
current_sequence(end+1) = toa(j);
current_time = toa(j);
last_index = j;
lost_count = lost_count + missed_count;
total_missed = total_missed + missed_count;
else
% 如果脉冲差距太大，可能不是同一个序列
if lost_count >= 2
break;
end
end
end
end
% 检查序列质量
if length(current_sequence) >= min_length
% 计算序列持续时间
duration = current_sequence(end) - current_sequence(1);
% 检查序列质量标准：
% 1. 序列持续时间至少为2倍PRI
% 2. 丢失脉冲比例不超过50%
if duration >= 2 * pri && (total_missed / length(current_sequence)) <= 0.5
if length(current_sequence) > best_length
best_sequence = current_sequence;
best_length = length(current_sequence);
end
end
end
end
% 验证序列长度
if best_length >= min_length
sequence = best_sequence;
% 从原始TOA中移除序列
[~, idx] = ismember(sequence, toa);
remaining = toa;
remaining(idx) = [];
end
end
% 设置参数
max_level = 3; % 可根据需要调整
k_threshold = 0.1;
tolerance = 0.02;
min_pulse_sequence = 5; % 提高最小序列长度要求
min_total_pulses = 5;
% 运行优化的CDIF算法
[sequences, pris, remaining] = optimized_CDIF('PDW2.txt', max_level, k_threshold, tolerance, min_pulse_sequence, min_total_pulses);
上述代码输出的序列中还是会存在相同pri被分到不同序列的情况
扫一扫
8166
到【灌水乐园】发言
