% si immagina di trasmettere (con modulazione BPSK su canale AWGN
% e decisioni Hard) messaggi di MESSAGE_LENGTH bit segmentati in
% parole di codice di Hamming o di Hamming esteso
%
% vogliamo confrontare:
%  - probabilità di errore sul bit in funzione di E_b/N_0
%  - throughput inteso come rate di consegna di messaggi completi
%    in funzione del rapporto segnale rumore sul canale
%
% risulta che:
%  - con decisioni Hard, i codici estesi hanno prestazioni
%    peggiori, sia in termini di BER (tentare di correggere
%    quando si sono rivelati due errori è più dannoso che altro
%    per il BER), sia in termini di throughput (i codici estesi
%    hanno rate minore)

function Hamming_demo
	
	%     Hamming      ExtHamming
	%
	% N-K |  N  |  K     N  | N-K
	% ----+-----+----- -----+-----
	%  3  |  7  |  4     8  |  4
	%  4  | 15  | 11    16  |  5
	%  5  | 31  | 26    32  |  6
	
	% confronto tra codici di Hamming e di Hamming estesi
	N_minus_K_Hamming = 3;
	
	% costruisci le parti non sistematiche della matrice di
	% controllo di parità del codice di Hamming e del
	% corrispondente codice di Hamming esteso
	P_Hamming = Hamming_Parity_Check(N_minus_K_Hamming);
	P_ExtHamming = ExtHamming_Parity_Check(P_Hamming);
	CL_ExtHamming = ExtHamming_CosetLeaders(P_ExtHamming);
	K = size(P_Hamming,1);
	R_Hamming = K / (K + size(P_Hamming,2));
	R_ExtHamming = K / (K + size(P_ExtHamming,2));
	
	% lunghezza del messaggio
	MESSAGE_LENGTH = 32; % 4 bytes
	FRAMES_PER_MESSAGE = ceil(MESSAGE_LENGTH/K);
	
	% rapporti segnale-rumore sul canale da simulare
	SNR_dB_min = -5;
	SNR_dB_max = 10;
	SNR_dB_step = 0.5;
	
	% durata della simulazione
	seconds_budget = inf;
	max_trials = 10000;
	
	% attendi 'frame_errors_threshold' eventi di errore per
	% fare la stima delle probabilità
	frame_errors_threshold = 300;
	
	% vettori per le figure
	x1_Hamming = [];
	x1_ExtHamming = [];
	y1_Hamming = [];
	y1_ExtHamming = [];
	x2 = [];
	y2_uncoded = [];
	y2_Hamming = [];
	y2_ExtHamming = [];
	y2_Shannon = [];
	
	% loop principale sul rapporto segnale rumore
	SNR_dB = SNR_dB_min;
	tic;
	while SNR_dB <= SNR_dB_max && toc < seconds_budget
		
		fprintf('SNR: %4g dB:',SNR_dB);
		
		% probabilità di errore di BPSK + AWGN + decisioni Hard
		SNR = 10^(SNR_dB/10);
		p = Q(sqrt(SNR));
		
		% contatori di eventi
		trials = 0;
		frame_errors_Hamming = 0;
		bit_errors_Hamming = 0;
		frame_errors_ExtHamming = 0;
		bit_errors_ExtHamming = 0;
		
		while (frame_errors_Hamming < frame_errors_threshold || ...
				frame_errors_ExtHamming < frame_errors_threshold) && ...
				trials <= max_trials && toc < seconds_budget
			
			% sorgente
			u = rand(1,K) < 0.5;
			
			% codifica
			c_Hamming = encode(u,P_Hamming);
			c_ExtHamming = encode(u,P_ExtHamming);
			
			% canale BSC che modella
			r_Hamming = BSC(p,c_Hamming);
			r_ExtHamming = BSC(p,c_ExtHamming);
			
			% decodifica
			u_est_Hamming = Hamming_decode(r_Hamming,P_Hamming);
			u_est_ExtHamming = ExtHamming_decode(r_ExtHamming,...
				P_ExtHamming,CL_ExtHamming);
			
			% conteggio eventi d'errore
			bit_errors = sum(u_est_Hamming ~= u);
			if bit_errors
				frame_errors_Hamming = ...
					frame_errors_Hamming + 1;
				bit_errors_Hamming = bit_errors_Hamming ...
					+ bit_errors;
			end
			bit_errors = sum(u_est_ExtHamming ~= u);
			if bit_errors
				frame_errors_ExtHamming = ...
					frame_errors_ExtHamming + 1;
				bit_errors_ExtHamming = bit_errors_ExtHamming ...
					+ bit_errors;
			end
			trials = trials + 1;
		end

		if frame_errors_Hamming >= frame_errors_threshold && ...
			frame_errors_ExtHamming >= frame_errors_threshold
			fprintf(' OK\n');
		elseif trials > max_trials
			fprintf(' Max trials reached (%d and %d errors / %d)\n',...
				frame_errors_Hamming,frame_errors_ExtHamming,...
				frame_errors_threshold);
		else
			fprintf(' Timeout (%d and %d errors / %d)\n',...
				frame_errors_Hamming,frame_errors_ExtHamming,...
				frame_errors_threshold);
		end
		
		% stima i rate di errore sul bit
		BER_Hamming = bit_errors_Hamming/(K*trials);
		BER_ExtHamming = bit_errors_ExtHamming/(K*trials);
		
		% a parità di rapporto segnale-rumore sul canale, il
		% rapporto energia di bit / densità spettrale monolatera
		% del rumore cambia col rate di codifica
		Eb_N0_dB_Hamming    = SNR_dB - 10*log10(2*R_Hamming);
		Eb_N0_dB_ExtHamming = SNR_dB - 10*log10(2*R_ExtHamming);
		
		% nella figura 1 rappresentiamo la probabilità di errore
		% sul bit in funzione di E_b/N_0
		x1_Hamming = [x1_Hamming,Eb_N0_dB_Hamming];
		x1_ExtHamming = [x1_ExtHamming,Eb_N0_dB_ExtHamming];
		y1_Hamming = [y1_Hamming,BER_Hamming];
		y1_ExtHamming = [y1_ExtHamming,BER_ExtHamming];
		
		% probabilità di errore sul frame
		FER_uncoded = 1-(1-p)^K;
		FER_Hamming = frame_errors_Hamming/trials;
		FER_ExtHamming = frame_errors_ExtHamming/trials;
		
		% nella figura 2 rappresentiamo il numero medio di bit
		% consegnati con successo (in un messaggio di lunghezza
		% MESSAGE_LENGTH) in funzione del rapporto segnale rumore
		% sul canale
		x2 = [x2,SNR_dB];
		y2_uncoded = [y2_uncoded,...
			(1-FER_uncoded)^FRAMES_PER_MESSAGE];
		y2_Hamming = [y2_Hamming,...
			R_Hamming*(1-FER_Hamming)^FRAMES_PER_MESSAGE];
		y2_ExtHamming = [y2_ExtHamming,...
			R_ExtHamming*(1-FER_ExtHamming)^FRAMES_PER_MESSAGE];
		y2_Shannon = [y2_Shannon,1+p*log2(p)+(1-p)*log2(1-p)];
		% y2_Shannon = [y2_Shannon,0.5*log2(1+SNR)];

		SNR_dB = SNR_dB + SNR_dB_step;
	end
	
	% le prestazioni del sistema non codificato sono prevedibili
	% teoricamente, quindi le calcoliamo qui al volo
	x1_uncoded = sort([x1_Hamming,x1_ExtHamming]);
	y1_uncoded = Q(sqrt(2*10.^(x1_uncoded/10)));
	
	subplot(2,1,1)
 	hold off
	semilogy(x1_uncoded,y1_uncoded)
	hold on
	semilogy(x1_Hamming,y1_Hamming,'g')
	semilogy(x1_ExtHamming,y1_ExtHamming,'r')
	grid on
	xlabel('E_b/N_0(dB)')
	ylabel('Bit Error Rate')
	legend('Uncoded','Hamming','Extended Hamming')
	
	subplot(2,1,2)
 	hold off
	plot(x2,y2_uncoded)
	hold on
	plot(x2,y2_Hamming,'g')
	plot(x2,y2_ExtHamming,'r')
	plot(x2,y2_Shannon,'k--')
	grid on
	xlabel('Channel Signal-to-Noise Ratio (dB)')
	ylabel('Throughput (bits/channel use)')
	legend('Uncoded','Hamming','Extended Hamming','Shannon')



% la funzione Q(x) = Pr{N(0,1) > x}
function val = Q(x)
	val = 0.5*erfc(x/sqrt(2));



% il canale binario simmetrico
function y = BSC(p,x)
	y = x;
	for i = 1:numel(x)
		if rand < p
			y(i) = 1 - x(i);
		end
	end



% costruisci la parte non sistematica della matrice dei controlli
% di parità del codice di Hamming
function P = Hamming_Parity_Check(N_minus_K)
	
	N = 2^N_minus_K - 1;
	K = N - N_minus_K;
	
	% le righe da considerare sono tutte e sole quelle con un
	% numero di uno maggiore o uguale a due
	P = zeros(K,N_minus_K);
	temp = dec2bin(1:N)-'0';
	next_row = 1;
	for n = 1:N
		if sum(temp(n,:)) >= 2
			P(next_row,:) = temp(n,:);
			next_row = next_row + 1;
		end
	end



% codifica sistematica
function c = encode(u,P)
	c = [u,mod(u*P,2)];



% calcolo della sindrome
function s = syndrome(r,P)
	K = size(P,1);
	s = mod(r(1:K)*P+r(K+1:end),2);



% decodifica del codice di Hamming
function u_est = Hamming_decode(r,P)
	[K, N_minus_K]=size(P);
	
	% calcola la sindrome
	s = syndrome(r,P);
	
	% decodifica a sindrome:
	u_est = r(1:K);
	switch sum(s)
		case 0
			% r è una parola di codice
		case 1
			% la sindrome indica un errore nei controlli di
			% parità, non è necessario cambiare nulla
		otherwise
			% cerca la sindrome tra le righe di P; l'indice della
			% riga è la posizione dell'errore
			for k = 1:K
				if sum(s == P(k,:)) == N_minus_K
					u_est(k) = 1 - r(k);
					break
				end
			end
	end



% la parte non sistematica della matrice dei controlli di parità
% del codice esteso si ottiene aggiungendo un controllo di parità
% semplice alle righe di P (+1 in quanto c'e' anche la parte
% sistematica)
function P = ExtHamming_Parity_Check(P)
	P = [P,mod(sum(P,2)+1,2)];



% costruzione della tabella dei coset leaders (di peso fino a
% due, il codice di Hamming esteso ha distanza quattro)
function CL = ExtHamming_CosetLeaders(P)
	[K, N_minus_K] = size(P);
	N = K + N_minus_K;
	
	H = [P; eye(N_minus_K)];
	CL = zeros(2^N_minus_K - 1, N);
	
	% loop sulle sindromi (i è il valore della
	% stringa di bit della sindrome)
	for i = 1:2^N_minus_K - 1
		s = dec2bin(i,N_minus_K)-'0';
		
		found = false;
		
		% prova a trovare la sindrome tra le
		% righe di H (coset leader di peso 1)
		for n = 1:N
			if sum(s == H(n,:)) == N_minus_K
				CL(i,n) = 1;
				found = true;
				break
			end
		end
		
		if found
			continue
		end
		
		% altrimenti, prova a trovare la sindrome
		% tra le combinazioni a due a due di
		% righe di H
		for n1 = 1:N-1
			for n2 = n1+1:N
				temp = mod(H(n1,:)+H(n2,:),2);
				if sum(s == temp) == N_minus_K
					CL(i,n1) = 1;
					CL(i,n2) = 1;
					found = true;
					break
				end
			end
			if found
				break
			end
		end
	end



% decodifica del codice di Hamming esteso (decidere una parola
% quando si è rivelato che gli errori sono DUE è probabilmente
% sbagliato, vedere la decodifica per codici prodotto)
function u_est = ExtHamming_decode(r,P,CL)
	[K, N_minus_K]=size(P);
	
	% calcola la sindrome
	s = syndrome(r,P);
	
	% decodifica sistematica
	u_est = r(1:K);
	
	switch sum(s)
		case {0,1}
			% r è una parola di codice oppure l'errore più
			% probabile è un errore singolo nei controlli di
			% parità
		otherwise
			% converti la sindrome in numero e correggi col
			% coset leader
			s = s * (2.^(N_minus_K-1:-1:0))';
			u_est = mod(u_est + CL(s,1:K),2);
	end
