www.gusucode.com > phased 案例源码 matlab代码程序 > phased/AnalyzeNoiseOnly2DCFARExample.m

    %% Set CFAR 2-D Threshold for Noise-Only Data
% This example shows how to set a CFAR threshold based upon a required
% probability of false alarm. The e
%%
% *Note:* You can replace each call to the function with the equivalent |step|
% syntax. For example, replace |myObject(x)| with |step(myObject,x)|.
%%
% Perform cell-averaging CFAR detection on a 41-by-41 matrix of cells
% containing Gaussian noise. Estimate the empirical pfa and compare it to
% the required pfa. To get a good estimate, perform this simulation on 1000
% similar matrices. First, set a threshold using the required pfa. In this
% case, there are no targets and the pfa can be estimated from the number
% of cells that exceed the threshold. Assume that the data is processed
% throught a square-law detector and that no pulse integration is
% performed. Use a training-cell band of 3 cells in width and 4 cells in
% height. Use a guard band of 3 cells in width and 2 cells in height to
% separate the cells under test (CUT) from the training cells. Specify a
% required pfa of 5.0e-4.
p = 5e-4;
rs = RandStream.create('mt19937ar','Seed',5);
N = 41;
ntrials = 1000;
detector = phased.CFARDetector2D('TrainingBandSize',[4,3], ...
    'ThresholdFactor','Auto','GuardBandSize',[2,3], ...
    'ProbabilityFalseAlarm',p,'Method','SOCA','ThresholdOutputPort',true);
%%
% Create a 41-by-41 image containing random complex data. Then, square the
% data to simulate a square-law detector.
x = 2/sqrt(2)*(randn(rs,N,N,ntrials) + 1i*randn(rs,N,N,ntrials));
x2 = abs(x).^2;
%%
% Process all the cells in each image. To do this, find the row and column
% of each CUT cell whose training region falls entirely within each image.
Ngc = detector.GuardBandSize(2);
Ngr = detector.GuardBandSize(1);
Ntc = detector.TrainingBandSize(2);
Ntr = detector.TrainingBandSize(1);
cutidx = [];
colstart = Ntc + Ngc + 1;
colend = N - ( Ntr + Ngr);
rowstart = Ntc + Ngc + 1;
rowend = N - ( Ntr + Ngr);
for m = colstart:colend
    for n = rowstart:rowend
        cutidx = [cutidx,[n;m]];
    end
end
ncutcells = size(cutidx,2);
%%
% Display the CUT cells.
cutimage = zeros(N,N);
for k = 1:ncutcells
    cutimage(cutidx(1,k),cutidx(2,k)) = 1;
end
imagesc(cutimage)
axis equal

%%
% Perform the detection on all CUT cells. Return the detection classification and the
% threshold used to classify the cell.
[dets,th] = detector(x2,cutidx);
%%
% Find and display an image with a false alarm for illustration.
di = [];
for k = 1:ntrials
    d = dets(:,k);
    if (any(d) > 0)
        di = [di,k];
    end
end
idx = di(1);
detimg = zeros(N,N);
for k = 1:ncutcells
    detimg(cutidx(1,k),cutidx(2,k)) = dets(k,idx);
end
imagesc(detimg)
axis equal
%%
% Compute the empirical pfa.
pfa = sum(dets(:))/ntrials/ncutcells
%%
% The empirical and specified pfa agree.

%%
% Display the average empirical threshold value over all images.
mean(th(:))

%%
% Compute the theoretical threshold factor for the required pfa.
threshfactor = npwgnthresh(p,1,'noncoherent');
threshfactor = 10^(threshfactor/10);
disp(threshfactor)
%%
% The theoretical threshold factor multiplied by the noise variance should
% agree with the measured threshold.
noisevar = mean(x2(:));
disp(threshfactor*noisevar);
%%
% The theoretical threshold and empirical threshold agree to within
% an acceptable difference.