% plots illustrating the 16ASK constellation in DVB-S2

gamma = 1/(1-sqrt(2-sqrt(3))); %sqrt(2/(2-sqrt(3)));
R1 = 1;
R2 = gamma*R1;

% signal energies
Es = R1^2*(1 + 3*gamma^2)/4;
Eb = Es/4;

ss = zeros(1, 16);

for nn = 1:4
    ss(nn) = R1*exp(1j*pi/4*(2*nn-1));
end

for nn = 5:16
    ss(nn) = R2*exp(1j*pi/12*(2*(nn-4)-1));
end

% show decision regions
figure(1)
voronoi(real(ss), imag(ss))
hold on
plot(ss, 'ro')
hold off
axis('equal')
title('Decision Regions')

% Find Pe
% a) there are foour relevant disance between signal points
d1 = sqrt(2)*R1;                       % distance between points on inner ring
d2 = (gamma-1)*R1;                     % distance between points on diagonal
d3 = sqrt(2-sqrt(3))*gamma*R1;         % dist b/w points on outer ring
d4 = sqrt(1+gamma^2-sqrt(3)*gamma)*R1; % dist from inner ring to off-diagonal outer

% b) translate to energy efficiency
eta1 = d1^2/Eb;
eta2 = d2^2/Eb;
eta3 = d3^2/Eb;
eta4 = d4^2/Eb;

eta_min = min([eta1, eta2, eta3, eta4]);

% c) compute conditional error probabilities
EbN0 = logspace(0, 1.5, 100);

Pe_inner = 2*Q(sqrt(eta1*EbN0/2)) + Q(sqrt(eta2*EbN0/2)) + 2*Q(sqrt(eta4*EbN0/2));
Pe_diag = Q(sqrt(eta2*EbN0/2)) + 2*Q(sqrt(eta3*EbN0/2));
Pe_off  = 2*Q(sqrt(eta3*EbN0/2)) + Q(sqrt(eta4*EbN0/2));

% d) average the conditional Pe
Pe = Pe_inner/4 + Pe_diag/4 + Pe_off/2;

figure(2)
semilogy(10*log10(EbN0), Pe, '-', ...
    10*log10(EbN0), 3*Q(sqrt(1.6*EbN0/2)), '--')
grid on
xlabel('E_b/N_0 (dB)')
ylabel('P_e')
legend('16APSK', 'QPSK')

title_str = sprintf('\\gamma = %4.2f, \\eta_{min} = %4.2f', gamma, eta_min);
title(title_str)