engee.clear() # 清除工作区
is_dinamic_plot = true # 选择图形类型：true(plotlyjs())-interactive,false(gr())-static
is_dinamic_plot ? plotlyjs() : gr()

let
    installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
    list_packages = ["JLD2","Hungarian"]
    for pack in list_packages
        pack in installed_packages || Pkg.add(pack)
    end
end

using LinearAlgebra,Hungarian,DataFrames
import JLD2

include("utils/plot_sig_and_spec.jl") # 建立示波器和频谱图的功能
include("utils/viz_scene.jl") # 构建雷达运行场景的功能
include("utils/calc_plot_rcs.jl") # 反向散射图构造函数
include("utils/sort_ang_mse.jl") # 排序轴承估计值的功能
include("utils/plot_spec_azel.jl") # 构造空间谱的功能
include("utils/calc_table_ang_analysis.jl") # 用于精度分析的表计算功能

# 基本雷达参数
fc = 1e9 # [Hz]，载波频率
c = physconst("LightSpeed") # [m/s]，信号的传播速度
lambda = c / fc; # 波长，m

println("载波频率：∞(FC*1E-9)GHz")
println("光速：＜c m/s")
println("波长：＜lambda m")

# ---信号参数 ---
bw = 100e6; # LFM频段，Hz
T = 10e-6; # 脉冲持续时间，s
Q = 2 # 脉冲连续性
prf = round(1 / (Q * T);sigdigits = 15); # 脉冲重复率（必要时舍入）
type_dir = "Up" # 频率变化方向["向上","向下"]
type_sweep = "Positive"  # 频率变化类型["正"，"对称"]
n_pulses = 2 # 累积脉冲数
n_groups_imps = 2 # 脉冲突发数

fs = 2*bw; # 采样率，Hz
NBand = 64; # 频率子带数

T_coh = n_pulses/prf # coregent积累的时间
ΔR = c/(2*bw) # 范围分辨率，m
ΔV = lambda/(2*T_coh) # 速度分辨率，m/s
Rmin = c*T/2 # 最小一位数范围，m
Rmax = c/(2*prf) # 最大一位数范围，m
Vmax = lambda*prf/4 # 最大个位数速度，m/s

# lfm信号发生器的系统对象的形成
waveform = EngeePhased.LinearFMWaveform( 
    SampleRate = fs, # 采样率，Hz
    PulseWidth = T, # 脉冲持续时间
    PRF = prf, # 脉冲重复率
    SweepBandwidth = bw,  # LCHM乐队
    SweepDirection = "Up", # 频率变化方向
    SweepInterval = "Positive", # 锯型，非对称或对称
    OutputFormat = "Pulses", # 输出格式类型[]
    NumPulses = n_pulses # 脉冲数
);

# LFM信号的计算
sig_lfm = waveform();

# 示波器和频谱图的构建
plot_sig_and_spec(sig_lfm;fs=fs,name_sig = "LCHM");

println("采样率：∞(round(fs*1e-6;sigdigits=5))MHz")
println("脉冲持续时间：△(round(T*1e6;sigdigits=5))微秒")
println("随访期：$(round(1e6/prf;sigdigits=5))mks")
println("重复率：△(round(prf*1e-3;sigdigits=5))kHz")
println("钻孔:$(圆(Q;sigdigits=5))")
println("相干累积时间：△(round(T_coh*1e3;sigdigits=5))ms")
println("每包脉冲数：$n_pulses")
println("包数：$n_groups_imps")
println("范围分辨率：△(round(ΔR;sigdigits=3))m") 
println("速度分辨率：△(round(ΔV;sigdigits=3))m/s") 
println("最大一位数范围：$(round(Rmax;sigdigits=3))m") 
println("最小一位数范围：$(round(Rmin;sigdigits=3))m") 
println("最大个位数速度:∞(round(Vmax;sigdigits=3))m/s") 

# 天线参数
d = lambda /2 # 元素之间的距离
freq_range = [fc-bw fc+bw] # 辐射的频率范围
Nx = 24 # X轴上的元素数（方位角分辨率）
Ny = 24 # Y轴上的元素数（按位置角度的分辨率）
Δaz = rad2deg(lambda/(d*(Nx-1))) # 方位分辨率，deg
Δel = rad2deg(lambda/(d*(Ny-1))) # 座椅角度分辨率，deg
# 元素总数N=Nx*Ny

# 具有各向同性天线元件
ant_elem = EngeePhased.IsotropicAntennaElement( 
    FrequencyRange = freq_range, # 频率范围
    BackBaffled = true # 去除天线底部的反向散射
) 

# 创建矩形大灯
array = EngeePhased.URA(
    Element = ant_elem, # 天线元件
    Size = [Nx Ny], # 网格大小
    ElementSpacing = [d d], # 元素之间的距离
);

pattern(array,fc) |> display;

println("元素之间的距离：$(round(d;sigdigits=3))m") 
println("矩形大灯尺寸:$(Nx)x$Ny[Nx,Ny]") 
println("辐射频率范围：∞(freq_range*1e-9)GHz")
println("方位角分辨率：△(round(Δaz;sigdigits=3))°") 
println("位置角分辨率：△(round(Δel;sigdigits=3))°") 

# 接收和发送路径的参数
Pt = 1e3 # W、峰值功率
Gtx = 20 # dB，发射天线的增益
LFtx = 1.5 # dB，传输路径中的损耗

Grx = 20 # dB，接收天线的增益
NF = 2 # dB，接收路径的噪声系数
LFrx = 1.5 # dB，接收路径中的损耗
TN = 290 # 噪音温度
NP = physconst("boltzmann") * bw * TN; # W、噪声功率

# 带有发射器的生物
transmitter = EngeePhased.Transmitter( 
    Gain = Gtx, # 发射机增益
    PeakPower = Pt, # 峰值功率
    LossFactor=LFtx # 传输路径中的损耗
)

# 具有宽带传输灯的设备
radiator = EngeePhased.WidebandRadiator(
    Sensor = array, # 天线元件的底部
    SampleRate = fs, # 采样率
    CarrierFrequency = fc, # 载波频率
    NumSubbands = NBand # 频率子带数
)

# 具有宽带接收头灯的设备
collector = EngeePhased.WidebandCollector(
    Sensor = array, # 天线元件的底部
    SampleRate = fs, # 采样率
    CarrierFrequency = fc, # 中央频率
    NumSubbands = NBand # 频率子带数
)

# 创建接收前置放大器(LNA)
receiver = EngeePhased.ReceiverPreamp(
    NoiseMethod="Noise power", # 噪音设定方法["噪音功率"、"噪音温度"]
    Gain = Grx, # 接收器增益，dB
    NoisePower = NP, # 噪音功率
    LossFactor=LFrx, # 接收道中的损失
    NoiseFigure = NF, # 路径中的噪声因子
    SampleRate = fs # 采样率
);

println("-----------------发射器-----------------------------")
println("发射机功率：△(round(Pt*1e-3;sigdigits=5))kW")
println("发射机增益：△(Round(Gtx;sigdigits=5))dB")
println("传输路径损耗：△(round(LFtx;sigdigits=5))dB")
println("--------------------------------------------------------")
println("-----------------接收器-----------------------------")
println("接收器增益：△(round(Gtx;sigdigits=5))dB")
println("接收路径中的损耗：△(round(LFtx;sigdigits=5))dB")
println("噪声温度：△(round(LFrx;sigdigits=5))dB")
println("噪声系数：△(round(nf;sigdigits=5))dB")
println("噪声温度：△(round(TN;sigdigits=5))°C")
println("噪声功率：△(round(NP;sigdigits=5))W")

az_scan = -20:0.5:20 # 方位扫描范围
el_scan = 0:0.5:30 # 按位置角度扫描范围
n_sig = 4 # 评分数（进球数）

algorithm = "MUSIC" # @param ["Beamscan","MVDR","MUSIC"]

 beamscan = EngeePhased.BeamscanEstimator2D( 
    SensorArray = array, # 天线
    OperatingFrequency = fc, # 载波频率
    AzimuthScanAngles = az_scan,  # 方位扫描范围
    ElevationScanAngles = el_scan,  # 按位置角度扫描范围
    DOAOutputPort = true, # 启用轴承估计值的输出
    NumSignals = n_sig # 评级数目
);

 mvdr = EngeePhased.MVDREstimator2D( 
    SensorArray = array, # 天线
    OperatingFrequency = fc, # 载波频率
    AzimuthScanAngles = az_scan,  # 方位扫描范围
    ElevationScanAngles = el_scan,  # 按位置角度扫描范围
    DOAOutputPort = true, # 启用轴承估计值的输出
    NumSignals = n_sig # 评级数目
);

music =  EngeePhased.MUSICEstimator2D(
    SensorArray = array, # 天线
    OperatingFrequency = fc, # 载波频率
    AzimuthScanAngles = az_scan,  # 方位扫描范围
    ElevationScanAngles = el_scan,  # 按位置角度扫描范围
    DOAOutputPort = true, # 启用轴承估计值的输出
    NumSignalsSource="Property",
    NumSignals = n_sig, # 评级数目
    ForwardBackwardAveraging=false # 计算正向和反向平均
);

doa_alg =  if algorithm == "Beamscan"
    println("算法：Beamscan")
    deepcopy(beamscan)
elseif algorithm == "MVDR"
    println("算法：MVDR");
    deepcopy(mvdr)
elseif algorithm == "MUSIC"
    println("算法：音乐");  
    deepcopy(music)
end;

# 分布环境参数
is_atmosphere = true # 大气对电磁波传播的影响的考虑
T_ch = 15 # 大气温度，°C
press_air = 101325.0 # 大气压力，PA
vap_density = 7.5 # 湿蒸汽密度
liq_density = 0.0 # 冷凝水密度
rain_rate = 0.0 # 降水强度

channel_fwd = EngeePhased.WidebandLOSChannel( 
    PropagationSpeed = c,
    SampleRate = fs,
    OperatingFrequency = fc,
    NumSubbands = NBand,
    TwoWayPropagation = false, # 考虑双向信号传播
    SpecifyAtmosphere = is_atmosphere,
    Temperature = T_ch,
    DryAirPressure = press_air,
    WaterVapourDensity = vap_density,
    LiquidWaterDensity = liq_density,
    RainRate = rain_rate
);

channel_bwd = deepcopy(channel_fwd); # 从目标到雷达的分配通道（相同）

data = JLD2.load("quad_rcs_pattern.jld2")

num_tgts = n_sig # 目标数量（设置相等数量的测向器评级）
rcs_az = data["az"] # 方位EPR记帐范围(-180:180),deg
rcs_el = data["el"] # EPR会计范围按座位角度(-90:90),deg
rcs_pattern = data["rcs"] # 目标门矩阵[181x361],m^2

type_tgt = "Nonfluctuating" # 目标类型["非波动","Swerling1","Swerling2","Swerling3","Swerling4"]

# 宽带SDR目标的形成
rcs_targets = EngeePhased.WidebandBackscatterRadarTarget( 
    AzimuthAngles = rcs_az, # 基于方位的EPR记帐范围
    ElevationAngles = rcs_el, # 按座位角度计算的EPR会计范围
    RCSPattern = rcs_pattern, # DOR目标
    PropagationSpeed = c, # 传播速度，m/s
    OperatingFrequency = fc, # 载波频率，Hz
    NumSubbands = NBand, # 频率子带数
    Model = type_tgt # 目标EPR的波动类型
);

text_title = "可视化和DOR目标"
k_norm  = 0.6 # 门尺寸缩放系数
calc_plot_rcs(rcs_pattern,rcs_az,rcs_el;title = text_title,k = k_norm)

# 目标移动模型
targets_pos = [
            [1000; 100; 200];; # 1目标
            [1200; -10; 200];; # 2目标
            [1100; -150; 150];; # 3目标
            [950; 50; 320];; # 4目标
] # 目标的初始位置，m
targets_vel = [
            [-130; 0; 0];; # 第1目标的速度（在平行于x轴的雷达上）
            [-150; 0; 0];; # 第二目标的速度（在平行于x轴的雷达上）
            [-200; 0; 0];; # 第3目标的速度（在平行于x轴的雷达上）
            [-230; 0; 0];; # 第4个目标的速度（在平行于x轴的雷达上）
] # 目标速度，m/s
type_traj = "Velocity" # 运动模型是均匀的["速度"，"加速度"，"自定义"]

targets_platform = EngeePhased.Platform(
    MotionModel = type_traj,
    InitialPosition = targets_pos, # 起始位置
    Velocity = targets_vel # 目标速度
)

# 雷达运动模型
radar_pos = [0; 0; 0] # 雷达位置
radar_vel = [0; 0; 0] # 雷达速度

radar_platform = EngeePhased.Platform(
    InitialPosition = radar_pos, # 雷达位置
    Velocity = radar_vel # 雷达速度
);

abs_dist,ang_true = rangeangle(targets_pos,radar_pos)

for i in 1:num_tgts
    println("----------------目标#$(i)----------------------------")
    println("方位角：ψ(round(ang_true[1,i];sigdigits=3))°，仰角：ψ(round(ang_true[2,i];sigdigits=3))°")
    println("雷达范围至目标号。$(i):$(round(abs_dist[i]*1e-3;sigdigits=5))km")
end

viz_scene(targets_pos, radar_pos) |> display;

println("-----------------雷达站-----------------------------")
println("载波频率：∞(round(fc*1e-9;sigdigits=3))GHz")
println("光速：＜c m/s")
println("波长:∞(round(lambda;sigdigits=3))m") 
println("频率子带数：△NBand m")
println("雷达位置:$radar_pos[Px,Py,Pz],m")
println("雷达速度:△radar_vel[Vx,Vyc,Vz],m/s")
println("--------------------------------------------------------")
println("-----------------发电机-----------------------------")
println("采样率：∞(round(fs*1e-6;sigdigits=5))MHz")
println("脉冲持续时间：△(round(T*1e6;sigdigits=5))微秒")
println("随访期：$(round(1e6/prf;sigdigits=5))mks")
println("重复率：△(round(prf*1e-3;sigdigits=5))kHz")
println("钻孔:$(圆(Q;sigdigits=5))")
println("相干累积时间：△(round(T_coh*1e3;sigdigits=5))ms")
println("每包脉冲数：$n_pulses")
println("包数：$n_groups_imps")
println("范围分辨率：△(round(ΔR;sigdigits=3))m") 
println("速度分辨率：△(round(ΔV;sigdigits=3))m/s") 
println("最大一位数范围：$(round(Rmax;sigdigits=3))m") 
println("最小一位数范围：$(round(Rmin;sigdigits=3))m") 
println("最大个位数速度:∞(round(Vmax;sigdigits=3))m/s") 
println("--------------------------------------------------------")
println("-----------------发射器-----------------------------")
println("发射机功率：△(round(Pt*1e-3;sigdigits=5))kW")
println("发射机增益：△(Round(Gtx;sigdigits=5))dB")
println("传输路径损耗：△(round(LFtx;sigdigits=5))dB")
println("--------------------------------------------------------")
println("-----------------接收器-----------------------------")
println("接收器增益：△(round(Gtx;sigdigits=5))dB")
println("接收路径中的损耗：△(round(LFtx;sigdigits=5))dB")
println("噪声温度：△(round(LFrx;sigdigits=5))dB")
println("噪声系数：△(round(nf;sigdigits=5))dB")
println("噪声温度：△(round(TN;sigdigits=5))°C")
println("噪声功率：△(round(NP;sigdigits=5))W")
println("--------------------------------------------------------")
println("-----------------天线-------------------------------")
println("元素之间的距离：$(round(d;sigdigits=3))m") 
println("矩形大灯尺寸:$(Nx)x$Ny[Nx,Ny]") 
println("辐射频率范围：∞(freq_range*1e-9)GHz")
println("方位角分辨率：△(round(Δaz;sigdigits=3))°") 
println("位置角分辨率：△(round(Δel;sigdigits=3))°") 
println("--------------------------------------------------------")
println("-----------------测向仪-----------------------------")
println("方位扫描范围：△(az_scan)deg.")
println("按位置角度扫描范围：△(el_scan)deg。")
println("评级数目：$n_sig")
println("--------------------------------------------------------")

# 输出信号和轴承估计的内存分配
data_y = zeros(ComplexF64,length(sig_lfm),Nx*Ny,n_groups_imps)
data_ang_true = zeros(2, num_tgts, n_groups_imps); # 目标的真正方向
data_ang_est = zeros(2, num_tgts, n_groups_imps); # 目标的真正方向
data_spec = zeros(length(el_scan),length(az_scan), n_groups_imps); # 测向器输出处视角的光谱数据

len_one_imps = round(Int,length(sig_lfm)/n_pulses) # 1脉冲的计数数
sig_lfm_one_imps = sig_lfm[1:len_one_imps] # 分配1个脉冲
sig_buffer_one_groups = similar(data_y[:,:,1]) # 用于累积脉冲串的缓冲器

println("场景的计算正在运行。..")
@inbounds for i in 1:n_groups_imps
    println("束数△(i)的计算。..")

    @inbounds for k in 1:n_pulses
        println("脉冲数△(k)的计算。..")
        # 以脉冲周期1/prf的增量更新雷达和目标的位置
        radar_pos, radar_vel = radar_platform(1/prf)
        tgts_pos, tgts_vel = targets_platform(1/prf)

        # 计算和维护目标轴承
        _,ang_doa = rangeangle(tgts_pos,radar_pos)
        data_ang_true[:,:,i] .= deepcopy(ang_doa)

        # 周三宽带信号发射
        sig_tx = transmitter(sig_lfm_one_imps) # 传送路径的通过
        sig_rad = radiator(sig_tx,ang_doa) # 向太空辐射信号

        # 从雷达到目标的传播
        sig_fwd = channel_fwd(sig_rad,radar_pos,tgts_pos,radar_vel,tgts_vel)
        # 目标的反思
        sig_refl = rcs_targets(sig_fwd,ang_doa)

        # 从目标到雷达的传播
        sig_ret = channel_bwd(sig_refl,tgts_pos,radar_pos,tgts_vel,radar_vel)

        # 接收道
        y = collector(sig_ret,ang_doa) # 接收灯
        sig_buffer_one_groups[(1:len_one_imps).+(k-1)*len_one_imps,:] .= receiver(y) # 基于接收器热噪声(LNA)的预放大
    end

    # 轴承计算
    data_spec[:,:,i],data_ang_est[:,:,i] = doa_alg(sig_buffer_one_groups); # 评估目标到达的范围和方向
    data_y[:,:,i] .= sig_buffer_one_groups; # 保存接收到的信号
end

println("场景计算完成！")

n_imp = 1 # 分析脉冲的数量
plot_size = (600,600) # 许可
plot_spec_azel(az_scan,el_scan,data_spec[:,:,n_imp];
    title = "测向器输出处的功率谱（算法△(algorithm))",
    size = plot_size
)

sort_ang_est =  sort_ang_mse(data_ang_true,data_ang_est) # 排序的角度估计数据

for k in 1:n_groups_imps
    println("---------------脉冲包号 $(k) -----------------------------")
    for i in 1:num_tgts
        println("   ----------------目标#$(i) ----------------------------")
        println("       目标：方位角=$(round(data_ang_true[1,i,k];sigdigits=4))°，仰角=$(round(data_ang_true[2,i,k];sigdigits=4))°");
        println("       估计：方位角=$(round(sort_ang_est[1,i,k];sigdigits=4))°，座角=$(round(sort_ang_est[2,i,k];sigdigits=4))°");
        println("       误差：方位角=$(round(abs(sort_ang_est[1,i,k]-data_ang_true[1,i,k]);sigdigits=4))°，座角=$(round(abs(sort_ang_est[2,i,k]-data_ang_true[2,i,k]);sigdigits=4))°");
    end
    println("最大误差：△(round(maximum(abs.（data_ang_true[:,:,k]。-sort_ang_est[:,:,k]));sigdigits=4))°")
    println("-------------------------------------------------------------")
end

tgtrng = maximum(abs_dist) # 目标范围
# SNR的计算，精度为3位
snr = round(radareqsnr(
	lambda, # 波长
	tgtrng, # 目标范围
	Pt, # 脉冲功率
	T, # 脉冲持续时间
	RCS = 0.01, # 目标的EPR（选择平均值）
	Gain = Grx, # 天线增益
	Loss = LFtx+LFrx,  # 路径中的损失
	Ts = TN # 噪音温度
)[1];sigdigits = 3)

println("当前信噪比：＜信噪比db")

y = deepcopy(data_y[:,:,end]) # 复制最后一个脉冲的信号数据
ang_trues_end =  data_ang_true[:,:,end]; # 复制目标的真实方向

noise_vec = collect(20:10:70) # 噪声功率的矢量相对于NP从20dB增加到70dB
snr_vec = snr .- noise_vec # SNR值的向量

println("SNR矢量：△(round.（snr_vec;sigdigits=3））dB")

# 为音乐不可比拟的频谱和角度估计分配内存
spec_music = zeros(length(el_scan),length(az_scan), length(noise_vec)) 
ang_music = zeros(2, num_tgts, length(noise_vec)) # 目标的真正方向

@inbounds for i in 1:length(noise_vec)
    println("SNR的估计计算△(round(snr_vec[i];sigdigits=3))dB")
    release!(music) # 将系统对象重置为初始状态
    noise_sig =  sqrt(NP*10^(noise_vec[i]/10))*randn(size(y)) # 噪声信号的新数据
    spec_music[:,:,i],ang_music[:,:,i] = music(y.+noise_sig); # 评估目标到达的范围和方向
end

fig_spec_music_vec = [] # 
@inbounds for i in 1:length(noise_vec)
    push!(fig_spec_music_vec,
        plot_spec_azel(az_scan,el_scan,spec_music[:,:,i];
            title = "音乐(SNR=$(round(snr_vec[i];sigdigits=3))db)"
        )
    )
end
plot(fig_spec_music_vec...)

sort_ang_music = similar(ang_music) # 排序的角度估计数据
[sort_ang_music[:,:,i] .= sort_ang_mse(ang_trues_end,ang_music[:,:,i]) for i in 1:length(snr_vec)]

# 建立一个表
calc_table_ang_analysis(ang_trues_end,sort_ang_music,snr_vec)

# 为MVDR分配不可比拟的频谱和角度估计的内存
spec_mvdr = zeros(length(el_scan),length(az_scan), length(noise_vec)) 
ang_mvdr = zeros(2, num_tgts, length(noise_vec)) # 目标的真正方向

@inbounds for i in 1:length(noise_vec)
    println("SNR的估计计算△(round(snr_vec[i];sigdigits=3))dB")
    release!(mvdr) # 将系统对象重置为初始状态
    noise_sig =  sqrt(NP*10^(noise_vec[i]/10))*randn(size(y)) # 噪声信号的新数据
    spec_mvdr[:,:,i],ang_mvdr[:,:,i] = mvdr(y.+noise_sig); # 评估目标到达的范围和方向
end

fig_spec_mvdr_vec = [] # 
@inbounds for i in 1:length(noise_vec)
    push!(fig_spec_mvdr_vec,
        plot_spec_azel(az_scan,el_scan,spec_mvdr[:,:,i];
            title = "MVDR(SNR=∞(round(snr_vec[i];sigdigits=3))db)"
        )
    )
end
plot(fig_spec_mvdr_vec...)

sort_ang_mvdr = similar(ang_mvdr) # 排序的角度估计数据
[sort_ang_mvdr[:,:,i] .= sort_ang_mse(ang_trues_end,ang_mvdr[:,:,i]) for i in 1:length(snr_vec)]

# 建立一个表
calc_table_ang_analysis(ang_trues_end,sort_ang_mvdr,snr_vec)

# 为波束扫描分配不可比拟的频谱和角度估计的内存
spec_beamscan = zeros(length(el_scan),length(az_scan), length(noise_vec)) 
ang_beamscan = zeros(2, num_tgts, length(noise_vec)) # 目标的真正方向

@inbounds for i in 1:length(noise_vec)
    println("SNR的估计计算△(round(snr_vec[i];sigdigits=3))dB")
    release!(beamscan) # 将系统对象重置为初始状态
    noise_sig =  sqrt(NP*10^(noise_vec[i]/10))*randn(size(y)) # 噪声信号的新数据
    spec_beamscan[:,:,i],ang_beamscan[:,:,i] = beamscan(y.+noise_sig); # 评估目标到达的范围和方向
end

fig_spec_beamscan_vec = [] # 
@inbounds for i in 1:length(noise_vec)
    push!(fig_spec_beamscan_vec,
        plot_spec_azel(az_scan,el_scan,spec_beamscan[:,:,i];
            title = "beamscan(SNR=∞(round(snr_vec[i];sigdigits=3))db)"
        )
    )
end
plot(fig_spec_beamscan_vec...)

sort_ang_beamscan = similar(ang_beamscan) # 排序的角度估计数据
[sort_ang_beamscan[:,:,i] .= sort_ang_mse(ang_trues_end,ang_beamscan[:,:,i]) for i in 1:length(snr_vec)]

# 建立一个表
calc_table_ang_analysis(ang_trues_end,sort_ang_beamscan,snr_vec)

plot_spec_2d = [];
for snr_i in eachindex(snr_vec)
    spec_beamscan_snr_0db_norm = maximum(spec_beamscan[:,:,snr_i];dims=1)[:]./maximum(spec_beamscan[:,:,snr_i])
    spec_mvdr_snr_0db_norm = maximum(spec_mvdr[:,:,snr_i];dims=1)[:]./maximum(spec_mvdr[:,:,snr_i])
    spec_music_snr_0db_norm = maximum(spec_music[:,:,snr_i];dims=1)[:]./maximum(spec_music[:,:,snr_i])

    fig = plot(az_scan,spec_beamscan_snr_0db_norm,lab="BeamScan",xlab="方位角，度",
        ylab="额定功率",lw=2,
        title="Snr=$(round(snr_vec[snr_i];sigdigits=3))dB的测向算法比较"
    )
    plot!(az_scan,spec_mvdr_snr_0db_norm,lab="MVDR",lw=2)
    plot!(az_scan,spec_music_snr_0db_norm,lab="MUSIC",lw=2)
    push!(plot_spec_2d,fig)
end
plot(plot_spec_2d...;layout = (length(plot_spec_2d),1),size = (600,300*length(plot_spec_2d)))