Radar Design Part I: From Power Budget Analysis to Dynamic Scenario Simulation

rng('default') % Set the random number generator for reproducible results
run("Part_I_EXPORT_Swerling1.m") % Script exported from Radar Designer

radarSensorClone = clone(radarSensorModel); % Save the original radar from the original app configuration
targetClone = target; % Save the fluctuating target for later analysis

run("Part_I_EXPORT_Swerling0.m")

run('Part_I_EXPORT_Swerling0_Noise.m')

run('Part_I_EXPORT_Swerling0_FalseAlarm.m')

rng('default') % Set the random number generator for reproducible results
target = targetClone; % Fluctuating target from the original export from the app
radarSensorModel = radarSensorClone; % Original radar configuration from the app
reset(radarSensorModel)
updateRate = .1; % (Hz)
scene = radarScenario('UpdateRate',updateRate);

% Place radar on a stationary platform
radarPos = [0,0,radar.Height]; % East-North-Up Cartesian Coordinates
radarOrientation = [0 0 0];  % (deg) Yaw, Pitch, Roll
radarPlat = platform(scene,'Position',radarPos,'Sensors',radarSensorModel,'Orientation',radarOrientation);

% Place a static target within the detectable range of the radar
targetPlat = platform(scene,'Position',191.6e3*[cosd(target.ElevationAngle) 0 sind(target.ElevationAngle)], ...
    'Signatures',{rcsSignature('Pattern',ones(2,2)*pow2db(target.RCS), ...
    'FluctuationModel',target.SwerlingCase)});

% Run Monte-Carlo Experiment to Test Detectability
T = 10000; % (s) Time of arrival of target at final position
n=1;
simDet = zeros(ceil(T*updateRate+1),1);
while n <= numel(simDet)
    dets = detect(scene);
    if ~isempty(dets)
        simDet(n) = 1;
    end
    n = n+1;
end
mean(simDet)

% Create Theater Plot and Plotters
[scene,tp] = helperConstructScene(radarSensorModel, radar, true);

rng('default') % Set the random number generator for reproducible results
reset(radarSensorModel) % Allows you to rerun the radar

detBuffer = {};
while advance(scene)
    [dets, config] = detect(scene);

    detBuffer = [detBuffer;dets]; %#ok<AGROW>

    % Is full scan complete?
    if config.IsScanDone
        break % yes
    end
end
statStaringDets = detBuffer; % Save results for final comparison
detPos = zeros(numel(detBuffer),3);
for m = 1:numel(detBuffer)
    detPos(m,:) = detBuffer{m}.Measurement(1:3).';
end
% Plot the detections on the theater plot
statDetPlotter = detectionPlotter(tp,'DisplayName','Statistical Detection',...
    'Marker','x','MarkerSize',20);
plotDetection(statDetPlotter,detPos)

rng('default') % Set the random number generator for reproducible results
release(radarSensorModel) % Allows us to modify the radar

rpm  = 6; % Revolutions per minute
scanrate = rpm*180/60;  % deg/s
updateRate = scanrate/radar.EffectiveAzRes; % Hz

radarSensorModel.ScanMode = "Electronic";
radarSensorModel.UpdateRate = updateRate;
radarSensorModel.FieldOfView = [radar.EffectiveAzRes,60]; % deg
radarSensorModel.ElectronicElevationLimits = [0 45];      % deg
radarSensorModel.ElectronicAzimuthLimits = [-90 90];      % deg

[scene,tp,FOVPlotter] = helperConstructScene(radarSensorModel, radar, true);

detBuffer = {};
% Run Simulation
while advance(scene)
    plotCoverage(FOVPlotter,coverageConfig(radarSensorModel))
    [dets, config] = detect(scene);

    detBuffer = [detBuffer;dets]; %#ok<AGROW>

    % Is full scan complete?
    if config.IsScanDone
        break % yes
    end
    pause(.001)
end
statScanningDets = detBuffer; % Save results for final comparison
% Plot the detections on the theater plot
detPos = zeros(numel(detBuffer),3);
for m = 1:numel(detBuffer)
    detPos(m,:) = detBuffer{m}.Measurement(1:3).';
end
statDetPlotter = detectionPlotter(tp,'DisplayName','Statistical Detection',...
    'Marker','x','MarkerSize',20);
plotDetection(statDetPlotter,detPos)

rgs = helperGetTargetRanges(scene); % (m) Ranges of targets

% Radar Designer App
radarSNRExpected = availableSNR + pow2db(radar.NumCoherentPulses) - radar.CustomLoss; % expected SNR (with processing gain) for non-fluctuating target
appSNR = interp1(target.ranges,radarSNRExpected,rgs);

% Staring Radar
statStaringSNRs = helperSortSNRs(statStaringDets);
statStaringSNRs = [statStaringSNRs; NaN];

% Scanning Radar
statScanningSNRs = helperSortSNRs(statScanningDets);
statScanningSNRs = [statScanningSNRs; NaN];

figure
bar(["Target 1","Target 2","Target 3","Target 4\newline(Undetected)"],[appSNR statStaringSNRs statScanningSNRs]')
yline(npwgnthresh(radar.Pfa),'LineWidth',2,'LineStyle','--')
ylabel('SNR (dB)')
title('Comparing Target SNR For Each Simulation Type')
grid
grid minor
legend('Predicted SNR', 'Simulated SNR (staring radar)','Simulated SNR (scanning radar)','SNR Detection Threshold')

% Save results for Part II results comparison 
data.results = [appSNR statStaringSNRs statScanningSNRs];
data.labels = ["Power-Level Prediction", "Measurement-Level Staring Radar","Measurement-Level Scanning Radar"];
save('AppAndStatResults.mat','-struct','data');

% Define targets in the scenario (same as in Part II of this example series)
function [tgt1, tgt2, tgt3, tgt4] = helperCreateTargets(scene)
    tgt1 = platform(scene, ...
        'ClassID',1, ... % set to 1 for targets
        'Position',[30e3 20e3 5e3], ...
        'Signatures',{rcsSignature(Pattern=5)}); % (dBsm)
    tgt2 = platform(scene, ...
        'ClassID',1, ...
        'Position',[95e3 -50e3 10e3], ...
        'Signatures',{rcsSignature(Pattern=5)}); % (dBsm)
    tgt3 = platform(scene, ...
        'ClassID',1, ...
        'Position',[125e3 30e3 12e3], ...
        'Signatures',{rcsSignature(Pattern=5)}); % (dBsm)
    tgt4 = platform(scene, ...
        'ClassID',1, ...
        'Position',[250e3 250e3 5e3], ...
        'Signatures',{rcsSignature(Pattern=5)}); % (dBsm)
end

% Construct radar scenario and dynamic plot
function [scene, tp, FOVPlotter] = helperConstructScene(radarSensorModel, radar, plotFlag)    
    
    scene = radarScenario('UpdateRate',0); % UpdateRate of 0 updates when new measurement is available
    
    % Define Targets
    [tgt1, tgt2, tgt3, tgt4] = helperCreateTargets(scene);
    
    % Define Radar Platform
    radarPos = [0,0,radar.Height]; % East-North-Up Cartesian Coordinates
    radarOrientation = [0 0 0];  % (deg) Yaw, Pitch, Roll
    % ClassID of the radar is 2 while the targets are set to 1
    radarPlat = platform(scene,'ClassID',2,'Position',radarPos,'Sensors',radarSensorModel,'Orientation',radarOrientation);
    
    tp = false;
    if plotFlag
        tp = theaterPlot('XLimits',[0 500e3],'YLimits',[-500e3 500e3],'ZLimits',[0,50e3],...
            'AxesUnits',["km","km","km"]);
        legend('Location','northeast','BackgroundAlpha',.6)
        title('Live Radar Simulation')
        view(157,41);grid on;
    
        % Plot Target Ground Truths
        gtplot = platformPlotter(tp,'DisplayName','Target Ground Truth',...
            'Marker','^','MarkerSize',8,'MarkerFaceColor','g');
        plotPlatform(gtplot,[tgt1.Position;tgt2.Position;tgt3.Position;tgt4.Position],{'Target1','Target2','Target3','Target4'});
    
        % Plot Radar
        radarplot = platformPlotter(tp,'DisplayName','Radar',...
            'Marker','s','MarkerSize',8,'MarkerFaceColor','b');
        plotPlatform(radarplot, radarPlat.Position,...
                radarPlat.Dimensions,quaternion(radarPlat.Orientation,'rotvecd'))
        
        % Plot Radar Field of View / Coverage
        FOVPlotter = coveragePlotter(tp,'DisplayName','Sensor Coverage');
        plotCoverage(FOVPlotter,coverageConfig(radarPlat.Sensors{1}))
    end
end

% Find the target platforms' ranges relative to radar
function rgs = helperGetTargetRanges(scene)
    plats = scene.Platforms;
    radarPlatIdx = cell2mat(cellfun(@(x) x.ClassID==2,plats,'UniformOutput',false));

    targetPosesVector = targetPoses(plats{radarPlatIdx}); % Relative to radar platform
    rgs = zeros(numel(plats)-1,1);

    for ii = 1:numel(targetPosesVector)
        rgs(ii) = norm(targetPosesVector(ii).Position);
    end
end

% Sort SNRs according to target indexes
function SNRVector = helperSortSNRs(dets)
    targetIdx = cellfun(@(x) x.ObjectAttributes{1,1}.TargetIndex, dets, 'UniformOutput', 1);
    [~,sortIdx] = sort(targetIdx,'ascend');
    SNRVector = cellfun(@(x) x.ObjectAttributes{1,1}.SNR, dets, 'UniformOutput', 1);
    SNRVector = SNRVector(sortIdx);
end