%BinaryComm
%
% Demonstration of Stark and Woods Example 1.6-1
% Brian D. Jeffs, 9/12/05

M = 10000;                          % number of symbols to average over
L = 50;                             % number of symbols to plot
sigma2 = input('Enter noise power (variance): ');

Px1 = input('Enter probability that X=1: ');

T = input('Enter detection threshold, T= ? ');

X = rand(M,1) < Px1;            %Generate binary sequence
N = sqrt(sigma2)*randn(M,1);    %Generate noise

W = X + N;
Y = W > T;                          %Threshold the noisy signal

plot(X(1:L),'o');
axis([0 L+1 -.5 1.5]);
title('Binary source signal, X');
xlabel('Symbol index, i.e. time sample n');

pause

plot(1:L,X(1:L),'o',1:L,W(1:L),'+');
axis([0 L+1 -.5 1.5]);
title('Signal with noise');
xlabel('Symbol index, i.e. time sample n');
legend('X','W');

pause

plot(1:L,X(1:L),'o',1:L,W(1:L),'+',1:L,Y(1:L),'d',[0,L+1],[T,T]);
axis([0 L+1 -.5 1.5]);
title('Thesholded received signal');
xlabel('Symbol index, i.e. time sample n');
legend('X','W','Y','T');

PY1gX1 = sum(Y.*X)/sum(X);
PY1gX0 = sum(Y.*~X)/sum(~X);
PY0gX0 = sum(~Y.*~X)/sum(~X);
PY0gX1 = sum(~Y.*X)/sum(X);
PY1X1 = PY1gX1*Px1;
PY1X0 = PY1gX0*(1-Px1);
PY0X0 = PY0gX0*(1-Px1);
PY0X1 = PY0gX1*Px1;

disp('Conditional probability computed from relative frequency:');
disp(['P[Y=1|X=1] = ',num2str(PY1gX1)]);
disp(['P[Y=1|X=0] = ',num2str(PY1gX0)]);
disp(['P[Y=0|X=1] = ',num2str(PY0gX1)]);
disp(['P[Y=0|X=0] = ',num2str(PY0gX0)]);
disp(' ');
disp(['P[Y=1,X=1] = ',num2str(PY1X1)]);
disp(['P[Y=1,X=0] = ',num2str(PY1X0)]);
disp(['P[Y=0,X=1] = ',num2str(PY0X1)]);
disp(['P[Y=0,X=0] = ',num2str(PY0X0)]);
disp(' ');
disp(['Total error prob. = P[Y=1,X=0] + P[Y=1,X=0] = ',num2str(PY1X0+PY0X1)]);