function hamming_demo

% 	decoder = 'Hard';
	decoder = 'Syndrome';
% 	decoder = 'Soft';
% 	decoder = 'SISO';
	
	% parametro del codice
	NmK = 3;
	N = 2^NmK-1;
	K = N-NmK;
	
	temp = de2bi(0:2^NmK-1);
	P = temp(sum(temp,2)>=2,:);
	G=[eye(K),P];
	Ht = [P;eye(N-K)];
	
	% tabella dei rappresentanti di laterale
	l_tab = gen_l_tab(Ht);
	
	% tabelle delle parole informative e di codice
	d_tab = de2bi(0:2^K-1);
	c_tab = mod(d_tab*G,2);
	
	% rapporto segnale rumore
	SNR_dB = 5;
	sigma_w = sqrt(10^(-SNR_dB/10));
	
	tic
	n_err = 0;
	n_loop = 0;
	n_ev = 0;
	n_ev_target = 3000;
	while n_ev < n_ev_target
		
		% sorgente
		d = d_tab(ceil(2^K*rand),:);
		
		% codifica
		c = mod(d*G,2);
		
		% mapping BPSK
		v = 1-2*c;
		
		% canale AWGN
		w = sigma_w*randn(1,N);
		r = v+w;

		switch decoder
			case 'Hard'
				
				% decisioni all'uscita del canale
				b = r<0;
				
				% calcola le distanze di Hamming
				dH = sum(mod(c_tab+repmat(b,2^K,1),2),2);
				
				% decisioni MV (Hard)
				[val,pos]=min(dH);
				d_est = d_tab(pos,:);
			
			case 'Syndrome'
				
				% decisioni all'uscita del canale
				b = r<0;
				
				% decodifica a sindrome
				s = mod(b*Ht,2);
				is = s * 2.^(0:NmK-1)'; % bi2de(s);
				c_est = mod(b+l_tab(1+is,:),2);
				
				% la parte sistematica
				d_est = c_est(1:K);
				
			case 'Soft'
				
				% calcola le distanze euclidee
				d2=sum(abs(repmat(r,2^K,1)-(1-2*c_tab)).^2,2);
				
				% decisioni MV (Soft)
				[val,pos]=min(d2);
				d_est = d_tab(pos,:);
				
			case 'SISO'
				
				% demapping LLR
				Lc = 2*r/sigma_w^2;
				
				% Hamming SISO
				Ld_a = zeros(1,K);
				Ld = hamming_siso(Ld_a,Lc,d_tab,c_tab);
				%Ld_e = Ld-Ld_a;
				
				% decisioni Hard
				d_est = Ld<0;
		end
		
		n_loop = n_loop + 1;
		err = sum(mod(d+d_est,2));
		n_err = n_err + err;
		if err
			n_ev = n_ev + 1;
		end
	end
	fprintf('SNR = %gdB - ER: %g (%d/%d) - BER: %g (%d/%d)\n',...
		SNR_dB,n_ev/n_loop,n_ev,n_loop,...
		n_err/(n_loop*K),n_err,n_loop*K);
	toc
	
end

% decodifica SISO
%
% TODO: implementare il metodo di Hartmann-Rudolph (basato sul calcolo
%       mediante la somma sulle parole del codice duale)
%
function Ld = hamming_siso(Ld_a,Lc,d_tab,c_tab)
	
	% implementare il calcolo degli LLR :)
	
	K = size(Ld_a,2);
	Ld = zeros(1,K);
	for bit = 1:K
		numeratore = 0;
		denominatore = 0;
		for i = 1:2^K
			if d_tab(i,bit) == 0
				numeratore = numeratore + exp((Lc*(0.5-c_tab(i,:))') + (Ld_a*(0.5-d_tab(i,:))'));
			elseif d_tab(i,bit) == 1
				denominatore = denominatore + exp((Lc*(0.5-c_tab(i,:))') + (Ld_a*(0.5-d_tab(i,:))'));
			end
		end
		Ld(bit) = log(numeratore/denominatore);
	end
end

% generazione della tabella dei rappresentanti di laterale
%
% TODO: generalizzare a codici lineari arbitrari (per ora funziona
%       solo per i codici di Hamming)
%
function l_tab = gen_l_tab(Ht)

	% inizializzazione
	[N,NmK] = size(Ht);
	l_tab = zeros(2^NmK,N);

	% flag per indicare che devono essere assegnate
	l_tab(2:end,1) = -1;
	
	% prova tutte le stringhe di peso 1
	for i1 = 1:N
		s = Ht(i1,:);
		is = bi2de(s);
		if l_tab(1+is,1) < 0
			% correggi l'inizializzazione
			l_tab(1+is,1) = 0;
			% assegna il bit del rappresentante
			l_tab(1+is,i1) = 1;
		end
	end
	
% 	% prova tutte le stringhe di peso 2
% 	for i1 = 1:N-1
% 		for i2 = i1+1:N
% 			s = mod(Ht(i1,:)+Ht(i2,:),2);
% 			is = bi2de(s);
% 			if l_tab(1+is,1) < 0
% 				% correggi l'inizializzazione
% 				l_tab(1+is,1) = 0;
% 				% assegna i bit del rappresentante
% 				l_tab(1+is,i1) = 1;
% 				l_tab(1+is,i2) = 1;
% 			end
% 		end
% 	end
end
