diff --git a/Code/Python/3.Robust and complexity test/Complexity_Test.py b/Code/Python/3.Robust and complexity test/Complexity_Test.py new file mode 100644 index 0000000..53facee --- /dev/null +++ b/Code/Python/3.Robust and complexity test/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/Code/Python/3.Robust and complexity test/robust_tsc_omp_cuda_showall.py b/Code/Python/3.Robust and complexity test/robust_tsc_omp_cuda_showall.py new file mode 100644 index 0000000..197cbed --- /dev/null +++ b/Code/Python/3.Robust and complexity test/robust_tsc_omp_cuda_showall.py @@ -0,0 +1,402 @@ + +import os +import math +import time +import numpy as np +import torch +import matplotlib.pyplot as plt + +# ============================================================ +# 0. 全局设置 +# ============================================================ +plt.rcParams['font.sans-serif'] = [ + 'Microsoft YaHei', 'SimHei', 'PingFang SC', + 'Noto Sans CJK SC', 'Arial Unicode MS', 'DejaVu Sans' +] +plt.rcParams['axes.unicode_minus'] = False +plt.rcParams['mathtext.fontset'] = 'stix' + +SAVE_DIR = "robust_tsc_omp_cuda_results" +os.makedirs(SAVE_DIR, exist_ok=True) + +device = torch.device("cuda" if torch.cuda.is_available() else "cpu") +torch.manual_seed(2024) +np.random.seed(2024) + +# ============================================================ +# 1. 参数区 +# 说明： +# - 本脚本默认使用 CUDA（若可用） +# - 为避免“卡死”，默认采用较稳妥的 Monte Carlo 与扫描步长 +# - 所有图都会生成并在脚本结束时统一以窗口形式显示 +# ============================================================ +ARRAY_TYPE = 'NULA' # 'ULA' 或 'NULA' +N_ant = 24 +N_RF = 4 +d_lambda = 0.5 + +target_angle = 5.0 +clutter_angle = -15.0 +noise_power = 1.0 +INR_dB = 30.0 +INR_linear = 10 ** (INR_dB / 10) + +L_snapshots = 100 +Monte_Carlo = 300 +SNR_dB_range = np.arange(-10, 22, 2) + +theta_prior = 5.0 +window_half_width = 5.0 +sector_min = theta_prior - window_half_width +sector_max = theta_prior + window_half_width + +# TSC-OMP 字典与测角搜索栅格 +tsc_grid_np = np.arange(sector_min, sector_max + 0.1, 0.1) +scan_grid_np = np.arange(sector_min, sector_max + 0.001, 0.05) # 比 0.01 更快 + +# 对角加载系数 +alpha_list = [0.00, 0.01, 0.03, 0.05, 0.10] + +# 是否保存图片 +SAVE_FIG = True + +# ============================================================ +# 2. 阵列位置 +# ============================================================ +def get_array_positions(array_type='ULA', N=24, d=0.5): + if array_type.upper() == 'ULA': + pos = np.arange(-(N - 1) / 2, (N - 1) / 2 + 1) * d + return pos.astype(np.float32) + elif array_type.upper() == 'NULA': + # 与论文一致的示意性 NULA + base = 0.7 + jump = 1.03 + pos = [] + x = 0.0 + for _ in range(12): + pos.append(x) + x += base + x += (jump - base) + for _ in range(12, 24): + pos.append(x) + x += base + pos = np.array(pos, dtype=np.float32) + pos = pos - np.mean(pos) + return pos + else: + raise ValueError("ARRAY_TYPE 只能是 'ULA' 或 'NULA'") + +array_pos_np = get_array_positions(ARRAY_TYPE, N_ant, d_lambda) +array_pos_t = torch.tensor(array_pos_np, dtype=torch.float32, device=device).view(-1, 1) + +# ============================================================ +# 3. 基础函数（CUDA） +# ============================================================ +def sync(): + if device.type == 'cuda': + torch.cuda.synchronize() + +def steering_vector(theta_deg): + theta = torch.as_tensor(theta_deg, dtype=torch.float32, device=device) + theta_rad = torch.deg2rad(theta) + phase = 2 * torch.pi * array_pos_t * torch.sin(theta_rad) + return torch.exp(1j * phase) + +def steering_matrix(theta_deg_array): + theta = torch.as_tensor(theta_deg_array, dtype=torch.float32, device=device).view(1, -1) + theta_rad = torch.deg2rad(theta) + phase = 2 * torch.pi * array_pos_t @ torch.sin(theta_rad) + return torch.exp(1j * phase) + +def build_dictionary(angle_grid_np, normalize=True): + A = steering_matrix(angle_grid_np) + if normalize: + A = A / math.sqrt(A.shape[0]) + return A + +def build_prior_broad(sector_min, sector_max, step=0.5): + angles = np.arange(sector_min, sector_max + 1e-12, step, dtype=np.float32) + A = steering_matrix(angles) + a = torch.sum(A, dim=1, keepdim=True) + a = a / torch.linalg.norm(a) + return a + +def extract_omp(w_opt, A_dic, angle_grid_np, N_RF): + """ + CUDA 版 OMP / TSC-OMP + """ + res = w_opt.clone() + F_RF = None + chosen_angles = [] + + for _ in range(N_RF): + proj = torch.abs(A_dic.mH @ res).flatten() + idx = int(torch.argmax(proj).item()) + atom = A_dic[:, idx:idx+1] + F_RF = atom if F_RF is None else torch.cat((F_RF, atom), dim=1) + chosen_angles.append(float(angle_grid_np[idx])) + res = w_opt - F_RF @ (torch.linalg.pinv(F_RF) @ w_opt) + + return F_RF, chosen_angles + +def build_prior_weight(target_angle, clutter_angle, sector_min, sector_max): + a_clutter = steering_vector(clutter_angle) + R_interf = noise_power * torch.eye(N_ant, dtype=torch.complex64, device=device) \ + + INR_linear * (a_clutter @ a_clutter.mH) + a_prior = build_prior_broad(sector_min, sector_max, step=0.5) + w_prior = torch.linalg.solve(R_interf, a_prior) + w_prior = w_prior / torch.linalg.norm(w_prior) + return w_prior + +def calc_hybrid_weight_loaded(F_RF, R_in, a_target, alpha_dl=0.05): + a_eff = F_RF.mH @ a_target + R_eff = F_RF.mH @ R_in @ F_RF + gamma = alpha_dl * torch.real(torch.trace(R_eff)) / R_eff.shape[0] + R_eff_loaded = R_eff + gamma * torch.eye(R_eff.shape[0], dtype=torch.complex64, device=device) + inv_R = torch.linalg.inv(R_eff_loaded) + F_BB = (inv_R @ a_eff) / (a_eff.mH @ inv_R @ a_eff) + w = F_RF @ F_BB + w = w / torch.linalg.norm(w) + cond_number = torch.linalg.cond(R_eff_loaded).real.item() + return w, R_eff_loaded, cond_number + +def robust_capon_spectrum_batch(F_RF, X_batch, A_scan, alpha_dl=0.05): + """ + X_batch: [MC, N_ant, L] + A_scan: [N_ant, Ns] + 返回： + P: [MC, Ns] + mean_cond: 平均条件数 + """ + Y = torch.matmul(F_RF.mH.unsqueeze(0), X_batch) # [MC, N_RF, L] + R_bb = torch.matmul(Y, Y.mH) / X_batch.shape[-1] # [MC, N_RF, N_RF] + + trace_R = torch.real(torch.diagonal(R_bb, dim1=-2, dim2=-1).sum(-1)) + gamma = alpha_dl * trace_R / F_RF.shape[1] + eye = torch.eye(F_RF.shape[1], dtype=torch.complex64, device=device).unsqueeze(0) + R_bb_loaded = R_bb + gamma[:, None, None] * eye + + inv_R = torch.linalg.inv(R_bb_loaded) # [MC, N_RF, N_RF] + A_eff = F_RF.mH @ A_scan # [N_RF, Ns] + + # P = 1 / |a^H invR a| + denom = torch.einsum('rn,brs,sn->bn', A_eff.conj(), inv_R, A_eff) + P = 1.0 / torch.clamp(torch.abs(denom), min=1e-12) + + cond_vals = torch.linalg.cond(R_bb_loaded).real + mean_cond = torch.mean(cond_vals).item() + return P, mean_cond + +def calc_sinr(w, a_target, a_clutter, sig_p): + R_interf = noise_power * torch.eye(N_ant, dtype=torch.complex64, device=device) \ + + INR_linear * (a_clutter @ a_clutter.mH) + num = sig_p * (torch.abs(w.mH @ a_target).item() ** 2) + den = torch.real(w.mH @ R_interf @ w).item() + return 10 * np.log10(max(num / max(den, 1e-12), 1e-12)) + +# ============================================================ +# 4. 固定 TSC-OMP 模拟矩阵 +# ============================================================ +print(f"当前设备: {device}") +sync() +tic = time.perf_counter() + +w_prior = build_prior_weight( + target_angle=target_angle, + clutter_angle=clutter_angle, + sector_min=sector_min, + sector_max=sector_max +) +A_dic_tsc = build_dictionary(tsc_grid_np, normalize=True) +F_RF_tsc, chosen_angles = extract_omp(w_prior, A_dic_tsc, tsc_grid_np, N_RF) + +a_target = steering_vector(target_angle) +a_clutter = steering_vector(clutter_angle) + +sync() +print(f"TSC-OMP 选中的模拟原子角度: {chosen_angles}") +print(f"固定模拟矩阵构建耗时: {time.perf_counter() - tic:.3f} s") + +# ============================================================ +# 5. 静态方向图：不同加载系数 +# ============================================================ +sig_power_static = 10 ** (0 / 10.0) +s_t = np.sqrt(sig_power_static / 2) * ( + torch.randn(1, L_snapshots, dtype=torch.float32, device=device) + + 1j * torch.randn(1, L_snapshots, dtype=torch.float32, device=device) +) +s_c = np.sqrt(INR_linear / 2) * ( + torch.randn(1, L_snapshots, dtype=torch.float32, device=device) + + 1j * torch.randn(1, L_snapshots, dtype=torch.float32, device=device) +) +n_t = np.sqrt(noise_power / 2) * ( + torch.randn(N_ant, L_snapshots, dtype=torch.float32, device=device) + + 1j * torch.randn(N_ant, L_snapshots, dtype=torch.float32, device=device) +) +X_static = a_target @ s_t + a_clutter @ s_c + n_t +R_in_static = (X_static @ X_static.mH) / L_snapshots + +theta_scan_np = np.arange(-60, 60.1, 0.1) +A_bp = steering_matrix(theta_scan_np) + +plt.figure(figsize=(10, 6)) +for alpha_dl in alpha_list: + w_hyb, _, _ = calc_hybrid_weight_loaded(F_RF_tsc, R_in_static, a_target, alpha_dl=alpha_dl) + BP = torch.abs(w_hyb.mH @ A_bp).flatten() ** 2 + BP = 10 * torch.log10(BP / torch.max(BP) + 1e-12) + plt.plot(theta_scan_np, BP.detach().cpu().numpy(), linewidth=2, label=f'α={alpha_dl:.2f}') + +plt.axvline(target_angle, color='g', linestyle='--', label=f'目标 {target_angle:.2f}°') +plt.axvline(clutter_angle, color='r', linestyle='--', label=f'杂波 {clutter_angle:.2f}°') +plt.xlabel('角度 θ (°)') +plt.ylabel('归一化增益 (dB)') +plt.title(f'不同对角加载系数下 TSC-OMP 波束方向图 ({ARRAY_TYPE})') +plt.grid(True, linestyle='--', alpha=0.5) +plt.legend() +plt.tight_layout() +if SAVE_FIG: + plt.savefig(os.path.join(SAVE_DIR, f'BP_loaded_compare_{ARRAY_TYPE}.png'), dpi=300) + +# ============================================================ +# 6. 输出 SINR / RMSE / 条件数 对比 +# ============================================================ +SINR_loaded = {alpha: [] for alpha in alpha_list} +RMSE_loaded = {alpha: [] for alpha in alpha_list} +COND_loaded = {alpha: [] for alpha in alpha_list} + +A_scan = steering_matrix(scan_grid_np) + +for snr_db in SNR_dB_range: + sync() + tic = time.perf_counter() + + sig_p = 10 ** (snr_db / 10.0) + + # 一次性批量生成 Monte Carlo 数据 + s_t = np.sqrt(sig_p / 2) * ( + torch.randn(Monte_Carlo, 1, L_snapshots, dtype=torch.float32, device=device) + + 1j * torch.randn(Monte_Carlo, 1, L_snapshots, dtype=torch.float32, device=device) + ) + s_c = np.sqrt(INR_linear / 2) * ( + torch.randn(Monte_Carlo, 1, L_snapshots, dtype=torch.float32, device=device) + + 1j * torch.randn(Monte_Carlo, 1, L_snapshots, dtype=torch.float32, device=device) + ) + n_t = np.sqrt(noise_power / 2) * ( + torch.randn(Monte_Carlo, N_ant, L_snapshots, dtype=torch.float32, device=device) + + 1j * torch.randn(Monte_Carlo, N_ant, L_snapshots, dtype=torch.float32, device=device) + ) + X_batch = a_target.unsqueeze(0) * s_t + a_clutter.unsqueeze(0) * s_c + n_t + R_in_batch = torch.matmul(X_batch, X_batch.mH) / L_snapshots + + for alpha_dl in alpha_list: + # --- SINR：使用每个 Monte Carlo 样本求权，最后取平均 --- + sinr_vals = [] + cond_eff_vals = [] + + for mc in range(Monte_Carlo): + w_hyb, _, cond_eff = calc_hybrid_weight_loaded( + F_RF_tsc, R_in_batch[mc], a_target, alpha_dl=alpha_dl + ) + sinr_vals.append(calc_sinr(w_hyb, a_target, a_clutter, sig_p)) + cond_eff_vals.append(cond_eff) + + SINR_loaded[alpha_dl].append(float(np.mean(sinr_vals))) + + # --- RMSE 与 R_bb 条件数：批量 Capon 谱 --- + P, mean_cond_bb = robust_capon_spectrum_batch( + F_RF_tsc, X_batch, A_scan, alpha_dl=alpha_dl + ) + est_idx = torch.argmax(P, dim=1) + est_angles = torch.as_tensor(scan_grid_np, dtype=torch.float32, device=device)[est_idx] + rmse = torch.sqrt(torch.mean((est_angles - target_angle) ** 2)).item() + + RMSE_loaded[alpha_dl].append(rmse) + COND_loaded[alpha_dl].append(mean_cond_bb) + + sync() + print(f"SNR={snr_db:>3} dB 完成，耗时 {time.perf_counter() - tic:.2f} s") + +# ============================================================ +# 7. 输出 SINR 图 +# ============================================================ +plt.figure(figsize=(10, 6)) +for alpha_dl in alpha_list: + plt.plot(SNR_dB_range, SINR_loaded[alpha_dl], marker='o', linewidth=2, label=f'α={alpha_dl:.2f}') +plt.xlabel('输入 SNR (dB)') +plt.ylabel('输出 SINR (dB)') +plt.title(f'不同对角加载系数下 TSC-OMP 输出 SINR 对比 ({ARRAY_TYPE})') +plt.grid(True, linestyle='--', alpha=0.5) +plt.legend() +plt.tight_layout() +if SAVE_FIG: + plt.savefig(os.path.join(SAVE_DIR, f'SINR_loaded_compare_{ARRAY_TYPE}.png'), dpi=300) + +# ============================================================ +# 8. RMSE 图 +# ============================================================ +plt.figure(figsize=(10, 6)) +for alpha_dl in alpha_list: + plt.semilogy(SNR_dB_range, RMSE_loaded[alpha_dl], marker='s', linewidth=2, label=f'α={alpha_dl:.2f}') +plt.xlabel('输入 SNR (dB)') +plt.ylabel('测角 RMSE (°)') +plt.title(f'不同对角加载系数下 TSC-OMP 测角 RMSE 对比 ({ARRAY_TYPE})') +plt.grid(True, which='both', linestyle='--', alpha=0.5) +plt.legend() +plt.tight_layout() +if SAVE_FIG: + plt.savefig(os.path.join(SAVE_DIR, f'RMSE_loaded_compare_{ARRAY_TYPE}.png'), dpi=300) + +# ============================================================ +# 9. 条件数图 +# ============================================================ +plt.figure(figsize=(10, 6)) +for alpha_dl in alpha_list: + plt.semilogy(SNR_dB_range, COND_loaded[alpha_dl], marker='d', linewidth=2, label=f'α={alpha_dl:.2f}') +plt.xlabel('输入 SNR (dB)') +plt.ylabel('低维协方差矩阵条件数') +plt.title(f'不同对角加载系数下低维协方差矩阵条件数 ({ARRAY_TYPE})') +plt.grid(True, which='both', linestyle='--', alpha=0.5) +plt.legend() +plt.tight_layout() +if SAVE_FIG: + plt.savefig(os.path.join(SAVE_DIR, f'COND_loaded_compare_{ARRAY_TYPE}.png'), dpi=300) + +# ============================================================ +# 10. 低SNR点 α 扫描汇总图 +# ============================================================ +low_snr_pick = 0 +alpha_arr = np.array(alpha_list, dtype=float) +rmse_pick = np.array([RMSE_loaded[a][np.where(SNR_dB_range == low_snr_pick)[0][0]] for a in alpha_list], dtype=float) +cond_pick = np.array([COND_loaded[a][np.where(SNR_dB_range == low_snr_pick)[0][0]] for a in alpha_list], dtype=float) + +fig = plt.figure(figsize=(10, 6)) +ax1 = fig.add_subplot(111) +ax1.plot(alpha_arr, rmse_pick, marker='o', linewidth=2, color='#0072BD', label='RMSE') +ax1.set_xlabel('对角加载系数 α') +ax1.set_ylabel('RMSE (°)', color='#0072BD') +ax1.tick_params(axis='y', labelcolor='#0072BD') +ax1.grid(True, linestyle='--', alpha=0.5) + +ax2 = ax1.twinx() +ax2.semilogy(alpha_arr, cond_pick, marker='s', linewidth=2, color='#D95319', label='条件数') +ax2.set_ylabel('条件数', color='#D95319') +ax2.tick_params(axis='y', labelcolor='#D95319') + +plt.title(f'低SNR={low_snr_pick} dB 下加载系数对 RMSE 与条件数的影响 ({ARRAY_TYPE})') +fig.tight_layout() +if SAVE_FIG: + plt.savefig(os.path.join(SAVE_DIR, f'alpha_sweep_lowSNR_{ARRAY_TYPE}.png'), dpi=300) + +# ============================================================ +# 11. 打印总结 +# ============================================================ +print("\n===== 鲁棒性提升实验完成 =====") +for alpha_dl in alpha_list: + print(f"α={alpha_dl:.2f} : " + f"平均SINR={np.mean(SINR_loaded[alpha_dl]):.3f} dB, " + f"平均RMSE={np.mean(RMSE_loaded[alpha_dl]):.4f}°, " + f"平均条件数={np.mean(COND_loaded[alpha_dl]):.3e}") + +print(f"\n结果图片已保存到: {SAVE_DIR}") +print("即将以窗口形式显示所有图片……") +plt.show() diff --git a/Code/Python/3.Robust and complexity test/verify_pseudo_peak_omp_tsc.py b/Code/Python/3.Robust and complexity test/verify_pseudo_peak_omp_tsc.py new file mode 100644 index 0000000..198b42f --- /dev/null +++ b/Code/Python/3.Robust and complexity test/verify_pseudo_peak_omp_tsc.py @@ -0,0 +1,501 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle +import os + +# ============================================================ +# 0. 全局绘图设置（尽量兼容中文） +# ============================================================ +plt.rcParams['font.sans-serif'] = [ + 'Microsoft YaHei', 'SimHei', 'PingFang SC', + 'Noto Sans CJK SC', 'Arial Unicode MS', 'DejaVu Sans' +] +plt.rcParams['axes.unicode_minus'] = False + +# ============================================================ +# 1. 用户参数区（你可以根据论文需要自行调整） +# ============================================================ +ARRAY_TYPE = 'ULA' # 'ULA' 或 'NULA' +N_ant = 24 # 阵元数 +N_RF = 4 # 射频链数 +d_lambda = 0.5 # 阵元间距（单位：波长），ULA 用 +target_angle = 5.0 # 目标角度 +clutter_angle = -15.0 # 强杂波角度 +sector_min = 0.0 # TSC-OMP 扇区下界 +sector_max = 10.0 # TSC-OMP 扇区上界 +clutter_band_half = 2.0 # 杂波邻域半宽（用于统计“杂波误选率”） +noise_power = 1.0 +SNR_dB = -15.0 # 故意设置得较低，便于显现伪峰 +INR_dB = 35.0 # 强干扰 +L_snapshots = 48 # 快拍数，故意设置得不高，便于显现伪峰 +diag_loading_ratio = 0.01 # 对角加载比例 +MC = 400 # Monte Carlo 次数 + +# 字典设置 +full_grid = np.arange(-90.0, 90.0 + 0.5, 0.5) # 传统OMP：全空间字典 +tsc_grid = np.arange(sector_min, sector_max + 0.1, 0.1) # TSC-OMP：局部字典 + +# 宽波束先验 +broad_step = 0.5 + +# 输出文件夹 +save_dir = 'pseudo_peak_results' +os.makedirs(save_dir, exist_ok=True) + + +# ============================================================ +# 2. 阵列位置定义 +# ============================================================ +def get_array_positions(array_type='ULA', N=24, d=0.5): + """ + 返回阵列位置（单位：波长） + """ + if array_type.upper() == 'ULA': + pos = np.arange(-(N - 1) / 2, (N - 1) / 2 + 1) * d + return pos.astype(float) + + elif array_type.upper() == 'NULA': + # 一个示例性的 NULA 阵列位置（非均匀、递增排序、近似关于0分布） + # 你也可以改成与你论文代码完全一致的 NULA 位置 + pos = np.array([ + -5.75, -5.10, -4.40, -3.85, -3.15, -2.60, + -2.05, -1.45, -0.95, -0.35, 0.20, 0.75, + 1.25, 1.85, 2.45, 3.05, 3.70, 4.30, + 4.95, 5.65, 6.35, 7.05, 7.85, 8.70 + ]) + pos = pos - np.mean(pos) + return pos.astype(float) + + else: + raise ValueError("ARRAY_TYPE 只能是 'ULA' 或 'NULA'") + + +array_pos = get_array_positions(ARRAY_TYPE, N_ant, d_lambda) + + +# ============================================================ +# 3. 基本函数 +# ============================================================ +def steering_vector(theta_deg, array_pos): + """ + 阵列导向矢量 a(theta) + theta_deg : 角度（度） + array_pos : 阵元位置（单位：波长） + """ + theta_rad = np.deg2rad(theta_deg) + phase = 2 * np.pi * array_pos * np.sin(theta_rad) + a = np.exp(1j * phase).reshape(-1, 1) + return a + + +def build_dictionary(angle_grid, array_pos, normalize=True): + """ + 构造字典矩阵 A = [a(theta1), a(theta2), ...] + """ + A = np.hstack([steering_vector(ang, array_pos) for ang in angle_grid]) + if normalize: + A = A / np.sqrt(A.shape[0]) + return A + + +def omp_trace(w_opt, A_dic, angle_grid, N_RF): + """ + 传统 OMP / TSC-OMP 的公共实现：只是在字典上不同 + + 返回： + F_RF : 最终模拟矩阵 + chosen_ang : 每一步选中的角度 + proj_hist : 每一步投影谱 |A^H r| + """ + res = w_opt.copy() + F_RF = np.zeros((A_dic.shape[0], 0), dtype=np.complex128) + chosen_ang = [] + proj_hist = [] + + for _ in range(N_RF): + proj = np.abs(A_dic.conj().T @ res).flatten() + proj_hist.append(proj.copy()) + + idx = np.argmax(proj) + chosen_ang.append(float(angle_grid[idx])) + + atom = A_dic[:, idx:idx+1] + F_RF = np.hstack((F_RF, atom)) + + coeff = np.linalg.pinv(F_RF) @ w_opt + res = w_opt - F_RF @ coeff + + return F_RF, np.array(chosen_ang), proj_hist + + +def build_broad_prior_vector(array_pos, sector_min, sector_max, step=0.5): + """ + 宽波束先验导向向量：对目标扇区做叠加 + """ + a_broad = np.zeros((len(array_pos), 1), dtype=np.complex128) + for th in np.arange(sector_min, sector_max + 1e-12, step): + a_broad += steering_vector(th, array_pos) + a_broad /= np.linalg.norm(a_broad) + return a_broad + + +def build_sample_prior(array_pos, + target_angle, + clutter_angle, + SNR_dB, + INR_dB, + noise_power, + L_snapshots, + diag_loading_ratio, + sector_min, + sector_max, + broad_step=0.5): + """ + 构造有限快拍样本协方差下的先验数字权值 w_prior + 这是让传统 OMP 更容易显现伪峰风险的关键 + """ + N_ant = len(array_pos) + sig_p = 10 ** (SNR_dB / 10.0) + INR_linear = 10 ** (INR_dB / 10.0) + + a_t = steering_vector(target_angle, array_pos) + a_c = steering_vector(clutter_angle, array_pos) + + s_t = np.sqrt(sig_p / 2) * ( + np.random.randn(1, L_snapshots) + 1j * np.random.randn(1, L_snapshots) + ) + s_c = np.sqrt(INR_linear / 2) * ( + np.random.randn(1, L_snapshots) + 1j * np.random.randn(1, L_snapshots) + ) + n_t = np.sqrt(noise_power / 2) * ( + np.random.randn(N_ant, L_snapshots) + 1j * np.random.randn(N_ant, L_snapshots) + ) + + X = a_t @ s_t + a_c @ s_c + n_t + R_hat = (X @ X.conj().T) / L_snapshots + + dl = diag_loading_ratio * np.real(np.trace(R_hat)) / N_ant + R_hat = R_hat + dl * np.eye(N_ant) + + a_broad = build_broad_prior_vector(array_pos, sector_min, sector_max, broad_step) + + # 用低快拍样本协方差构造“先验数字权值” + w_prior = np.linalg.solve(R_hat, a_broad) + w_prior /= np.linalg.norm(w_prior) + + return w_prior + + +def count_stats(chosen_angles, sector_min, sector_max, clutter_angle, clutter_band_half): + """ + 统计： + P_hit : 目标扇区命中率 + P_clutter : 杂波误选率 + P_out : 扇区外误选率 + """ + chosen_angles = np.array(chosen_angles) + + hit = np.sum((chosen_angles >= sector_min) & (chosen_angles <= sector_max)) + clutter = np.sum((chosen_angles >= clutter_angle - clutter_band_half) & + (chosen_angles <= clutter_angle + clutter_band_half)) + out = np.sum((chosen_angles < sector_min) | (chosen_angles > sector_max)) + + return hit, clutter, out + + +def compute_badness_score(proj, angle_grid, sector_min, sector_max, clutter_angle, clutter_band_half): + """ + 用于自动挑选一个“传统OMP伪峰最明显”的代表性实验样本 + """ + proj = np.array(proj).flatten() + + target_mask = (angle_grid >= sector_min) & (angle_grid <= sector_max) + clutter_mask = (angle_grid >= clutter_angle - clutter_band_half) & (angle_grid <= clutter_angle + clutter_band_half) + out_mask = ~target_mask + + target_peak = np.max(proj[target_mask]) if np.any(target_mask) else 1e-12 + clutter_peak = np.max(proj[clutter_mask]) if np.any(clutter_mask) else 0.0 + out_peak = np.max(proj[out_mask]) if np.any(out_mask) else 0.0 + + score = max(clutter_peak / (target_peak + 1e-12), + out_peak / (target_peak + 1e-12)) + return score + + +# ============================================================ +# 4. 构造字典 +# ============================================================ +A_dic_full = build_dictionary(full_grid, array_pos, normalize=True) +A_dic_tsc = build_dictionary(tsc_grid, array_pos, normalize=True) + +# ============================================================ +# 5. Monte Carlo 验证 +# ============================================================ +all_angles_full = [] +all_angles_tsc = [] + +full_hit_total = 0 +full_clutter_total = 0 +full_out_total = 0 + +tsc_hit_total = 0 +tsc_clutter_total = 0 +tsc_out_total = 0 + +# 保存一个“代表性最差样本”，用于画第一步投影谱 +best_bad_score = -1.0 +rep_full_proj0 = None +rep_tsc_proj0 = None +rep_full_chosen = None +rep_tsc_chosen = None + +for mc in range(MC): + w_prior = build_sample_prior( + array_pos=array_pos, + target_angle=target_angle, + clutter_angle=clutter_angle, + SNR_dB=SNR_dB, + INR_dB=INR_dB, + noise_power=noise_power, + L_snapshots=L_snapshots, + diag_loading_ratio=diag_loading_ratio, + sector_min=sector_min, + sector_max=sector_max, + broad_step=broad_step + ) + + # 传统OMP（全空间字典） + _, chosen_full, proj_hist_full = omp_trace(w_prior, A_dic_full, full_grid, N_RF) + + # TSC-OMP（目标扇区局部字典） + _, chosen_tsc, proj_hist_tsc = omp_trace(w_prior, A_dic_tsc, tsc_grid, N_RF) + + all_angles_full.extend(chosen_full.tolist()) + all_angles_tsc.extend(chosen_tsc.tolist()) + + h, c, o = count_stats(chosen_full, sector_min, sector_max, clutter_angle, clutter_band_half) + full_hit_total += h + full_clutter_total += c + full_out_total += o + + h, c, o = count_stats(chosen_tsc, sector_min, sector_max, clutter_angle, clutter_band_half) + tsc_hit_total += h + tsc_clutter_total += c + tsc_out_total += o + + # 选一个“传统OMP最容易看到伪峰”的样本 + cur_score = compute_badness_score(proj_hist_full[0], full_grid, + sector_min, sector_max, + clutter_angle, clutter_band_half) + # 优先考虑“传统OMP至少有一个原子落在扇区外或杂波邻域”的样本 + bad_event = np.any((chosen_full < sector_min) | (chosen_full > sector_max)) or \ + np.any((chosen_full >= clutter_angle - clutter_band_half) & + (chosen_full <= clutter_angle + clutter_band_half)) + + if bad_event and cur_score > best_bad_score: + best_bad_score = cur_score + rep_full_proj0 = proj_hist_full[0] + rep_tsc_proj0 = proj_hist_tsc[0] + rep_full_chosen = chosen_full.copy() + rep_tsc_chosen = chosen_tsc.copy() + +# 如果一个“明确坏样本”都没抓到，就退而选投影谱最坏的那一个 +if rep_full_proj0 is None: + # 再跑一遍，抓 score 最高的 + best_bad_score = -1.0 + for mc in range(MC): + w_prior = build_sample_prior( + array_pos=array_pos, + target_angle=target_angle, + clutter_angle=clutter_angle, + SNR_dB=SNR_dB, + INR_dB=INR_dB, + noise_power=noise_power, + L_snapshots=L_snapshots, + diag_loading_ratio=diag_loading_ratio, + sector_min=sector_min, + sector_max=sector_max, + broad_step=broad_step + ) + _, chosen_full, proj_hist_full = omp_trace(w_prior, A_dic_full, full_grid, N_RF) + _, chosen_tsc, proj_hist_tsc = omp_trace(w_prior, A_dic_tsc, tsc_grid, N_RF) + + cur_score = compute_badness_score(proj_hist_full[0], full_grid, + sector_min, sector_max, + clutter_angle, clutter_band_half) + if cur_score > best_bad_score: + best_bad_score = cur_score + rep_full_proj0 = proj_hist_full[0] + rep_tsc_proj0 = proj_hist_tsc[0] + rep_full_chosen = chosen_full.copy() + rep_tsc_chosen = chosen_tsc.copy() + +# ============================================================ +# 6. 统计量计算 +# ============================================================ +total_selections = MC * N_RF + +P_hit_full = full_hit_total / total_selections +P_clutter_full = full_clutter_total / total_selections +P_out_full = full_out_total / total_selections + +P_hit_tsc = tsc_hit_total / total_selections +P_clutter_tsc = tsc_clutter_total / total_selections +P_out_tsc = tsc_out_total / total_selections + +print("==================================================") +print(f"阵列类型: {ARRAY_TYPE}") +print(f"SNR = {SNR_dB:.1f} dB, INR = {INR_dB:.1f} dB, 快拍数 = {L_snapshots}, MC = {MC}") +print("--------------------------------------------------") +print("传统 OMP（全空间字典）:") +print(f"目标扇区命中率 P_hit = {P_hit_full:.4f}") +print(f"杂波误选率 P_clutter = {P_clutter_full:.4f}") +print(f"扇区外误选率 P_out = {P_out_full:.4f}") +print("--------------------------------------------------") +print("TSC-OMP（目标扇区约束字典）:") +print(f"目标扇区命中率 P_hit = {P_hit_tsc:.4f}") +print(f"杂波误选率 P_clutter = {P_clutter_tsc:.4f}") +print(f"扇区外误选率 P_out = {P_out_tsc:.4f}") +print("==================================================") + +# 保存统计结果到 txt +with open(os.path.join(save_dir, 'pseudo_peak_statistics.txt'), 'w', encoding='utf-8') as f: + f.write(f"阵列类型: {ARRAY_TYPE}\n") + f.write(f"SNR = {SNR_dB:.1f} dB, INR = {INR_dB:.1f} dB, 快拍数 = {L_snapshots}, MC = {MC}\n\n") + f.write("传统 OMP（全空间字典）:\n") + f.write(f"P_hit = {P_hit_full:.6f}\n") + f.write(f"P_clutter = {P_clutter_full:.6f}\n") + f.write(f"P_out = {P_out_full:.6f}\n\n") + f.write("TSC-OMP（目标扇区约束字典）:\n") + f.write(f"P_hit = {P_hit_tsc:.6f}\n") + f.write(f"P_clutter = {P_clutter_tsc:.6f}\n") + f.write(f"P_out = {P_out_tsc:.6f}\n") + +# ============================================================ +# 7. 画图1：第一步投影谱对比（最适合论文说明“伪峰”） +# ============================================================ +fig1, axes = plt.subplots(2, 1, figsize=(10, 9)) + +# 上图：传统OMP第一步投影谱 +ax = axes[0] +ax.plot(full_grid, rep_full_proj0, 'b-', linewidth=1.8, label='传统OMP第一步投影谱') +ax.axvline(target_angle, color='g', linestyle='--', linewidth=1.5, label=f'目标 {target_angle:.2f}°') +ax.axvline(clutter_angle, color='r', linestyle='--', linewidth=1.5, label=f'杂波 {clutter_angle:.2f}°') + +# 标出目标扇区 +ymax = 1.05 * np.max(rep_full_proj0) +rect = Rectangle((sector_min, 0), sector_max - sector_min, ymax, + facecolor='green', alpha=0.10, edgecolor=None) +ax.add_patch(rect) +ax.text((sector_min + sector_max) / 2, 0.93 * ymax, '目标扇区', color='green', + ha='center', va='top', fontsize=11) + +# 标出杂波邻域 +rect2 = Rectangle((clutter_angle - clutter_band_half, 0), 2 * clutter_band_half, ymax, + facecolor='red', alpha=0.10, edgecolor=None) +ax.add_patch(rect2) +ax.text(clutter_angle, 0.80 * ymax, '杂波邻域', color='red', + ha='center', va='top', fontsize=11) + +# 标出传统OMP选中的角度 +for idx, ang in enumerate(rep_full_chosen): + val = rep_full_proj0[np.argmin(np.abs(full_grid - ang))] + ax.plot(ang, val, 'ko') + ax.text(ang, val + 0.02 * ymax, f'{idx+1}:{ang:.1f}°', fontsize=9, ha='center') + +ax.set_title('传统OMP第一步投影谱', fontsize=13, fontweight='bold') +ax.set_xlabel('角度 (°)') +ax.set_ylabel(r'$|a^H(\theta)r^{(0)}|$') +ax.grid(True, linestyle='--', alpha=0.5) +ax.set_xlim([-90, 90]) +ax.legend(loc='upper right', fontsize=10) + +# 下图：TSC-OMP第一步投影谱 +ax = axes[1] +ax.plot(tsc_grid, rep_tsc_proj0, 'm-', linewidth=1.8, label='TSC-OMP第一步投影谱') +ax.axvline(target_angle, color='g', linestyle='--', linewidth=1.5, label=f'目标 {target_angle:.2f}°') + +for idx, ang in enumerate(rep_tsc_chosen): + val = rep_tsc_proj0[np.argmin(np.abs(tsc_grid - ang))] + ax.plot(ang, val, 'ko') + ax.text(ang, val + 0.02 * np.max(rep_tsc_proj0), f'{idx+1}:{ang:.1f}°', fontsize=9, ha='center') + +ax.set_title('TSC-OMP第一步投影谱', fontsize=13, fontweight='bold') +ax.set_xlabel('角度 (°)') +ax.set_ylabel(r'$|a^H(\theta)r^{(0)}|$') +ax.grid(True, linestyle='--', alpha=0.5) +ax.set_xlim([sector_min, sector_max]) +ax.legend(loc='upper right', fontsize=10) + +plt.tight_layout() +fig1.savefig(os.path.join(save_dir, f'Fig1_projection_compare_{ARRAY_TYPE}.png'), dpi=300, bbox_inches='tight') + +# ============================================================ +# 8. 画图2：被选原子角度分布直方图 +# ============================================================ +fig2, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=False) + +# 传统OMP +ax = axes[0] +bins_full = np.arange(-90, 90 + 1, 1.0) +ax.hist(all_angles_full, bins=bins_full, color='steelblue', edgecolor='black', alpha=0.85) +ax.axvline(target_angle, color='g', linestyle='--', linewidth=1.5, label='目标') +ax.axvline(clutter_angle, color='r', linestyle='--', linewidth=1.5, label='杂波') +ax.axvspan(sector_min, sector_max, color='green', alpha=0.10, label='目标扇区') +ax.set_title('传统OMP所选原子角度分布', fontsize=13, fontweight='bold') +ax.set_xlabel('被选原子角度 (°)') +ax.set_ylabel('出现次数') +ax.grid(True, linestyle='--', alpha=0.5) +ax.legend(fontsize=10) + +# TSC-OMP +ax = axes[1] +bins_tsc = np.arange(sector_min, sector_max + 0.2, 0.2) +ax.hist(all_angles_tsc, bins=bins_tsc, color='orchid', edgecolor='black', alpha=0.85) +ax.axvline(target_angle, color='g', linestyle='--', linewidth=1.5, label='目标') +ax.axvspan(sector_min, sector_max, color='green', alpha=0.10, label='目标扇区') +ax.set_title('TSC-OMP所选原子角度分布', fontsize=13, fontweight='bold') +ax.set_xlabel('被选原子角度 (°)') +ax.set_ylabel('出现次数') +ax.grid(True, linestyle='--', alpha=0.5) +ax.legend(fontsize=10) + +plt.tight_layout() +fig2.savefig(os.path.join(save_dir, f'Fig2_hist_compare_{ARRAY_TYPE}.png'), dpi=300, bbox_inches='tight') + +# ============================================================ +# 9. 画图3：伪峰统计指标柱状图 +# ============================================================ +fig3, ax = plt.subplots(figsize=(9, 6)) + +metrics = ['目标扇区命中率', '杂波误选率', '扇区外误选率'] +omp_values = [P_hit_full, P_clutter_full, P_out_full] +tsc_values = [P_hit_tsc, P_clutter_tsc, P_out_tsc] + +x = np.arange(len(metrics)) +width = 0.34 + +bars1 = ax.bar(x - width/2, omp_values, width, label='传统OMP', color='steelblue') +bars2 = ax.bar(x + width/2, tsc_values, width, label='TSC-OMP', color='orchid') + +for bars in [bars1, bars2]: + for b in bars: + h = b.get_height() + ax.text(b.get_x() + b.get_width()/2, h + 0.01, f'{h:.3f}', + ha='center', va='bottom', fontsize=10) + +ax.set_xticks(x) +ax.set_xticklabels(metrics, fontsize=11) +ax.set_ylim([0, 1.08]) +ax.set_ylabel('概率 / 比例', fontsize=12) +ax.set_title('传统OMP与TSC-OMP伪峰统计指标对比', fontsize=13, fontweight='bold') +ax.grid(True, axis='y', linestyle='--', alpha=0.5) +ax.legend(fontsize=11) + +plt.tight_layout() +fig3.savefig(os.path.join(save_dir, f'Fig3_stat_compare_{ARRAY_TYPE}.png'), dpi=300, bbox_inches='tight') + +# ============================================================ +# 10. 显示图像 +# ============================================================ +plt.show() \ No newline at end of file