背景简介
采用“M/N”检测,即 N N N次独立检测, M M M次检测到目标即确定目标存在。对于虚警,只需要将假设修改为目标不存在,其他一致。
两假设如下:
H
0
:
z
(
k
)
=
n
(
k
)
H
1
:
z
(
k
)
=
A
+
n
(
k
)
k
=
0
,
1
,
⋯
,
N
−
1
\begin{matrix}H_0:&z(k)=n(k)\\H_1:&z(k)=A+n(k)\\ \end{matrix}\quad k=0,1,\cdots,N-1
H0:H1:z(k)=n(k)z(k)=A+n(k)k=0,1,⋯,N−1
其中,
n
(
k
)
n(k)
n(k)服从参量
σ
\sigma
σ已知的瑞利分布,其PDF如下
f
(
z
)
=
z
σ
2
⋅
e
−
z
2
2
σ
2
f(z)=\frac{z}{\sigma^2}\cdot e^{-\frac{z^2}{2\sigma^2}}
f(z)=σ2z⋅e−2σ2z2
A
A
A为确定性常数。
虚警、检测概率
综上,假设
A
>
0
A>0
A>0、虚警概率为
P
F
A
P_{FA}
PFA,可以写出单次检测的门限
γ
\gamma
γ和检测概率
P
D
P_D
PD如下:
γ
=
−
2
σ
2
⋅
ln
P
F
A
P
D
=
e
−
(
A
−
γ
)
2
2
σ
2
\gamma=\sqrt{-2\sigma^2\cdot\ln P_{FA}}\\ P_D=e^{-\frac{(A-\gamma)^2}{2\sigma^2}}
γ=−2σ2⋅lnPFAPD=e−2σ2(A−γ)2
对于 A < 0 A<0 A<0的情况,可以推到类似结论,只是在上述结论的基础上减一取负。
对于“1/N”检验,显然其总检测概率和虚警概率如下:
P
D
′
=
1
−
(
1
−
P
D
)
N
P
F
A
′
=
1
−
(
1
−
P
F
A
)
N
P_D^\prime=1-(1-P_D)^N\\P_{FA}^\prime=1-(1-P_{FA})^N
PD′=1−(1−PD)NPFA′=1−(1−PFA)N
MATLAB仿真
可以通过MATLAB仿真验证上述结论。设 σ = 1 \sigma=1 σ=1、 A = 1 A=1 A=1、 N = 10 N=10 N=10,假设单次虚警概率 P F A = 0.01 P_{FA}=0.01 PFA=0.01。直接给出MATLAB代码。
主程序:
%% clear
close all;
clearvars;
clc;
%% variables
pfa = 0.01; % 虚警
sigma = 1; % 标准差,是sigma不是sigma^2
A = 1; % A必须大于0
N = 10;
M = 1;
simulationTimes = 10000;
threshold = sqrt(-2 * sigma^2 * log(pfa));
pd = exp(-1 * (A - threshold)^2 / (2 * sigma^2));
%% draw histogram
noise = RayleighDistRandArr(sigma,[1,N * simulationTimes]);
% simulationTimes倍
signal = A + RayleighDistRandArr(sigma,[1,N * simulationTimes]);
f1 = figure(1);
subplot(2,1,1);
histogram(noise);
title("噪声的PDF");
xlim([0,10]);
subplot(2,1,2);
histogram(signal);
title("信号的PDF");
xlim([0,10]);
sgtitle("噪声和信号的概率分布统计直方图");
% close(f1);
%% monte-carlo simulation
noiseDetections = zeros([1,simulationTimes]);
for mcT = 1:simulationTimes
noise = RayleighDistRandArr(sigma,[1,N]);
detection = (noise >= threshold);
sumDetection = (sum(detection,'all') >= M); % 1/N检测的总检测
noiseDetections(1,mcT) = sumDetection;
end
pfa_simu = sum(noiseDetections,'all') / simulationTimes;
% 仿真的虚警率
disp(strcat("理论单次虚警率 ",num2str(pfa)," ,仿真的总虚警率 ",num2str(pfa_simu)));
signalDetections = zeros([1,simulationTimes]);
for mcT = 1:simulationTimes
signal = A + RayleighDistRandArr(sigma,[1,N]);
detection = (signal >= threshold);
sumDetection = (sum(detection,'all') >= M); % 1/N检测的总检测
signalDetections(1,mcT) = sumDetection;
end
pd_simu = sum(signalDetections,'all') / simulationTimes;
% 仿真的检测率
disp(strcat("理论单次检测率 ",num2str(pd)," ,仿真的总检测率 ",num2str(pd_simu)));
RayleighDistRandArr.m
function value = RayleighDistRandArr(sigma,arrayShape)
%RayleighDistRand 生成服从瑞利分布的随机数
% value = RayleighDistRand(sigma) 生成一个服从瑞利分布的随机数,
% 其中,瑞利分布内的参数为sigma,sigma必须大于0,且为实数。
%
% value = RayleighDistRand(sigma,arrayShape) 根据arrayShape
% 指定的大小生成服从瑞利分布的随机数矩阵。
%
% Copyright 小涛29 2022
if nargin < 1
error(message('必须输入sigma'));
else
if nargin >= 2
% 输入两个数,生成arrayShape形状的矩阵
normRandValue = rand(arrayShape);
else
% 输入一个数,内部生成一个1x1的normRandValue值
normRandValue = rand;
end
end
sigma = real(sigma);
if (any(normRandValue < 0,'all') && any(normRandValue > 1,'all'))
% 随机值不是[0,1]内的随机分布,报错
error(message('输入随机数必须是[0,1]内的均匀分布'));
end
if (sigma < 0)
% sigma必须大于0
error(message('sigma必须大于0'));
end
value = sqrt(-2 * log(normRandValue)) * sigma;
% 将均匀分布转换为瑞利分布
end