Вычисление вероятности одноцелевого трека при отсутствии ложных предупреждений

Бернулли-Сум

Сумма Бернулли позволяет быстро и легко анализировать эффективность в случае одной цели при отсутствии ложных предупреждений. Вероятность обнаружения дорожки $$P_{dt}$$ может быть определено в терминах вероятности обнаружения приемника ${\mathit{P}}_{\mathrm{d}}$ для периода окна, заданного как

${\mathit{T}}_{\mathit{w}}=\mathit{NT}$,

где T - основной период дискретизации, и N представляет количество возможностей для обнаружения.

Для логики порога подтверждения M-of-N, вероятность целевой дорожки ${\mathit{P}}_{\mathrm{dt}}$ определяется как

${\mathit{P}}_{\mathrm{dt}}={\sum }_{\mathit{i}=\mathit{M}}^{\mathit{N}}\mathit{C}\left(\mathit{N},\mathit{i}\right){\mathit{P}}_{\mathrm{d}}^{\mathit{i}}{\left(1-{\mathit{P}}_{\mathrm{d}}\right)}^{\left(\mathit{N}-\mathit{i}\right)}$,

где

$\mathit{C}\left(\mathit{N},\mathit{i}\right)=\frac{\mathit{N}!}{\left(\mathit{N}-\mathit{i}\right)!\mathit{i}!}$.

Логика порога подтверждения, обозначенная как M-of-N или M/N, является одноэтапной логикой, где дорожка должна ассоциироваться с обнаружением, также известным как попадание, по меньшей мере, M раз из N последовательных образов. Для примера рассмотрим логику 2 из 3. На рисунке ниже твердый желтый цвет представляет удар, который может быть либо от цели, либо от ложного предупреждения. Узорчатые синие блоки представляют промахи. Представленные ниже случаи не являются исчерпывающими, но, как можно видеть, случаи 2 и 3 удовлетворяют порогу, а случай 1 - нет.

Исследуйте вероятность целевого трека от вероятности обнаружения для порога подтверждения 2/3 с помощью метода суммы Бернулли. Обратите внимание, как быстро выполняется расчет. Выполните вычисление суммы Бернулли, принимая, что для инициализации дорожки требуется одно попадание.

% Define probabilities for analysis 
Pd = linspace(0,1,100).';

% Define confirmation threshold M/N
M = 2; % Number of hits
N = 3; % Number of observations or opportunities

% Perform Bernoulli sum assuming 1 hit is required to initialize a track
tic 
PdtBernoulli = helperBernoulliSum(Pd,M,N);
elapsedTime = toc;
helperUpdate('Bernoulli',elapsedTime);

Bernoulli calculation completed. Total computation time is 0.0185 seconds.

% Plot the probability of detection versus the probability of target track 
hAxes = helperPlot(Pd,PdtBernoulli,'M/N = 2/3','P_D','P_{DT}', ...
    sprintf('Bernoulli Sum\nProbability of Target Track in the Absence of False Alarms'));

% Set desired probability of target track
yline(hAxes,0.9,'--','LineWidth',2,'DisplayName','Desired P_{DT}'); 

Принимая, что была желательна вероятность целевого трека 0,9, вышеописанный график информирует аналитика о том, что требуется вероятность обнаружения около 0,7.

Вычисление вероятностей для цели в загромождении

Значения, полученные в результате вычисления суммы Бернулли, полезны в быстром анализе, но в целом не являются показательными для реальных окружений слежения, где на целевые дорожки влияет наличие ложных предупреждений. Рассмотрим сценарий, где цель работает в присутствии загромождения.

Принимая, что ложные аварийные сигналы происходят на основе per-look per-cell, вероятность ложных аварийных сигналов в шлюзе слежения зависит от количества камер в шлюзе. Предположим, что существует три типа событий,

Количество камер в воротах зависит от истории событий и порядка, в котором происходят события. Эти факторы диктуют последовательность роста ворот трекера.

Метод суммы Бернулли принимает, что нет ложных предупреждений и что порядок обнаружений не имеет значения. Таким образом, при использовании суммы Бернулли в цели в загроможденном сценарии, это приводит к чрезмерно оптимистическим результатам.

Одним из подходов, который используется для анализа таких сценариев, является оценка каждой возможной последовательности дорожек и определение, какие последовательности удовлетворяют логике порога подтверждения. Этот подход грубой силы создания марковской цепи, как правило, слишком интенсивен в вычислительном отношении.

Другой подход заключается в использовании анализа типа Монте-Карло. Вместо того, чтобы генерировать полную цепь Маркова, симуляция Монте-Карло вручную генерирует случайные последовательности N событий. Порог подтверждения применяется к каждой последовательности, и статистика агрегируется. Метод Монте-Карло основан на законе больших чисел, поэтому эффективность улучшается с увеличениями количества итераций. Анализ Монте-Карло хорошо поддается параллелизации, но в случае небольших вероятностей ложного предупреждения количество итераций может стать быстро несостоятельным. Таким образом, мы должны рассмотреть альтернативные методы, чтобы быстро вычислить трековые меры вероятности.

Алгоритм истории общих ворот (CGH)

Алгоритм Common Gate History значительно сокращает время расчетов и требования к памяти. Алгоритм предотвращает необходимость ручной генерации последовательностей, как в случае анализа Монте-Карло, который может быть дорогим для событий с низкой вероятностью.

Алгоритм начинается с предположения, что существует три типа треков, которые могут содержать:

Целевая дорожка определяется как любая дорожка, которая содержит, по меньшей мере, одно целевое обнаружение и удовлетворяет порогу подтверждения M/N. Таким образом, типы дорожек 1 и 2 рассматриваются как целевые дорожки, в то время как 3 рассматривается как ложная дорожка.

Учитывая ранее определенные типы дорожек и с логикой M/N порога подтверждения, предлагается статистическая величина для идентификации состояния дорожки. Состояние дорожки определяется как

где ${\omega }_{\mathit{l}}$ количество временных шагов с момента последнего обнаружения (целевое или ложное предупреждение), ${\omega }_{\mathit{lt}}$количество временных шагов с момента последнего обнаружения цели, и $\lambda$ - общее количество обнаружений (целевые объекты или ложные аварийные сигналы). Когда алгоритм протекает, вектор состояния дорожки развивается согласно марковской цепи.

Алгоритм принимает, что дорожка может быть запущена, задав только два типа событий:

После инициирования дорожки следующие четыре типа событий продолжают дорожку:

Вероятность дорожки в взгляде m умножается на вероятность события, которое продолжает дорожку в взгляде $\mathit{m}+1$. Затем дорожки суммируются путем добавления вероятностей дорожки файлов дорожки с общим вектором истории ворот. Это держит количество трековых состояний в марковской цепи в разумных пределах.

Предположения алгоритма Common Gate History следующие:

Алгоритм Common Gate History может использоваться, чтобы вычислить все четыре типа вероятностей:

Вероятность целевой дорожки и вероятность ложного пути формируют базис рабочей характеристики трекера (TOC). (TOC) дополняет рабочую характеристику приемника (ROC), которая обычно используется в анализе и предсказании эффективности приемников. Объединение ROC и TOC предоставляет сквозной инструмент системного анализа.

Вычислите и постройте график ТОС с помощью кривых ROC из rocsnr в качестве входов. Примите отношение сигнал/шум (ОСШ) 8 дБ. Продолжите использовать логику порога подтверждения 2/3, которая использовалась в примере суммы Бернулли. Используйте toccgh встроенный трекер. Обратите внимание, как быстро расчет, несмотря на низкие вероятности ложных предупреждений.

% Receiver operating characteristic (ROC)
snrdB = 8; % SNR (dB) 
[Pd,Pfa] = rocsnr(snrdB,'MaxPfa',1e-3,'MinPfa',1e-12,'NumPoints',20);

% Plot ROC
helperPlotLog(Pfa,Pd,snrdB, ...
    'Probability of False Alarm (P_{FA})', ...
    'Probability of Detection (P_D)', ...
    'Receiver Operating Characteristic');

% Common gate history algorithm
tic
[PdtCGH,PftCGH] = toccgh(Pd,Pfa,'ConfirmationThreshold',[M N]);
elapsedTime = toc;
helperUpdate('Common Gate History',elapsedTime);

Common Gate History calculation completed. Total computation time is 0.5575 seconds.

% Plot Common Gate History results 
hAxes = helperPlotLog(PftCGH,PdtCGH,'CGH','P_{FT}','P_{DT}', ...
    'Tracker Operating Characteristic (TOC) Curve');

Алгоритм Common Gate History позволяет оценивать эффективность трекера, подобную анализу Монте-Карло, но с приемлемым временем расчета, несмотря на низковероятностные события. Алгоритм Common Gate History, таким образом, позволяет проводить высокоуровневое исследование и выбирать опции перед более интенсивными, подробными симуляциями.

Использование CGH с пользовательскими трекерами

Рассмотрите разработку трекера для автомобильного приложения. Задайте пользовательский одномерный трекер почти постоянной скорости (NCV) с помощью trackingKF. Примите, что частота обновления $$\Delta t$$ составляет 1 секунду. Примите, что матрица переходов имеет вид

и технологический шум имеет вид

$\mathit{Q}=\left[\begin{array}{cc}\frac{\Delta {\mathit{t}}^4}{4}& \frac{\Delta {\mathit{t}}^3}{2}\\ \frac{\Delta {\mathit{t}}^3}{2}& \Delta {\mathit{t}}^2\end{array}\right]\mathit{q}$,

где $\mathit{q}$ - коэффициент настройки, заданный как

$\mathit{q}={\mathit{a}}_{\max }^2\Delta \mathit{t}$.

Вход ${\mathit{a}}_{\max }^2$- это максимальное ожидаемое целевое ускорение. Предположим, что максимальное ускорение составляет 4 м/с$${}^2$$ ожидается для транспортных средств.

% Define the state transition matrix
dt = 1; % Update rate (sec)
A = [1 dt; 0 1];

% Define the process noise 
Q = [dt^4/4 dt^3/2; dt^3/2 dt^2];

% Tune the process noise
amax = 4; % Maximum target acceleration (m/s^2)
q = amax^2*dt; 

% Update the process noise 
Q = Q.*q; 

% Initialize the Kalman filter 
trkfilt = trackingKF('MotionModel','Custom', ...
    'StateCovariance', [0 0; 0 0], ...
    'StateTransitionModel',A, ...
    'ProcessNoise',Q, ...
    'MeasurementNoise',0, ...
    'MeasurementModel',[1 0]);

Эллипс ошибки используется для моделирования неопределенности отслеживания. Из этого эллипса неопределенности может быть вычислена последовательность роста затвора.

1-$\sigma$ значения эллипса ошибки вычисляются как квадратный корень собственных значений $\lambda$ предсказанного состояния ковариации ${\mathit{P}}_{\mathit{k}+1|\mathit{k}}$:

$$[\sqrt{{\lambda }_1},\sqrt{{\lambda }_2}]=\sqrt{eig(P_{k+1|k})}$$.

Площадь эллипса ошибки затем вычисляется как

$\mathrm{Error}\mathrm{Ellipse}\mathrm{Area}=\pi \sqrt{\lambda {}_1}\sqrt{\lambda {}_2}$.

Площадь интервалов вычисляется как

$\mathrm{Bin}\mathrm{Area}=\left(\Delta \mathrm{Range}\right)\left(\Delta \mathrm{Range}\mathrm{Rate}\right)$.

Наконец, размер ворот в интервалах

$\mathrm{Gate}\mathrm{Size}=\frac{\mathrm{Error}\mathrm{Ellipse}\mathrm{Area}}{\mathrm{Bin}\mathrm{Area}}$.

Таким образом, размер затвора зависит от трекера, последовательности событий и разрешения интервалов. Вычислите последовательность роста затвора, принимая порог N подтверждения, равный 3. Предположим, что разрешения области значений и области значений для автомобильного радара равны 1 м и 1 м/с, соответственно.

% Calculate gate growth sequence
res = [1, 1]; % Bin resolutions [range (m), range-rate (m/s)]
gs = helperCalculateGateSize(N,trkfilt,res) 

gs = 1×3

     1    51   124

% Common gate history algorithm
tic
[PdtCGHcustom,PftCGHcustom] = toccgh(Pd,Pfa,'ConfirmationThreshold',[M N],'GateGrowthSequence',gs);
elapsedTime = toc;
helperUpdate('Common Gate History',elapsedTime);

Common Gate History calculation completed. Total computation time is 0.1840 seconds.

% Add plot to previous plot
helperAddPlotLog(hAxes,PftCGHcustom,PdtCGHcustom,'CGH with Custom Gate Growth Sequence');

Оценка эффективности трекера для автомобильной радиолокационной системы

Вероятность ложного предупреждения и вероятность требований целевого трека

При помощи ROC и TOC в сочетании системный аналитик может выбрать рабочую точку детектора, которая удовлетворяет общим системным требованиям. Рассмотрим корпус автомобильного радара. Из-за особенностей приложения желательно, чтобы ложные предупреждения оставались очень низкими событиями вероятности. Кроме того, вероятность целевого трека должна быть высокой в целях безопасности. Примите во внимание следующие два требования.

Вычислите кривые ROC для ОСШ, равные 6, 8, 10 и 12 дБ.

% Calculate ROC curves
snrdB = [6 8 10 12]; % SNR (dB)
[Pd,Pfa] = rocsnr(snrdB,'MaxPfa',1e-3,'MinPfa',1e-10,'NumPoints',10);

Используйте вычисленные кривые ROC в качестве входов для toccgh функция для генерации связанных ТОС-кривых. Используйте тот же порог подтверждения и ранее сгенерированную последовательность роста ворот.

% Generate TOC curves
tic
toccgh(Pd,Pfa,'ConfirmationThreshold',[M N],'GateGrowthSequence',gs);
elapsedTime = toc;
helperUpdate('Common Gate History',elapsedTime);

Common Gate History calculation completed. Total computation time is 0.6742 seconds.

% Requirement 1: Probability of false alarms must be less than 1e-6
hAxesROC = subplot(2,1,1);
xlim([1e-10 1e-2])
ylim([0 1])
reqPfa = 1e-6; 
helperColorZonesReqPfa(hAxesROC,reqPfa)
legend(hAxesROC,'Location','eastoutside')

% Requirement 2: Probability of target track must be equal to 0.9 or above
hAxesTOC = subplot(2,1,2);
xlim([1e-14 1e-4])
ylim([0 1])
reqPdt = 0.9; 
helperColorZonesReqPdt(hAxesTOC,reqPdt)

Требование 1 диктует, что только точки с 1 по 6 на кривых ROC могут быть включены для последующего анализа, поскольку они удовлетворяют вероятности ниже, чем ${10}^{-6}$. Точки с 7 по 10 не отвечают этому требованию.

Учитывая требование 2, можно видеть, что кривая ROC, соответствующая ОСШ на 6 дБ, не удовлетворяет второму требованию ни в одной точке. Единственными кривыми, которые нужно продолжить, являются кривые 8, 10 и 12 дБ. Требование 2 удовлетворяется только для точки 10 на кривой 8 дБ, точек 9 и 10 на кривой 10 дБ и точек 5-10 на кривой 12 дБ.

Объединяя требования 1 и 2, существует только две точки анализа, которые удовлетворяют обоим требованиям: Точки 5 и 6 на кривой 12 дБ. Точка 5 соответствует вероятности целевой дорожки 0,90 и вероятности ложной дорожки $1.27\times {10}^{-13}$. Перемещая их назад на кривую ROC, это соответствует вероятности обнаружения 0,68 и вероятности ложного предупреждения $1.29\times {10}^{-7}$. Точно так же точка 6 соответствует вероятности целевой дорожки 0,96, вероятность ложной дорожки $1.64\times {10}^{-12}$, вероятность обнаружения 0,80 и вероятность ложного предупреждения $7.74\times {10}^{-7}$. В этом случае выбирается рабочая точка 6. Это компромисс улучшенной вероятности целевого трека за счет несколько более высокой, но разумной вероятности ложного трека.

Алгоритм Common Gate History позволяет оценить ожидаемое количество ложных треков на основе количества целей, ожидаемых в окружении и количества камер в радиолокационных данных. Ожидаемое количество ложных треков ${\mathit{E}}_{\mathrm{ft}}$ вычисляется как

${\mathit{E}}_{\mathrm{ft}}={\mathit{P}}_{\mathrm{ft},\mathrm{nt}}{\mathit{N}}_{\mathrm{c}}+{\mathit{P}}_{\mathrm{ft}}{\mathit{N}}_{\mathrm{t}}$,

где ${\mathit{P}}_{\mathrm{ft},\mathrm{nt}}$ - вероятность ложного следа при отсутствии целей, ${\mathit{N}}_{\mathrm{c}}$ количество камер, ${\mathit{P}}_{\mathrm{ft}}$ - вероятность ложного следа в присутствии целей, и ${\mathit{N}}_{\mathrm{t}}$ количество целевых объектов.

Рассмотрим окружение, в которой количество целей, как ожидается, будет равным 10, а количество камер -

% Calculate expected number of false tracks using toccgh
numCells = 1e5;         % Number of cells in radar data
numTargets = 10;        % Number of targets in scenario
selectedPd = Pd(6,4) ;  % Selected probability of detection
selectedPfa = Pfa(6);   % Selected probability of false alarm
[Pdt,Pft,Eft] = toccgh(selectedPd,selectedPfa, ...
    'ConfirmationThreshold',[M N],'GateGrowthSequence',gs, ...
    'NumCells',numCells,'NumTargets',10);

% Output results
helperPrintTrackProbabilities(Pdt,Pft,Eft);

Probability of Target Track in Presence of False Alarms = 0.9581
Probability of False Track in the Presence of Targets = 1.6410e-12
Expected Number of False Tracks = 5

Таким образом, исходя из системных параметров, аналитик мог ожидать около пяти ложных треков.

Анализ порогов подтверждения

Рассмотрим тот же случай проекта автомобильного радара, но исследуйте эффект порогов подтверждения 2/4, 3/4 и 4/4. Предположим следующее:

Сначала вычислите кривую ROC с помощью rocpfa функция.

% Calculate ROC curve assuming a probability of false alarm of 1e-6
Pfa = 1e-6;
numPoints = 20;
[Pd,snrdB] = rocpfa(Pfa,'NumPoints',numPoints,'MaxSNR',15);

Обновите последовательность роста ворот из-за большего количества наблюдений.

% Update the gate growth sequence
N = 4; 
trkfilt = trackingKF('MotionModel','Custom', ...
    'StateCovariance', [0 0; 0 0], ...
    'StateTransitionModel',A, ...
    'ProcessNoise',Q, ...
    'MeasurementNoise',0, ...
    'MeasurementModel',[1 0]);
gs = helperCalculateGateSize(N,trkfilt,res)

gs = 1×4

     1    51   124   225

Вычислите TOC, учитывая кривые ROC в качестве входов. Анализируйте пороги подтверждения, равные 2/4, 3/4 и 4/4.

% Calculate TOC
cp = [2 4; 3 4; 4 4];
numCp = size(cp,1);
PdtMat = zeros(numPoints,numCp);
PftMat = zeros(numPoints,numCp);
EftMat = zeros(numPoints,numCp);
for ii = 1:numCp
    [PdtMat(:,ii),PftMat(:,ii),EftMat(:,ii)] = toccgh(Pd.',Pfa, ...
        'ConfirmationThreshold',cp(ii,:),'GateGrowthSequence',gs, ...
        'NumCells',numCells,'NumTargets',10);
end

% Plot ROC and TOC
reqPdt = 0.9; 
helperPlotROCTOC(reqPdt,Pfa,Pd,snrdB,PdtMat,cp);

helperPrintReqValues(reqPdt,Pd,snrdB,PdtMat,EftMat,cp);

Confirmation Threshold            = 2/4
Required Probability of Detection = 0.55
Required SNR (dB)                 = 10.76
Expected Number of False Tracks   = 18

Confirmation Threshold            = 3/4
Required Probability of Detection = 0.81
Required SNR (dB)                 = 12.03
Expected Number of False Tracks   = 1

Confirmation Threshold            = 4/4
Required Probability of Detection = 0.97
Required SNR (dB)                 = 13.37
Expected Number of False Tracks   = 1

Рассматривая результаты, видно, что чем более жесткий порог подтверждения, тем выше необходимый ОСШ. Однако, более строгие пороги подтверждения приводят к улучшению количества ложных треков.

Сводные данные

В оценке эффективности трекера существует четыре типа вероятностей, которые часто желательны:

Для вычисления 1 могут использоваться суммы Бернулли. Однако для получения других вероятностей необходимо использовать другой метод. В то время как анализ Монте-Карло может использоваться для расчета последних трех типов вероятностей, вычислительные ресурсы и необходимое время могут стать несостоятельными, что особенно верно для событий с низкой вероятностью. Алгоритм Common Gate History (CGH) может использоваться, чтобы вычислить все четыре величины и значительно уменьшает вычислительные ресурсы, необходимые.

Алгоритм Common Gate History может использоваться, чтобы сгенерировать рабочую характеристику трекера (TOC). TOC соответствует эксплуатационной характеристики приемника (ROC) и предоставляет средство для оценки общей производительности системы. Кривые TOC и ROC могут использоваться множеством способов, таких как определение:

Наконец, алгоритм CGH позволяет вычислять ожидаемое количество ложных треков, что предлагает представление о емкости системы. Ожидаемое количество ложных дорожек может использоваться, чтобы установить вычислительную нагрузку и помочь с решениями, связанными с оборудованием и обработкой.

Ссылки

Вспомогательные функции

function Pcnf = helperBernoulliSum(Pd,Mc,Nc)
% Calculate simple Bernoulli sum. Use the start TOC logic, which assumes
% that there is already one hit that initializes the track. 

% Update M and N for probability of deletion
Nd = Nc - 1;      % Need one hit to start counting. Assume first hit initializes track. 
Md = Nc - Mc + 1; % Need this many misses to delete

% Bernoulli sum. Probability of deletion calculation.
ii = Md:Nd;
C = arrayfun(@(k) nchoosek(Nd,k),ii);
P = (1 - Pd);
Pdel = sum(C.*P(:).^ii.*(1 - P(:)).^(Nd - ii),2);

% Probability of confirmation
Pcnf = 1 - Pdel;
end

function helperUpdate(calculationType,elapsedTime)
% Output elapsed time 
fprintf('%s calculation completed. Total computation time is %.4f seconds.\n', ...
        calculationType,elapsedTime);
end

function varargout = helperPlot(x,y,displayName,xAxisName,yAxisName,titleName,varargin)
% Create a plot with logarithmic scaling on the x-axis

% Create a figure
figure('Name',titleName)
hAxes = gca;

% Plot data
plot(hAxes,x,y,'LineWidth',2,'DisplayName',displayName,varargin{:})
hold(hAxes,'on')
grid(hAxes,'on')

% Update axes
hAxes.Title.String = titleName;
hAxes.XLabel.String = xAxisName;
hAxes.YLabel.String = yAxisName;

% Make sure legend is on and in best location
legend(hAxes,'Location','Best')

% Set axes as optional output 
if nargout == 1
    varargout{1} = hAxes;
end
end

function varargout = helperPlotLog(x,y,displayName,xAxisName,yAxisName,titleName,varargin)
% Create a plot with logarithmic scaling on the x-axis

% Create a figure
figure('Name',titleName)
hAxes = gca;

% Plot data
numCol = size(y,2);
for ii = 1:numCol
    idxX = min(ii,size(x,2));
    hLine = semilogx(hAxes,x(:,idxX),y(:,ii),'LineWidth',2,varargin{:});
    if ischar(displayName)
        hLine.DisplayName = displayName;
    else
        hLine.DisplayName = sprintf('SNR (dB) = %.2f',displayName(ii));
    end
    hold on
end
grid on

% Update axes
hAxes.Title.String = titleName;
hAxes.XLabel.String = xAxisName;
hAxes.YLabel.String = yAxisName;

% Make sure legend is on and in best location 
legend(hAxes,'Location','Best')

% Set axes as optional output 
if nargout == 1
    varargout{1} = hAxes;
end
end

function helperAddPlotLog(hAxes,x,y,displayName,varargin)
% Add an additional plot to the axes hAxes with logarithmic scaling on the
% x-axis

% Plot data
hold on
hLine = semilogx(hAxes,x,y,'LineWidth',2,varargin{:});
hLine.DisplayName = displayName;
end

function gs = helperCalculateGateSize(N,trkfilt,res)
% Calculate a gate growth sequence in bins

% Initialize tracker gate growth sequence
gs = zeros(1,N); % Gate growth sequence

% Calculate gate growth sequence by projecting state uncertainty using
% linear approximations.
areaBin = prod(res(:),1);
for n = 1:N
    [~,Ppred] = predict(trkfilt); % Predict
    
    % Calculate the products of the 1-sigma values
    E = eig(Ppred);
    E(E<0) = 0; % Remove negative values
    sigma1Prod = sqrt(prod(E(:),1));
    
    % Calculate error ellipse area
    areaErrorEllipse = pi*sigma1Prod; % Area of ellipse = pi*a*b
    
    % Translate to bins
    gs(n) = max(ceil(areaErrorEllipse/areaBin),1);
end
end

function helperColorZonesReqPfa(hAxes,req)
% Plot color zones for requirement type 1 

% Vertical requirement line
xline(req,'--','DisplayName','Requirement',...
    'HitTest','off'); 

% Get axes limits
xlims = get(hAxes,'XLim');
ylims = get(hAxes,'YLim');

% Green box
pos = [xlims(1) ylims(1) req ylims(2)];
x = [pos(1) pos(1) pos(3) pos(3) pos(1)];
y = [pos(1) pos(4) pos(4) pos(1) pos(1)];
hP = patch(hAxes,x,y,[0.4660 0.6740 0.1880], ...
    'FaceAlpha',0.3,'EdgeColor','none','DisplayName','Requirement Met');
uistack(hP,'bottom');

% Red box
pos = [req ylims(1) xlims(2) ylims(2)];
x = [pos(1) pos(1) pos(3) pos(3) pos(1)];
y = [pos(1) pos(4) pos(4) pos(1) pos(1)];
hP = patch(hAxes,x,y,[0.6350 0.0780 0.1840], ...
    'FaceAlpha',0.3,'EdgeColor','none','DisplayName','Requirement Not Met',...
    'HitTest','off');
uistack(hP,'bottom');
end

function helperColorZonesReqPdt(hAxes,req)
% Plot color zones for requirement type 2 

% Horizontal requirement line
yline(req,'--','DisplayName','Requirement',...
    'HitTest','off')

% Get axes limits
xlims = get(hAxes,'XLim');
ylims = get(hAxes,'YLim');

% Green box
pos = [xlims(1) req xlims(2) ylims(2)];
x = [pos(1) pos(1) pos(3) pos(3) pos(1)];
y = [pos(1) pos(4) pos(4) pos(1) pos(1)];
hP = patch(hAxes,x,y,[0.4660 0.6740 0.1880], ...
    'FaceAlpha',0.3,'EdgeColor','none','DisplayName','Requirement Met',...
    'HitTest','off');
uistack(hP,'bottom');

% Red box
pos = [xlims(1) req xlims(2) req];
x = [pos(1) pos(1) pos(3) pos(3) pos(1)];
y = [pos(1) pos(4) pos(4) pos(1) pos(1)];
hP = patch(hAxes,x,y,[0.6350 0.0780 0.1840], ...
    'FaceAlpha',0.3,'EdgeColor','none','DisplayName','Requirement Not Met',...
    'HitTest','off');
uistack(hP,'bottom');
end

function helperPrintTrackProbabilities(Pdt,Pft,Eft)
% Print out results 

fprintf('Probability of Target Track in Presence of False Alarms = %.4f\n',Pdt)
fprintf('Probability of False Track in the Presence of Targets = %.4e\n',Pft)
fprintf('Expected Number of False Tracks = %d\n',Eft)
end

function helperPlotROCTOC(reqPdt,Pfa,Pd,snrdB,PdtMat,cp)
% Plot ROC/TOC 

% Plot ROC curves
figure
hAxesROC = subplot(2,1,1);
plot(hAxesROC,snrdB,Pd,'-o')
title(hAxesROC,'Receiver Operating Characteristic (ROC)')
xlabel(hAxesROC,'SNR (dB)')
ylabel(hAxesROC,'P_D')
grid(hAxesROC,'on')
legend(hAxesROC,sprintf('%.1e',Pfa),'Location','Best')

% Plot SNR versus probability of target track
hAxesTOC = subplot(2,1,2);
numCp = size(cp,1);
for ii = 1:numCp
    plot(hAxesTOC,snrdB,PdtMat(:,ii),'-o', ...
        'DisplayName',sprintf('%d/%d',cp(ii,1),cp(ii,2)))
    hold(hAxesTOC,'on')
end
title(hAxesTOC,'SNR versus P_{DT}')
xlabel(hAxesTOC,'SNR (dB)')
ylabel(hAxesTOC,'P_{DT}')
grid(hAxesTOC,'on')
legend(hAxesTOC,'Location','Best')

% Label points
colorVec = get(hAxesROC,'ColorOrder'); 
numSnr = numel(snrdB);
textArray = arrayfun(@(x) sprintf('  %d',x),1:numSnr,'UniformOutput',false).';
xPosROC = snrdB;
colorFont = brighten(colorVec,-0.75);
numColors = size(colorVec,1);
idxC = mod(1:numCp,numColors); % Use only available default colors
idxC(idxC == 0) = numColors; % Do not let color index equal 0

% Label points ROC
yPosROC = Pd;
text(hAxesROC,xPosROC,yPosROC,textArray,'FontSize',6,'Color',colorFont(1,:),'Clipping','on')

% Label points TOC
xPosTOC = snrdB;
for ii = 1:numCp
    yPosTOC = PdtMat(:,ii);
    text(hAxesTOC,xPosTOC,yPosTOC,textArray,'FontSize',6,'Color',colorFont(idxC(ii),:),'Clipping','on')
end

% Add requirement zone color blocks
helperColorZonesReqPdt(hAxesTOC,reqPdt);
end

function helperPrintReqValues(reqPdt,Pd,snrdB,PdtMat,EftMat,cp)
% Output information about required values given a required probability of
% target track requirement

% Get values 
numCp = size(PdtMat,2);
reqPd = zeros(1,numCp);
reqSNRdB = zeros(1,numCp); 
expEft = zeros(1,numCp);
for ii = 1:numCp
    reqPd(ii) = interp1(PdtMat(:,ii),Pd,reqPdt);
    reqSNRdB(ii) = interp1(PdtMat(:,ii),snrdB,reqPdt);
    expEft(ii) = interp1(PdtMat(:,ii),EftMat(:,ii),reqPdt);
end

 % Display required probability of detection, SNR, and expected false
 % tracks
 for ii = 1:numCp
     fprintf('Confirmation Threshold            = %d/%d\n',cp(ii,1),cp(ii,2));
     fprintf('Required Probability of Detection = %.2f\n',reqPd(ii));
     fprintf('Required SNR (dB)                 = %.2f\n',reqSNRdB(ii));
     fprintf('Expected Number of False Tracks   = %d\n\n',expEft(ii));
 end
end