641 lines
25 KiB
Python
641 lines
25 KiB
Python
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) |