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

Бернуллиевая сумма

Бернуллиевая сумма допускает быстрый и легкий анализ эффективности в случае единой цели в отсутствие ложных предупреждений. Вероятность обнаружения дорожки $$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 требуется.

Вычисление вероятностей для цели в помехе

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

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

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

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

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

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

Алгоритм Общей истории логического элемента (CGH)

Общий алгоритм Истории Логического элемента значительно уменьшает время вычисления и требования к памяти. Алгоритм предотвращает потребность в ручной генерации последовательностей как в случае анализа Монте-Карло, который может быть дорогостоящим для событий низкой вероятности.

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

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

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

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

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

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

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

Предположения об Общем алгоритме Истории Логического элемента следующие:

Общий алгоритм Истории Логического элемента может использоваться, чтобы вычислить все четыре типа вероятности:

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

Вычислите и постройте TOC с помощью кривых ROC от rocsnr как входные параметры. Примите отношение сигнал-шум (SNR) 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');

Общий алгоритм Истории Логического элемента разрешает оценку эффективности средства отслеживания, похожей на анализ Монте-Карло, но с приемлемым временем вычисления несмотря на события низкой вероятности. Общий алгоритм Истории Логического элемента таким образом разрешает высокоуровневое расследование и выбор опций до более интенсивных, подробных симуляций.

Используя 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 для SNRs, равного 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 функция, чтобы сгенерировать связанные кривые TOC. Используйте тот же порог подтверждения и пропустите последовательность роста, ранее сгенерированную.

% 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 выбрана. Это - компромисс улучшенной вероятности целевой дорожки за счет немного более высокой, но разумной вероятности ложной дорожки.

Общий алгоритм Истории Логического элемента разрешает оценку ожидаемого количества ложных дорожек на основе количества целей, ожидаемых в среде и количестве ячеек в радарных данных. Ожидаемое количество ложных дорожек ${\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, могут использоваться Бернуллиевые суммы. Однако, чтобы получить другие вероятности, различный метод должен использоваться. В то время как анализ Монте-Карло может использоваться для расчета последних трех типов вероятностей, вычислительные ресурсы и требуемое время могут стать ненадежными, который особенно верен для событий с низкой вероятностью. Алгоритм Общей истории логического элемента (CGH) может использоваться, чтобы вычислить все четыре количества и значительно уменьшает вычислительные необходимые ресурсы.

Общий алгоритм Истории Логического элемента может использоваться, чтобы сгенерировать Рабочую характеристику средства отслеживания (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