% Simulation demonstrating diversity gain
%
% (c) 2005 Klaus Witrisal
% Technical University Graz, Austria
% Signal Processing and Speech Communication Lab
% witrisal@tugraz.at; www.spsc.tugraz.at

N_d = 1000;	% number of data symbols sent over each channel
N_ch = 1000;	% number of channel simulations
E_s = 1;	% symbol energy
N_0 = 10.^(0:-.2:-3);	% vector of noise powers

% parameters for QPSK
N_e_av = 2;			% 2 nearest neighbors
d_min = sqrt(2);	% min distance for dominant error event


rho_v = E_s./N_0;		% average SISO SNR


%%%%%%%%%%%%%%% theoretical BER over AWGN channel

% note: Q(x) = 0.5*erfc(x/sqrt(2))
SER_AWGN = N_e_av*0.5*erfc(sqrt(rho_v*d_min^2/2)/sqrt(2));

figure(1); clf
semilogy(10*log10(rho_v),SER_AWGN)
xlabel 'E_s/N_0 [dB]'
ylabel 'symbol error rate (SER)'
axis([0 30 1e-20 1])
hold on

disp 'approx. exact SER computed by Q-function for AWGN channel'
pause

% Chernoff bound
SER_AWGN_cb = 1/sqrt(2*pi)*N_e_av*exp(-rho_v*d_min^2/4);

semilogy(10*log10(rho_v),SER_AWGN_cb,'--')

disp 'upper bound on SER computed by Chernoff-bound for AWGN channel'
pause


%%%%%%%%%%%%%%% simulation of BER over AWGN channel
N = 10000;	% nb. of symbols for AWGN simulation

% transmitted QPSK symbols 
s = sqrt(E_s/2) * (sign(randn(1,N)) + j*sign(randn(1,N)));

% loop over SNR
SER_AWGN_sim = [];
for sig_n2 = N_0

	rho = E_s/sig_n2

	% noise vector at variance sig_n2
	n = sqrt(sig_n2/2)*(randn(1,N) + j*randn(1,N));
	
	% received signal
	r = s + n;
	
	% determine SER
	SER = 1- mean((sign(real(r)) == sign(real(s))) & ... 
		(sign(imag(r)) == sign(imag(s))));
		
	SER_AWGN_sim = [SER_AWGN_sim SER];
	
	
	if rho == 10
		figure(2), clf
		plot(r,'.g')
		hold on
		plot(r.*((sign(real(r)) ~= sign(real(s))) | ... 
		(sign(imag(r)) ~= sign(imag(s)))),'r.')
		plot(s,'o')
		axis equal
		grid on
		%pause
	end

end

figure(1)	
semilogy(10*log10(rho_v),SER_AWGN_sim,'o')
axis([0 30 1e-8 1])

disp 'simulated SER for AWGN channel'
pause


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% for SISO Rayleigh channel
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%% theoretical SER over Rayleigh channel

% exact equation for BER (see Proakis, (14-3-7))
BER_Rayl = 0.5*(1-sqrt(rho_v./(1+rho_v)));

% convert to approx. eqn. for SER
SER_Rayl = N_e_av*0.5*(1-sqrt(rho_v*d_min^2/4./(1+rho_v*d_min^2/4)));

figure(1)	
semilogy(10*log10(rho_v),SER_Rayl,'r-')
axis([0 30 1e-5 1])

disp 'analytical SER for Rayleigh channel'
pause

% Chernoff bound
SER_Rayl_cb = 1/sqrt(2*pi)*N_e_av./(1+rho_v*d_min^2/4);

semilogy(10*log10(rho_v),SER_Rayl_cb,'r--')

disp 'upper bound on SER computed by Chernoff-bound for Rayleigh channel'
pause




%%%%%%%%%%%%%%% simulation of SER over Rayleigh channel

% transmitted QPSK symbols 
s = sqrt(E_s/2) * (sign(randn(1,N_d)) + j*sign(randn(1,N_d)));

% loop over SNR
SER_Rayl_sim = [];
for sig_n2 = N_0

	rho = E_s/sig_n2		% average SISO SNR
	
	% noise vector at variance sig_n2
	n = sqrt(sig_n2/2)*(randn(1,N_d) + j*randn(1,N_d));

	%loop over channel realizations
	SER = 0;
	for i=1:N_ch;

		% channel transfer coefficient (complex Gaussian)
		h = sqrt(1/2)*(randn(1) + j*randn(1));

		% received signal
		y = h*s + n;
		
		% channel "equalization": 
		z = exp(-j*angle(h))*y;
	
		% determine SER
		S_err = ((sign(real(z)) ~= sign(real(s))) | ... 
		(sign(imag(z)) ~= sign(imag(s))));
		SER = SER + mean(S_err);
		
		
		% visualize received symbols
		if (rho == 10) & (i<10)
			figure(2), clf
			plot(y,'.g')
			hold on
			plot(real(y.*S_err),imag(y.*S_err),'r.')
			plot(h*s,'o')
			axis([-3 3 -3 3])
			axis equal
			grid on
			pause
		end

	end
		
	SER_Rayl_sim = [SER_Rayl_sim SER/N_ch];
	
	
end

figure(1)	
semilogy(10*log10(rho_v),SER_Rayl_sim,'ro')
axis([0 30 1e-5 1])

disp 'simulated SER for Rayleigh channel'
pause


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% for SIMO Rayleigh channel; M = 2 diversity branches
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%% analytical SER over M = 2 diversity branches Rayleigh channel

M = 2;	% number of diversity branches

% Chernoff bound
SER_M2_cb = 1/sqrt(2*pi)*N_e_av*(rho_v*d_min^2/4/M).^(-M);

semilogy(10*log10(rho_v),SER_M2_cb,'g--','color',[0 .6 0])

disp 'upper bound on SER computed by Chernoff-bound for Rayleigh channel; M = 2 diversity branches'
pause




%%%%%%%%%%%%%%% simulation of SER over M = 2 diversity branches Rayleigh channel

% transmitted QPSK symbols 
s = sqrt(E_s/2/M) * (sign(randn(1,N_d)) + j*sign(randn(1,N_d)));

% loop over SNR
SER_M2_sim = [];
for sig_n2 = N_0

	rho = E_s/sig_n2		% average SISO SNR
	
	% noise vector at variance sig_n2
	n = sqrt(sig_n2/2)*(randn(M,N_d) + j*randn(M,N_d));

	%loop over channel realizations
	SER = 0;
	for i=1:N_ch;

		% channel transfer coefficient (complex Gaussian)
		h = sqrt(1/2)*(randn(M,1) + j*randn(M,1));

		% received signal
		y = h*s + n;
		
		% maximum ratio combining 
		z = h'*y;
	
		% determine SER
		S_err = ((sign(real(z)) ~= sign(real(s))) | ... 
		(sign(imag(z)) ~= sign(imag(s))));
		SER = SER + mean(S_err);
		
		
		% visualize received symbols
		if (rho == 10) & (i<10)
			figure(2), clf
			for m=1:M
				subplot(1,M,m)
				plot(y(m,:),'.g')
				hold on
				plot(real(y(m,:).*S_err),imag(y(m,:).*S_err),'r.')
				plot(h(m)*s,'o')
				axis([-2 2 -2 2])
				axis equal
				title(['diversity branch ' int2str(m)])
				grid on
			end
			figure(3), clf
			plot(z,'.g')
			hold on
			plot(real(z.*S_err),imag(z.*S_err),'r.')
			%plot(h(m)*s,'o')
			axis([-3 3 -3 3])
			axis equal
			grid on
			title 'combined channels'
			
			pause
		end

	end
		
	SER_M2_sim = [SER_M2_sim SER/N_ch];
	
	
end

figure(1)	
semilogy(10*log10(rho_v),SER_M2_sim,'go','color',[0 .6 0])
axis([0 30 1e-5 1])

disp 'simulated SER for Rayleigh channel with diversity; M = 2'
pause


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% for SIMO Rayleigh channel; M = 4 diversity branches
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%% analytical SER over M = 4 diversity branches Rayleigh channel

M = 4;	% number of diversity branches

% Chernoff bound
SER_M4_cb = 1/sqrt(2*pi)*N_e_av*(rho_v*d_min^2/4/M).^(-M);

semilogy(10*log10(rho_v),SER_M4_cb,'m--')

disp 'upper bound on SER computed by Chernoff-bound for Rayleigh channel; M = 4 diversity branches'
pause




%%%%%%%%%%%%%%% simulation of SER over M = 2 diversity branches Rayleigh channel

% transmitted QPSK symbols 
s = sqrt(E_s/2/M) * (sign(randn(1,N_d)) + j*sign(randn(1,N_d)));

% loop over SNR
SER_M4_sim = [];
for sig_n2 = N_0

	rho = E_s/sig_n2		% average SISO SNR
	
	% noise vector at variance sig_n2
	n = sqrt(sig_n2/2)*(randn(M,N_d) + j*randn(M,N_d));

	%loop over channel realizations
	SER = 0;
	for i=1:N_ch;

		% channel transfer coefficient (complex Gaussian)
		h = sqrt(1/2)*(randn(M,1) + j*randn(M,1));

		% received signal
		y = h*s + n;
		
		% maximum ratio combining 
		z = h'*y;
	
		% determine SER
		S_err = ((sign(real(z)) ~= sign(real(s))) | ... 
		(sign(imag(z)) ~= sign(imag(s))));
		SER = SER + mean(S_err);
		
		
		% visualize received symbols
		if (rho == 10) & (i<10)
			figure(2), clf
			for m=1:M
				subplot(1,M,m)
				plot(y(m,:),'.g')
				hold on
				plot(real(y(m,:).*S_err),imag(y(m,:).*S_err),'r.')
				plot(h(m)*s,'o')
				axis([-2 2 -2 2])
				axis equal
				title(['diversity branch ' int2str(m)])
				grid on
			end
			figure(3), clf
			plot(z,'.g')
			hold on
			plot(real(z.*S_err),imag(z.*S_err),'r.')
			%plot(h(m)*s,'o')
			axis([-3 3 -3 3])
			axis equal
			grid on
			title 'combined channels'
			
			pause
		end

	end
		
	SER_M4_sim = [SER_M4_sim SER/N_ch];
	
	
end

figure(1)	
semilogy(10*log10(rho_v),SER_M4_sim,'mo')
axis([0 30 1e-5 1])

disp 'simulated SER for Rayleigh channel with diversity; M = 4'
legend('AWGN exact','AWGN bound','AWGN simul','SISO; exact','SISO; bound','SISO; simul','M = 2; bound','M = 2; simul','M = 4; bound','M = 4; simul',3)

pause

