% Various routines that simulate or compute aspects of digital communication systems
%
% Simulate matched filter receiver
figure(1);
N = 100;
noise_amp = 3;
signal_set = 'bpsk'; % either 'bpsk', 'ask', or 'fsk'
bits = ['1', '0', '1', '0', '0', '1'];
%
if strcmp(signal_set, 'bpsk')
	signal1 = ones(1,N);
	signal0 = -signal1;
elseif strcmp(signal_set,'fsk')
	signal1 = sqrt(2)*sin(2*pi*2*[0:N-1]/N);
	signal0 = sqrt(2)*sin(2*pi*3*[0:N-1]/N);
elseif strcmp(signal_set, 'ask');
	signal1 = ones(1,N);
	signal0 = zeros(1,N);
else
	perror(sprintf('Unknown signal set %s\n',signal_set));
end
color0='r';color1='b';
x = []; xcolor = [];
for n=1:length(bits)
	x=[x eval(strcat('signal',bits(n)))];
	xcolor = [xcolor eval(strcat('color',bits(n)))];
end
% Send signal through white noise channel
r = x + noise_amp*randn(1,length(x));
% Run matched filters
y1=filter(signal1(N:-1:1),1,r);
y0=filter(signal0(N:-1:1),1,r);
% Graphics
subplot(211)
t=[0:length(r)-1];
plot(t,r);hold on
a = axis;
xp=x*(0.75*max(abs([a(3) a(4)])/max(x)));
for n=1:length(bits)
	plot(t((n-1)*N+1:n*N),xp((n-1)*N+1:n*N),[xcolor(n) '--']);
	h = text((n-1)*N+N/2,max(xp),bits(n));
	set(h,'fontsize',16);set(h,'color',xcolor(n));
end
for n=N*[1:length(bits)],h=line([n n],a(3:4));set(h,'linestyle','--');end
h=title('Received signal');set(h,'fontsize',18);
hold off
subplot(212)
plot(t,y0,color0,t,y1,color1)
a = axis;
for n=1:length(bits)
	if y1(n*N)>= y0(n*N)
		h = text(n*N-10,.75*a(4),'1');
		set(h,'fontsize',16);set(h,'color',color1);
		if bits(n) == '0'
			set(h,'fontweight','bold');
		end
	else 
		h = text(n*N-10,.75*a(4),'0');
		set(h,'fontsize',16);set(h,'color',color0);
		if bits(n) == '1'
			set(h,'fontweight','bold');
		end
	end
end
for n=N*[1:length(bits)],h=line([n n],a(3:4));set(h,'linestyle','--');end
h=title('Matched Filter Output');set(h,'fontsize',18);
%
% Compute Pr[e] curves
%
figure(2)
snrdb = [-10:.5:12];
snr = 10.^(snrdb/10);
p_bpsk = Q(sqrt(2*snr));
p_fsk = Q(sqrt(snr));
h=semilogy(snrdb,p_bpsk,snrdb,p_fsk,'r--');grid;axis([-10 12 10^(-8) 1])
set(gca,'fontsize',18);
h=xlabel('Signal-to-Noise Ratio (dB)');set(h,'fontsize',18);
h=ylabel('Bit Error Probability');set(h,'fontsize',18);
legend('BPSK','FSK');
%
% Error Correction (Repetition Code)
%
figure(3)
snrdb = [-10:.5:12];
snr = 10.^(snrdb/10);
p_bpsk = Q(sqrt(2*snr));
pb_bpsk = p_bpsk;
pb_rep = 1-(1-p_bpsk).^3-3*p_bpsk.*(1-p_bpsk).^2;
p_31 = Q(sqrt(2*snr/3)); % Same transmission transmission time
pb_31 = p_31;
semilogy(snrdb,pb_bpsk,'b-',snrdb,pb_31,'r--',snrdb,pb_rep,'k-.');
grid;axis([-10 12 10^(-15) 1])
set(gca,'fontsize',18);
h=xlabel('Signal-to-Noise Ratio (dB)');set(h,'fontsize',18);
h=ylabel('Block Error Probability');set(h,'fontsize',18);
legend('No ECC','Same time/data bit','Length 3 Repetition Code')
title('Repetition Code (3,1)')
%
% Hamming Code Performance Curves
%
figure(4)
snrdb = [-10:.5:12];
snr = 10.^(snrdb/10);
p_bpsk = Q(sqrt(2*snr));
p_31 = Q(sqrt(2*1/3*snr));
p_74 = Q(sqrt(2*4/7*snr));
p_1511 = Q(sqrt(2*11/15*snr));
p_3126 = Q(sqrt(2*26/31*snr));
pb_uc1 = p_bpsk;
pb_uc4 = 1-(1-p_bpsk).^4;
pb_uc11 = 1-(1-p_bpsk).^11;
pb_uc26 = 1-(1-p_bpsk).^26;
pb_hamming31 = 1-(1-p_31).^3-3*p_31.*(1-p_31).^2;
pb_hamming74 = 1-(1-p_74).^7-7*p_74.*(1-p_74).^6;
pb_hamming1511 = 1-(1-p_1511).^15-15*p_1511.*(1-p_1511).^10;
pb_hamming3126 = 1-(1-p_3126).^31-31*p_3126.*(1-p_3126).^30;

clf
subplot(2,2,1);
semilogy(snrdb,pb_uc1,'k-',snrdb,pb_hamming31,'k-.');
grid;axis([-10 12 10^(-10) 1]);
h=ylabel('Block Error Probability'); set(h,'fontsize',18);
h = title('(3,1) Hamming Code'); set(h,'fontsize',18);
set(gca,'fontsize',18);

subplot(2,2,2);
semilogy(snrdb,pb_uc4,'b-',snrdb,pb_hamming74,'b-.');
grid;axis([-10 12 10^(-10) 1]);
h = title('(7,4) Hamming Code'); set(h,'fontsize',18);
set(gca,'fontsize',18);

subplot(2,2,3);
semilogy(snrdb,pb_uc11,'r-',snrdb,pb_hamming1511,'r-.');
grid;axis([-10 12 10^(-10) 1]);
h=ylabel('Block Error Probability'); set(h,'fontsize',18);
h=xlabel('Signal-to-Noise Ratio (dB)'); set(h,'fontsize',18);
h = title('(15,11) Hamming Code'); set(h,'fontsize',18);
set(gca,'fontsize',18);

subplot(2,2,4);
semilogy(snrdb,pb_uc26,'g-',snrdb,pb_hamming3126,'g-.');
grid;axis([-10 12 10^(-10) 1]);
h=xlabel('Signal-to-Noise Ratio (dB)');% set(h,'fontsize',18);
h = title('(31,26) Hamming Code'); set(h,'fontsize',18);
set(gca,'fontsize',18);
legend('Uncoded','Hamming Code',3);
%
% Capacity Calculations
%
figure(5)
p=[0:.01:.5];
C=1+p.*log2(p)+(1-p).*log2(1-p);
C(1) = 1;
subplot(211)
plot(p,C)
set(gca,'fontsize',18);grid
h = xlabel('Error Probability');set(h,'fontsize',18);
h = ylabel('Capacity (bits)');set(h,'fontsize',18);
subplot(212)
C=1+p_bpsk.*log2(p_bpsk)+(1-p_bpsk).*log2(1-p_bpsk);
plot(snrdb,C);
set(gca,'fontsize',18);grid;legend('Capacity using BPSK')
h = xlabel('Signal-to-Noise Ratio (dB)');set(h,'fontsize',18);
h = ylabel('Capacity (bits)');set(h,'fontsize',18);