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

小涛29

已于 2022-05-30 16:11:47 修改

阅读量350

点赞数 2

文章标签：算法概率论 matlab

于 2022-05-25 11:02:19 首次发布

本文链接：https://blog.csdn.net/lovelyed/article/details/124961725

版权

背景简介

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

两假设如下：
$\begin{matrix}H_0:&z(k)=n(k)\\H_1:&z(k)=A+n(k)\\ \end{matrix}\quad 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$ 如下：
$\gamma=\sqrt{-2\sigma^2\cdot\ln P_{FA}}\\ P_D=e^{-\frac{(A-\gamma)^2}{2\sigma^2}}$

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

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

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