calcmofncfarthreshold.m

Compute Classical CFAR binary detection threshold for radar detection under additive Gaussian white

function [T,M,P1] = calcMofNCFARthreshold(nTest, nRef, N, thePFA)
% calculates classical 2-threshold radar detector CFAR thresholds 
% for 0-mean unknown, but constant variance AWGN for N sequential tests
% that is, it is a 2 threshold CFAR detector.
%
% classical Neyman-Pearson detection threshold for radar detection
% under additive Gaussian white noise criterion and specifid false alarm
% probability under the binary detector architecture
%
% the primary thresholds T are calculated such that under the
% noise-only condition for N trials:
%              z = (Ptest > T * Pref)                  Eq(1)
%              prob(sum(z) >= M | N trials) <= pfa     Eq(2)
%
% where Ptest is the sum of normed-squares of the nTest cells (total power)
% Pref is the sum of normed-squares of the nRef cells (total power)
% pfa is the required maximum type-I (aka false alarm prob) error
% M is the secondary threshold
% N is the number of trials
%
%   Independent identically distributed random variables IID-RV.
%   Each cell in nTest and nRef sets is the square of a Gaussian RV of 
%   0-mean and unknown but constant variance. This detector is variance 
%   independent and thus said to be CFAR. 
%   The test scenario looks like this:
%
%    trial 1             trial 2        ....    trial N
%  nTest    nRef       nTest    nRef         nTest    nRef      N trials
% ------------------------------------------------------------------------
% [-----] [-------]   [-----] [-------] ... [-----] [-------]   squared    
% [-----] [-------]   [-----] [-------] ... [-----] [-------] Gaussian RVs
%   ...      ...        ...      ...        ...      ...       
% [-----] [-------]   [-----] [-------] ... [-----] [-------] summed 
%   sum      sum        sum      sum          sum      sum
%    |        |          |        |            |        |
%    V        V          V        V            V        V
%      test 1              test 2                test N        | detector
% Ptest > Pref*T     Ptest > Pref*T   ...   Ptest > Pref*T     |   #1
%     0 OR 1      +       0 OR 1    + ... +       0 OR 1       |
%
%                    note that min sum = 0, max sum = N
%
%     < ----- detector # 2--------------------->
%     final sum >= M ??? [Yes] -> detection.    | the probability noise
%                        [No ] -> no detection  | alone triggers detector
%                                               | #2 is thePFA.
%
% ------------------------ algorithm parameters ----------------------
% inputs: nTest  - the number of noise-only test cells
%         nRef   - the number of noise-only reference cells
%         thePfA - probability noise only exceeds primary threshold at
%                  at least M times out of N trials and M=1:N
%         N      - number of trials
%
% outputs: M - the secondary thresholds = [1:N];
%          T - the primary thresholds such that the ratio of test cell's
%              power to reference power exceeds T at least M times out of N
%         P1 - probability of a single event (prob of threshold 1 crossing)
%
% see also calcCFARthreshold
% http://www.mathworks.com/matlabcentral/fx_files/22773
%
% demo: calcMofNCFARthreshold; % no inputs
%
% ***** this function requires the matlab statistics toolbox for both 
%       threshold calculations and the random noise in the monte carlo test
%
% michaelB brost. as usual, fully sharable under GPLv3
%

% demo similar to demo of calcCFARthreshold
if(nargin == 0),
nTest  = 100;
nRef   = 250;
thePFA = 1e-3;
N = 15; 
[T, M, P1] = calcMofNCFARthreshold(nTest, nRef, N, thePFA);

% simulate at most 4 T-values (no need to test everything) 
doTest(nTest, nRef, thePFA, T);
clear('T');
return;
end

% secondary thresholds are calculated here
M = [1:N]';

% primary thresholds are preallocated
T = zeros(size(M));

% prob of single event thresholds are preallocated
P1 = zeros(size(M));

% primary thresholds (T) are calculated. secondary thresholds (M) are free
for k1=1:length(M)
% single event probability for a specified second threshold
P1(k1) = betainv(thePFA, M(k1), N+1-M(k1));

% primary threshold depends upon secondary threshold and number of ref and test cells
T(k1) = finv(1 - P1(k1), nTest, nRef) * nTest / nRef;
end

return

% -------------------------------------------------------------------------------- %
% testing - brute force monte-carlo

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

% for brute force MC, take at most 500,000 samples. to save memory, coerce to
% single precision
nPt = min(1e5, ceil(50/thePfa));
nIter = 15;
nHit = zeros(nIter, 1);

wHnd = waitbar(0, '');


% only test at most 4 thresholds
if(length(T) > 4),
thresholdIndex = round([1, length(T)/3, 2*length(T)/3, length(T)]);
else,
thresholdIndex = 1:length(T);
end

N = length(T);

% overwrite M and T for testing
M = thresholdIndex;
T = T(thresholdIndex);

for t=1:length(T),

nHit = zeros(nIter, 1);
for k1=1:nIter,

mnCnt = zeros(nPt, 1);

for k2=1:N

% use fast f-distributed random number generator from statistics toolbox
% this is equivalent to the ratio of nTest summed squared normals to
% nRef summed squared normals scaled by the inverse of (nRef/nTest)
noiseRatio = single(frnd(nTest, nRef, [nPt, 1])) .* (nTest/nRef);

% primary threshold
thisT = single(T(t));

% find primary threshold crossings
index = find(noiseRatio >= thisT);

% count the number of primary crossings
mnCnt(index) = mnCnt(index) + 1;

end

waitbar(k1/nIter, wHnd, ...
    sprintf('PFA simulation iteration %d of %d for threshold %d', k1, nIter, t));

% secondary threshold applied to primary threshold crossing counts
index = find(mnCnt >= M(t));

hitCnt(k1) = length(index);
end

% 2-threshold false alarm probability 
pHit = hitCnt ./ nPt;

% make figure 50 pct wider than usual to accomodate the title
figure;
figSiz = get(gcf, 'position');
set(gcf, 'position', [figSiz(1)-figSiz(3)/4, figSiz(2), figSiz(3)*1.5, figSiz(4)]);

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

% calc mean and variance of our pfa vector
pfaAvg = mean(pHit);
pfaStd = sqrt(var(pHit));

% plot stuff.......
title(sprintf(...
'[T,M] threshold pair(%d of %d) pfa: %5.3e, 1-sigma: %5.3e (design: %5.2e)\n', ...
thresholdIndex(t), 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');

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

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


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


end

close(wHnd);