统计信号处理小作业——瑞利分布噪声中确定性直流信号的检测

背景简介

采用“M/N”检测，即$N$次独立检测，$M$次检测到目标即确定目标存在。对于虚警，只需要将假设修改为目标不存在，其他一致。

两假设如下：
$$\begin{array}{cc}{H}_0:& z(k)=n(k)\\ H_1:& z(k)=A+n(k)\end{array}k=0,1,\cdots ,N-1$$
其中，$n(k)$服从参量$\sigma$已知的瑞利分布，其PDF如下
$$f(z)=\frac{z}{{\sigma }^2}\cdot e^{-\frac{z^2}{2{\sigma }^2}}$$
$A$为确定性常数。

虚警、检测概率

综上，假设$A>0$、虚警概率为$P_{FA}$，可以写出单次检测的门限$\gamma$和检测概率$P_D$如下：
$$\begin{gathered} \gamma =\sqrt{-2{\sigma }^2\cdot \ln P_{FA}} \\ {P}_D=e^{-\frac{(A-\gamma {)}^2}{2{\sigma }^2}} \end{gathered}$$

对于$A<0$的情况，可以推到类似结论，只是在上述结论的基础上减一取负。

对于“1/N”检验，显然其总检测概率和虚警概率如下：
$$\begin{gathered} P_D^{\prime }=1-(1-P_D{)}^N \\ {P}_{FA}^{\prime }=1-(1-P_{FA}{)}^N \end{gathered}$$

MATLAB仿真

可以通过MATLAB仿真验证上述结论。设$\sigma =1$、$A=1$、$N=10$，假设单次虚警概率$P_{FA}=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