calccfarthreshold.m

利用matlab编写的横虚警检测程序

function T = calcCFARthreshold(nTest, nRef, thePFA)
% calculates classical radar CFAR threshold for 0-mean AWGN
%
% classical Neyman-Pearson detection threshold for radar detection
% under additive Gaussian white noise criterion and specifid false alarm
% probability.
%
% the threshold T is calculated such that under the noise-only condition
%
% prob(Ptest > T * Pref) <= pfa                        Eq(1)
%
% where Ptest is the sum of normed-squares of the nTest cells
%       Pref  is the sum of normed-squares of the nRef cells
%       pfa is the required maximum type-I (aka false alarm prob) error
%
%
% inputs: nTest  - the number of noise-only test cells
%         nRef   - the number of noise-only reference cells
%         thePfA - maximum probability such that 
%                  prob(Ptest/Pref > T) <= pfa         Eq(2)
%                  This is equivalent to Eq(1) 
%
% demo:  calcCFARthreshold;  % no inputs
%
% the classical CFAR processor compares the total power in nTest cells
% to the power in nRef cells. The threshold T is selected so that the 
% ratio of Eq(2) is satisfied. it is assumed that the both sets' cells
% contain only thermal (AWGN) noise. If this ratio exceeds T, it assumed a 
% radar target exists. The probability that only noise exceeds T is given 
% by the threshold T (usually on the order of 10^-4 to 10^-7). 
%
% michaelB brost. as usual, fully sharable under GPLv3
%

% demo
if(nargin == 0),
    nTest = 10;
    nRef  = 25;
    thePFA = 1e-4;
    T = calcCFARthreshold(nTest, nRef, thePFA);
    doTest(nTest, nRef, thePFA, T);
    clear('T');
    return;
end

% requires the statistics toolbox
T = finv(1 - thePFA, nTest, nRef) * nTest / nRef;

return

function doTest(nTest, nRef, thePfa, T)
% simple test of threshold calc

nPt   = min(1e5, ceil(50/thePfa));
nIter = 15;
nHit  = zeros(nIter, 1);

wHnd = waitbar(0, '');
for k1=1:nIter
    % test set noise power
    testPower = sum(randn(nPt, nTest).^2, 2);

    % reference set noise power
    refPower = sum(randn(nPt, nRef).^2, 2);
    
    % test
    index = find(testPower >= (T .* refPower));
    
    % count up contacts
    nHit(k1) = length(index);
    
    waitbar(k1/nIter, wHnd, sprintf('PFA simulation iteration %d of %d', k1, nIter));
   
end

close(wHnd);

pHit = nHit / nPt;

figure;

stem(pHit);
ylabel('P_F_A');
xlabel('runs');

pfaAvg = mean(pHit);
pfaStd = sqrt(var(pHit));
title(sprintf(...
    'calculated average pfa: %5.3e, 1-sigma: %5.3e (design: %5.2e)\n', ...
    pfaAvg, pfaStd, thePfa));

hold on;
lHnd = line([1, nIter], [thePfa, thePfa]);
set(lHnd, 'color', 'k');

lHnd = line([1, nIter], [pfaAvg, pfaAvg]);
set(lHnd, 'color', 'b');

lHnd = line([1, nIter], [thePfa, thePfa] - pfaStd);
set(lHnd, 'color', 'r');

lHnd = line([1, nIter], [thePfa, thePfa] + pfaStd);
set(lHnd, 'color', 'r');

legend('simulated PFA', 'specified PFA', 'mean PFA', '1 \sigma limits', ...
'location', 'south');

set(gca, 'xlim', [0, nIter+1]);