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

2026-05-10 14:35:38 +08:00
parent 623735fe94
commit a203889de0
1 changed files with 535 additions and 0 deletions
@@ -0,0 +1,535 @@
+# -*- coding: utf-8 -*-
+"""
+TSC-OMP 局限性演示 GIF 生成脚本（顺滑增强版）
+------------------------------------------------
+生成 3 个 GIF：
+1) gif_prior_error_smooth.gif      先验误差敏感性
+2) gif_window_width_smooth.gif     波门宽度敏感性
+3) gif_clutter_approach_smooth.gif 强干扰逼近导致性能退化
+
+特点：
+- 帧率更高（默认 8 fps）
+- 参数扫描更密，动画更顺滑
+- 独立可运行，不依赖其他脚本
+"""
+
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation, PillowWriter
+
+# =========================================================
+# 0. 全局绘图设置
+# =========================================================
+plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei', 'Arial Unicode MS']
+plt.rcParams['axes.unicode_minus'] = False
+
+# =========================================================
+# 1. 系统参数
+# =========================================================
+N_ant = 24
+N_RF = 4
+d_lambda = 0.5
+
+target_angle = 5.04
+clutter_angle_default = -15.0
+
+noise_power = 1.0
+INR_dB = 30.0
+INR_linear = 10 ** (INR_dB / 10)
+
+L_snapshots = 100
+theta_scan = np.arange(-90, 90.01, 0.1)
+
+# Monte Carlo 参数（如果电脑慢，可以把 30 改成 20）
+Monte_Carlo = 30
+SNR_dB = 0.0  # GIF 中用于 RMSE 演示的默认 SNR
+
+# =========================================================
+# 2. GIF 平滑参数（核心改进）
+# =========================================================
+GIF_FPS = 8          # 帧率：越大越顺滑，但文件也更大
+GIF_INTERVAL = 80    # 播放间隔（毫秒）
+FIG_DPI = 110
+
+# 线条样式
+LINE_WIDTH = 2.6
+MARKER_SIZE = 7
+
+# 输出目录
+SAVE_DIR = "gif_outputs"
+os.makedirs(SAVE_DIR, exist_ok=True)
+
+# 随机种子
+GLOBAL_SEED = 2025
+
+# =========================================================
+# 3. 阵列导向矢量（ULA）
+# =========================================================
+def gen_a(theta_deg, N=N_ant, d=d_lambda):
+    """
+    ULA 导向矢量
+    """
+    theta_rad = np.deg2rad(theta_deg)
+    n_idx = np.arange(-(N - 1) / 2, (N - 1) / 2 + 1)
+    return np.exp(1j * 2 * np.pi * d * n_idx * np.sin(theta_rad)).reshape(-1, 1)
+
+# =========================================================
+# 4. 字典构造
+# =========================================================
+def build_dictionary(angle_grid):
+    A = np.zeros((N_ant, len(angle_grid)), dtype=complex)
+    for i, th in enumerate(angle_grid):
+        A[:, i:i+1] = gen_a(th) / np.sqrt(N_ant)
+    return A
+
+# =========================================================
+# 5. OMP 选择射频矩阵
+# =========================================================
+def omp_select(w_target, A_dic, Nrf):
+    """
+    基础 OMP：根据目标权向量 w_target 在字典 A_dic 中贪婪选择 Nrf 个原子
+    """
+    residual = w_target.copy()
+    selected = []
+    F_RF = np.zeros((A_dic.shape[0], 0), dtype=complex)
+
+    for _ in range(Nrf):
+        proj = A_dic.conj().T @ residual
+        idx = int(np.argmax(np.abs(proj)))
+        selected.append(idx)
+        F_RF = A_dic[:, selected]
+        residual = w_target - F_RF @ (np.linalg.pinv(F_RF) @ w_target)
+
+    return F_RF
+
+# =========================================================
+# 6. 计算静态先验数字权值（MVDR）
+# =========================================================
+def calc_prior_weight(a_target, a_clutter, sig_power_static=1.0):
+    """
+    用静态协方差设计数字先验权值
+    """
+    R_static = (noise_power * np.eye(N_ant) +
+                sig_power_static * (a_target @ a_target.conj().T) +
+                INR_linear * (a_clutter @ a_clutter.conj().T))
+    dl = 0.01 * np.trace(R_static) / N_ant
+    R_static_dl = R_static + dl * np.eye(N_ant)
+
+    inv_R = np.linalg.inv(R_static_dl)
+    w_prior = inv_R @ a_target / (a_target.conj().T @ inv_R @ a_target)
+    w_prior = w_prior / np.linalg.norm(w_prior)
+    return w_prior, R_static
+
+# =========================================================
+# 7. 构造 TSC-OMP 混合波束
+# =========================================================
+def build_tsc_structure(target_ang, clutter_ang, theta_prior, window_half_width):
+    """
+    根据目标角、干扰角、先验中心、波门半宽，构造 TSC-OMP 混合波束
+    """
+    a_target = gen_a(target_ang)
+    a_clutter = gen_a(clutter_ang)
+
+    # 局部约束字典（TSC）
+    dic_angles = np.arange(theta_prior - window_half_width,
+                           theta_prior + window_half_width + 1e-9,
+                           0.1)
+    A_dic = build_dictionary(dic_angles)
+
+    # 静态先验权值
+    w_prior, R_static = calc_prior_weight(a_target, a_clutter, sig_power_static=1.0)
+
+    # OMP 选模拟 RF
+    F_RF = omp_select(w_prior, A_dic, N_RF)
+
+    # 基带 MVDR
+    a_eff = F_RF.conj().T @ a_target
+    R_eff = F_RF.conj().T @ R_static @ F_RF
+    dl = 0.01 * np.trace(R_eff) / N_RF
+    R_eff_dl = R_eff + dl * np.eye(N_RF)
+
+    inv_R_eff = np.linalg.inv(R_eff_dl)
+    F_BB = inv_R_eff @ a_eff / (a_eff.conj().T @ inv_R_eff @ a_eff)
+
+    w_hybrid = F_RF @ F_BB
+    w_hybrid = w_hybrid / np.linalg.norm(w_hybrid)
+
+    return {
+        "a_target": a_target,
+        "a_clutter": a_clutter,
+        "F_RF": F_RF,
+        "w_tsc": w_hybrid
+    }
+
+# =========================================================
+# 8. 方向图
+# =========================================================
+def calc_beam_pattern(w, scan_grid):
+    bp = np.zeros(len(scan_grid))
+    for i, th in enumerate(scan_grid):
+        a = gen_a(th)
+        bp[i] = np.abs(a.conj().T @ w).item() ** 2
+    bp_db = 10 * np.log10(bp / np.max(bp) + 1e-12)
+    return bp_db
+
+# =========================================================
+# 9. 输出 SINR
+# =========================================================
+def calc_output_sinr(w, a_target, a_clutter, sig_power):
+    R_i_n = noise_power * np.eye(N_ant) + INR_linear * (a_clutter @ a_clutter.conj().T)
+    numerator = sig_power * np.abs(w.conj().T @ a_target).item() ** 2
+    denominator = np.real((w.conj().T @ R_i_n @ w).item())
+    return 10 * np.log10(numerator / (denominator + 1e-12) + 1e-12)
+
+# =========================================================
+# 10. 单次 Monte Carlo 下的 TSC-OMP 测角 RMSE
+# =========================================================
+def estimate_rmse_tsc(target_ang, clutter_ang, theta_prior, window_half_width,
+                      snr_db, monte_carlo=Monte_Carlo, seed_offset=0):
+    """
+    用低维 Capon 光谱做测角，返回 RMSE
+    为了兼顾速度和稳定性，这里做了适度的快速化处理
+    """
+    rng = np.random.default_rng(GLOBAL_SEED + seed_offset)
+
+    structure = build_tsc_structure(target_ang, clutter_ang, theta_prior, window_half_width)
+    F_RF = structure["F_RF"]
+    a_target = structure["a_target"]
+    a_clutter = structure["a_clutter"]
+
+    sig_power = 10 ** (snr_db / 10)
+
+    # 搜索网格只在波门内搜索
+    scan_grid = np.arange(theta_prior - window_half_width,
+                          theta_prior + window_half_width + 1e-9,
+                          0.05)
+    A_scan = np.hstack([gen_a(th) for th in scan_grid])
+    A_eff = F_RF.conj().T @ A_scan
+
+    err_sq = []
+
+    for _ in range(monte_carlo):
+        s_t = np.sqrt(sig_power / 2) * (rng.standard_normal((1, L_snapshots)) +
+                                        1j * rng.standard_normal((1, L_snapshots)))
+        s_c = np.sqrt(INR_linear / 2) * (rng.standard_normal((1, L_snapshots)) +
+                                         1j * rng.standard_normal((1, L_snapshots)))
+        n_t = np.sqrt(noise_power / 2) * (rng.standard_normal((N_ant, L_snapshots)) +
+                                          1j * rng.standard_normal((N_ant, L_snapshots)))
+
+        X = a_target @ s_t + a_clutter @ s_c + n_t
+        Y = F_RF.conj().T @ X
+
+        R_bb = (Y @ Y.conj().T) / L_snapshots
+
+        # 对角加载（提升稳定性）
+        gamma = 0.05 * np.trace(R_bb) / N_RF
+        R_bb_dl = R_bb + gamma * np.eye(N_RF)
+
+        inv_R = np.linalg.inv(R_bb_dl)
+
+        # Capon 谱：P(theta)=1/(a^H inv(R) a)
+        temp = inv_R @ A_eff
+        denom = np.sum(np.conj(A_eff) * temp, axis=0)
+        spectrum = 1.0 / (np.abs(denom) + 1e-12)
+
+        est_ang = scan_grid[np.argmax(spectrum)]
+        err_sq.append((est_ang - target_ang) ** 2)
+
+    return float(np.sqrt(np.mean(err_sq)))
+
+# =========================================================
+# 11. 一些辅助绘图函数
+# =========================================================
+def beautify_axes(ax):
+    ax.grid(True, linestyle='--', alpha=0.5)
+    ax.tick_params(labelsize=11)
+
+# =========================================================
+# 12. GIF 1：先验误差敏感性
+# =========================================================
+def make_prior_error_gif():
+    print("正在生成 GIF 1：先验误差敏感性 ...")
+
+    # 帧更密，动画更顺滑
+    prior_errors = np.arange(-8.0, 8.01, 0.5)
+
+    patterns = []
+    sinrs = []
+    rmses = []
+
+    for i, e in enumerate(prior_errors):
+        theta_prior = round(target_angle) + e
+        structure = build_tsc_structure(target_angle, clutter_angle_default,
+                                        theta_prior, 5.0)
+        bp = calc_beam_pattern(structure["w_tsc"], theta_scan)
+        sinr = calc_output_sinr(structure["w_tsc"],
+                                structure["a_target"],
+                                structure["a_clutter"],
+                                10 ** (SNR_dB / 10))
+        rmse = estimate_rmse_tsc(target_angle, clutter_angle_default,
+                                 theta_prior, 5.0, SNR_dB,
+                                 monte_carlo=Monte_Carlo,
+                                 seed_offset=1000 + i)
+
+        patterns.append(bp)
+        sinrs.append(sinr)
+        rmses.append(rmse)
+
+    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
+    ax1, ax2 = axes
+
+    line_bp, = ax1.plot([], [], lw=LINE_WIDTH, color="#0072BD", label="TSC-OMP方向图")
+    ax1.axvline(target_angle, color='g', linestyle='--', lw=2, label="真实目标")
+    prior_line = ax1.axvline(target_angle, color='m', linestyle='--', lw=2, label="先验中心")
+    ax1.set_xlim([-20, 20])
+    ax1.set_ylim([-70, 2])
+    ax1.set_xlabel("角度 θ (°)", fontsize=13)
+    ax1.set_ylabel("归一化增益 (dB)", fontsize=13)
+    ax1.set_title("先验误差增大时的性能退化", fontsize=15, fontweight='bold')
+    beautify_axes(ax1)
+    ax1.legend(fontsize=11)
+
+    line_rmse, = ax2.plot([], [], color="#D95319", lw=LINE_WIDTH,
+                          marker='o', markersize=MARKER_SIZE, label="RMSE")
+    current_pt, = ax2.plot([], [], 'ro', markersize=9)
+    ax2.set_xlim(prior_errors[0], prior_errors[-1])
+    ax2.set_ylim(0, max(rmses) * 1.25 + 0.1)
+    ax2.set_xlabel("先验角误差 (°)", fontsize=13)
+    ax2.set_ylabel("测角 RMSE (°)", fontsize=13)
+    ax2.set_title("先验误差越大，RMSE 越容易升高", fontsize=15, fontweight='bold')
+    beautify_axes(ax2)
+    ax2.legend(fontsize=11)
+
+    info = ax1.text(0.02, 0.96, '', transform=ax1.transAxes, va='top',
+                    bbox=dict(boxstyle='round', facecolor='white', alpha=0.88), fontsize=11)
+
+    def update(frame):
+        e = prior_errors[frame]
+        theta_prior = round(target_angle) + e
+
+        line_bp.set_data(theta_scan, patterns[frame])
+        prior_line.set_xdata([theta_prior, theta_prior])
+
+        line_rmse.set_data(prior_errors[:frame+1], rmses[:frame+1])
+        current_pt.set_data([e], [rmses[frame]])
+
+        info.set_text(
+            f"先验角误差 = {e:+.2f}°\n"
+            f"先验中心 = {theta_prior:.2f}°\n"
+            f"输出 SINR = {sinrs[frame]:.2f} dB\n"
+            f"RMSE = {rmses[frame]:.3f}°"
+        )
+
+        return line_bp, prior_line, line_rmse, current_pt, info
+
+    ani = FuncAnimation(fig, update, frames=len(prior_errors),
+                        interval=GIF_INTERVAL, blit=False, repeat=True)
+
+    save_path = os.path.join(SAVE_DIR, "gif_prior_error_smooth.gif")
+    ani.save(save_path, writer=PillowWriter(fps=GIF_FPS), dpi=FIG_DPI)
+    plt.close(fig)
+
+    print(f"GIF 1 保存完成：{save_path}")
+
+# =========================================================
+# 13. GIF 2：波门宽度敏感性
+# =========================================================
+def make_window_width_gif():
+    print("正在生成 GIF 2：波门宽度敏感性 ...")
+
+    widths = np.arange(2.0, 10.01, 0.5)
+
+    patterns = []
+    sinrs = []
+    rmses = []
+
+    for i, w in enumerate(widths):
+        theta_prior = round(target_angle)
+        structure = build_tsc_structure(target_angle, clutter_angle_default,
+                                        theta_prior, w)
+        bp = calc_beam_pattern(structure["w_tsc"], theta_scan)
+        sinr = calc_output_sinr(structure["w_tsc"],
+                                structure["a_target"],
+                                structure["a_clutter"],
+                                10 ** (SNR_dB / 10))
+        rmse = estimate_rmse_tsc(target_angle, clutter_angle_default,
+                                 theta_prior, w, SNR_dB,
+                                 monte_carlo=Monte_Carlo,
+                                 seed_offset=2000 + i)
+
+        patterns.append(bp)
+        sinrs.append(sinr)
+        rmses.append(rmse)
+
+    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
+    ax1, ax2 = axes
+
+    line_bp, = ax1.plot([], [], lw=LINE_WIDTH, color="#0072BD", label="TSC-OMP方向图")
+    ax1.axvline(target_angle, color='g', linestyle='--', lw=2, label="真实目标")
+    gate_left = ax1.axvline(target_angle - widths[0], color='m', linestyle=':', lw=2)
+    gate_right = ax1.axvline(target_angle + widths[0], color='m', linestyle=':', lw=2, label="波门边界")
+    ax1.set_xlim([-20, 20])
+    ax1.set_ylim([-70, 2])
+    ax1.set_xlabel("角度 θ (°)", fontsize=13)
+    ax1.set_ylabel("归一化增益 (dB)", fontsize=13)
+    ax1.set_title("波门宽度变化对方向图的影响", fontsize=15, fontweight='bold')
+    beautify_axes(ax1)
+    ax1.legend(fontsize=11)
+
+    line_rmse, = ax2.plot([], [], color="#D95319", lw=LINE_WIDTH,
+                          marker='o', markersize=MARKER_SIZE, label="RMSE")
+    current_pt, = ax2.plot([], [], 'ro', markersize=9)
+    ax2.set_xlim(widths[0], widths[-1])
+    ax2.set_ylim(0, max(rmses) * 1.25 + 0.1)
+    ax2.set_xlabel("波门半宽 (°)", fontsize=13)
+    ax2.set_ylabel("测角 RMSE (°)", fontsize=13)
+    ax2.set_title("波门过窄或过宽都会带来性能损失", fontsize=15, fontweight='bold')
+    beautify_axes(ax2)
+    ax2.legend(fontsize=11)
+
+    info = ax1.text(0.02, 0.96, '', transform=ax1.transAxes, va='top',
+                    bbox=dict(boxstyle='round', facecolor='white', alpha=0.88), fontsize=11)
+
+    def update(frame):
+        w = widths[frame]
+
+        line_bp.set_data(theta_scan, patterns[frame])
+        gate_left.set_xdata([target_angle - w, target_angle - w])
+        gate_right.set_xdata([target_angle + w, target_angle + w])
+
+        line_rmse.set_data(widths[:frame+1], rmses[:frame+1])
+        current_pt.set_data([w], [rmses[frame]])
+
+        info.set_text(
+            f"波门半宽 = ±{w:.2f}°\n"
+            f"输出 SINR = {sinrs[frame]:.2f} dB\n"
+            f"RMSE = {rmses[frame]:.3f}°"
+        )
+
+        return line_bp, gate_left, gate_right, line_rmse, current_pt, info
+
+    ani = FuncAnimation(fig, update, frames=len(widths),
+                        interval=GIF_INTERVAL, blit=False, repeat=True)
+
+    save_path = os.path.join(SAVE_DIR, "gif_window_width_smooth.gif")
+    ani.save(save_path, writer=PillowWriter(fps=GIF_FPS), dpi=FIG_DPI)
+    plt.close(fig)
+
+    print(f"GIF 2 保存完成：{save_path}")
+
+# =========================================================
+# 14. GIF 3：强干扰逼近目标
+# =========================================================
+def make_clutter_approach_gif():
+    print("正在生成 GIF 3：强干扰逼近目标时的性能退化 ...")
+
+    # 更密的干扰角扫描，让动画更顺滑
+    clutter_list = np.arange(-25.0, -1.99, 0.5)
+
+    patterns = []
+    sinrs = []
+    rmses = []
+    separations = []
+
+    for i, cang in enumerate(clutter_list):
+        theta_prior = round(target_angle)
+        structure = build_tsc_structure(target_angle, cang,
+                                        theta_prior, 5.0)
+        bp = calc_beam_pattern(structure["w_tsc"], theta_scan)
+        sinr = calc_output_sinr(structure["w_tsc"],
+                                structure["a_target"],
+                                structure["a_clutter"],
+                                10 ** (SNR_dB / 10))
+        rmse = estimate_rmse_tsc(target_angle, cang,
+                                 theta_prior, 5.0, SNR_dB,
+                                 monte_carlo=Monte_Carlo,
+                                 seed_offset=3000 + i)
+
+        patterns.append(bp)
+        sinrs.append(sinr)
+        rmses.append(rmse)
+        separations.append(abs(target_angle - cang))
+
+    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
+    ax1, ax2 = axes
+
+    line_bp, = ax1.plot([], [], lw=LINE_WIDTH, color="#0072BD", label="TSC-OMP方向图")
+    ax1.axvline(target_angle, color='g', linestyle='--', lw=2, label="目标")
+    clutter_line = ax1.axvline(clutter_list[0], color='r', linestyle='--', lw=2, label="强干扰")
+    ax1.set_xlim([-40, 20])
+    ax1.set_ylim([-70, 2])
+    ax1.set_xlabel("角度 θ (°)", fontsize=13)
+    ax1.set_ylabel("归一化增益 (dB)", fontsize=13)
+    ax1.set_title("强干扰逐步逼近目标时的方向图退化", fontsize=15, fontweight='bold')
+    beautify_axes(ax1)
+    ax1.legend(fontsize=11)
+
+    line_rmse, = ax2.plot([], [], color="#D95319", lw=LINE_WIDTH,
+                          marker='o', markersize=MARKER_SIZE, label="RMSE")
+    current_pt, = ax2.plot([], [], 'ro', markersize=9)
+    ax2.set_xlim(max(separations), min(separations))
+    ax2.invert_xaxis()
+    ax2.set_ylim(0, max(rmses) * 1.25 + 0.1)
+    ax2.set_xlabel("目标-干扰角距 |θt-θc| (°)", fontsize=13)
+    ax2.set_ylabel("测角 RMSE (°)", fontsize=13)
+    ax2.set_title("角距减小后，可分性下降，RMSE 增大", fontsize=15, fontweight='bold')
+    beautify_axes(ax2)
+    ax2.legend(fontsize=11)
+
+    info = ax1.text(0.02, 0.96, '', transform=ax1.transAxes, va='top',
+                    bbox=dict(boxstyle='round', facecolor='white', alpha=0.88), fontsize=11)
+
+    def update(frame):
+        cang = clutter_list[frame]
+        sep = separations[frame]
+
+        line_bp.set_data(theta_scan, patterns[frame])
+        clutter_line.set_xdata([cang, cang])
+
+        line_rmse.set_data(separations[:frame+1], rmses[:frame+1])
+        current_pt.set_data([sep], [rmses[frame]])
+
+        info.set_text(
+            f"干扰角度 = {cang:.2f}°\n"
+            f"目标-干扰角距 = {sep:.2f}°\n"
+            f"输出 SINR = {sinrs[frame]:.2f} dB\n"
+            f"RMSE = {rmses[frame]:.3f}°"
+        )
+
+        return line_bp, clutter_line, line_rmse, current_pt, info
+
+    ani = FuncAnimation(fig, update, frames=len(clutter_list),
+                        interval=GIF_INTERVAL, blit=False, repeat=True)
+
+    save_path = os.path.join(SAVE_DIR, "gif_clutter_approach_smooth.gif")
+    ani.save(save_path, writer=PillowWriter(fps=GIF_FPS), dpi=FIG_DPI)
+    plt.close(fig)
+
+    print(f"GIF 3 保存完成：{save_path}")
+
+# =========================================================
+# 15. 主函数
+# =========================================================
+def main():
+    print("====================================================")
+    print("开始生成 TSC-OMP 局限性分析 GIF（顺滑增强版）")
+    print("输出目录：", os.path.abspath(SAVE_DIR))
+    print(f"GIF 帧率：{GIF_FPS} fps")
+    print(f"动画间隔：{GIF_INTERVAL} ms")
+    print("====================================================")
+
+    make_prior_error_gif()
+    make_window_width_gif()
+    make_clutter_approach_gif()
+
+    print("====================================================")
+    print("全部 GIF 生成完毕！")
+    print("生成文件：")
+    print("1) gif_prior_error_smooth.gif")
+    print("2) gif_window_width_smooth.gif")
+    print("3) gif_clutter_approach_smooth.gif")
+    print("====================================================")
+
+if __name__ == "__main__":
+    main()