diff --git a/Graduation Design/Complexity_Test.py b/Graduation Design/Complexity_Test.py new file mode 100644 index 0000000..53facee --- /dev/null +++ b/Graduation Design/Complexity_Test.py @@ -0,0 +1,641 @@ +import os +import time +import math +import numpy as np +import torch +import matplotlib.pyplot as plt + +# ========================================================= +# 0. 全局设置 +# ========================================================= +plt.rcParams['axes.unicode_minus'] = False +plt.rcParams['font.sans-serif'] = ['Microsoft YaHei'] +plt.rcParams['mathtext.fontset'] = 'stix' + +device = torch.device("cuda" if torch.cuda.is_available() else "cpu") +torch.manual_seed(2024) +np.random.seed(2024) + +print(f"使用设备: {device}") + +SAVE_DIR = "analysis_results_ula_nula" +os.makedirs(SAVE_DIR, exist_ok=True) + +ARRAY_TYPES = ["ULA", "NULA"] + +# ========================================================= +# 1. 基本参数 +# ========================================================= +N_ant_default = 24 +N_RF_default = 4 +target_angle_default = 5.0 +clutter_angle_default = -15.0 +noise_power = 1.0 +INR_linear = 10 ** (30 / 10) +L_snapshots = 100 + +# 分析配置 +SNR_dB_test = 10 +Monte_Carlo_test = 80 +SCAN_STEP = 0.02 + +# 颜色 +method_colors = { + "Analog": "#77AC30", + "Digital": "#000000", + "DFT": "#4DBEEE", + "OMP": "#D95319", + "TSC-OMP": "#0072BD", +} +array_colors = { + "ULA": "#7E2F8E", + "NULA": "#EDB120", +} + +# ========================================================= +# 2. 阵列几何 +# ========================================================= +def get_positions_lambda(array_type: str, N_ant: int): + """ + 返回阵列位置,单位统一为“波长 λ” + - ULA: 半波长均匀线阵 + - NULA: + * N_ant == 24 时,严格复刻你现有主线脚本中的非均匀阵列物理位置 + * 否则使用一个可复现的广义非均匀阵列,用于复杂度规模扫描 + """ + array_type = array_type.upper() + + if array_type == "ULA": + pos = 0.5 * np.arange(-(N_ant - 1) / 2, (N_ant - 1) / 2 + 1) + return pos.astype(np.float32) + + # 精确复刻你当前 NULA 主线 + if array_type == "NULA" and N_ant == 24: + fc = 4.9e9 + lam = 3e8 / fc + a = np.arange(0, 11 * 43 + 1, 43) + b = np.array([11 * 43 + 63]) + c = 11 * 43 + 63 + np.arange(43, 11 * 43 + 1, 43) + delta_d = 1e-3 * np.concatenate((a, b, c)) # meter + pos = delta_d / lam + return pos.astype(np.float32) + + # 广义 NULA,用于 N_ant != 24 的复杂度缩放分析 + base = 0.5 * np.arange(N_ant) + jitter = 0.08 * np.sin(np.linspace(0, 3 * np.pi, N_ant)) + pos = base + jitter + pos = pos - np.mean(pos) + return pos.astype(np.float32) + + +def gen_a(theta_deg, N_ant, array_type="NULA"): + """ + 导向矢量,统一使用 exp(-j 2π p sinθ) + p 为以 λ 为单位的阵元位置 + """ + pos = get_positions_lambda(array_type, N_ant) + pos_t = torch.tensor(pos, dtype=torch.float32, device=device).unsqueeze(1) + theta_rad = torch.deg2rad(torch.tensor(theta_deg, dtype=torch.float32, device=device)) + return torch.exp(-1j * 2 * torch.pi * pos_t * torch.sin(theta_rad)) + + +def build_dictionary(angle_grid, N_ant, array_type="NULA"): + return torch.cat( + [gen_a(float(th), N_ant, array_type) / math.sqrt(N_ant) for th in angle_grid], + dim=1 + ) + + +def sync_if_needed(): + if device.type == "cuda": + torch.cuda.synchronize() + + +# ========================================================= +# 3. 五种架构模块 +# ========================================================= +def extract_omp(w_opt, A_dic, N_RF): + """ + 传统 OMP / TSC-OMP 公用提取函数 + """ + res = w_opt.clone() + F_RF = None + for _ in range(N_RF): + proj = A_dic.mH @ res + idx = torch.argmax(torch.abs(proj)) + atom = A_dic[:, idx:idx+1] + F_RF = atom if F_RF is None else torch.cat((F_RF, atom), dim=1) + res = w_opt - F_RF @ (torch.linalg.pinv(F_RF) @ w_opt) + return F_RF + + +def build_dft_rf(theta_prior, window_half_width, N_ant, N_RF, array_type="NULA"): + angles_dft = torch.linspace( + theta_prior - window_half_width, + theta_prior + window_half_width, + N_RF, + device=device + ) + return torch.cat( + [gen_a(float(th), N_ant, array_type) / math.sqrt(N_ant) for th in angles_dft], + dim=1 + ) + + +def make_prior_broad(theta_prior, window_half_width, N_ant, array_type="NULA"): + a_prior_broad = torch.zeros((N_ant, 1), dtype=torch.complex64, device=device) + for tb in torch.arange( + theta_prior - window_half_width, + theta_prior + window_half_width + 0.1, + 0.5, + device=device + ): + a_prior_broad += gen_a(float(tb), N_ant, array_type) + a_prior_broad /= torch.norm(a_prior_broad) + return a_prior_broad + + +def calc_hybrid_weight(F_RF, a_target, R_in, N_RF): + a_eff = F_RF.mH @ a_target + R_eff = F_RF.mH @ R_in @ F_RF + R_eff = R_eff + (0.05 * torch.trace(R_eff).real / N_RF) * torch.eye(N_RF, device=device) + inv_R_eff = torch.linalg.inv(R_eff) + F_BB = (inv_R_eff @ a_eff) / (a_eff.mH @ inv_R_eff @ a_eff) + w = F_RF @ F_BB + return w / torch.norm(w) + + +def build_all_structures(true_target_angle, clutter_angle, N_ant, N_RF, + theta_prior, window_half_width, array_type="NULA"): + """ + 构建五种架构所需的静态权值/射频矩阵 + """ + a_target = gen_a(true_target_angle, N_ant, array_type) + a_clutter = gen_a(clutter_angle, N_ant, array_type) + + R_interf = noise_power * torch.eye(N_ant, device=device) + INR_linear * (a_clutter @ a_clutter.mH) + + # 用于 OMP / TSC-OMP 的先验宽波束 + a_prior_broad = make_prior_broad(theta_prior, window_half_width, N_ant, array_type) + w_digital_prior = torch.linalg.inv(R_interf) @ a_prior_broad + w_digital_prior /= torch.norm(w_digital_prior) + + full_grid = np.arange(-90, 90.5, 0.5) + tsc_grid = np.arange( + theta_prior - window_half_width, + theta_prior + window_half_width + 0.1, + 0.1 + ) + + A_dic_full = build_dictionary(full_grid, N_ant, array_type) + A_dic_tsc = build_dictionary(tsc_grid, N_ant, array_type) + + F_RF_dft = build_dft_rf(theta_prior, window_half_width, N_ant, N_RF, array_type) + F_RF_omp = extract_omp(w_digital_prior, A_dic_full, N_RF) + F_RF_tsc = extract_omp(w_digital_prior, A_dic_tsc, N_RF) + + # 静态协方差,用于纯数字与混合权值展示 + sig_power_static = 10 ** (0 / 10) + R_static = (noise_power * torch.eye(N_ant, device=device) + + sig_power_static * (a_target @ a_target.mH) + + INR_linear * (a_clutter @ a_clutter.mH)) + R_static_dl = R_static + (0.01 * torch.trace(R_static).real / N_ant) * torch.eye(N_ant, device=device) + + # 纯模拟 + w_analog = a_target / torch.norm(a_target) + + # 纯数字 MVDR + inv_R = torch.linalg.inv(R_static_dl) + w_digital = (inv_R @ a_target) / (a_target.mH @ inv_R @ a_target) + w_digital = w_digital / torch.norm(w_digital) + + return { + "a_target": a_target, + "a_clutter": a_clutter, + "w_analog": w_analog, + "w_digital": w_digital, + "F_RF_dft": F_RF_dft, + "F_RF_omp": F_RF_omp, + "F_RF_tsc": F_RF_tsc, + } + + +# ========================================================= +# 4. 统一 RMSE 评估 +# ========================================================= +def estimate_rmse_all_methods(true_target_angle, clutter_angle, N_ant, N_RF, + snr_db, theta_prior, window_half_width, + array_type="NULA", mc=80, scan_step=0.02): + """ + 返回: + analog_rmse, digital_rmse, dft_rmse, omp_rmse, tsc_rmse + """ + sig_power = 10 ** (snr_db / 10) + + structures = build_all_structures( + true_target_angle, clutter_angle, N_ant, N_RF, + theta_prior, window_half_width, array_type + ) + + a_target = structures["a_target"] + a_clutter = structures["a_clutter"] + w_analog = structures["w_analog"] + w_digital = structures["w_digital"] + F_RF_dft = structures["F_RF_dft"] + F_RF_omp = structures["F_RF_omp"] + F_RF_tsc = structures["F_RF_tsc"] + + scan_grid = torch.arange( + theta_prior - window_half_width, + theta_prior + window_half_width + 0.001, + scan_step, + device=device + ) + A_scan = torch.cat([gen_a(float(th), N_ant, array_type) for th in scan_grid], dim=1) + + s_t = np.sqrt(sig_power / 2) * ( + torch.randn(mc, 1, L_snapshots, device=device) + + 1j * torch.randn(mc, 1, L_snapshots, device=device) + ) + s_c = np.sqrt(INR_linear / 2) * ( + torch.randn(mc, 1, L_snapshots, device=device) + + 1j * torch.randn(mc, 1, L_snapshots, device=device) + ) + n_t = np.sqrt(noise_power / 2) * ( + torch.randn(mc, N_ant, L_snapshots, device=device) + + 1j * torch.randn(mc, N_ant, L_snapshots, device=device) + ) + X = a_target.unsqueeze(0) * s_t + a_clutter.unsqueeze(0) * s_c + n_t + + # -------- 纯模拟:Bartlett -------- + P_ana = torch.abs(torch.einsum('si,bij,js->bs', A_scan.conj().T, torch.bmm(X, X.mH) / L_snapshots, A_scan)) + ang_ana = scan_grid[torch.argmax(P_ana, dim=1)] + + # -------- 纯数字 MVDR:Capon -------- + R_in = torch.bmm(X, X.mH) / L_snapshots + trace_R_in = torch.diagonal(R_in, dim1=-2, dim2=-1).sum(-1).real + R_in_dl = R_in + (0.05 * trace_R_in / N_ant).view(-1, 1, 1) * torch.eye(N_ant, device=device).unsqueeze(0) + inv_R_in = torch.linalg.inv(R_in_dl) + P_dig = 1.0 / torch.abs(torch.einsum('si,bij,js->bs', A_scan.conj().T, inv_R_in, A_scan)) + ang_dig = scan_grid[torch.argmax(P_dig, dim=1)] + + # -------- 混合三种 -------- + def process_hybrid(F_RF): + A_eff = F_RF.mH @ A_scan + Y = F_RF.mH.unsqueeze(0) @ X + R_hyb = torch.bmm(Y, Y.mH) / L_snapshots + trace_R = torch.diagonal(R_hyb, dim1=-2, dim2=-1).sum(-1).real + R_hyb = R_hyb + (0.05 * trace_R / N_RF).view(-1, 1, 1) * torch.eye(N_RF, device=device).unsqueeze(0) + inv_R_hyb = torch.linalg.inv(R_hyb) + P = 1.0 / torch.abs(torch.einsum('si,bij,js->bs', A_eff.conj().T, inv_R_hyb, A_eff)) + return scan_grid[torch.argmax(P, dim=1)] + + ang_dft = process_hybrid(F_RF_dft) + ang_omp = process_hybrid(F_RF_omp) + ang_tsc = process_hybrid(F_RF_tsc) + + def rmse(x): + return torch.sqrt(torch.mean((x - true_target_angle) ** 2)).item() + + return rmse(ang_ana), rmse(ang_dig), rmse(ang_dft), rmse(ang_omp), rmse(ang_tsc) + + +# ========================================================= +# 5. 复杂度分析 +# ========================================================= +def complexity_proxy(N_ant, N_RF, D_full, D_tsc): + """ + 五种架构的复杂度代理量 + """ + c_ana = N_ant + c_dig = N_ant ** 3 + c_dft = N_ant * N_RF + N_RF ** 3 + c_omp = N_RF * N_ant * D_full + N_RF ** 4 + c_tsc = N_RF * N_ant * D_tsc + N_RF ** 4 + return c_ana, c_dig, c_dft, c_omp, c_tsc + + +def run_complexity_analysis(): + print("\n========== 开始复杂度分析(ULA + NULA) ==========") + + target_angle = target_angle_default + clutter_angle = clutter_angle_default + N_RF = N_RF_default + window_half_width = 5.0 + repeat_times = 20 + + N_ant_list = [8, 12, 16, 20, 24, 28, 32] + results_runtime = {} + results_proxy = {} + + for array_type in ARRAY_TYPES: + print(f"\n--- 复杂度分析: {array_type} ---") + t_ana, t_dig, t_dft, t_omp, t_tsc = [], [], [], [], [] + c_ana_list, c_dig_list, c_dft_list, c_omp_list, c_tsc_list = [], [], [], [], [] + + for N_ant in N_ant_list: + a_target = gen_a(target_angle, N_ant, array_type) + a_clutter = gen_a(clutter_angle, N_ant, array_type) + + theta_prior = round(target_angle) + full_grid = np.arange(-90, 90.5, 0.5) + tsc_grid = np.arange(theta_prior - window_half_width, + theta_prior + window_half_width + 0.1, + 0.1) + + A_dic_full = build_dictionary(full_grid, N_ant, array_type) + A_dic_tsc = build_dictionary(tsc_grid, N_ant, array_type) + + R_interf = noise_power * torch.eye(N_ant, device=device) + INR_linear * (a_clutter @ a_clutter.mH) + a_prior_broad = make_prior_broad(theta_prior, window_half_width, N_ant, array_type) + + # 纯模拟 + sync_if_needed() + tic = time.perf_counter() + for _ in range(repeat_times): + _ = a_target / torch.norm(a_target) + sync_if_needed() + t_ana.append((time.perf_counter() - tic) / repeat_times * 1e3) + + # 纯数字 MVDR + sync_if_needed() + tic = time.perf_counter() + for _ in range(repeat_times): + R_static = (noise_power * torch.eye(N_ant, device=device) + + (a_target @ a_target.mH) + + INR_linear * (a_clutter @ a_clutter.mH)) + R_static_dl = R_static + (0.01 * torch.trace(R_static).real / N_ant) * torch.eye(N_ant, device=device) + inv_R = torch.linalg.inv(R_static_dl) + _ = (inv_R @ a_target) / (a_target.mH @ inv_R @ a_target) + sync_if_needed() + t_dig.append((time.perf_counter() - tic) / repeat_times * 1e3) + + # DFT + sync_if_needed() + tic = time.perf_counter() + for _ in range(repeat_times): + F_RF_dft = build_dft_rf(theta_prior, window_half_width, N_ant, N_RF, array_type) + R_static = (noise_power * torch.eye(N_ant, device=device) + + (a_target @ a_target.mH) + + INR_linear * (a_clutter @ a_clutter.mH)) + _ = calc_hybrid_weight(F_RF_dft, a_target, R_static, N_RF) + sync_if_needed() + t_dft.append((time.perf_counter() - tic) / repeat_times * 1e3) + + # OMP/TSC + w_digital_prior = torch.linalg.inv(R_interf) @ a_prior_broad + w_digital_prior /= torch.norm(w_digital_prior) + + sync_if_needed() + tic = time.perf_counter() + for _ in range(repeat_times): + _ = extract_omp(w_digital_prior, A_dic_full, N_RF) + sync_if_needed() + t_omp.append((time.perf_counter() - tic) / repeat_times * 1e3) + + sync_if_needed() + tic = time.perf_counter() + for _ in range(repeat_times): + _ = extract_omp(w_digital_prior, A_dic_tsc, N_RF) + sync_if_needed() + t_tsc.append((time.perf_counter() - tic) / repeat_times * 1e3) + + c_ana, c_dig, c_dft, c_omp, c_tsc = complexity_proxy( + N_ant, N_RF, len(full_grid), len(tsc_grid) + ) + c_ana_list.append(c_ana) + c_dig_list.append(c_dig) + c_dft_list.append(c_dft) + c_omp_list.append(c_omp) + c_tsc_list.append(c_tsc) + + results_runtime[array_type] = { + "N_ant": N_ant_list, + "Analog": t_ana, + "Digital": t_dig, + "DFT": t_dft, + "OMP": t_omp, + "TSC-OMP": t_tsc, + } + results_proxy[array_type] = { + "N_ant": N_ant_list, + "Analog": c_ana_list, + "Digital": c_dig_list, + "DFT": c_dft_list, + "OMP": c_omp_list, + "TSC-OMP": c_tsc_list, + } + + # 每种阵列单独出图 + plt.figure(figsize=(10, 6)) + for method in ["Analog", "Digital", "DFT", "OMP", "TSC-OMP"]: + plt.plot(N_ant_list, results_runtime[array_type][method], linewidth=2, label=method, color=method_colors[method]) + plt.grid(True, ls='--', alpha=0.6) + plt.xlabel('阵元数 N_ant') + plt.ylabel('平均运行时间 (ms)') + plt.title(f'复杂度分析:运行时间随阵元数变化 ({array_type})') + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(SAVE_DIR, f'complexity_runtime_vs_ant_{array_type}.png'), dpi=300) + plt.show() + + plt.figure(figsize=(10, 6)) + for method in ["Analog", "Digital", "DFT", "OMP", "TSC-OMP"]: + plt.plot(N_ant_list, results_proxy[array_type][method], linewidth=2, label=method, color=method_colors[method]) + plt.grid(True, ls='--', alpha=0.6) + plt.xlabel('阵元数 N_ant') + plt.ylabel('复杂度代理量') + plt.title(f'复杂度分析:复杂度代理量随阵元数变化 ({array_type})') + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(SAVE_DIR, f'complexity_proxy_vs_ant_{array_type}.png'), dpi=300) + plt.show() + + # ULA vs NULA 对照图(固定 24 阵元) + idx_24 = N_ant_list.index(24) + methods = ["Analog", "Digital", "DFT", "OMP", "TSC-OMP"] + ula_vals = [results_runtime["ULA"][m][idx_24] for m in methods] + nula_vals = [results_runtime["NULA"][m][idx_24] for m in methods] + + x = np.arange(len(methods)) + width = 0.35 + + plt.figure(figsize=(10, 6)) + plt.bar(x - width/2, ula_vals, width, label='ULA', color=array_colors["ULA"]) + plt.bar(x + width/2, nula_vals, width, label='NULA', color=array_colors["NULA"]) + plt.xticks(x, methods) + plt.ylabel('平均运行时间 (ms)') + plt.title('复杂度分析:24阵元下 ULA 与 NULA 运行时间对比') + plt.grid(True, axis='y', ls='--', alpha=0.5) + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(SAVE_DIR, 'complexity_runtime_compare_ULA_NULA_24ant.png'), dpi=300) + plt.show() + + print("复杂度分析完成。") + return results_runtime, results_proxy + + +# ========================================================= +# 6. 局限性分析 +# ========================================================= +def run_limitation_analysis(): + print("\n========== 开始局限性分析(ULA + NULA) ==========") + + N_ant = N_ant_default + N_RF = N_RF_default + true_target_angle = target_angle_default + clutter_angle = clutter_angle_default + snr_db = SNR_dB_test + + limitation_results = {} + + for array_type in ARRAY_TYPES: + print(f"\n--- 局限性分析: {array_type} ---") + limitation_results[array_type] = {} + + # 6.1 先验误差敏感性 + prior_errors = np.arange(-8, 9, 1) + curves_prior = {m: [] for m in ["Analog", "Digital", "DFT", "OMP", "TSC-OMP"]} + + for e in prior_errors: + theta_prior = round(true_target_angle) + e + r = estimate_rmse_all_methods( + true_target_angle, clutter_angle, N_ant, N_RF, + snr_db, theta_prior, 5.0, array_type, + mc=Monte_Carlo_test, scan_step=SCAN_STEP + ) + for m, v in zip(["Analog", "Digital", "DFT", "OMP", "TSC-OMP"], r): + curves_prior[m].append(v) + + limitation_results[array_type]["prior_errors"] = prior_errors + limitation_results[array_type]["prior_curves"] = curves_prior + + plt.figure(figsize=(10, 6)) + for method in ["Analog", "Digital", "DFT", "OMP", "TSC-OMP"]: + plt.semilogy(prior_errors, curves_prior[method], linewidth=2, label=method, color=method_colors[method]) + plt.grid(True, which='both', ls='--', alpha=0.6) + plt.xlabel('先验角误差 (°)') + plt.ylabel('RMSE (°)') + plt.title(f'局限性分析:先验误差敏感性 ({array_type}, SNR={snr_db}dB)') + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(SAVE_DIR, f'limitation_prior_error_{array_type}.png'), dpi=300) + plt.show() + + # 6.2 波门宽度敏感性 + widths = [2, 3, 4, 5, 6, 8, 10] + curves_width = {m: [] for m in ["Analog", "Digital", "DFT", "OMP", "TSC-OMP"]} + + for w in widths: + theta_prior = round(true_target_angle) + r = estimate_rmse_all_methods( + true_target_angle, clutter_angle, N_ant, N_RF, + snr_db, theta_prior, float(w), array_type, + mc=Monte_Carlo_test, scan_step=SCAN_STEP + ) + for m, v in zip(["Analog", "Digital", "DFT", "OMP", "TSC-OMP"], r): + curves_width[m].append(v) + + limitation_results[array_type]["widths"] = widths + limitation_results[array_type]["width_curves"] = curves_width + + plt.figure(figsize=(10, 6)) + for method in ["Analog", "Digital", "DFT", "OMP", "TSC-OMP"]: + plt.semilogy(widths, curves_width[method], linewidth=2, label=method, color=method_colors[method]) + plt.grid(True, which='both', ls='--', alpha=0.6) + plt.xlabel('波门半宽 (°)') + plt.ylabel('RMSE (°)') + plt.title(f'局限性分析:波门宽度敏感性 ({array_type}, SNR={snr_db}dB)') + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(SAVE_DIR, f'limitation_window_width_{array_type}.png'), dpi=300) + plt.show() + + # 6.3 杂波逼近目标时的性能退化 + clutter_list = np.array([-25, -20, -15, -10, -8, -6, -4, -2]) + sep_deg = np.abs(clutter_list - true_target_angle) + curves_sep = {m: [] for m in ["Analog", "Digital", "DFT", "OMP", "TSC-OMP"]} + + for cang in clutter_list: + theta_prior = round(true_target_angle) + r = estimate_rmse_all_methods( + true_target_angle, float(cang), N_ant, N_RF, + snr_db, theta_prior, 5.0, array_type, + mc=Monte_Carlo_test, scan_step=SCAN_STEP + ) + for m, v in zip(["Analog", "Digital", "DFT", "OMP", "TSC-OMP"], r): + curves_sep[m].append(v) + + limitation_results[array_type]["sep_deg"] = sep_deg + limitation_results[array_type]["sep_curves"] = curves_sep + + plt.figure(figsize=(10, 6)) + for method in ["Analog", "Digital", "DFT", "OMP", "TSC-OMP"]: + plt.semilogy(sep_deg, curves_sep[method], linewidth=2, label=method, color=method_colors[method]) + plt.grid(True, which='both', ls='--', alpha=0.6) + plt.xlabel('目标-杂波角距 |θ_t-θ_c| (°)') + plt.ylabel('RMSE (°)') + plt.title(f'局限性分析:杂波逼近目标时的性能退化 ({array_type})') + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(SAVE_DIR, f'limitation_clutter_separation_{array_type}.png'), dpi=300) + plt.show() + + # ---------------- ULA vs NULA 对照图(重点方法可读性更高) ---------------- + # 先验误差:比较 TSC-OMP + plt.figure(figsize=(10, 6)) + for array_type in ARRAY_TYPES: + x = limitation_results[array_type]["prior_errors"] + y = limitation_results[array_type]["prior_curves"]["TSC-OMP"] + plt.semilogy(x, y, linewidth=2.5, label=f'TSC-OMP ({array_type})', color=array_colors[array_type]) + plt.grid(True, which='both', ls='--', alpha=0.6) + plt.xlabel('先验角误差 (°)') + plt.ylabel('RMSE (°)') + plt.title('ULA 与 NULA 对照:TSC-OMP 先验误差敏感性') + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(SAVE_DIR, 'compare_prior_error_TSC_ULA_NULA.png'), dpi=300) + plt.show() + + # 波门宽度:比较 TSC-OMP + plt.figure(figsize=(10, 6)) + for array_type in ARRAY_TYPES: + x = limitation_results[array_type]["widths"] + y = limitation_results[array_type]["width_curves"]["TSC-OMP"] + plt.semilogy(x, y, linewidth=2.5, label=f'TSC-OMP ({array_type})', color=array_colors[array_type]) + plt.grid(True, which='both', ls='--', alpha=0.6) + plt.xlabel('波门半宽 (°)') + plt.ylabel('RMSE (°)') + plt.title('ULA 与 NULA 对照:TSC-OMP 波门宽度敏感性') + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(SAVE_DIR, 'compare_window_width_TSC_ULA_NULA.png'), dpi=300) + plt.show() + + # 杂波逼近:比较 TSC-OMP + plt.figure(figsize=(10, 6)) + for array_type in ARRAY_TYPES: + x = limitation_results[array_type]["sep_deg"] + y = limitation_results[array_type]["sep_curves"]["TSC-OMP"] + plt.semilogy(x, y, linewidth=2.5, label=f'TSC-OMP ({array_type})', color=array_colors[array_type]) + plt.grid(True, which='both', ls='--', alpha=0.6) + plt.xlabel('目标-杂波角距 |θ_t-θ_c| (°)') + plt.ylabel('RMSE (°)') + plt.title('ULA 与 NULA 对照:TSC-OMP 杂波邻近敏感性') + plt.legend() + plt.tight_layout() + plt.savefig(os.path.join(SAVE_DIR, 'compare_clutter_sep_TSC_ULA_NULA.png'), dpi=300) + plt.show() + + print("局限性分析完成。") + return limitation_results + + +# ========================================================= +# 7. 主程序 +# ========================================================= +if __name__ == "__main__": + runtime_results, proxy_results = run_complexity_analysis() + limitation_results = run_limitation_analysis() + print("\n全部分析完成,结果保存在:", SAVE_DIR) \ No newline at end of file diff --git a/Graduation Design/RMSE_Test_ula.py b/Graduation Design/RMSE_Test_ula.py new file mode 100644 index 0000000..9cb50ac --- /dev/null +++ b/Graduation Design/RMSE_Test_ula.py @@ -0,0 +1,146 @@ +import torch +import numpy as np +import matplotlib.pyplot as plt +import time + +# ========================================== +# 0. 全局绘图参数与 CUDA 引擎初始化 +# ========================================== +plt.rcParams['axes.unicode_minus'] = False +plt.rcParams['font.sans-serif'] = ['Microsoft YaHei'] +plt.rcParams['mathtext.fontset'] = 'stix' + +device = torch.device("cuda" if torch.cuda.is_available() else "cpu") +torch.manual_seed(2024) + +N_ant, N_RF, d_lambda = 24, 4, 0.5 +clutter_angle = -15.0 +noise_power = 1.0 +INR_linear = 10**(30/10) + +def gen_a(theta_deg): + theta_rad = torch.deg2rad(torch.tensor(theta_deg, dtype=torch.float32, device=device)) + n_idx = torch.arange(-(N_ant-1)/2, (N_ant-1)/2 + 0.1, device=device) + return torch.exp(1j * 2 * torch.pi * d_lambda * n_idx * torch.sin(theta_rad)).unsqueeze(1) + +a_clutter = gen_a(clutter_angle) +R_interference = noise_power * torch.eye(N_ant, device=device) + INR_linear * (a_clutter @ a_clutter.mH) + +dic_angles_full = torch.arange(-90, 90.5, 0.5, device=device) +A_dic_full = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in dic_angles_full], dim=1) + +angle_ranges = [np.arange(0, 11), np.arange(11, 21), np.arange(21, 31), np.arange(31, 41)] +snr_test_list = [0, 5, 10, 15, 20] +Monte_Carlo = 100 +L_snapshots = 100 +colors = ['#0072BD', '#D95319', '#EDB120', '#7E2F8E', '#77AC30'] + +def run_omp(w_opt, A_dic): + res = w_opt.clone() + F_RF_list = [] + for k in range(N_RF): + idx = torch.argmax(torch.abs(A_dic.mH @ res)) + F_RF_list.append(A_dic[:, idx:idx+1]) + F_RF_mat = torch.cat(F_RF_list, dim=1) + res = w_opt - F_RF_mat @ (torch.linalg.pinv(F_RF_mat) @ w_opt) + return F_RF_mat + +print(f"启动 CUDA 并行引擎,对 4 大扇区进行极速扫雷测试 (MC={Monte_Carlo})...") +start_time = time.time() + +# 取消了这里原本的 fig, axs = plt.subplots(2, 2) +# 改为在循环内部动态创建单图 + +for fig_idx, test_angles in enumerate(angle_ranges): + print(f"正在计算扇区 {test_angles[0]}° ~ {test_angles[-1]}° ...") + + # 【修改点 1】:为每一个扇区单独生成一张独立画布 + fig, ax = plt.subplots(figsize=(9, 7)) + + RMSE_DFT_All = np.zeros((len(snr_test_list), len(test_angles))) + RMSE_FULL_All = np.zeros((len(snr_test_list), len(test_angles))) + RMSE_TSC_All = np.zeros((len(snr_test_list), len(test_angles))) + + for snr_idx, snr_dB in enumerate(snr_test_list): + sig_p = 10**(snr_dB/10) + + for a_idx, t_angle in enumerate(test_angles): + true_target_angle = t_angle + 0.04 + theta_prior = round(true_target_angle) + window = 5 + + dic_tsc = torch.arange(theta_prior - window, theta_prior + window + 0.1, 0.1, device=device) + A_dic_tsc = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in dic_tsc], dim=1) + + angles_dft = torch.linspace(theta_prior - window, theta_prior + window, N_RF, device=device) + F_RF_dft = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in angles_dft], dim=1) + + a_prior_broad = torch.zeros((N_ant, 1), dtype=torch.complex64, device=device) + for tb in torch.arange(theta_prior - window, theta_prior + window + 0.1, 0.5, device=device): + a_prior_broad += gen_a(tb) + a_prior_broad /= torch.norm(a_prior_broad) + + w_dig_prior = torch.linalg.inv(R_interference) @ a_prior_broad + w_dig_prior /= torch.norm(w_dig_prior) + + F_RF_full = run_omp(w_dig_prior, A_dic_full) + F_RF_tsc = run_omp(w_dig_prior, A_dic_tsc) + + scan_grid = torch.arange(theta_prior - window, theta_prior + window + 0.01, 0.01, device=device) + A_scan = torch.cat([gen_a(th) for th in scan_grid], dim=1) + + A_eff_dft = F_RF_dft.mH @ A_scan + A_eff_full = F_RF_full.mH @ A_scan + A_eff_tsc = F_RF_tsc.mH @ A_scan + + a_target = gen_a(true_target_angle) + s_t = np.sqrt(sig_p/2) * (torch.randn(Monte_Carlo, 1, L_snapshots, device=device) + 1j * torch.randn(Monte_Carlo, 1, L_snapshots, device=device)) + s_c = np.sqrt(INR_linear/2) * (torch.randn(Monte_Carlo, 1, L_snapshots, device=device) + 1j * torch.randn(Monte_Carlo, 1, L_snapshots, device=device)) + n_t = np.sqrt(noise_power/2) * (torch.randn(Monte_Carlo, N_ant, L_snapshots, device=device) + 1j * torch.randn(Monte_Carlo, N_ant, L_snapshots, device=device)) + + X = a_target.unsqueeze(0) * s_t + a_clutter.unsqueeze(0) * s_c + n_t + + def process_hybrid(F_RF, A_eff): + Y = F_RF.mH.unsqueeze(0) @ X + R = (Y @ Y.mH) / L_snapshots + R += (0.05 * torch.diagonal(R, dim1=-2, dim2=-1).sum(-1).view(-1, 1, 1) / N_RF) * torch.eye(N_RF, device=device).unsqueeze(0) + inv_R = torch.linalg.inv(R) + P = 1.0 / torch.abs(torch.einsum('si, bij, js -> bs', A_eff.conj().T, inv_R, A_eff)) + return np.sqrt(np.mean(((scan_grid[torch.argmax(P, dim=1)] - true_target_angle)**2).cpu().numpy())) + + RMSE_DFT_All[snr_idx, a_idx] = process_hybrid(F_RF_dft, A_eff_dft) + RMSE_FULL_All[snr_idx, a_idx] = process_hybrid(F_RF_full, A_eff_full) + RMSE_TSC_All[snr_idx, a_idx] = process_hybrid(F_RF_tsc, A_eff_tsc) + + c_str = colors[snr_idx] + lbl_suffix = f" (SNR={snr_dB}dB)" + ax.semilogy(test_angles, RMSE_DFT_All[snr_idx, :], '--^', color=c_str, linewidth=1.5, label='DFT'+lbl_suffix) + ax.semilogy(test_angles, RMSE_FULL_All[snr_idx, :], ':x', color=c_str, linewidth=2, label='传统OMP'+lbl_suffix) + ax.semilogy(test_angles, RMSE_TSC_All[snr_idx, :], '-o', color=c_str, linewidth=2.5, label='TSC-OMP'+lbl_suffix) + + ax.set_xlim([test_angles[0], test_angles[-1]]) + ax.set_xticks(test_angles) + ax.grid(True, which="both", ls="--", alpha=0.5) + ax.set_ylim([1e-3, 20]) + ax.set_xlabel('目标真实出现角度 $\\theta (\\circ)$', fontsize=12, fontweight='bold') + ax.set_ylabel(r'测角 RMSE $(^\circ)$', fontsize=12, fontweight='bold') + ax.set_title(f'宽波束先验引导下,扇区 ({test_angles[0]}°~{test_angles[-1]}°) 鲁棒性评估', fontsize=13, fontweight='bold') + + + # 【修改点 2】:独立设置每张图的图例。 + plt.tight_layout() + fig.subplots_adjust(bottom=0.25) + + # 【关键调整】:将 bbox_to_anchor 的 Y 坐标从 -0.15 提上来,改为 -0.08 + # 这样图例就会离 X 轴更近,避开底部操作条 + ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.10), ncol=3, fontsize=10) + + # 【修改点 3】:自动保存 + file_name = f'RMSE_Sector_{test_angles[0]}_{test_angles[-1]}.png' + plt.savefig(file_name, dpi=300, bbox_inches='tight') + print(f" 已保存独立图表: {file_name}") + +print(f"全扇区计算完毕!总耗时: {time.time() - start_time:.4f} 秒。") + +# 一次性在屏幕上弹出生成的4张独立图表 +plt.show() \ No newline at end of file diff --git a/Graduation Design/main_ula.py b/Graduation Design/main_ula.py new file mode 100644 index 0000000..e358ed2 --- /dev/null +++ b/Graduation Design/main_ula.py @@ -0,0 +1,280 @@ +import torch +import numpy as np +import matplotlib.pyplot as plt +import time + +# ========================================== +# 0. 全局设置与 CUDA 引擎 +# ========================================== +plt.rcParams['axes.unicode_minus'] = False +plt.rcParams['font.sans-serif'] = ['Microsoft YaHei'] +device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + +# 锁定随机数种子 +torch.manual_seed(2024) +np.random.seed(2024) + +# 颜色设置 +c_ana = '#77AC30' +c_dig = '#000000' +c_dft = '#4DBEEE' +c_omp = '#D95319' +c_tsc = '#0072BD' + +# ========================================== +# 1. 物理层参数初始化 (标准 ULA) +# ========================================== +N_ant, N_RF = 24, 4 #24根天线,4个射频链路 +#target_angle = float(input("请输入目标位置角度(°)(范围0-20°):")) #目标位置 +#clutter_angle = float(input("请输入干扰位置角度(°)(范围-3°~-20°):")) #干扰位置 +target_angle = 5.0 +clutter_angle = -15.0 +L_snapshots = 100 +noise_power = 1.0 +INR_linear = 10**(30/10) + +d_lambda = 0.5 # 经典的半波长间距 + +def gen_a(theta_deg): + """标准的均匀线阵 (ULA) 导向矢量""" + theta_rad = torch.deg2rad(torch.tensor(theta_deg, dtype=torch.float32, device=device)) + # 以阵列中心为参考点 + n_idx = torch.arange(-(N_ant-1)/2, (N_ant-1)/2 + 0.1, device=device) + return torch.exp(1j * 2 * torch.pi * d_lambda * n_idx * torch.sin(theta_rad)).unsqueeze(1) + +a_target = gen_a(target_angle) +a_clutter = gen_a(clutter_angle) + +# ========================================== +# 2. 动态生成目标约束扇区与字典 +# ========================================== +theta_coarse_prior = round(target_angle) +window_half_width = 5.0 + +dic_angles_constrained = torch.arange(theta_coarse_prior - window_half_width, theta_coarse_prior + window_half_width + 0.1, 0.1, device=device) +A_dic_constrained = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in dic_angles_constrained], dim=1) + +dic_angles_full = torch.arange(-90, 90.5, 0.5, device=device) +A_dic_full = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in dic_angles_full], dim=1) + +angles_dft = torch.linspace(theta_coarse_prior - window_half_width, theta_coarse_prior + window_half_width, N_RF, device=device) +F_RF_dft = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in angles_dft], dim=1) + +# ========================================== +# 3. 提取防发散的慢时间先验射频矩阵 F_RF +# ========================================== +print('正在进行慢时间 (Slow-Time) 模拟射频网络硬件配置...') +a_prior_broad = torch.zeros((N_ant, 1), dtype=torch.complex64, device=device) +for theta_b in torch.arange(theta_coarse_prior - window_half_width, theta_coarse_prior + window_half_width + 0.1, 0.5, device=device): + a_prior_broad += gen_a(theta_b) +a_prior_broad /= torch.norm(a_prior_broad) + +R_interference_prior = noise_power * torch.eye(N_ant, device=device) + INR_linear * (a_clutter @ a_clutter.mH) +w_digital_prior = torch.linalg.inv(R_interference_prior) @ a_prior_broad +w_digital_prior /= torch.norm(w_digital_prior) + +def extract_omp(w_opt, A_dic): + res = w_opt.clone() + F_RF_mat = None + for k in range(N_RF): + proj = A_dic.mH @ res + idx = torch.argmax(torch.abs(proj)) + new_atom = A_dic[:, idx:idx+1] + F_RF_mat = new_atom if F_RF_mat is None else torch.cat((F_RF_mat, new_atom), dim=1) + res = w_digital_prior - F_RF_mat @ (torch.linalg.pinv(F_RF_mat) @ w_digital_prior) + return F_RF_mat + +F_RF_omp_full_static = extract_omp(w_digital_prior, A_dic_full) +F_RF_omp_static = extract_omp(w_digital_prior, A_dic_constrained) + +# ========================================== +# 4. 计算波束方向图与理论 SINR +# ========================================== +print('正在计算静态波束方向图与理论 SINR...') +w_analog = a_target / np.sqrt(N_ant) + +def calc_hybrid_weight(F_RF, R_in_stat): + a_eff = F_RF.mH @ a_target + R_eff = F_RF.mH @ R_in_stat @ F_RF + R_eff += (0.05 * torch.trace(R_eff).real / N_RF) * torch.eye(N_RF, device=device) + inv_R_eff = torch.linalg.inv(R_eff) + F_BB = (inv_R_eff @ a_eff) / (a_eff.mH @ inv_R_eff @ a_eff) + w_hyb = F_RF @ F_BB + return w_hyb / torch.norm(w_hyb) + +sig_power_static = 10**(0/10) +# 直接使用无穷大快拍下的理论理想协方差矩阵,消除所有随机扰动,展示极限压制能力 +R_in_static = noise_power * torch.eye(N_ant, device=device) + \ + sig_power_static * (a_target @ a_target.mH) + \ + INR_linear * (a_clutter @ a_clutter.mH) + +w_hybrid_dft = calc_hybrid_weight(F_RF_dft, R_in_static) +w_hybrid_omp_full = calc_hybrid_weight(F_RF_omp_full_static, R_in_static) +w_hybrid_omp = calc_hybrid_weight(F_RF_omp_static, R_in_static) + +R_in_static_dl = R_in_static + (0.01 * torch.trace(R_in_static).real / N_ant) * torch.eye(N_ant, device=device) +inv_R_stat = torch.linalg.inv(R_in_static_dl) +w_digital_ideal = (inv_R_stat @ a_target) / (a_target.mH @ inv_R_stat @ a_target) +w_digital_ideal /= torch.norm(w_digital_ideal) + +theta_scan = torch.arange(-90, 90.1, 0.1, device=device) +A_scan_full = torch.cat([gen_a(th) for th in theta_scan], dim=1) + +def get_bp(w): return 10 * torch.log10((torch.abs(A_scan_full.mH @ w)**2).squeeze() / torch.max(torch.abs(A_scan_full.mH @ w)**2)) +BP_Ana = get_bp(w_analog).cpu().numpy() +BP_Dig = get_bp(w_digital_ideal).cpu().numpy() +BP_DFT = get_bp(w_hybrid_dft).cpu().numpy() +BP_OMP_FULL = get_bp(w_hybrid_omp_full).cpu().numpy() +BP_OMP = get_bp(w_hybrid_omp).cpu().numpy() +theta_scan_np = theta_scan.cpu().numpy() + +SNR_dB_range = np.arange(-10, 22, 2) +SINR_Ana, SINR_Dig, SINR_DFT, SINR_OMP_FULL, SINR_OMP = [], [], [], [], [] +for snr in SNR_dB_range: + sig_p = 10**(snr/10) + def calc_sinr(w): return (sig_p * torch.abs(w.mH @ a_target)**2) / torch.real(w.mH @ R_interference_prior @ w) + SINR_Ana.append(10 * torch.log10(calc_sinr(w_analog)).item()) + SINR_Dig.append(10 * torch.log10(calc_sinr(w_digital_ideal)).item()) + SINR_DFT.append(10 * torch.log10(calc_sinr(w_hybrid_dft)).item()) + SINR_OMP_FULL.append(10 * torch.log10(calc_sinr(w_hybrid_omp_full)).item()) + SINR_OMP.append(10 * torch.log10(calc_sinr(w_hybrid_omp)).item()) + +# ========================================== +# 5. 基于张量 Batch 操作的极致蒙特卡洛评估 +# ========================================== +print('正在启动 CUDA 高维张量批量计算 (Batched Evaluation)...') +Monte_Carlo = 500 +scan_grid = torch.arange(theta_coarse_prior - window_half_width, theta_coarse_prior + window_half_width + 0.01, 0.01, device=device) +A_scan = torch.cat([gen_a(th) for th in scan_grid], dim=1) + +A_eff_dft = F_RF_dft.mH @ A_scan +A_eff_omp_full = F_RF_omp_full_static.mH @ A_scan +A_eff_omp = F_RF_omp_static.mH @ A_scan + +RMSE_Ana, RMSE_Dig, RMSE_DFT, RMSE_OMP_FULL, RMSE_OMP = [], [], [], [], [] +start_time = time.time() + +for snr in SNR_dB_range: + sig_power = 10**(snr/10) + s_t = np.sqrt(sig_power/2) * (torch.randn(Monte_Carlo, 1, L_snapshots, device=device) + 1j*torch.randn(Monte_Carlo, 1, L_snapshots, device=device)) + s_c = np.sqrt(INR_linear/2) * (torch.randn(Monte_Carlo, 1, L_snapshots, device=device) + 1j*torch.randn(Monte_Carlo, 1, L_snapshots, device=device)) + n_t = np.sqrt(noise_power/2) * (torch.randn(Monte_Carlo, N_ant, L_snapshots, device=device) + 1j*torch.randn(Monte_Carlo, N_ant, L_snapshots, device=device)) + + X_mc = a_target.unsqueeze(0) * s_t + a_clutter.unsqueeze(0) * s_c + n_t + R_in = torch.bmm(X_mc, X_mc.mH) / L_snapshots + trace_R_in = torch.diagonal(R_in, dim1=-2, dim2=-1).sum(-1).real + R_in_dl = R_in + (0.05 * trace_R_in / N_ant).view(-1, 1, 1) * torch.eye(N_ant, device=device).unsqueeze(0) + inv_R_in_dl = torch.linalg.inv(R_in_dl) + + P_ana = torch.einsum('si, bij, js -> bs', A_scan.conj().T, R_in, A_scan).abs() + P_dig = torch.einsum('si, bij, js -> bs', A_scan.conj().T, inv_R_in_dl, A_scan).abs() + + def process_hybrid_batched(F_RF, A_eff_hyb): + Y = F_RF.mH.unsqueeze(0) @ X_mc + R_hyb = torch.bmm(Y, Y.mH) / L_snapshots + trace_R = torch.diagonal(R_hyb, dim1=-2, dim2=-1).sum(-1).real + R_hyb += (0.05 * trace_R / N_RF).view(-1, 1, 1) * torch.eye(N_RF, device=device).unsqueeze(0) + inv_R_hyb = torch.linalg.inv(R_hyb) + P_hyb = torch.einsum('si, bij, js -> bs', A_eff_hyb.conj().T, inv_R_hyb, A_eff_hyb).abs() + return scan_grid[torch.argmin(P_hyb, dim=1)] + + ang_ana = scan_grid[torch.argmax(P_ana, dim=1)] + ang_dig = scan_grid[torch.argmin(P_dig, dim=1)] + ang_dft = process_hybrid_batched(F_RF_dft, A_eff_dft) + ang_omp_full = process_hybrid_batched(F_RF_omp_full_static, A_eff_omp_full) + ang_omp = process_hybrid_batched(F_RF_omp_static, A_eff_omp) + + RMSE_Ana.append(torch.sqrt(torch.mean((ang_ana - target_angle)**2)).item()) + RMSE_Dig.append(torch.sqrt(torch.mean((ang_dig - target_angle)**2)).item()) + RMSE_DFT.append(torch.sqrt(torch.mean((ang_dft - target_angle)**2)).item()) + RMSE_OMP_FULL.append(torch.sqrt(torch.mean((ang_omp_full - target_angle)**2)).item()) + RMSE_OMP.append(torch.sqrt(torch.mean((ang_omp - target_angle)**2)).item()) + +print(f"CUDA 并行测角完毕!耗时: {time.time() - start_time:.4f} 秒。") + +# ========================================== +# 6. 绘图:单图输出与汇总输出 +# ========================================== +bp_data = [BP_Ana, BP_Dig, BP_DFT, BP_OMP_FULL, BP_OMP] +bp_labels = ['纯模拟架构', '纯数字架构(MVDR)', 'DFT混合架构', '传统OMP混合架构', 'TSC-OMP混合架构'] +bp_colors = [c_ana, c_dig, c_dft, c_omp, c_tsc] +bp_ls = ['-.', '-', '--', ':', '-'] + +# 6.1 分别绘制 5 张独立的波束方向图 +for i in range(5): + fig, ax = plt.subplots(figsize=(8, 6)) + ax.plot(theta_scan_np, bp_data[i], color=bp_colors[i], linestyle=bp_ls[i], linewidth=2.5, label='主波束') + ax.axvline(x=target_angle, color='g', linestyle='--', linewidth=2, label=f'机动目标 ({target_angle}°)') + ax.axvline(x=clutter_angle, color='r', linestyle='--', linewidth=2, label=f'强杂波 ({clutter_angle}°)') + ax.set_ylim([-60, 0]) + ax.set_xlim([-60, 60]) + ax.grid(True, linestyle='--', alpha=0.6) + # 标题更新为 ULA + ax.set_title(f'波束方向图 (ULA) - {bp_labels[i]}', fontsize=14, fontweight='bold') + ax.set_xlabel(r'角度 $\theta (^\circ)$', fontsize=12) + ax.set_ylabel('归一化增益 (dB)', fontsize=12) + + plt.tight_layout() + fig.subplots_adjust(bottom=0.22) + ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12), ncol=3, fontsize=10) + + file_name = f'BP_Single_{i+1}_ULA_{bp_labels[i].split("(")[0]}.png' + plt.savefig(file_name, dpi=300, bbox_inches='tight') + +# 6.2 绘制五大架构综合对比波束图 +fig_all, ax_all = plt.subplots(figsize=(10, 6)) +for i in range(5): + ax_all.plot(theta_scan_np, bp_data[i], color=bp_colors[i], linestyle=bp_ls[i], linewidth=2.5, label=bp_labels[i].split("(")[0]) +ax_all.axvline(x=target_angle, color='g', linestyle='--', linewidth=2) +ax_all.axvline(x=clutter_angle, color='r', linestyle='--', linewidth=2) +ax_all.text(target_angle + 1, -5, f'机动目标 ({target_angle}°)', color='g', fontweight='bold') +ax_all.text(clutter_angle - 1, -5, f'强杂波 ({clutter_angle}°)', color='r', fontweight='bold', ha='right') +ax_all.set_ylim([-60, 0]) +ax_all.set_xlim([-60, 60]) +ax_all.grid(True, linestyle='--', alpha=0.6) +ax_all.set_title('波束方向图综合对比 (ULA)', fontsize=14, fontweight='bold') +ax_all.set_xlabel(r'角度 $\theta (^\circ)$', fontsize=12) +ax_all.set_ylabel('归一化增益 (dB)', fontsize=12) + +plt.tight_layout() +fig_all.subplots_adjust(bottom=0.22) +ax_all.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12), ncol=5, fontsize=10) +plt.savefig('BP_All_Combined_ULA.png', dpi=300, bbox_inches='tight') + +# 6.3 绘制 SINR 对比图 +fig_sinr, ax_sinr = plt.subplots(figsize=(10, 6)) +ax_sinr.plot(SNR_dB_range, SINR_Ana, '-v', color=c_ana, linewidth=2, label='纯模拟') +ax_sinr.plot(SNR_dB_range, SINR_Dig, '-o', color=c_dig, linewidth=2, label='纯数字(上限)') +ax_sinr.plot(SNR_dB_range, SINR_DFT, '-d', color=c_dft, linewidth=2, label='DFT') +ax_sinr.plot(SNR_dB_range, SINR_OMP_FULL, '-*', color=c_omp, linewidth=2, label='传统OMP') +ax_sinr.plot(SNR_dB_range, SINR_OMP, '-s', color=c_tsc, linewidth=2.5, label='TSC-OMP(所提)') +ax_sinr.grid(True, linestyle='--', alpha=0.6) +ax_sinr.set_title(f'全架构输出 SINR 理论逼近测试 (ULA - 目标 {target_angle}°)', fontsize=14, fontweight='bold') +ax_sinr.set_xlabel('输入 SNR (dB)', fontsize=12) +ax_sinr.set_ylabel('输出 SINR (dB)', fontsize=12) + +plt.tight_layout() +fig_sinr.subplots_adjust(bottom=0.22) +ax_sinr.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12), ncol=5, fontsize=10) +plt.savefig('SINR_All_Combined_ULA.png', dpi=300, bbox_inches='tight') + +# 6.4 绘制 RMSE 对比图 +fig_rmse, ax_rmse = plt.subplots(figsize=(10, 6)) +ax_rmse.semilogy(SNR_dB_range, RMSE_Ana, '-.v', color=c_ana, linewidth=2, label='纯模拟') +ax_rmse.semilogy(SNR_dB_range, RMSE_Dig, '-o', color=c_dig, linewidth=2, label='纯数字(上限)') +ax_rmse.semilogy(SNR_dB_range, RMSE_DFT, '--d', color=c_dft, linewidth=2, label='DFT') +ax_rmse.semilogy(SNR_dB_range, RMSE_OMP_FULL, ':*', color=c_omp, linewidth=2.5, label='传统OMP') +ax_rmse.semilogy(SNR_dB_range, RMSE_OMP, '-s', color=c_tsc, linewidth=2.5, label='TSC-OMP (所提)') +ax_rmse.grid(True, which="both", ls="--", alpha=0.6) +ax_rmse.set_ylim([1e-3, 20]) +ax_rmse.set_title(f'全架构实战恶劣环境 DOA 测角精度评估 (ULA)', fontsize=14, fontweight='bold') +ax_rmse.set_xlabel('信噪比 SNR (dB)', fontsize=12) +ax_rmse.set_ylabel(r'测角 RMSE $(^\circ)$', fontsize=12) + +plt.tight_layout() +fig_rmse.subplots_adjust(bottom=0.22) +ax_rmse.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12), ncol=5, fontsize=10) +plt.savefig('RMSE_All_Combined_ULA.png', dpi=300, bbox_inches='tight') + +plt.show() +print('计算完毕,单图与综合图表均已独立生成。') \ No newline at end of file