From a203889de0e117eae8add1f96a7105f9ac1fffdd Mon Sep 17 00:00:00 2001 From: david1465833828 Date: Sun, 10 May 2026 14:35:38 +0800 Subject: [PATCH] =?UTF-8?q?=E4=B8=8A=E4=BC=A0=E6=96=87=E4=BB=B6=E8=87=B3?= =?UTF-8?q?=E3=80=8CCode/Python/3.Robust=20and=20complexity=20test?= =?UTF-8?q?=E3=80=8D?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../gif_generator.py | 535 ++++++++++++++++++ 1 file changed, 535 insertions(+) create mode 100644 Code/Python/3.Robust and complexity test/gif_generator.py diff --git a/Code/Python/3.Robust and complexity test/gif_generator.py b/Code/Python/3.Robust and complexity test/gif_generator.py new file mode 100644 index 0000000..2fa600a --- /dev/null +++ b/Code/Python/3.Robust and complexity test/gif_generator.py @@ -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() \ No newline at end of file