clear all;
Eb = 1;
Eb_NO_dBs = -10:2:10;
nb_symbols = 100000;

modulation = 'qam';
modulation_length = 4;

SER = zeros(1,length(Eb_NO_dBs));
SER_theory = zeros(1,length(Eb_NO_dBs)); 

iter = 1;
for Eb_N0_dB = Eb_NO_dBs
    Eb_N0 = 10^(Eb_N0_dB/10);
    
    N0 = Eb/Eb_N0;
    
    %% Modulations
    Es = Eb*log2(modulation_length);

    if modulation=='pam'
        symbols_index = floor(rand(1,nb_symbols)*modulation_length);
        x = 2*symbols_index-modulation_length+1;
        x = x*sqrt(Es*3/(modulation_length^2-1));
    elseif modulation=='qam'
        symbols_index_real = floor(rand(1,nb_symbols)*sqrt(modulation_length));
        symbols_index_imag = floor(rand(1,nb_symbols)*sqrt(modulation_length));
        x = (2*symbols_index_real-sqrt(modulation_length)+1) ...
            + sqrt(-1)*(2*symbols_index_imag-sqrt(modulation_length)+1);
        x = x*sqrt(Es*3/2*1/(modulation_length-1));
    end

    %% Transmission, ML decision and Symbol Error Rate

    z = sqrt(N0)/sqrt(2)*(randn(1,nb_symbols)+sqrt(-1)*randn(1,nb_symbols));
    y = x + z;
    
    if modulation=='pam'
        denorm_y = (y/sqrt(Es*3/(modulation_length^2-1))+modulation_length-1)/2;
        hat_x = max(min(round(real(denorm_y)),modulation_length-1),0);

        SER(iter) = sum(abs(hat_x-symbols_index)>.25)/nb_symbols;
        SER_theory(iter) = 2*(modulation_length-1)/modulation_length*qfunc(sqrt(6*log2(modulation_length)/(modulation_length^2-1)*Eb_N0));
        
    elseif modulation=='qam'
        denorm_y = (y/sqrt(Es*3/2*1/(modulation_length-1))+sqrt(modulation_length)-1+sqrt(-1)*(sqrt(modulation_length)-1))/2;
        hat_x_real = max(min(round(real(denorm_y)),modulation_length-1),0);
        hat_x_imag = max(min(round(imag(denorm_y)),modulation_length-1),0);
        hat_x = hat_x_real+sqrt(-1)*hat_x_imag;
        
        SER(iter) = 1-sum((abs(hat_x_real-symbols_index_real)<.25).*(abs(hat_x_imag-symbols_index_imag)<.25))/nb_symbols;
        SER_theory(iter) = 1 - (1-2*(sqrt(modulation_length)-1)/sqrt(modulation_length)*qfunc(sqrt(3*log2(modulation_length)/(modulation_length-1)*Eb_N0)))^2;
    end
    
    iter = iter + 1;
end

%% Comparative plots
figure;
hold on;
grid on;
plot(Eb_NO_dBs,log10(SER),'b');
plot(Eb_NO_dBs,log10(SER_theory),'r');