上传文件至「Code/Python/3.Robust and complexity test」

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)

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()

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()