clear all; close all;

P = 10;
N = 1;

y = 3*sqrt(P+N)*linspace(-1,1);

p0 = 1/sqrt(2*pi)*1/sqrt(N)*exp(-y.^2/(2*N));
p1 = 1/sqrt(2*pi)*1/sqrt(P+N)*exp(-y.^2/(2*(P+N)));
plot(y,p0); hold on;
plot(y,p1,'r');

% Compute threshold
yth = sqrt((P+N)*N/P*log((P+N)/N));
plot(yth*[1 1],[0 1/sqrt(2*pi*N)],'k--');
plot(-yth*[1 1],[0 1/sqrt(2*pi*N)],'k--');

% SNR
snr = linspace(0.1,10);
p_err1 = 2*(1-normcdf(sqrt(1+snr.^-1).*log(1+snr)));
p_err2 = normcdf(sqrt(snr.^-1).*log(1+snr),0,1)-...
         normcdf(-sqrt(snr.^-1).*log(1+snr),0,1);
figure; 
p_err = 1/2*(p_err1+p_err2);
plot(snr,p_err);
xlabel('SNR');
ylabel('Probability of error');