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

Бернулли Сум

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

${\mathit{}}_{\mathit{Tw}}\mathit{=}$NT,

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

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

${\mathit{}}_{\mathrm{}}{}_{\mathit{}\mathit{}}^{\mathit{}}\mathit{Pdt=\sum i=MNC}\mathit{(}N\mathit{,}{\mathit{i}}_{\mathrm{)}}^{\mathit{}}{\mathrm{}Pdi{\mathit{}}_{\mathrm{(}}}^{\mathit{}1-Pd\mathit{)}}$(N-i),

где

$\mathit{C}\mathit{(}N\mathit{,}i\frac{\mathit{)}}{=\mathit{}N\mathit{!}(\mathit{}}$N-i)! 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)

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

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

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

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

где ${}_{\mathit{startl}}$ - количество временных шагов с момента последнего обнаружения (цель или ложная тревога), ${}_{\mathit{startlt}}$- количество временных шагов с момента последнего обнаружения цели, λ $$- общее число обнаружений (цели или ложные тревоги). По мере продвижения алгоритма вектор состояния дорожки развивается в соответствии с цепью Маркова.

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

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

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

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

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

Вероятность целевой дорожки и вероятность ложной дорожки составляют основу рабочей характеристики трекера (TOC). (TOC) дополняет рабочую характеристику приемника (ROC), которая обычно используется при анализе и прогнозировании производительности приемников. Объединение ROC и 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');

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

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

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

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

$\mathit{Q}=\begin{array}{cc}\frac{{\mathit{[}}^{\mathrm{}}}{\mathrm{}}& \frac{{\mathit{}}^{\mathrm{}}}{\mathrm{\Delta t44\Delta t3}}\\ \frac{{\mathit{}}^{\mathrm{}}}{\mathrm{}}& {\mathit{}}^{\mathrm{}}\end{array}\mathit{2\Delta t32\Delta t2}$] q,

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

$\mathit{q}{\mathit{=}}_{\mathrm{}}^{\mathrm{amax}}\mathit{}$2Δt.

Входное ${\mathit{}}_{\mathrm{amax}}^{\mathrm{}}$2 - это максимальное ожидаемое целевое ускорение. Предположим, что для $${}^{\mathrm{}}$$транспортных средств ожидается максимальное ускорение 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]);

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

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

$$[\sqrt{{}_{\mathrm{\lambda \; 1}}},\sqrt{{}_{\mathrm{}}}\lambda \; 2]\sqrt{={}_{\mathrm{eig}\mathrm{(}\mathrm{Pk}}}$$+ 1 | k).

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

$\mathrm{}\mathrm{}\mathrm{Площадь\; эллипса\; ошибки}=\sqrt{\lambda {}_{\mathrm{\lambda }}}\sqrt{1{}_{\mathrm{\lambda }}}$2.

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

$\mathrm{}\mathrm{Площадь\; ячейки}=\mathrm{(}\Delta Range)\mathrm{(}\mathrm{\Delta Range}$Rate).

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

$\mathrm{}\mathrm{Размер\; литника}\frac{\mathrm{=}\mathrm{Ошибка}\mathrm{эллипса}}{\mathrm{}\mathrm{Область\; корзины}}$.

Таким образом, размер затвора зависит от трекера, последовательности событий и разрешения ячеек. Вычислите последовательность роста затвора, предполагая, что порог подтверждения 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 для SNR, равных 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 могут быть включены для дальнейшего анализа, поскольку они удовлетворяют вероятности ниже ${\mathrm{}}^{\mathrm{10-6}}$. Пункты 7-10 не соответствуют этому требованию.

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

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

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

${\mathit{}}_{\mathrm{Eft}}{\mathit{=}}_{\mathrm{}\mathrm{Pft,}}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{ntNc}}{\mathit{+}}_{\mathrm{}}$PftNt,

где ${\mathit{}}_{\mathrm{Pft},\mathrm{}}$nt - вероятность ложной дорожки при отсутствии целей, ${\mathit{}}_{\mathrm{}}$Nc - число ячеек, ${\mathit{}}_{\mathrm{}}$Pft - вероятность ложной дорожки при наличии целей, а ${\mathit{}}_{\mathrm{}}$Nt - число целей.

Рассмотрим среду, в которой число целей должно быть равно 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

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

Резюме

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

Для вычисления 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