clear all
close all
clc

%%
% Summary:
%   Script used for empirically determine the ROC curve of an
%   acquisition system.
%   A single correlator output is modeled as a complex Gaussian random
%   variable. Under the H1 hyphotesis, the correlator mean is determined
%   from the input C/N0 and the coherent integration time. The variance of
%   the real and imaginary parts of thecorrelator outputs is normalized to 
%   be unity.
%
% Author:
%   Daniele Borio
%%

CN0 = 35;   % dB-Hz
Tc = 0.001; % ms -coherent integration time
SNRp = CN0 + 10 * log10( 2 * Tc );  % post-integration SNR
snrp = 10^( SNRp / 10 );    % express the SNR in linear unit
sigma2 = 1;                 % variance of the in-phase and quadrature components

lambda = snrp * sigma2 ; % non-centrality parameters
N = 10000;   % Number of simulated random variables

% Simulate the correlator outputs
CH0 = sqrt(sigma2) * ( randn(1, N) + j * randn( 1, N ) );   % Under CH0
CH1 = CH0 + sqrt(lambda) * exp( j * 2 * pi * rand(1, N));   % Under CH1 
                                                      % The signal
                                                      % component is
                                                      % affected by a
                                                      % random phase
% Compute the decision variable
DH0 = abs( CH0 ).^2;
DH1 = abs( CH1 ).^2;

% Determine the decision threshold
Tmax = max( DH1 );
Tmin = 0;
Nt = N / 10;    % Number of threshold points
T = linspace( Tmin, Tmax, Nt );

% Compute the false alarm and detection probabilities
Pfa = zeros( size(T) );
Pd = zeros( size(T) );

for ii = 1:Nt,
    Pfa( ii ) = length( find( DH0 > T(ii) ) ) / N;
    Pd( ii ) = length( find( DH1 > T(ii) ) ) / N;
end

% Plot the ROC
loglog( Pfa, Pd, 'O','LineWidth', 2 )   % use log-log scale
grid on
axis([1/N 1 1/N 1])
set( gca, 'FontSize', 16 )
xlabel('False Alarm Probability')
ylabel('Detection Probability')

% Theoretical curves
Pfa = exp(-T/2);
Pd = marcumq( sqrt(lambda / sigma2 ), sqrt(T / sigma2 ) );
hold on
plot(Pfa, Pd, 'r', 'LineWidth', 2)