function [all_sequences, all_pri, remaining_toa] = optimized_CDIF(filename, max_level, k_threshold, 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;
% 主循环:处理所有序列
while length(remaining_toa) >= min_total_pulses
[sequence, pri, level, remaining_toa] = extract_next_sequence(remaining_toa, PRI_max, k_threshold, ...
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);
else
break;
end
end
% 显示最终结果
fprintf('\n===== 分析完成 =====\n');
fprintf('提取序列总数: %d\n', length(all_sequences));
fprintf('剩余脉冲数: %d\n', length(remaining_toa));
% 可视化最终结果
plot_final_results(all_sequences, remaining_toa);
end
function [sequence, pri, level, remaining_toa] = extract_next_sequence(toa, PRI_max, k_threshold, tolerance, …
min_pulse_sequence, min_total_pulses, max_level)
sequence = [];
pri = [];
level = [];
current_level = 1;
% CDIF分析循环(不移除图形绘制)
while current_level <= max_level
cumulative_hist = zeros(1, ceil(PRI_max));
current_T_total = max(toa) - min(toa);
% 构建累积直方图
for lvl = 1:current_level
if length(toa) > lvl
dtoa = toa(lvl+1:end) - toa(1:end-lvl);
valid_indices = (dtoa >= 1) & (dtoa <= PRI_max);
dtoa = dtoa(valid_indices);
for i = 1:length(dtoa)
idx = min(ceil(dtoa(i)), length(cumulative_hist));
cumulative_hist(idx) = cumulative_hist(idx) + 1;
end
end
end
% 计算门限并寻找候选PRI
thresholds = k_threshold * current_T_total ./ (1:length(cumulative_hist));
candidate_pri = [];
for pri_val = 1:length(cumulative_hist)
if cumulative_hist(pri_val) > thresholds(pri_val)
% 跳过过小值(降噪)
if pri_val < 10, continue; end
% 检查二倍值
double_pri = 2 * pri_val;
if double_pri <= length(cumulative_hist) && cumulative_hist(double_pri) > thresholds(double_pri)
candidate_pri(end+1) = pri_val;
else
candidate_pri(end+1) = pri_val;
end
end
end
% 尝试提取序列
if ~isempty(candidate_pri)
candidate_pri = sort(candidate_pri);
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
current_level = current_level + 1;
end
end
function plot_final_results(all_sequences, remaining_toa)
% 创建新图形
figure(‘Position’, [100, 100, 1200, 800], ‘Name’, ‘CDIF脉冲分选结果’);
hold on;
grid on;
title(‘CDIF脉冲分选结果’);
xlabel(‘到达时间 (us)’);
ylabel(‘脉冲序列’);
% 确定颜色方案
colors = lines(length(all_sequences));
% 确定Y轴位置范围
max_y = max(1, length(all_sequences)); % 确保至少为1
y_min = -0.5;
y_max = max_y + 0.5;
% 为每个序列创建不同的Y轴位置
for i = 1:length(all_sequences)
% 计算当前序列的平均PRI
sequence_toa = sort(all_sequences{i});
sequence_diffs = diff(sequence_toa);
avg_pri = mean(sequence_diffs(sequence_diffs > 0));
% 为序列中的每个脉冲分配相同的Y值 (i)
y_value = i * ones(size(all_sequences{i}));
% 绘制序列脉冲
scatter(all_sequences{i}, y_value, 80, colors(i,:), 'filled', 'o');
% 标记序列号和PRI值
min_time = min(all_sequences{i});
text(min_time - max(all_sequences{i})*0.02, i, ...
sprintf('序列%d (PRI=%.1fus)', i, avg_pri), ...
'VerticalAlignment', 'middle', 'FontWeight', 'bold', ...
'BackgroundColor', 'white', 'Margin', 1);
% 添加序列趋势线
plot(sequence_toa, i * ones(size(sequence_toa)), ...
'Color', colors(i,:), 'LineWidth', 1.5);
end
% 绘制剩余脉冲(位于Y=0位置)
if ~isempty(remaining_toa)
y_remaining = zeros(size(remaining_toa));
scatter(remaining_toa, y_remaining, 50, [0.5 0.5 0.5], 'x', 'LineWidth', 1.5);
text(min(remaining_toa) - max(all_sequences{end})*0.02, 0, '未分类脉冲', ...
'Color', [0.5 0.5 0.5], 'FontWeight', 'bold', ...
'VerticalAlignment', 'middle', 'BackgroundColor', 'white');
end
% 设置Y轴范围和标签
yticks(0:max_y);
ytick_labels = arrayfun(@(i) {sprintf('序列%d', i)}, 1:max_y);
if ~isempty(remaining_toa)
ytick_labels{1} = '未分类脉冲';
end
yticklabels(ytick_labels);
% 添加网格和美观设置
set(gca, 'YGrid', 'on', 'XGrid', 'on', 'GridLineStyle', '--');
% 添加图例
legends = arrayfun(@(i) sprintf('序列%d', i), 1:length(all_sequences), 'UniformOutput', false);
if ~isempty(remaining_toa)
legends = [{'未分类脉冲'}, legends]; % 确保未分类在首位
end
legend(legends, 'Location', 'eastoutside', 'FontSize', 10);
% 调整视图范围
min_x = min(cellfun(@min, all_sequences));
max_x = max(cellfun(@max, all_sequences));
if ~isempty(remaining_toa)
min_x = min(min_x, min(remaining_toa));
max_x = max(max_x, max(remaining_toa));
end
x_range = max_x - min_x;
xlim([min_x - x_range * 0.05, max_x + x_range * 0.05]);
% 自动调整Y轴范围
ylim([y_min y_max]);
hold off;
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(10, length(toa))
current_sequence = toa(start_idx);
current_time = toa(start_idx);
last_index = start_idx;
lost_count = 0;
% 向前搜索
for j = (start_idx+1):length(toa)
expected_time = current_time + pri;
% 检查当前脉冲是否在期望时间附近
if toa(j) >= expected_time - tolerance_val && ...
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
% 检查是否在下一个期望时间附近
next_expected = expected_time + pri;
if toa(j) >= next_expected - tolerance_val && ...
toa(j) <= next_expected + tolerance_val
% 接受丢失一个脉冲
current_sequence(end+1) = toa(j);
current_time = toa(j);
last_index = j;
lost_count = 1;
% 检查是否在更远的期望时间附近
elseif toa(j) > next_expected + tolerance_val
% 检查是否在第二个丢失脉冲的位置
next_next_expected = next_expected + pri;
if toa(j) >= next_next_expected - tolerance_val && ...
toa(j) <= next_next_expected + tolerance_val
% 接受丢失两个脉冲
current_sequence(end+1) = toa(j);
current_time = toa(j);
last_index = j;
lost_count = 2;
else
% 如果脉冲差距太大,可能不是同一个序列
if lost_count >= 2
break;
end
end
end
end
end
% 检查序列长度
if length(current_sequence) > best_length
best_sequence = current_sequence;
best_length = length(current_sequence);
end
end
% 验证序列长度
if best_length >= min_length
sequence = best_sequence;
% 从原始TOA中移除序列
[~, idx] = ismember(sequence, toa);
remaining = toa;
remaining(idx) = [];
end
end
% 设置参数
max_level = 2;
k_threshold = 0.2;
tolerance = 0.05;
min_pulse_sequence = 3;
min_total_pulses = 5;
% 运行优化的CDIF算法
[sequences, pris, remaining] = optimized_CDIF(‘PDW1.txt’, max_level, k_threshold, tolerance, min_pulse_sequence, min_total_pulses);
将上述代码中绘图函数,改为绘制结果的直方图以及门限函数,横轴为PRI纵轴为统计数