上传文件至「Graduation Design」

Graduation Design/RMSE_Test_nula.py

@@ -0,0 +1,156 @@
+import torch
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+
+# ==========================================
+# 0. 全局绘图参数与 CUDA 引擎初始化
+# ==========================================
+plt.rcParams['axes.unicode_minus'] = False
+plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
+plt.rcParams['mathtext.fontset'] = 'stix' 
+
+device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+torch.manual_seed(2024)
+
+N_ant, N_RF = 24, 4
+clutter_angle = -15.0
+noise_power = 1.0
+INR_linear = 10**(30/10)
+
+# 统一使用非均匀线阵 (NULA) 物理间隔
+fc = 4.9e9
+lam = 3e8 / fc
+an_a = np.arange(0, 11*43 + 1, 43)
+an_b = np.array([11*43 + 63])
+an_c = 11*43 + 63 + np.arange(43, 11*43 + 1, 43)
+delta_d = 1e-3 * np.concatenate((an_a, an_b, an_c))  # Shape: (24,)
+delta_d_tensor = torch.tensor(delta_d, dtype=torch.float32, device=device).unsqueeze(1)
+
+def gen_a(theta_deg):
+    theta_rad = torch.deg2rad(torch.tensor(theta_deg, dtype=torch.float32, device=device))
+    # 注意这里统一使用 NULA 的物理模型和负指数
+    return torch.exp(-1j * 2 * torch.pi * delta_d_tensor * torch.sin(theta_rad) / lam)
+
+a_clutter = gen_a(clutter_angle)
+R_interference = noise_power * torch.eye(N_ant, device=device) + INR_linear * (a_clutter @ a_clutter.mH)
+
+dic_angles_full = torch.arange(-90, 90.5, 0.5, device=device)
+A_dic_full = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in dic_angles_full], dim=1)
+
+
+angle_ranges = [np.arange(0, 11), np.arange(11, 21), np.arange(21, 31), np.arange(31, 41)]
+snr_test_list = [0, 5, 10, 15, 20]
+Monte_Carlo = 100
+L_snapshots = 100
+colors = ['#0072BD', '#D95319', '#EDB120', '#7E2F8E', '#77AC30']
+
+def run_omp(w_opt, A_dic):
+    res = w_opt.clone()
+    F_RF_list = []
+    for k in range(N_RF):
+        idx = torch.argmax(torch.abs(A_dic.mH @ res))
+        F_RF_list.append(A_dic[:, idx:idx+1])
+        F_RF_mat = torch.cat(F_RF_list, dim=1)
+        res = w_opt - F_RF_mat @ (torch.linalg.pinv(F_RF_mat) @ w_opt)
+    return F_RF_mat
+
+print(f"启动 CUDA 并行引擎，对 4 大扇区进行极速扫雷测试 (MC={Monte_Carlo})...")
+start_time = time.time()
+
+# 取消了这里原本的 fig, axs = plt.subplots(2, 2)
+# 改为在循环内部动态创建单图
+
+for fig_idx, test_angles in enumerate(angle_ranges):
+    print(f"正在计算扇区 {test_angles[0]}° ~ {test_angles[-1]}° ...")
+    
+    # 【修改点 1】：为每一个扇区单独生成一张独立画布
+    fig, ax = plt.subplots(figsize=(9, 7))
+    
+    RMSE_DFT_All = np.zeros((len(snr_test_list), len(test_angles)))
+    RMSE_FULL_All = np.zeros((len(snr_test_list), len(test_angles)))
+    RMSE_TSC_All = np.zeros((len(snr_test_list), len(test_angles)))
+    
+    for snr_idx, snr_dB in enumerate(snr_test_list):
+        sig_p = 10**(snr_dB/10)
+        
+        for a_idx, t_angle in enumerate(test_angles):
+            true_target_angle = t_angle + 0.04
+            theta_prior = round(true_target_angle)
+            window = 5
+            
+            dic_tsc = torch.arange(theta_prior - window, theta_prior + window + 0.1, 0.1, device=device)
+            A_dic_tsc = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in dic_tsc], dim=1)
+            
+            angles_dft = torch.linspace(theta_prior - window, theta_prior + window, N_RF, device=device)
+            F_RF_dft = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in angles_dft], dim=1)
+            
+            a_prior_broad = torch.zeros((N_ant, 1), dtype=torch.complex64, device=device)
+            for tb in torch.arange(theta_prior - window, theta_prior + window + 0.1, 0.5, device=device):
+                a_prior_broad += gen_a(tb)
+            a_prior_broad /= torch.norm(a_prior_broad)
+            
+            w_dig_prior = torch.linalg.inv(R_interference) @ a_prior_broad
+            w_dig_prior /= torch.norm(w_dig_prior)
+            
+            F_RF_full = run_omp(w_dig_prior, A_dic_full)
+            F_RF_tsc = run_omp(w_dig_prior, A_dic_tsc)
+            
+            scan_grid = torch.arange(theta_prior - window, theta_prior + window + 0.01, 0.01, device=device)
+            A_scan = torch.cat([gen_a(th) for th in scan_grid], dim=1)
+            
+            A_eff_dft = F_RF_dft.mH @ A_scan
+            A_eff_full = F_RF_full.mH @ A_scan
+            A_eff_tsc = F_RF_tsc.mH @ A_scan
+            
+            a_target = gen_a(true_target_angle)
+            s_t = np.sqrt(sig_p/2) * (torch.randn(Monte_Carlo, 1, L_snapshots, device=device) + 1j * torch.randn(Monte_Carlo, 1, L_snapshots, device=device))
+            s_c = np.sqrt(INR_linear/2) * (torch.randn(Monte_Carlo, 1, L_snapshots, device=device) + 1j * torch.randn(Monte_Carlo, 1, L_snapshots, device=device))
+            n_t = np.sqrt(noise_power/2) * (torch.randn(Monte_Carlo, N_ant, L_snapshots, device=device) + 1j * torch.randn(Monte_Carlo, N_ant, L_snapshots, device=device))
+            
+            X = a_target.unsqueeze(0) * s_t + a_clutter.unsqueeze(0) * s_c + n_t
+            
+            def process_hybrid(F_RF, A_eff):
+                Y = F_RF.mH.unsqueeze(0) @ X
+                R = (Y @ Y.mH) / L_snapshots
+                R += (0.05 * torch.diagonal(R, dim1=-2, dim2=-1).sum(-1).view(-1, 1, 1) / N_RF) * torch.eye(N_RF, device=device).unsqueeze(0)
+                inv_R = torch.linalg.inv(R)
+                P = 1.0 / torch.abs(torch.einsum('si, bij, js -> bs', A_eff.conj().T, inv_R, A_eff))
+                return np.sqrt(np.mean(((scan_grid[torch.argmax(P, dim=1)] - true_target_angle)**2).cpu().numpy()))
+
+            RMSE_DFT_All[snr_idx, a_idx] = process_hybrid(F_RF_dft, A_eff_dft)
+            RMSE_FULL_All[snr_idx, a_idx] = process_hybrid(F_RF_full, A_eff_full)
+            RMSE_TSC_All[snr_idx, a_idx] = process_hybrid(F_RF_tsc, A_eff_tsc)
+            
+        c_str = colors[snr_idx]
+        lbl_suffix = f" (SNR={snr_dB}dB)"
+        ax.semilogy(test_angles, RMSE_DFT_All[snr_idx, :], '--^', color=c_str, linewidth=1.5, label='DFT'+lbl_suffix)
+        ax.semilogy(test_angles, RMSE_FULL_All[snr_idx, :], ':x', color=c_str, linewidth=2, label='传统OMP'+lbl_suffix)
+        ax.semilogy(test_angles, RMSE_TSC_All[snr_idx, :], '-o', color=c_str, linewidth=2.5, label='TSC-OMP'+lbl_suffix)
+
+    ax.set_xlim([test_angles[0], test_angles[-1]])
+    ax.set_xticks(test_angles)
+    ax.grid(True, which="both", ls="--", alpha=0.5)
+    ax.set_ylim([1e-3, 20])
+    ax.set_xlabel('目标真实出现角度 $\\theta (\\circ)$', fontsize=12, fontweight='bold')
+    ax.set_ylabel(r'测角 RMSE $(^\circ)$', fontsize=12, fontweight='bold')
+    ax.set_title(f'宽波束先验引导下，扇区 ({test_angles[0]}°~{test_angles[-1]}°) 鲁棒性评估', fontsize=13, fontweight='bold')
+
+
+    # 【修改点 2】：独立设置每张图的图例。
+    plt.tight_layout() 
+    fig.subplots_adjust(bottom=0.25) 
+    
+    # 【关键调整】：将 bbox_to_anchor 的 Y 坐标从 -0.15 提上来，改为 -0.08 
+    # 这样图例就会离 X 轴更近，避开底部操作条
+    ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.10), ncol=3, fontsize=10)
+
+    # 【修改点 3】：自动保存
+    file_name = f'RMSE_Sector_{test_angles[0]}_{test_angles[-1]}.png'
+    plt.savefig(file_name, dpi=300, bbox_inches='tight')
+    print(f"  已保存独立图表: {file_name}")
+
+print(f"全扇区计算完毕！总耗时: {time.time() - start_time:.4f} 秒。")
+
+# 一次性在屏幕上弹出生成的4张独立图表
+plt.show()

Graduation Design/TSC-DL-Animation.py

@@ -0,0 +1,148 @@
+import torch
+import torch.nn as nn
+import torch.optim as optim
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+import time
+
+# ==========================================
+# 0. 全局设置
+# ==========================================
+plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
+plt.rcParams['axes.unicode_minus'] = False
+device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+print(f"========== 启动 AI 动态追踪生成器 ==========\n使用设备: {device}")
+
+N_ant, N_RF, d_lambda = 24, 4, 0.5
+noise_power = 1.0
+INR_linear = 10**(30/10) # 30dB 强杂波
+
+def gen_a_batch(theta_deg_tensor):
+    theta_rad = torch.deg2rad(theta_deg_tensor)
+    n_idx = torch.arange(-(N_ant-1)/2, (N_ant-1)/2 + 0.1, device=device).unsqueeze(0)
+    phases = 2 * torch.pi * d_lambda * n_idx * torch.sin(theta_rad.unsqueeze(1))
+    return torch.exp(1j * phases).unsqueeze(2)
+
+# ==========================================
+# 1. 广域追踪网络 (彻底解决相位死锁)
+# ==========================================
+class DynamicHBFNet(nn.Module):
+    def __init__(self):
+        super().__init__()
+        self.net = nn.Sequential(
+            nn.Linear(48, 512), nn.GELU(),
+            nn.Linear(512, 512), nn.GELU(),
+            nn.Linear(512, N_ant * N_RF * 2) 
+        )
+
+    def forward(self, w_opt):
+        w_flat = w_opt.squeeze(2)
+        x = torch.cat([w_flat.real, w_flat.imag], dim=1)
+        out = self.net(x).view(-1, N_ant, N_RF, 2)
+        complex_out = out[..., 0] + 1j * out[..., 1]
+        # 严格恒模约束
+        F_RF = complex_out / (torch.abs(complex_out) + 1e-8) / np.sqrt(N_ant)
+        return F_RF
+
+# ==========================================
+# 2. 广域动态特训
+# ==========================================
+model = DynamicHBFNet().to(device)
+optimizer = optim.Adam(model.parameters(), lr=0.002)
+
+Batch_Size, Epochs = 500, 10000      
+
+print("正在对 AI 进行大范围目标与杂波抗干扰特训...")
+for epoch in range(Epochs):
+    model.train()
+    optimizer.zero_grad()
+    
+    # 【动态大扫掠特训】：目标 0~40°，杂波 -20~-3°
+    theta_t = torch.rand(Batch_Size, device=device) * 40.0
+    theta_c = -20.0 + torch.rand(Batch_Size, device=device) * 17.0
+    
+    a_t, a_c = gen_a_batch(theta_t), gen_a_batch(theta_c)
+    
+    R_interf = noise_power * torch.eye(N_ant, device=device).unsqueeze(0) + INR_linear * (a_c @ a_c.mH)
+    w_opt = torch.linalg.pinv(R_interf) @ a_t
+    w_opt = w_opt / torch.norm(w_opt, dim=1, keepdim=True)
+    
+    F_RF = model(w_opt)
+    F_BB = torch.linalg.pinv(F_RF.mH @ F_RF + 1e-3*torch.eye(N_RF, device=device).unsqueeze(0)) @ (F_RF.mH @ w_opt)
+    w_ai = F_RF @ F_BB
+    w_ai = w_ai / torch.norm(w_ai, dim=1, keepdim=True)
+    
+    sim = torch.abs(torch.sum(w_opt.conj() * w_ai, dim=(1,2)))**2
+    loss_sim = torch.mean(1.0 - sim)
+    loss_null = torch.mean(torch.abs(torch.sum(a_c.conj() * w_ai, dim=(1,2)))**2)
+    
+    loss = loss_sim + 150.0 * loss_null 
+    loss.backward()
+    optimizer.step()
+
+print("特训完成！正在生成实时追踪 GIF 动图...")
+
+# ==========================================
+# 3. 生成精美动态 GIF
+# ==========================================
+model.eval()
+theta_scan = torch.arange(-60, 60.1, 0.5, device=device) # 限制扫描角度让图像更聚焦
+A_scan = gen_a_batch(theta_scan).squeeze(2).T 
+
+fig, ax = plt.subplots(figsize=(10, 6))
+ax.set_xlim(-60, 60)
+ax.set_ylim(-60, 0)
+ax.set_xlabel('角度 θ (°)', fontsize=12, fontweight='bold')
+ax.set_ylabel('归一化增益 (dB)', fontsize=12, fontweight='bold')
+ax.grid(True, ls='--')
+
+line_opt, = ax.plot([], [], 'k-', lw=2, alpha=0.5, label='纯数字最优下限')
+line_ai,  = ax.plot([], [], color='#D95319', ls='--', lw=2.5, label='AI 瞬时感知混合波束')
+vline_t   = ax.axvline(0, color='g', ls='-.', lw=2, label='机动目标轨迹')
+vline_c   = ax.axvline(0, color='r', ls='-.', lw=2, label='强杂波轨迹')
+
+title_text = ax.set_title('', fontsize=14, fontweight='bold')
+ax.legend(loc='lower center', fontsize=10)
+
+frames_count = 60 # 60帧平滑过渡
+
+def update(frame):
+    # 目标从 0 匀速扫到 40
+    t_ang = 0.0 + frame * (40.0 / (frames_count - 1))
+    # 杂波从 -3 匀速扫到 -20
+    c_ang = -3.0 - frame * (17.0 / (frames_count - 1))
+    
+    a_t = gen_a_batch(torch.tensor([t_ang], device=device))
+    a_c = gen_a_batch(torch.tensor([c_ang], device=device))
+    R_interf = noise_power * torch.eye(N_ant, device=device).unsqueeze(0) + INR_linear * (a_c @ a_c.mH)
+    
+    w_opt = torch.linalg.pinv(R_interf) @ a_t
+    w_opt /= torch.norm(w_opt, dim=1, keepdim=True)
+    
+    # AI 仅凭一次前向传播，瞬间生成波束！
+    with torch.no_grad():
+        F_RF = model(w_opt)
+        F_BB = torch.linalg.pinv(F_RF.mH @ F_RF + 1e-3*torch.eye(N_RF, device=device).unsqueeze(0)) @ (F_RF.mH @ w_opt)
+        w_ai = F_RF @ F_BB
+        w_ai /= torch.norm(w_ai, dim=1, keepdim=True)
+        
+    gain_opt = torch.abs(torch.sum(w_opt.conj() * A_scan.unsqueeze(0), dim=1))**2
+    bp_opt = 10 * torch.log10(gain_opt).cpu().numpy().flatten()
+    bp_opt -= np.max(bp_opt)
+    
+    gain_ai = torch.abs(torch.sum(w_ai.conj() * A_scan.unsqueeze(0), dim=1))**2
+    bp_ai = 10 * torch.log10(gain_ai).cpu().numpy().flatten()
+    bp_ai -= np.max(bp_ai)
+    
+    line_opt.set_data(theta_scan.cpu().numpy(), bp_opt)
+    line_ai.set_data(theta_scan.cpu().numpy(), bp_ai)
+    vline_t.set_xdata([t_ang, t_ang])
+    vline_c.set_xdata([c_ang, c_ang])
+    
+    title_text.set_text(f'AI 端到端波束实时追踪: 目标 {t_ang:.1f}° | 强杂波 {c_ang:.1f}°')
+    return line_opt, line_ai, vline_t, vline_c, title_text
+
+# 渲染动图
+anim = animation.FuncAnimation(fig, update, frames=frames_count, interval=100, blit=False)
+anim.save('E:\\Coding\\graduation design\\MIMO_Beamforming_Simulation_Platform\\final\\python\\Deeplearning\\AI_Beam_Tracking.gif', writer='pillow', fps=10)

Graduation Design/TSC-DL.py

@@ -0,0 +1,153 @@
+import torch
+import torch.nn as nn
+import torch.optim as optim
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+
+# ==========================================
+# 0. 全局设置 (对标你的 -35° 强杂波挑战)
+# ==========================================
+plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
+plt.rcParams['axes.unicode_minus'] = False
+device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+print(f"========== 启动 TSC-DL 扇区聚焦引擎 ==========\n使用设备: {device}")
+
+N_ant, N_RF, d_lambda = 24, 4, 0.5
+noise_power = 1.0
+INR_linear = 10**(30/10) # 30dB 强杂波
+
+#核心测试场景
+test_t_ang, test_c_ang = 15.0, -25.0 
+
+def gen_a_batch(theta_deg_tensor):
+    theta_rad = torch.deg2rad(theta_deg_tensor)
+    n_idx = torch.arange(-(N_ant-1)/2, (N_ant-1)/2 + 0.1, device=device).unsqueeze(0)
+    phases = 2 * torch.pi * d_lambda * n_idx * torch.sin(theta_rad.unsqueeze(1))
+    return torch.exp(1j * phases).unsqueeze(2)
+
+# ==========================================
+# 1. 复数域正交网络 (彻底消灭相位梯度死锁)
+# ==========================================
+class SectorHBFNet(nn.Module):
+    def __init__(self):
+        super().__init__()
+        self.net = nn.Sequential(
+            nn.Linear(48, 512), nn.GELU(),
+            nn.Linear(512, 512), nn.GELU(),
+            # 【核心架构升级】：不直接输出相位！输出复数的实部和虚部 (2个通道)
+            nn.Linear(512, N_ant * N_RF * 2) 
+        )
+
+    def forward(self, w_opt):
+        w_flat = w_opt.squeeze(2)
+        x = torch.cat([w_flat.real, w_flat.imag], dim=1)
+        out = self.net(x).view(-1, N_ant, N_RF, 2)
+        
+        complex_out = out[..., 0] + 1j * out[..., 1]
+        # 【物理强制约束】：将复数除以自身的模长，绝对精准地满足恒模约束，且梯度无比丝滑！
+        F_RF = complex_out / (torch.abs(complex_out) + 1e-8) / np.sqrt(N_ant)
+        return F_RF
+
+# ==========================================
+# 2. TSC-DL 扇区聚焦训练
+# ==========================================
+model = SectorHBFNet().to(device)
+optimizer = optim.Adam(model.parameters(), lr=0.002)
+
+Batch_Size, Epochs = 500, 1500   #训练次数 
+loss_history = []
+
+print("启动 TSC-DL 训练")
+start_time = time.time()
+for epoch in range(Epochs):
+    model.train()
+    optimizer.zero_grad()
+    
+    # 【与 TSC-OMP 呼应的核心逻辑】：扇区聚焦训练！
+    # 不再全空域瞎猜，而是在目标和杂波的局部战术扇区（±3°）内生成海量数据进行特训
+    theta_t = test_t_ang + (torch.rand(Batch_Size, device=device) - 0.5) * 6.0
+    theta_c = test_c_ang + (torch.rand(Batch_Size, device=device) - 0.5) * 6.0
+    
+    a_t, a_c = gen_a_batch(theta_t), gen_a_batch(theta_c)
+    
+    # 构建纯净的基准协方差矩阵
+    R_interf = noise_power * torch.eye(N_ant, device=device).unsqueeze(0) + INR_linear * (a_c @ a_c.mH)
+    w_opt = torch.linalg.pinv(R_interf) @ a_t
+    w_opt = w_opt / torch.norm(w_opt, dim=1, keepdim=True)
+    
+    F_RF = model(w_opt)
+    F_BB = torch.linalg.pinv(F_RF.mH @ F_RF + 1e-3*torch.eye(N_RF, device=device).unsqueeze(0)) @ (F_RF.mH @ w_opt)
+    w_ai = F_RF @ F_BB
+    w_ai = w_ai / torch.norm(w_ai, dim=1, keepdim=True)
+    
+    # 波束相似度
+    sim = torch.abs(torch.sum(w_opt.conj() * w_ai, dim=(1,2)))**2
+    loss_sim = torch.mean(1.0 - sim)
+    # 杂波方向绝对抑制
+    loss_null = torch.mean(torch.abs(torch.sum(a_c.conj() * w_ai, dim=(1,2)))**2)
+    
+    # 【重拳出击】：给予杂波方向 150 倍的极限死亡惩罚！逼迫 AI 挖出深坑！
+    loss = loss_sim + 150.0 * loss_null 
+    loss.backward()
+    optimizer.step()
+    loss_history.append(loss.item())
+
+print(f"训练完成！耗时: {time.time() - start_time:.2f} 秒\n")
+
+# ==========================================
+# 3. 静态成果检验：生成完美的波束方向图
+# ==========================================
+print(f"正在生成目标 {test_t_ang}°，强杂波 {test_c_ang}° 的终极抗干扰方向图...")
+theta_scan = torch.arange(-90, 90.1, 0.1, device=device)
+A_scan = gen_a_batch(theta_scan).squeeze(2).T 
+
+a_t_stat = gen_a_batch(torch.tensor([test_t_ang], device=device))
+a_c_stat = gen_a_batch(torch.tensor([test_c_ang], device=device))
+R_interf_stat = noise_power * torch.eye(N_ant, device=device).unsqueeze(0) + INR_linear * (a_c_stat @ a_c_stat.mH)
+
+with torch.no_grad():
+    model.eval()
+    w_opt_stat = torch.linalg.pinv(R_interf_stat) @ a_t_stat
+    w_opt_stat /= torch.norm(w_opt_stat, dim=1, keepdim=True)
+    
+    F_RF_stat = model(w_opt_stat)
+    F_BB_stat = torch.linalg.pinv(F_RF_stat.mH @ F_RF_stat + 1e-3*torch.eye(N_RF, device=device).unsqueeze(0)) @ (F_RF_stat.mH @ w_opt_stat)
+    w_ai_stat = F_RF_stat @ F_BB_stat
+    w_ai_stat /= torch.norm(w_ai_stat, dim=1, keepdim=True)
+    
+    # 极简且精准的点乘测向
+    gain_dig = torch.abs(torch.sum(w_opt_stat.conj() * A_scan.unsqueeze(0), dim=1))**2
+    BP_dig = 10 * torch.log10(gain_dig).cpu().numpy().flatten()
+    BP_dig -= np.max(BP_dig)
+    
+    gain_ai = torch.abs(torch.sum(w_ai_stat.conj() * A_scan.unsqueeze(0), dim=1))**2
+    BP_ai = 10 * torch.log10(gain_ai).cpu().numpy().flatten()
+    BP_ai -= np.max(BP_ai)
+
+# ==========================================
+# 4. 绘制终极展示双图
+# ==========================================
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
+
+ax1.plot(loss_history, color='#0072BD', linewidth=2)
+ax1.set_title('扇区聚焦网络 (TSC-DL) 极限优化损失曲线', fontweight='bold', fontsize=12)
+ax1.set_xlabel('迭代轮数 (Epochs)')
+ax1.set_ylabel('联合 Loss')
+ax1.grid(True, ls='--')
+
+ax2.plot(theta_scan.cpu().numpy(), BP_dig, color='#000000', linestyle='-', linewidth=2, alpha=0.5, label='纯数字理想波束')
+ax2.plot(theta_scan.cpu().numpy(), BP_ai, color='#D95319', linestyle='--', linewidth=2.5, label='扇区聚焦 AI 混合波束')
+ax2.axvline(test_t_ang, color='g', linestyle='-.', label=f'机动目标 ({test_t_ang}°)')
+ax2.axvline(test_c_ang, color='r', linestyle='-.', label=f'强杂波 ({test_c_ang}°)')
+ax2.set_ylim([-60, 0])
+ax2.set_xlim([-60, 60])
+ax2.set_title(f'深度学习端到端波束响应 (目标 {test_t_ang}° 杂波 {test_c_ang}° )', fontweight='bold', fontsize=12)
+ax2.set_xlabel('角度 θ (°)')
+ax2.set_ylabel('归一化增益 (dB)')
+ax2.grid(True, ls='--')
+ax2.legend(loc='lower center', fontsize=10)
+
+plt.tight_layout()
+plt.savefig('DL_Sector_Focus_Final.png', dpi=300)
+plt.show()

Graduation Design/main_nula.py

@@ -0,0 +1,292 @@
+import torch
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+
+# ==========================================
+# 0. 全局设置与 CUDA 引擎
+# ==========================================
+plt.rcParams['axes.unicode_minus'] = False
+plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
+device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+
+# 锁定随机数种子，对标 MATLAB 的 rng(2024)
+torch.manual_seed(2024)
+np.random.seed(2024)
+
+# 颜色设置对标 MATLAB
+c_ana = '#77AC30'
+c_dig = '#000000'
+c_dft = '#4DBEEE'
+c_omp = '#D95319'
+c_tsc = '#0072BD'
+
+# ==========================================
+# 1. 物理层参数初始化 (对标 setpar.m)
+# ==========================================
+N_ant, N_RF = 24, 4
+target_angle = 5.0
+clutter_angle = -15.0
+L_snapshots = 100
+noise_power = 1.0
+INR_linear = 10**(30/10)
+
+fc = 4.9e9
+lam = 3e8 / fc
+# 完全复刻 MATLAB 中的非均匀阵列 (NULA) 物理间隔
+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))  # Shape: (24,)
+delta_d_tensor = torch.tensor(delta_d, dtype=torch.float32, device=device).unsqueeze(1) # Shape: (24, 1)
+
+def gen_a(theta_deg):
+    """对标 gen_a.m, 注意此处的 -1j 以及物理位置 delta_d"""
+    theta_rad = torch.deg2rad(torch.tensor(theta_deg, dtype=torch.float32, device=device))
+    return torch.exp(-1j * 2 * torch.pi * delta_d_tensor * torch.sin(theta_rad) / lam)
+
+a_target = gen_a(target_angle)
+a_clutter = gen_a(clutter_angle)
+
+# ==========================================
+# 2. 动态生成目标约束扇区与字典
+# ==========================================
+theta_coarse_prior = round(target_angle)
+window_half_width = 5.0
+
+dic_angles_constrained = torch.arange(theta_coarse_prior - window_half_width, theta_coarse_prior + window_half_width + 0.1, 0.1, device=device)
+A_dic_constrained = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in dic_angles_constrained], dim=1)
+
+dic_angles_full = torch.arange(-90, 90.5, 0.5, device=device)
+A_dic_full = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in dic_angles_full], dim=1)
+
+angles_dft = torch.linspace(theta_coarse_prior - window_half_width, theta_coarse_prior + window_half_width, N_RF, device=device)
+F_RF_dft = torch.cat([gen_a(th) / np.sqrt(N_ant) for th in angles_dft], dim=1)
+
+# ==========================================
+# 3. 提取防发散的慢时间先验射频矩阵 F_RF
+# ==========================================
+print('正在进行慢时间 (Slow-Time) 模拟射频网络硬件配置...')
+a_prior_broad = torch.zeros((N_ant, 1), dtype=torch.complex64, device=device)
+for theta_b in torch.arange(theta_coarse_prior - window_half_width, theta_coarse_prior + window_half_width + 0.1, 0.5, device=device):
+    a_prior_broad += gen_a(theta_b)
+a_prior_broad /= torch.norm(a_prior_broad)
+
+R_interference_prior = noise_power * torch.eye(N_ant, device=device) + INR_linear * (a_clutter @ a_clutter.mH)
+w_digital_prior = torch.linalg.inv(R_interference_prior) @ a_prior_broad
+w_digital_prior /= torch.norm(w_digital_prior)
+
+def extract_omp(w_opt, A_dic):
+    res = w_opt.clone()
+    F_RF_mat = None
+    for k in range(N_RF):
+        proj = A_dic.mH @ res
+        idx = torch.argmax(torch.abs(proj))
+        new_atom = A_dic[:, idx:idx+1]
+        F_RF_mat = new_atom if F_RF_mat is None else torch.cat((F_RF_mat, new_atom), dim=1)
+        res = w_digital_prior - F_RF_mat @ (torch.linalg.pinv(F_RF_mat) @ w_digital_prior)
+    return F_RF_mat
+
+F_RF_omp_full_static = extract_omp(w_digital_prior, A_dic_full)
+F_RF_omp_static = extract_omp(w_digital_prior, A_dic_constrained)
+
+# ==========================================
+# 4. 计算波束方向图与理论 SINR
+# ==========================================
+print('正在计算静态波束方向图与理论 SINR...')
+w_analog = a_target / np.sqrt(N_ant)
+
+def calc_hybrid_weight(F_RF, R_in_stat):
+    a_eff = F_RF.mH @ a_target
+    R_eff = F_RF.mH @ R_in_stat @ F_RF
+    R_eff += (0.05 * torch.trace(R_eff).real / N_RF) * torch.eye(N_RF, device=device)
+    inv_R_eff = torch.linalg.inv(R_eff)
+    F_BB = (inv_R_eff @ a_eff) / (a_eff.mH @ inv_R_eff @ a_eff)
+    w_hyb = F_RF @ F_BB
+    return w_hyb / torch.norm(w_hyb)
+
+
+#随机数
+#sig_power_static = 10**(0/10)
+#s_t_st = np.sqrt(sig_power_static/2) * (torch.randn(1, L_snapshots, device=device) + 1j*torch.randn(1, L_snapshots, device=device))
+#s_c_st = np.sqrt(INR_linear/2) * (torch.randn(1, L_snapshots, device=device) + 1j*torch.randn(1, L_snapshots, device=device))
+#n_t_st = np.sqrt(noise_power/2) * (torch.randn(N_ant, L_snapshots, device=device) + 1j*torch.randn(N_ant, L_snapshots, device=device))
+#X_static = a_target @ s_t_st + a_clutter @ s_c_st + n_t_st
+#R_in_static = (X_static @ X_static.mH) / L_snapshots
+
+sig_power_static = 10**(0/10)
+# 直接使用无穷大快拍下的理论理想协方差矩阵，消除所有随机扰动
+R_in_static = noise_power * torch.eye(N_ant, device=device) + \
+              sig_power_static * (a_target @ a_target.mH) + \
+              INR_linear * (a_clutter @ a_clutter.mH)
+
+w_hybrid_dft = calc_hybrid_weight(F_RF_dft, R_in_static)
+w_hybrid_omp_full = calc_hybrid_weight(F_RF_omp_full_static, R_in_static)
+w_hybrid_omp = calc_hybrid_weight(F_RF_omp_static, R_in_static)
+
+R_in_static_dl = R_in_static + (0.01 * torch.trace(R_in_static).real / N_ant) * torch.eye(N_ant, device=device)
+inv_R_stat = torch.linalg.inv(R_in_static_dl)
+w_digital_ideal = (inv_R_stat @ a_target) / (a_target.mH @ inv_R_stat @ a_target)
+w_digital_ideal /= torch.norm(w_digital_ideal)
+
+theta_scan = torch.arange(-90, 90.1, 0.1, device=device)
+A_scan_full = torch.cat([gen_a(th) for th in theta_scan], dim=1)
+
+def get_bp(w): return 10 * torch.log10((torch.abs(A_scan_full.mH @ w)**2).squeeze() / torch.max(torch.abs(A_scan_full.mH @ w)**2))
+BP_Ana = get_bp(w_analog).cpu().numpy()
+BP_Dig = get_bp(w_digital_ideal).cpu().numpy()
+BP_DFT = get_bp(w_hybrid_dft).cpu().numpy()
+BP_OMP_FULL = get_bp(w_hybrid_omp_full).cpu().numpy()
+BP_OMP = get_bp(w_hybrid_omp).cpu().numpy()
+theta_scan_np = theta_scan.cpu().numpy()
+
+SNR_dB_range = np.arange(-10, 22, 2)
+SINR_Ana, SINR_Dig, SINR_DFT, SINR_OMP_FULL, SINR_OMP = [], [], [], [], []
+for snr in SNR_dB_range:
+    sig_p = 10**(snr/10)
+    def calc_sinr(w): return (sig_p * torch.abs(w.mH @ a_target)**2) / torch.real(w.mH @ R_interference_prior @ w)
+    SINR_Ana.append(10 * torch.log10(calc_sinr(w_analog)).item())
+    SINR_Dig.append(10 * torch.log10(calc_sinr(w_digital_ideal)).item())
+    SINR_DFT.append(10 * torch.log10(calc_sinr(w_hybrid_dft)).item())
+    SINR_OMP_FULL.append(10 * torch.log10(calc_sinr(w_hybrid_omp_full)).item())
+    SINR_OMP.append(10 * torch.log10(calc_sinr(w_hybrid_omp)).item())
+
+# ==========================================
+# 5. 基于张量 Batch 操作的极致蒙特卡洛评估
+# ==========================================
+print('正在启动 CUDA 高维张量批量计算 (Batched Evaluation)...')
+Monte_Carlo = 500
+scan_grid = torch.arange(theta_coarse_prior - window_half_width, theta_coarse_prior + window_half_width + 0.01, 0.01, device=device)
+A_scan = torch.cat([gen_a(th) for th in scan_grid], dim=1) 
+
+A_eff_dft = F_RF_dft.mH @ A_scan               
+A_eff_omp_full = F_RF_omp_full_static.mH @ A_scan
+A_eff_omp = F_RF_omp_static.mH @ A_scan
+
+RMSE_Ana, RMSE_Dig, RMSE_DFT, RMSE_OMP_FULL, RMSE_OMP = [], [], [], [], []
+start_time = time.time()
+
+for snr in SNR_dB_range:
+    sig_power = 10**(snr/10)
+    s_t = np.sqrt(sig_power/2) * (torch.randn(Monte_Carlo, 1, L_snapshots, device=device) + 1j*torch.randn(Monte_Carlo, 1, L_snapshots, device=device))
+    s_c = np.sqrt(INR_linear/2) * (torch.randn(Monte_Carlo, 1, L_snapshots, device=device) + 1j*torch.randn(Monte_Carlo, 1, L_snapshots, device=device))
+    n_t = np.sqrt(noise_power/2) * (torch.randn(Monte_Carlo, N_ant, L_snapshots, device=device) + 1j*torch.randn(Monte_Carlo, N_ant, L_snapshots, device=device))
+    
+    X_mc = a_target.unsqueeze(0) * s_t + a_clutter.unsqueeze(0) * s_c + n_t
+    R_in = torch.bmm(X_mc, X_mc.mH) / L_snapshots
+    trace_R_in = torch.diagonal(R_in, dim1=-2, dim2=-1).sum(-1).real
+    R_in_dl = R_in + (0.05 * trace_R_in / N_ant).view(-1, 1, 1) * torch.eye(N_ant, device=device).unsqueeze(0)
+    inv_R_in_dl = torch.linalg.inv(R_in_dl)
+    
+    P_ana = torch.einsum('si, bij, js -> bs', A_scan.conj().T, R_in, A_scan).abs()
+    P_dig = torch.einsum('si, bij, js -> bs', A_scan.conj().T, inv_R_in_dl, A_scan).abs()
+    
+    def process_hybrid_batched(F_RF, A_eff_hyb):
+        Y = F_RF.mH.unsqueeze(0) @ X_mc
+        R_hyb = torch.bmm(Y, Y.mH) / L_snapshots
+        trace_R = torch.diagonal(R_hyb, dim1=-2, dim2=-1).sum(-1).real
+        R_hyb += (0.05 * trace_R / N_RF).view(-1, 1, 1) * torch.eye(N_RF, device=device).unsqueeze(0)
+        inv_R_hyb = torch.linalg.inv(R_hyb)
+        P_hyb = torch.einsum('si, bij, js -> bs', A_eff_hyb.conj().T, inv_R_hyb, A_eff_hyb).abs()
+        return scan_grid[torch.argmin(P_hyb, dim=1)] 
+
+    ang_ana = scan_grid[torch.argmax(P_ana, dim=1)]
+    ang_dig = scan_grid[torch.argmin(P_dig, dim=1)]
+    ang_dft = process_hybrid_batched(F_RF_dft, A_eff_dft)
+    ang_omp_full = process_hybrid_batched(F_RF_omp_full_static, A_eff_omp_full)
+    ang_omp = process_hybrid_batched(F_RF_omp_static, A_eff_omp)
+
+    RMSE_Ana.append(torch.sqrt(torch.mean((ang_ana - target_angle)**2)).item())
+    RMSE_Dig.append(torch.sqrt(torch.mean((ang_dig - target_angle)**2)).item())
+    RMSE_DFT.append(torch.sqrt(torch.mean((ang_dft - target_angle)**2)).item())
+    RMSE_OMP_FULL.append(torch.sqrt(torch.mean((ang_omp_full - target_angle)**2)).item())
+    RMSE_OMP.append(torch.sqrt(torch.mean((ang_omp - target_angle)**2)).item())
+
+print(f"CUDA 并行测角完毕！耗时: {time.time() - start_time:.4f} 秒。")
+
+# ==========================================
+# 6. 绘图：单图输出与汇总输出
+# ==========================================
+bp_data = [BP_Ana, BP_Dig, BP_DFT, BP_OMP_FULL, BP_OMP]
+bp_labels = ['纯模拟架构', '纯数字架构(性能上限)', 'DFT混合架构', '传统OMP混合架构', 'TSC-OMP混合架构(所提)']
+bp_colors = [c_ana, c_dig, c_dft, c_omp, c_tsc]
+bp_ls = ['-.', '-', '--', ':', '-']
+
+# 6.1 分别绘制 5 张独立的波束方向图
+for i in range(5):
+    fig, ax = plt.subplots(figsize=(8, 6))
+    ax.plot(theta_scan_np, bp_data[i], color=bp_colors[i], linestyle=bp_ls[i], linewidth=2.5, label='主波束')
+    ax.axvline(x=target_angle, color='g', linestyle='--', linewidth=2, label=f'机动目标 ({target_angle}°)')
+    ax.axvline(x=clutter_angle, color='r', linestyle='--', linewidth=2, label=f'强杂波 ({clutter_angle}°)')
+    ax.set_ylim([-60, 0])
+    ax.set_xlim([-60, 60])
+    ax.grid(True, linestyle='--', alpha=0.6)
+    ax.set_title(f'波束方向图 - {bp_labels[i]}', fontsize=14, fontweight='bold')
+    ax.set_xlabel(r'角度 $\theta (^\circ)$', fontsize=12)
+    ax.set_ylabel('归一化增益 (dB)', fontsize=12)
+    
+    # 调整图例到下方
+    plt.tight_layout()
+    fig.subplots_adjust(bottom=0.22)
+    ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12), ncol=3, fontsize=10)
+    
+    file_name = f'BP_Single_{i+1}_{bp_labels[i].split("(")[0]}.png'
+    plt.savefig(file_name, dpi=300, bbox_inches='tight')
+
+# 6.2 绘制五大架构综合对比波束图
+fig_all, ax_all = plt.subplots(figsize=(10, 6))
+for i in range(5):
+    ax_all.plot(theta_scan_np, bp_data[i], color=bp_colors[i], linestyle=bp_ls[i], linewidth=2.5, label=bp_labels[i].split("(")[0])
+ax_all.axvline(x=target_angle, color='g', linestyle='--', linewidth=2)
+ax_all.axvline(x=clutter_angle, color='r', linestyle='--', linewidth=2)
+ax_all.text(target_angle + 1, -5, f'机动目标 ({target_angle}°)', color='g', fontweight='bold')
+ax_all.text(clutter_angle - 1, -5, f'强杂波 ({clutter_angle}°)', color='r', fontweight='bold', ha='right')
+ax_all.set_ylim([-60, 0])
+ax_all.set_xlim([-60, 60])
+ax_all.grid(True, linestyle='--', alpha=0.6)
+ax_all.set_title('五大架构波束方向图综合对比', fontsize=14, fontweight='bold')
+ax_all.set_xlabel(r'角度 $\theta (^\circ)$', fontsize=12)
+ax_all.set_ylabel('归一化增益 (dB)', fontsize=12)
+
+plt.tight_layout()
+fig_all.subplots_adjust(bottom=0.22)
+ax_all.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12), ncol=5, fontsize=10)
+plt.savefig('BP_All_Combined.png', dpi=300, bbox_inches='tight')
+
+# 6.3 绘制 SINR 对比图
+fig_sinr, ax_sinr = plt.subplots(figsize=(10, 6))
+ax_sinr.plot(SNR_dB_range, SINR_Ana, '-v', color=c_ana, linewidth=2, label='纯模拟')
+ax_sinr.plot(SNR_dB_range, SINR_Dig, '-o', color=c_dig, linewidth=2, label='纯数字(上限)')
+ax_sinr.plot(SNR_dB_range, SINR_DFT, '-d', color=c_dft, linewidth=2, label='DFT')
+ax_sinr.plot(SNR_dB_range, SINR_OMP_FULL, '-*', color=c_omp, linewidth=2, label='传统OMP')
+ax_sinr.plot(SNR_dB_range, SINR_OMP, '-s', color=c_tsc, linewidth=2.5, label='TSC-OMP(所提)')
+ax_sinr.grid(True, linestyle='--', alpha=0.6)
+ax_sinr.set_title(f'全架构输出 SINR 理论逼近测试 (目标 {target_angle}°)', fontsize=14, fontweight='bold')
+ax_sinr.set_xlabel('输入 SNR (dB)', fontsize=12)
+ax_sinr.set_ylabel('输出 SINR (dB)', fontsize=12)
+
+plt.tight_layout()
+fig_sinr.subplots_adjust(bottom=0.22)
+ax_sinr.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12), ncol=5, fontsize=10)
+plt.savefig('SINR_All_Combined.png', dpi=300, bbox_inches='tight')
+
+# 6.4 绘制 RMSE 对比图
+fig_rmse, ax_rmse = plt.subplots(figsize=(10, 6))
+ax_rmse.semilogy(SNR_dB_range, RMSE_Ana, '-.v', color=c_ana, linewidth=2, label='纯模拟')
+ax_rmse.semilogy(SNR_dB_range, RMSE_Dig, '-o', color=c_dig, linewidth=2, label='纯数字(上限)')
+ax_rmse.semilogy(SNR_dB_range, RMSE_DFT, '--d', color=c_dft, linewidth=2, label='DFT')
+ax_rmse.semilogy(SNR_dB_range, RMSE_OMP_FULL, ':*', color=c_omp, linewidth=2.5, label='传统OMP')
+ax_rmse.semilogy(SNR_dB_range, RMSE_OMP, '-s', color=c_tsc, linewidth=2.5, label='TSC-OMP (所提)')
+ax_rmse.grid(True, which="both", ls="--", alpha=0.6)
+ax_rmse.set_ylim([1e-3, 20])
+ax_rmse.set_title(f'全架构实战恶劣环境 DOA 测角精度评估 (机动目标 {target_angle}°)', fontsize=14, fontweight='bold')
+ax_rmse.set_xlabel('信噪比 SNR (dB)', fontsize=12)
+ax_rmse.set_ylabel(r'测角 RMSE $(^\circ)$', fontsize=12)
+
+plt.tight_layout()
+fig_rmse.subplots_adjust(bottom=0.22)
+ax_rmse.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12), ncol=5, fontsize=10)
+plt.savefig('RMSE_All_Combined.png', dpi=300, bbox_inches='tight')
+
+plt.show()
+print('计算完毕，单图与综合图表均已独立生成。')