% Try:
%------
% N=10
% s2=5
%------
% N=1
% s2=0.5
%------
% First N=60 s2=200 in Figure 1
% and then N=60 s2=10 in figure 2
% Compare with theoretical result


N=60;
s2=10;
figure(2)
hold off
snr=10*log10(1/s2);
Niter=500000;


s=ones(N,1);
s=s./sqrt(s'*s);
s=s*sqrt(N);
x=randn(N,Niter)*sqrt(s2);
H0=s'*x;

x=s*ones(1,Niter)+randn(N,Niter)*sqrt(s2);
H1=s'*x;


T=linspace(min(H0)-1,max(H1)+1,100);
P=[];
for t=1:length(T)
    Pfa=length(find(H0>T(t)))/Niter;
    Pd=length(find(H1>T(t)))/Niter;
    P=[P [Pd;Pfa]];
end
plot(P(2,:),P(1,:),'b','LineWidth',2)
hold on;


x=randn(N,Niter)*sqrt(s2);
H0=sum(x.^2,1);
x=randn(N,Niter)*sqrt(1+s2);
H1=sum(x.^2,1);

T=linspace(min(H0)-1,max(H1)+1,100);
P=[];
for t=1:length(T)
    Pfa=length(find(H0>T(t)))/Niter;
    Pd=length(find(H1>T(t)))/Niter;
    P=[P [Pd;Pfa]];
end

plot(P(2,:),P(1,:),'g','LineWidth',2)

x=randn(N,Niter)*sqrt(s2);
H0=sum(x.^2,1);

x=s*ones(1,Niter)+randn(N,Niter)*sqrt(s2);
H1=sum(x.^2,1);

T=linspace(min(H0)-1,max(H1)+1,100);
P=[];
for t=1:length(T)
    Pfa=length(find(H0>T(t)))/Niter;
    Pd=length(find(H1>T(t)))/Niter;
    P=[P [Pd;Pfa]];
end

plot(P(2,:),P(1,:),'r','LineWidth',2)




xlabel('P_F_A','FontSize',18,'FontWeight','bold')
ylabel('P_D','FontSize',18,'FontWeight','bold')
grid on
legend({'Chapter 4 - known signal','Chapter 5 - random signal','Chapter 7 - unknown signal'},'FontSize',18,'FontWeight','bold','Location','southeast')
title([ '\fontsize{18}N = ',num2str(N),',    A^2/\sigma^2 =',num2str(snr),' dB'])