蒙特卡罗ROC模拟
本例展示了如何使用蒙特卡罗模拟生成雷达系统的接收机工作特征(ROC)曲线。接收机的工作特性决定了当目标不存在(假警报)时,系统在拒绝大的伪信号值的同时检测目标的能力。探测系统将通过将接收到的信号值与预设的阈值进行比较来宣布目标的存在或不存在。被发现的概率(Pd)是指当目标实际存在时,瞬时信号值大于阈值的概率。虚报的概率(Pfa)为目标不存在时信号值大于阈值的概率。在这种情况下,信号是由噪声引起的,其性质取决于噪声统计量。蒙特卡罗模拟产生大量的雷达回波,无论是否有目标存在。仿真计算Pd而且Pfa是通过计算每种情况下超过阈值的信号值的比例。
ROC曲线Pd作为函数Pfa.ROC曲线的形状取决于信号的接收信噪比。如果到达的信号信噪比是已知的,那么ROC曲线显示系统在方面的表现如何Pd而且Pfa.如果你指定Pd而且Pfa,然后可以确定需要多少功率才能达到这个要求。
你可以使用这个函数rocsnr
计算理论ROC曲线。本例显示了由蒙特卡罗模拟单天线雷达系统生成的ROC曲线,并将该曲线与理论曲线进行了比较。
指定雷达要求
设置期望的检测概率为0.9,虚警概率为 .设置雷达最大距离为4000米,距离分辨率为50米。设定实际目标距离为3000米。将目标雷达截面设置为1.5平方米,工作频率设置为10ghz。所有的计算都在基带中进行。
C = physconst(“光速”);Pd = 0.9;Pfa = 1e-6;Max_range = 4000;Target_range = 300;Range_res = 50;Tgt_rcs = 1.5;Fc = 10e9;Lambda = c/fc;
任何模拟计算Pfa而且pd需要处理许多信号。为了保持较低的内存需求,可以以脉冲块的形式处理信号。设置要处理的脉冲数为45000,设置每个块的大小为10000。
Npulse = 45000;Npulsebuffsize = 10000;
选择波形和信号参数
利用脉冲距离分辨率计算波形脉冲带宽。从最大范围计算脉冲重复频率。因为信号是基带,所以将采样频率设置为带宽的两倍。根据脉冲带宽计算脉冲持续时间。
Pulse_bw = c/(2*range_res);PRF = c/(2*max_range);Fs = 2*pulse_bw;脉冲时长= 10/pulse_bw;波形=相控。LinearFMWaveform (“脉冲宽度”pulse_duration,...“SampleRate”fs,“SweepBandwidth”,...pulse_bw,脉冲重复频率的脉冲重复频率);
实现一个特定的目标Pd而且Pfa要求在目标反射信号后,有足够的信号功率到达接收机。利用Albersheim方程计算实现指定虚警概率和检测概率所需的最小信噪比。
Snr_min = albersheim(pd,pfa);
为了达到这个信噪比,必须向目标传输足够的能量。利用雷达方程估计峰值发射功率,peak_power
,使目标在3000米范围内达到指定的信噪比(以分贝为单位)。接收到的信号也取决于目标雷达截面(RCS)。假设遵循一个非波动模型(转向0)。设置雷达具有相同的发射和接收增益20 dB。给出了雷达方程
Txrx_gain = 20;Peak_power = ((4*pi)^3*noisepow(1/pulse_duration)*target_range^4*...db2pow (snr_min)) / (db2pow (2 * txrx_gain) * tgt_rcs *λ^ 2)
Peak_power = 293.1827
设置发送系统对象
创建组成模拟传输部分的System对象:雷达平台、天线、发射机和散热器。
天线平台=相控。平台(...“InitialPosition”, (0;0;0),...“速度”, (0;0;0]);天线=相控。IsotropicAntennaElement (...“FrequencyRange”[5 e9 15 e9]);发射机=相控。发射机(...“获得”txrx_gain,...“PeakPower”peak_power,...“InUseOutputPort”,真正的);散热器=阶段性。散热器(...“传感器”、天线、...“OperatingFrequency”、fc);
设置目标系统对象
创建一个目标系统对象™,对应于具有非零目标横截面的实际反射目标。这个目标的反射将模拟实际的雷达回波。为了计算假警报,创建第二个雷达横截面为零的目标系统对象。这个目标的反射除噪声外为零。
目标{1}=阶段性。RadarTarget (...“MeanRCS”tgt_rcs,...“OperatingFrequency”、fc);目标平台{1}=阶段性。平台(...“InitialPosition”, (target_range;0;0]);目标{2}=阶段性。RadarTarget (...“MeanRCS”0,...“OperatingFrequency”、fc);目标平台{2}=阶段性。平台(...“InitialPosition”, (target_range;0;0]);
设置自由空间传播系统对象
建模从雷达到目标的传播环境。
信道{1}=相控。空闲空间(...“SampleRate”fs,...“TwoWayPropagation”,真的,...“OperatingFrequency”、fc);信道{2}=相控。空闲空间(...“SampleRate”fs,...“TwoWayPropagation”,真的,...“OperatingFrequency”、fc);
设置接收系统对象
属性指定噪声NoiseMethod
财产噪声温度的
和ReferenceTemperature
属性为290 K。
收集器=阶段性。收集器(...“传感器”、天线、...“OperatingFrequency”、fc);接收器=阶段性。ReceiverPreamp (...“获得”txrx_gain,...“NoiseMethod”,噪声温度的,...“ReferenceTemperature”, 290.0,...“NoiseFigure”0,...“SampleRate”fs,...“EnableInputPort”,真正的);接收器。SeedSource =“属性”;接收器。Seed = 2010;
指定快速网格
快速时间网格是一个脉冲重复时间间隔内的时间样本集。每个样本对应一个范围箱。
Fast_time_grid = unigrid(0,1/fs,1/prf,“()”);Rangebins = c*fast_time_grid/2;
从波形创建传输脉冲
创建要传输的波形。
Wavfrm =波形();
创建包含发射天线增益的发射信号。
[sigtrans,tx_status] =发射器(wavfrm);
从波形系统对象创建匹配的过滤器系数。然后创建匹配的筛选器System对象。
MFCoeff = getMatchedFilter(波形);matchingdelay = size(MFCoeff,1) - 1;过滤器=阶段性。MatchedFilter (...“系数”MFCoeff,...“GainOutputPort”、假);
计算目标范围箱
计算目标范围,然后计算范围bin数组的索引。因为目标和雷达是静止的,所以在整个模拟循环中使用相同的位置和速度值。您可以假设整个模拟的范围bin索引是恒定的。
ant_pos = antenaplatform . initialposition;ant_vel =天线平台。速度;tgt_pos =目标平台{1}.InitialPosition;tgt_vel =目标平台{1}.Velocity;[tgt_rng,tgt_ang] = rangeangle(tgt_pos,ant_pos);Rangeidx = val2ind(tgt_rng,rangebins(2)-rangebins(1),rangebins(1));
脉冲回路
创建一个信号处理循环。每个步骤都通过执行System对象来完成。循环对脉冲进行两次处理,一次用于目标存在的条件,一次用于目标不存在的条件。
利用向太空辐射信号
分阶段。散热器
.将信号传播到目标并返回到使用的天线
分阶段。空闲空间
.反映来自目标使用的信号
分阶段。目标
.接收反射信号在天线使用
分阶段。收集器
.将接收信号通过使用的接收放大器传递
分阶段。ReceiverPreamp
.这一步还将随机噪声添加到信号中。匹配滤波器放大信号使用
分阶段。MatchedFilter
.将匹配的过滤器输出存储在目标范围bin索引处,以供进一步分析。
rcv_pulse = 0(长度(sigtrans),Npulsebuffsize);h1 = 0 (Npulse,1);h0 = 0 (Npulse,1);Nbuff = floor(Npulse/Npulsebuffsize);Nrem = Npulse - Nbuff*Npulsebuffsize;为N = 1:2H1和H0假设Trsig =散热器(sigtrans,tgt_ang);Trsig =通道{n}(Trsig,...ant_pos tgt_pos,...ant_vel tgt_vel);Rcvsig =目标{n}(trsig);Rcvsig = collector(Rcvsig,tgt_ang);为k = 1:Nbuff为m = 1:Npulsebuffsize rcv_pulse (:,m) = receiver(rcvsig,~(tx_status>0));结束rcv_pulse = filter(rcv_pulse);rcv_脉冲=缓冲区(rcv_脉冲(匹配延迟+1:结束),大小(rcv_脉冲,1));如果n == 1 h1((1:Npulsebuffsize) + (k-1)*Npulsebuffsize) = rcv_pulse (rangeidx,:).';其他的h0((1:Npulsebuffsize) + (k-1)*Npulsebuffsize) = rcv_pulse (rangeidx,:).';结束结束如果(Nrem > 0)为m = 1:Nrem rcv_pulse (:,m) = receiver(rcvsig,~(tx_status>0));结束rcv_pulse = filter(rcv_pulse);rcv_脉冲=缓冲区(rcv_脉冲(匹配延迟+1:结束),大小(rcv_脉冲,1));如果n == 1 h1((1:Nrem) + Nbuff*Npulsebuffsize) = rcv_pulse (rangeidx,1:Nrem).';其他的h0((1:Nrem) + Nbuff*Npulsebuffsize) = rcv_pulse (rangeidx,1:Nrem).';结束结束结束
创建匹配过滤器输出的直方图
计算目标存在和目标不存在的回报的直方图。使用100个箱子来粗略估计信号值的分布。设置直方图值的范围从最小的信号到最大的信号。
H1a = abs(h1);H0a = abs(h0);Thresh_low = min([h1a;h0a]);Thresh_hi = max([h1a;h0a]);Nbins = 100;Binedges = linspace(thresh_low,thresh_hi,nbins);图直方图(h0a,binedges) hold住在直方图(h1a binedges)从标题(“目标-缺失Vs目标-存在直方图”)传说(没有目标的,“目标”)
模拟与理论比较Pd而且Pfa
来计算Pd而且Pfa,计算目标缺席回报和目标出现回报超过给定阈值的实例数。这组阈值的粒度比前面模拟中用于创建直方图的箱子更细。然后,通过脉冲数将这些计数归一化,得到概率的估计值。向量sim_pfa
是模拟虚警概率作为阈值的函数,打
.向量sim_pd
是模拟的检测概率,也是阈值的函数。接收端设置阈值,以确定目标是否存在。上面的直方图表明,最佳阈值在1.8左右。
Nbins = 1000;Thresh_steps = linspace(thresh_low,thresh_hi,nbins);Sim_pd = 0 (1,nbins);Sim_pfa = 0 (1,nbins);为K = 1:nbins thresh = thresh_steps(K);Sim_pd (k) = sum(h1a >= thresh);Sim_pfa (k) = sum(h0a >= thresh);结束sim_pd = sim_pd/Npulse;sim_pfa = sim_pfa/Npulse;
要绘制实验ROC曲线,必须反转Pfa曲线,以便绘制Pd反对Pfa.你可以求倒数Pfa只有当你能表达的时候才会曲线Pfa的严格单调递减函数打
.来表达Pfa这样,可以找到所有数组下标Pfa是常数除以相邻指标。对象中删除这些值Pd而且Pfa数组。
Pfa_diff = diff(sim_pfa);Idx = (pfa_diff = 0);Sim_pfa (idx) = [];Sim_pd (idx) = [];
将最小的Pfa限制为 .
Minpfa = 1e-6;N = sum(sim_pfa >= minpfa);sim_pfa = fliplr(sim_pfa(1:N)).';sim_pd = fliplr(sim_pd(1:N)).';
计算理论Pfa而且Pd最小的值Pfa为1。然后绘制理论图Pfa曲线。
[theor_pd,theor_pfa] = rocsnr(snr_min, snr_min)“SignalType”,...“NonfluctuatingNoncoherent”,...“MinPfa”minpfa,“NumPoints”N“NumPulses”1);semilogx (theor_pfa theor_pd)在semilogx (sim_pfa sim_pd,“r”。)标题(“模拟和理论ROC曲线”)包含(Pfa的) ylabel (“Pd”网格)在传奇(“理论”,“模拟”,“位置”,“本身”)
利用一百万次脉冲改进模拟
在前面的模拟中,Pd低值Pfa不要沿着平滑的曲线下落,甚至不要向下延伸到指定的操作状态。原因是在非常低的水平Pfa水平,极少,如果有的话,样本超过阈值。在低时生成曲线Pfa的倒数顺序,你必须使用大量的样本Pfa.这种类型的模拟需要很长时间。下面的曲线使用了一百万次脉冲而不是45,000次。要运行此模拟,请重复示例,但要设置Npulse
到1000000。