上传文件至「Graduation Design Matlab」
This commit is contained in:
@@ -0,0 +1,121 @@
|
||||
%% 全空域四扇区动态先验鲁棒性测试 (多核并行加速版)
|
||||
clear; clc; close all;
|
||||
rng(2024);
|
||||
|
||||
par = setpar();
|
||||
N_ant = par.M; N_RF = par.Nrf;
|
||||
clutter_angle = -15; a_clutter = gen_a(clutter_angle, par);
|
||||
INR_linear = 10^(30/10); noise_power = 1;
|
||||
R_interference = noise_power * eye(N_ant) + INR_linear * (a_clutter * a_clutter');
|
||||
|
||||
disp('正在构建全空域字典 (传统OMP使用)...');
|
||||
dic_angles_full = -90:0.5:90; A_dic_full = zeros(N_ant, length(dic_angles_full));
|
||||
for i = 1:length(dic_angles_full), A_dic_full(:,i) = gen_a(dic_angles_full(i), par) / sqrt(N_ant); end
|
||||
|
||||
angle_ranges = {0:1:10, 11:1:20, 21:1:30, 31:1:40};
|
||||
snr_test_list = [0, 5, 10, 15, 20];
|
||||
Monte_Carlo = 500;
|
||||
L_snapshots = 100;
|
||||
color_array = {'#0072BD', '#D95319', '#EDB120', '#7E2F8E', '#77AC30'};
|
||||
|
||||
disp('正在启动多核并行引擎 (parfor) 进行全空域极速测试...');
|
||||
for fig_idx = 1:4
|
||||
test_angles = angle_ranges{fig_idx};
|
||||
fprintf('\n>>> 正在进行第 %d 张图的并行仿真:扇区 %d° ~ %d° ...\n', fig_idx, min(test_angles), max(test_angles));
|
||||
|
||||
RMSE_DFT_All = zeros(length(snr_test_list), length(test_angles));
|
||||
RMSE_OMP_FULL_All = zeros(length(snr_test_list), length(test_angles));
|
||||
RMSE_OMP_All = zeros(length(snr_test_list), length(test_angles));
|
||||
|
||||
figure('Name', sprintf('扇区 %d°~%d° 测试', min(test_angles), max(test_angles)), 'Color', 'w', 'Position', [50+fig_idx*50, 50+fig_idx*50, 850, 600]);
|
||||
hold on; grid on;
|
||||
|
||||
for snr_idx = 1:length(snr_test_list)
|
||||
snr_test = snr_test_list(snr_idx);
|
||||
sig_power = 10^(snr_test/10);
|
||||
|
||||
err_dft = zeros(length(test_angles), Monte_Carlo);
|
||||
err_omp_full = zeros(length(test_angles), Monte_Carlo);
|
||||
err_omp = zeros(length(test_angles), Monte_Carlo);
|
||||
|
||||
for a_idx = 1:length(test_angles)
|
||||
true_target_angle = test_angles(a_idx) + 0.04;
|
||||
theta_prior = round(true_target_angle); window = 5;
|
||||
|
||||
dic_tsc = (theta_prior - window) : 0.1 : (theta_prior + window);
|
||||
A_dic_tsc = zeros(N_ant, length(dic_tsc));
|
||||
for i = 1:length(dic_tsc), A_dic_tsc(:,i) = gen_a(dic_tsc(i), par) / sqrt(N_ant); end
|
||||
|
||||
angles_dft = linspace(theta_prior - window, theta_prior + window, N_RF);
|
||||
F_RF_dft = zeros(N_ant, N_RF);
|
||||
for i = 1:N_RF, F_RF_dft(:,i) = gen_a(angles_dft(i), par) / sqrt(N_ant); end
|
||||
|
||||
a_prior_broad = zeros(N_ant, 1);
|
||||
for theta_b = (theta_prior - window) : 0.5 : (theta_prior + window)
|
||||
a_prior_broad = a_prior_broad + gen_a(theta_b, par);
|
||||
end
|
||||
a_prior_broad = a_prior_broad / norm(a_prior_broad);
|
||||
w_dig_prior = inv(R_interference) * a_prior_broad; w_dig_prior = w_dig_prior / norm(w_dig_prior);
|
||||
|
||||
res_full = w_dig_prior; F_RF_omp_full = [];
|
||||
for k = 1:N_RF
|
||||
proj = A_dic_full' * res_full; [~, midx] = max(abs(proj));
|
||||
F_RF_omp_full = [F_RF_omp_full, A_dic_full(:, midx)]; res_full = w_dig_prior - F_RF_omp_full * (pinv(F_RF_omp_full)*w_dig_prior);
|
||||
end
|
||||
res_tsc = w_dig_prior; F_RF_omp = [];
|
||||
for k = 1:N_RF
|
||||
proj = A_dic_tsc' * res_tsc; [~, midx] = max(abs(proj));
|
||||
F_RF_omp = [F_RF_omp, A_dic_tsc(:, midx)]; res_tsc = w_dig_prior - F_RF_omp * (pinv(F_RF_omp)*w_dig_prior);
|
||||
end
|
||||
|
||||
scan_grid = (theta_prior - window) : 0.01 : (theta_prior + window);
|
||||
|
||||
% 为 parfor 准备的切片临时变量
|
||||
err_dft_tmp = zeros(1, Monte_Carlo); err_omp_full_tmp = zeros(1, Monte_Carlo); err_omp_tmp = zeros(1, Monte_Carlo);
|
||||
|
||||
% 【核心提速】:parfor 多核并行循环!
|
||||
parfor mc = 1:Monte_Carlo
|
||||
s_t = sqrt(sig_power/2)*(randn(1,L_snapshots) + 1j*randn(1,L_snapshots));
|
||||
s_c = sqrt(INR_linear/2)*(randn(1,L_snapshots) + 1j*randn(1,L_snapshots));
|
||||
n_t = sqrt(noise_power/2)*(randn(N_ant,L_snapshots) + 1j*randn(N_ant,L_snapshots));
|
||||
X = gen_a(true_target_angle, par)*s_t + a_clutter*s_c + n_t;
|
||||
|
||||
Y_dft = F_RF_dft' * X; Y_omp_full = F_RF_omp_full' * X; Y_omp = F_RF_omp' * X;
|
||||
R_dft = (Y_dft*Y_dft')/L_snapshots; R_dft = R_dft + (0.05*trace(R_dft)/N_RF)*eye(N_RF);
|
||||
R_omp_full = (Y_omp_full*Y_omp_full')/L_snapshots; R_omp_full = R_omp_full + (0.05*trace(R_omp_full)/N_RF)*eye(N_RF);
|
||||
R_omp = (Y_omp*Y_omp')/L_snapshots; R_omp = R_omp + (0.05*trace(R_omp)/N_RF)*eye(N_RF);
|
||||
|
||||
inv_dft = inv(R_dft); inv_full = inv(R_omp_full); inv_tsc = inv(R_omp);
|
||||
P_dft = zeros(1, length(scan_grid)); P_full = zeros(1, length(scan_grid)); P_tsc = zeros(1, length(scan_grid));
|
||||
|
||||
for s_idx = 1:length(scan_grid)
|
||||
a_test = gen_a(scan_grid(s_idx), par);
|
||||
P_dft(s_idx) = 1/abs((F_RF_dft'*a_test)'*inv_dft*(F_RF_dft'*a_test));
|
||||
P_full(s_idx) = 1/abs((F_RF_omp_full'*a_test)'*inv_full*(F_RF_omp_full'*a_test));
|
||||
P_tsc(s_idx) = 1/abs((F_RF_omp'*a_test)'*inv_tsc*(F_RF_omp'*a_test));
|
||||
end
|
||||
|
||||
[~, m1] = max(P_dft); err_dft_tmp(mc) = (scan_grid(m1) - true_target_angle)^2;
|
||||
[~, m2] = max(P_full); err_omp_full_tmp(mc) = (scan_grid(m2) - true_target_angle)^2;
|
||||
[~, m3] = max(P_tsc); err_omp_tmp(mc) = (scan_grid(m3) - true_target_angle)^2;
|
||||
end
|
||||
err_dft(a_idx, :) = err_dft_tmp; err_omp_full(a_idx, :) = err_omp_full_tmp; err_omp(a_idx, :) = err_omp_tmp;
|
||||
end
|
||||
RMSE_DFT_All(snr_idx, :) = sqrt(mean(err_dft, 2))'; RMSE_OMP_FULL_All(snr_idx, :) = sqrt(mean(err_omp_full, 2))'; RMSE_OMP_All(snr_idx, :) = sqrt(mean(err_omp, 2))';
|
||||
|
||||
c_str = color_array{snr_idx};
|
||||
semilogy(test_angles, RMSE_DFT_All(snr_idx, :), '--^', 'Color', c_str, 'LineWidth', 1.5, 'MarkerSize', 5);
|
||||
semilogy(test_angles, RMSE_OMP_FULL_All(snr_idx, :), ':x', 'Color', c_str, 'LineWidth', 2, 'MarkerSize', 7);
|
||||
semilogy(test_angles, RMSE_OMP_All(snr_idx, :), '-o', 'Color', c_str, 'LineWidth', 2.5, 'MarkerSize', 6);
|
||||
end
|
||||
xlim([min(test_angles), max(test_angles)]); xticks(min(test_angles):1:max(test_angles));
|
||||
set(gca, 'YScale', 'log', 'YMinorGrid', 'on', 'XMinorGrid', 'off', 'FontSize', 12); ylim([1e-3, 20]); yticks(10.^(-3:1:2));
|
||||
xlabel('目标真实出现角度 \theta (\circ)', 'FontSize', 14, 'FontWeight', 'bold'); ylabel('测角均方根误差 RMSE (\circ)', 'FontSize', 14, 'FontWeight', 'bold'); title(sprintf('宽波束先验引导下,扇区 (%d°~%d°) 动态波门鲁棒性评估', min(test_angles), max(test_angles)), 'FontSize', 15, 'FontWeight', 'bold');
|
||||
|
||||
legend_str = {};
|
||||
for snr_idx = 1:length(snr_test_list)
|
||||
legend_str{end+1} = sprintf('SNR=%ddB, DFT', snr_test_list(snr_idx)); legend_str{end+1} = sprintf('SNR=%ddB, 传统OMP(发散)', snr_test_list(snr_idx)); legend_str{end+1} = sprintf('SNR=%ddB, TSC-OMP(稳健)', snr_test_list(snr_idx));
|
||||
end
|
||||
lgd = legend(legend_str, 'Location', 'southoutside', 'FontSize', 10); set(lgd, 'NumColumns', 3);
|
||||
end
|
||||
disp('全部四组扇区极速测试完毕!');
|
||||
@@ -0,0 +1,174 @@
|
||||
%% 大规模阵列 MIMO 数模混合波束成形综合仿真平台 (6核全自动并行加速版)
|
||||
clear; clc; close all;
|
||||
rng(2024); % 锁定随机数种子
|
||||
|
||||
%% 1. 系统参数与环境初始化
|
||||
par = setpar();
|
||||
N_ant = par.M; N_RF = par.Nrf;
|
||||
target_angle = 5;
|
||||
clutter_angle = -15;
|
||||
theta_scan = -90:0.1:90;
|
||||
|
||||
a_target = gen_a(target_angle, par);
|
||||
a_clutter = gen_a(clutter_angle, par);
|
||||
L_snapshots = 100; noise_power = 1; INR_linear = 10^(30/10);
|
||||
|
||||
c_ana = '#77AC30'; c_dig = '#000000'; c_dft = '#4DBEEE'; c_omp = '#D95319'; c_tsc = '#0072BD';
|
||||
|
||||
%% 2. 动态生成目标约束扇区与字典
|
||||
theta_coarse_prior = round(target_angle); window_half_width = 5;
|
||||
dic_angles_constrained = (theta_coarse_prior - window_half_width) : 0.1 : (theta_coarse_prior + window_half_width);
|
||||
A_dic_constrained = zeros(N_ant, length(dic_angles_constrained));
|
||||
for i = 1:length(dic_angles_constrained)
|
||||
A_dic_constrained(:,i) = gen_a(dic_angles_constrained(i), par) / sqrt(N_ant);
|
||||
end
|
||||
dic_angles_full = -90:0.5:90;
|
||||
A_dic_full = zeros(N_ant, length(dic_angles_full));
|
||||
for i = 1:length(dic_angles_full)
|
||||
A_dic_full(:,i) = gen_a(dic_angles_full(i), par) / sqrt(N_ant);
|
||||
end
|
||||
F_RF_dft = zeros(N_ant, N_RF);
|
||||
angles_dft = linspace(theta_coarse_prior - window_half_width, theta_coarse_prior + window_half_width, N_RF);
|
||||
for i = 1:N_RF
|
||||
F_RF_dft(:,i) = gen_a(angles_dft(i), par)/sqrt(N_ant);
|
||||
end
|
||||
|
||||
%% 3. 提取防发散的慢时间先验射频矩阵 F_RF
|
||||
disp('正在进行慢时间 (Slow-Time) 模拟射频网络硬件配置...');
|
||||
a_prior_broad = zeros(N_ant, 1);
|
||||
for theta_b = (theta_coarse_prior - window_half_width) : 0.5 : (theta_coarse_prior + window_half_width)
|
||||
a_prior_broad = a_prior_broad + gen_a(theta_b, par);
|
||||
end
|
||||
a_prior_broad = a_prior_broad / norm(a_prior_broad);
|
||||
R_interference_prior = noise_power * eye(N_ant) + INR_linear * (a_clutter * a_clutter');
|
||||
w_digital_prior = inv(R_interference_prior) * a_prior_broad;
|
||||
w_digital_prior = w_digital_prior / norm(w_digital_prior);
|
||||
|
||||
res_full = w_digital_prior; F_RF_omp_full_static = [];
|
||||
for k=1:N_RF
|
||||
proj = A_dic_full' * res_full; [~, midx] = max(abs(proj));
|
||||
F_RF_omp_full_static = [F_RF_omp_full_static, A_dic_full(:, midx)];
|
||||
res_full = w_digital_prior - F_RF_omp_full_static * (pinv(F_RF_omp_full_static)*w_digital_prior);
|
||||
end
|
||||
res_tsc = w_digital_prior; F_RF_omp_static = [];
|
||||
for k=1:N_RF
|
||||
proj = A_dic_constrained' * res_tsc; [~, midx] = max(abs(proj));
|
||||
F_RF_omp_static = [F_RF_omp_static, A_dic_constrained(:, midx)];
|
||||
res_tsc = w_digital_prior - F_RF_omp_static * (pinv(F_RF_omp_static)*w_digital_prior);
|
||||
end
|
||||
|
||||
%% 4. 计算波束方向图并分别出图
|
||||
disp('正在计算并绘制静态波束方向图...');
|
||||
w_analog = a_target / sqrt(N_ant);
|
||||
sig_power_static = 10^(0/10);
|
||||
s_t_static = sqrt(sig_power_static/2)*(randn(1,L_snapshots)+1j*randn(1,L_snapshots));
|
||||
s_c_static = sqrt(INR_linear/2)*(randn(1,L_snapshots)+1j*randn(1,L_snapshots));
|
||||
n_t_static = sqrt(noise_power/2)*(randn(N_ant,L_snapshots)+1j*randn(N_ant,L_snapshots));
|
||||
X_static = a_target*s_t_static + a_clutter*s_c_static + n_t_static;
|
||||
R_in_static = (X_static*X_static')/L_snapshots;
|
||||
|
||||
a_eff_dft = F_RF_dft'*a_target;
|
||||
R_eff_dft = F_RF_dft'*R_in_static*F_RF_dft; R_eff_dft = R_eff_dft + (0.05*trace(R_eff_dft)/N_RF)*eye(N_RF);
|
||||
F_BB_dft = (inv(R_eff_dft)*a_eff_dft)/(a_eff_dft'*inv(R_eff_dft)*a_eff_dft); w_hybrid_dft = F_RF_dft * F_BB_dft; w_hybrid_dft = w_hybrid_dft / norm(w_hybrid_dft);
|
||||
a_eff_full = F_RF_omp_full_static'*a_target;
|
||||
R_eff_full = F_RF_omp_full_static'*R_in_static*F_RF_omp_full_static; R_eff_full = R_eff_full + (0.05*trace(R_eff_full)/N_RF)*eye(N_RF);
|
||||
F_BB_full = (inv(R_eff_full)*a_eff_full)/(a_eff_full'*inv(R_eff_full)*a_eff_full); w_hybrid_omp_full = F_RF_omp_full_static * F_BB_full; w_hybrid_omp_full = w_hybrid_omp_full / norm(w_hybrid_omp_full);
|
||||
a_eff_tsc = F_RF_omp_static'*a_target;
|
||||
R_eff_tsc = F_RF_omp_static'*R_in_static*F_RF_omp_static; R_eff_tsc = R_eff_tsc + (0.05*trace(R_eff_tsc)/N_RF)*eye(N_RF);
|
||||
F_BB_tsc = (inv(R_eff_tsc)*a_eff_tsc)/(a_eff_tsc'*inv(R_eff_tsc)*a_eff_tsc); w_hybrid_omp = F_RF_omp_static * F_BB_tsc; w_hybrid_omp = w_hybrid_omp / norm(w_hybrid_omp);
|
||||
R_in_static_dl = R_in_static + (0.01*trace(R_in_static)/N_ant)*eye(N_ant);
|
||||
w_digital_ideal = (inv(R_in_static_dl)*a_target)/(a_target'*inv(R_in_static_dl)*a_target); w_digital_ideal = w_digital_ideal / norm(w_digital_ideal);
|
||||
|
||||
BP_Analog = zeros(1,length(theta_scan)); BP_Digital = zeros(1,length(theta_scan)); BP_DFT = zeros(1,length(theta_scan)); BP_OMP_FULL = zeros(1,length(theta_scan)); BP_OMP = zeros(1,length(theta_scan));
|
||||
for p = 1:length(theta_scan)
|
||||
a_scan = gen_a(theta_scan(p), par);
|
||||
BP_Analog(p) = abs(a_scan'*w_analog)^2; BP_Digital(p) = abs(a_scan'*w_digital_ideal)^2;
|
||||
BP_DFT(p) = abs(a_scan'*w_hybrid_dft)^2; BP_OMP_FULL(p) = abs(a_scan'*w_hybrid_omp_full)^2; BP_OMP(p) = abs(a_scan'*w_hybrid_omp)^2;
|
||||
end
|
||||
BP_Analog_dB = 10*log10(BP_Analog/max(BP_Analog)); BP_Digital_dB = 10*log10(BP_Digital/max(BP_Digital)); BP_DFT_dB = 10*log10(BP_DFT/max(BP_DFT)); BP_OMP_FULL_dB = 10*log10(BP_OMP_FULL/max(BP_OMP_FULL)); BP_OMP_dB = 10*log10(BP_OMP/max(BP_OMP));
|
||||
|
||||
label_target = sprintf('机动目标 (%.2f°)', target_angle); label_clutter = sprintf('强杂波 (%.2f°)', clutter_angle);
|
||||
bp_labels = {'纯模拟架构', '纯数字架构(性能上限)', 'DFT混合架构', '传统OMP混合架构', 'TSC-OMP混合架构(所提)'};
|
||||
bp_data = {BP_Analog_dB, BP_Digital_dB, BP_DFT_dB, BP_OMP_FULL_dB, BP_OMP_dB}; bp_colors = {c_ana, c_dig, c_dft, c_omp, c_tsc};
|
||||
for fig_i = 1:5
|
||||
figure('Name', ['波束图-', bp_labels{fig_i}], 'Position', [100+fig_i*20, 100+fig_i*20, 600, 400]);
|
||||
plot(theta_scan, bp_data{fig_i}, 'Color', bp_colors{fig_i}, 'LineWidth', 2.5); hold on;
|
||||
xline(target_angle,'g--','LineWidth',1.5,'Label',label_target,'LabelVerticalAlignment','top','LabelHorizontalAlignment','right','FontSize',11,'FontWeight','bold');
|
||||
xline(clutter_angle,'r--','LineWidth',1.5,'Label',label_clutter,'LabelVerticalAlignment','top','LabelHorizontalAlignment','left','FontSize',11,'FontWeight','bold');
|
||||
grid on; ylim([-60, 0]); xlim([-60, 60]); title(sprintf('波束方向图 - %s', bp_labels{fig_i})); xlabel('角度 (\circ)'); ylabel('归一化增益 (dB)'); legend('主波束', 'Location', 'southwest');
|
||||
end
|
||||
figure('Name', '波束图-五大架构汇总', 'Position', [250, 250, 700, 500]);
|
||||
plot(theta_scan, BP_Analog_dB, '-.', 'Color', c_ana, 'LineWidth', 2); hold on; plot(theta_scan, BP_Digital_dB, '-', 'Color', c_dig, 'LineWidth', 2); plot(theta_scan, BP_DFT_dB, '--', 'Color', c_dft, 'LineWidth', 2); plot(theta_scan, BP_OMP_FULL_dB, ':', 'Color', c_omp, 'LineWidth', 2.5); plot(theta_scan, BP_OMP_dB, '-', 'Color', c_tsc, 'LineWidth', 2.5);
|
||||
xline(target_angle,'g--','LineWidth',2,'Label',label_target,'LabelVerticalAlignment','top','LabelHorizontalAlignment','right','FontSize',11,'FontWeight','bold'); xline(clutter_angle,'r--','LineWidth',2,'Label',label_clutter,'LabelVerticalAlignment','top','LabelHorizontalAlignment','left','FontSize',11,'FontWeight','bold');
|
||||
grid on; ylim([-60, 0]); xlim([-60, 60]); title('五大架构波束方向图综合对比'); xlabel('角度 (\circ)'); ylabel('归一化增益 (dB)'); legend('纯模拟', '纯数字', 'DFT', '传统OMP', 'TSC-OMP(所提)', 'Location', 'southwest');
|
||||
|
||||
%% 5. 理论 SINR 对比
|
||||
disp('正在计算 SINR...');
|
||||
SNR_dB_range = -10:2:20;
|
||||
SINR_Ana = zeros(1,length(SNR_dB_range)); SINR_Dig = zeros(1,length(SNR_dB_range)); SINR_DFT = zeros(1,length(SNR_dB_range)); SINR_OMP_FULL = zeros(1,length(SNR_dB_range)); SINR_OMP = zeros(1,length(SNR_dB_range));
|
||||
for i = 1:length(SNR_dB_range)
|
||||
sig_power = 10^(SNR_dB_range(i)/10);
|
||||
SINR_Ana(i) = (sig_power*abs(w_analog'*a_target)^2) / real(w_analog'*R_interference_prior*w_analog);
|
||||
SINR_Dig(i) = (sig_power*abs(w_digital_ideal'*a_target)^2) / real(w_digital_ideal'*R_interference_prior*w_digital_ideal);
|
||||
SINR_DFT(i) = (sig_power*abs(w_hybrid_dft'*a_target)^2) / real(w_hybrid_dft'*R_interference_prior*w_hybrid_dft);
|
||||
SINR_OMP_FULL(i) = (sig_power*abs(w_hybrid_omp_full'*a_target)^2) / real(w_hybrid_omp_full'*R_interference_prior*w_hybrid_omp_full);
|
||||
SINR_OMP(i) = (sig_power*abs(w_hybrid_omp'*a_target)^2) / real(w_hybrid_omp'*R_interference_prior*w_hybrid_omp);
|
||||
end
|
||||
figure('Name', 'SINR理论性能', 'Position', [300, 300, 600, 450]);
|
||||
plot(SNR_dB_range,10*log10(SINR_Ana),'-v','Color',c_ana,'LineWidth',2); hold on; plot(SNR_dB_range,10*log10(SINR_Dig),'-o','Color',c_dig,'LineWidth',2); plot(SNR_dB_range,10*log10(SINR_DFT),'-d','Color',c_dft,'LineWidth',2); plot(SNR_dB_range,10*log10(SINR_OMP_FULL),'-*','Color',c_omp,'LineWidth',2); plot(SNR_dB_range,10*log10(SINR_OMP),'-s','Color',c_tsc,'LineWidth',2);
|
||||
grid on; legend('纯模拟','纯数字(上限)','DFT','传统OMP','TSC-OMP(所提)','Location','northwest'); xlabel('输入 SNR (dB)'); ylabel('输出 SINR (dB)'); title(sprintf('全架构输出 SINR 理论逼近测试 (目标 %.2f°)', target_angle));
|
||||
|
||||
%% 6. 【核心优化:开启多核并行计算 parfor】
|
||||
disp('正在调用多核并行引擎 (parfor) 加速蒙特卡洛评估...');
|
||||
Monte_Carlo = 500;
|
||||
RMSE_Ana = zeros(1,length(SNR_dB_range)); RMSE_Dig = zeros(1,length(SNR_dB_range)); RMSE_DFT = zeros(1,length(SNR_dB_range)); RMSE_OMP_FULL = zeros(1,length(SNR_dB_range)); RMSE_OMP = zeros(1,length(SNR_dB_range));
|
||||
scan_grid = (theta_coarse_prior - window_half_width) : 0.01 : (theta_coarse_prior + window_half_width);
|
||||
|
||||
for i = 1:length(SNR_dB_range)
|
||||
sig_power = 10^(SNR_dB_range(i)/10);
|
||||
err_ana = zeros(1,Monte_Carlo); err_dig = zeros(1,Monte_Carlo); err_dft = zeros(1,Monte_Carlo); err_omp_full = zeros(1,Monte_Carlo); err_omp = zeros(1,Monte_Carlo);
|
||||
|
||||
% 替换为 parfor,自动分配到 6 个 Worker 极速计算!
|
||||
parfor mc = 1:Monte_Carlo
|
||||
s_t = sqrt(sig_power/2)*(randn(1,L_snapshots)+1j*randn(1,L_snapshots));
|
||||
s_c = sqrt(INR_linear/2)*(randn(1,L_snapshots)+1j*randn(1,L_snapshots));
|
||||
n_t = sqrt(noise_power/2)*(randn(N_ant,L_snapshots)+1j*randn(N_ant,L_snapshots));
|
||||
X_mc = a_target*s_t + a_clutter*s_c + n_t;
|
||||
|
||||
R_in = (X_mc*X_mc')/L_snapshots;
|
||||
inv_R_in_dl = inv(R_in + (0.05*trace(R_in)/N_ant)*eye(N_ant));
|
||||
|
||||
Y_dft = F_RF_dft' * X_mc; Y_omp_full = F_RF_omp_full_static' * X_mc; Y_omp = F_RF_omp_static' * X_mc;
|
||||
R_dft = (Y_dft*Y_dft')/L_snapshots; R_omp_full = (Y_omp_full*Y_omp_full')/L_snapshots; R_omp = (Y_omp*Y_omp')/L_snapshots;
|
||||
|
||||
R_dft = R_dft + (0.05*trace(R_dft)/N_RF)*eye(N_RF); R_omp_full = R_omp_full + (0.05*trace(R_omp_full)/N_RF)*eye(N_RF); R_omp = R_omp + (0.05*trace(R_omp)/N_RF)*eye(N_RF);
|
||||
inv_R_dft = inv(R_dft); inv_R_omp_full = inv(R_omp_full); inv_R_omp = inv(R_omp);
|
||||
|
||||
P_ana = zeros(1,length(scan_grid)); P_dig = zeros(1,length(scan_grid)); P_dft = zeros(1,length(scan_grid)); P_omp_full = zeros(1,length(scan_grid)); P_omp = zeros(1,length(scan_grid));
|
||||
|
||||
for a_idx = 1:length(scan_grid)
|
||||
a_test = gen_a(scan_grid(a_idx), par);
|
||||
P_ana(a_idx) = abs(a_test' * R_in * a_test); P_dig(a_idx) = 1/abs(a_test' * inv_R_in_dl * a_test);
|
||||
P_dft(a_idx) = 1/abs((F_RF_dft'*a_test)'*inv_R_dft*(F_RF_dft'*a_test));
|
||||
P_omp_full(a_idx) = 1/abs((F_RF_omp_full_static'*a_test)'*inv_R_omp_full*(F_RF_omp_full_static'*a_test));
|
||||
P_omp(a_idx) = 1/abs((F_RF_omp_static'*a_test)'*inv_R_omp*(F_RF_omp_static'*a_test));
|
||||
end
|
||||
|
||||
[~,id1] = max(P_ana); [~,id2] = max(P_dig); [~,id3] = max(P_dft); [~,id4] = max(P_omp_full); [~,id5] = max(P_omp);
|
||||
err_ana(mc) = (scan_grid(id1)-target_angle)^2; err_dig(mc) = (scan_grid(id2)-target_angle)^2;
|
||||
err_dft(mc) = (scan_grid(id3)-target_angle)^2; err_omp_full(mc) = (scan_grid(id4)-target_angle)^2; err_omp(mc) = (scan_grid(id5)-target_angle)^2;
|
||||
end
|
||||
RMSE_Ana(i) = sqrt(mean(err_ana)); RMSE_Dig(i) = sqrt(mean(err_dig)); RMSE_DFT(i) = sqrt(mean(err_dft)); RMSE_OMP_FULL(i) = sqrt(mean(err_omp_full)); RMSE_OMP(i) = sqrt(mean(err_omp));
|
||||
end
|
||||
|
||||
%% 7. 绘制分别与汇总的 RMSE 对比图
|
||||
rmse_data = {RMSE_Ana, RMSE_Dig, RMSE_DFT, RMSE_OMP_FULL, RMSE_OMP};
|
||||
for fig_i = 1:5
|
||||
figure('Name', ['RMSE-', bp_labels{fig_i}], 'Position', [150+fig_i*20, 150+fig_i*20, 600, 400]);
|
||||
semilogy(SNR_dB_range, rmse_data{fig_i}, '-o', 'Color', bp_colors{fig_i}, 'LineWidth', 2.5);
|
||||
grid on; ylim([1e-3, 20]); yticks(10.^(-3:1:2)); title(sprintf('DOA 测角精度 - %s (目标 %.2f°)', bp_labels{fig_i}, target_angle)); xlabel('信噪比 SNR (dB)'); ylabel('测角 RMSE (\circ)');
|
||||
end
|
||||
figure('Name', 'RMSE-五大架构汇总', 'Position', [400, 400, 700, 500]);
|
||||
semilogy(SNR_dB_range, RMSE_Ana, '-.v', 'Color', c_ana, 'LineWidth', 2); hold on; semilogy(SNR_dB_range, RMSE_Dig, '-o', 'Color', c_dig, 'LineWidth', 2); semilogy(SNR_dB_range, RMSE_DFT, '--d', 'Color', c_dft, 'LineWidth', 2); semilogy(SNR_dB_range, RMSE_OMP_FULL, ':*', 'Color', c_omp, 'LineWidth', 2.5); semilogy(SNR_dB_range, RMSE_OMP, '-s', 'Color', c_tsc, 'LineWidth', 2.5);
|
||||
grid on; ylim([1e-3, 20]); yticks(10.^(-3:1:2)); title(sprintf('全架构实战恶劣环境 DOA 测角精度评估 (机动目标 %.2f°)', target_angle)); xlabel('信噪比 SNR (dB)'); ylabel('测角 RMSE (\circ)'); legend('纯模拟 (极低分辨率)','纯数字 (性能上限)','DFT','传统OMP (极易发散)','TSC-OMP (所提稳健算法)','Location','southwest');
|
||||
disp('并行测试完成!');
|
||||
@@ -0,0 +1,29 @@
|
||||
function [ par ] = setpar
|
||||
|
||||
par.M = 24; % number of TX/RX antennas
|
||||
par.Nrf = 4; % number of RF chains
|
||||
par.Ns = 1; % number of data streams
|
||||
|
||||
par.fc = 4.9e9; % system bandwidth
|
||||
|
||||
par.Et = 1; % power
|
||||
par.Er = 1;
|
||||
%%----------------------------------ULAs-----------------------------------
|
||||
|
||||
par.labmda = 3e8/par.fc;
|
||||
a = [0:43:11*43]';
|
||||
b = 11*43 + 63;
|
||||
c = 11*43 + 63 + [43:43:11*43]';
|
||||
par.delta_d = 1e-3*[a;b;c];
|
||||
|
||||
phi = [-355, -240, -275, -160, -195, -80, -115, 0].';
|
||||
par.pos = exp(1j*phi);
|
||||
% PS network
|
||||
par.K1 = [1 1 1 1; 0 0 0 1; 1 1 1 1];
|
||||
par.K2 = [1 1; 0 1];
|
||||
|
||||
%%--------------------discrete frequency space-----------------------------
|
||||
par.diff_p = 0.05; % total discreted points
|
||||
par.theta = round(-90:par.diff_p:90, 6);
|
||||
|
||||
end
|
||||
Reference in New Issue
Block a user