%  Title       : QAM Constellation plots for 4, 16, 64
%  Author      : goutam.ramamurthy.06<AT>student.lth.se
%  History     : 2009-01-30      created
%  References  :
%  Description :
%
%        1         2         3         4         5         6         7         
%23456789012345678901234567890123456789012345678901234567890123456789012345
%
%

close all;
clear all;
clc;

%%%%%%%%%% GLOBAL SYSTEM PARAMETERS 
graycode_on = 1;

%--Case-Specific Simulation Parameters
SIMU{1}.name = '4QAM';
SIMU{1}.M_ary = 4;
SIMU{1}.legend = ' Coherent 4QAM AWGN Simulation';
SIMU{1}.linestyle = '-gd'; % linestyle: green Squares
SIMU{1}.title ='Ps(Sym Err) curve for 4QAM modulation';

SIMU{2}.name = '16QAM';
SIMU{2}.M_ary = 16;
SIMU{2}.legend = ' Coherent 16QAM AWGN Simulation';
SIMU{2}.linestyle = 'rd-'; % linestyle: red diamonds
SIMU{2}.title ='Ps(Sym Err) curve for 16QAM modulation';

SIMU{3}.name = '64QAM';
SIMU{3}.M_ary = 64;
SIMU{3}.legend = ' Coherent 64QAM AWGN Simulation';
SIMU{3}.linestyle = '-md'; % linestyle: maroon diamonds
SIMU{3}.title ='Ps(Sym Err) curve for 64QAM modulation';

%--Common Parameters 
SNRdB = 0:2:24;       % Es/N0 range for simulation

%~~~~~~~~~~~~~~

tic;
for iSIMU = 1:length(SIMU)     % iterating through Simulations
    
    M_ary = SIMU{iSIMU}.M_ary; % number of constellation points
    N_bits = 8e4;
    Es = (M_ary-1)*2/3;          % Avg. Ener per QAM sym, norm =sqrt(Es)
    d_minsq = 3*log2(M_ary)/(M_ary-1); % for M-QAM
    k = log2(M_ary);
    
    m = (1:sqrt(M_ary)/2);        % In rect QAMs, for M=2^k & k ~ {even I}
    a_MQAM = [-(2*m-1) (2*m-1)];  % alphabets
    
    
    for iSNRdB = 1:length(SNRdB) % iterating through SNRs
        fprintf('\nEs/N0: %d dB\n', SNRdB(iSNRdB));
        SNR = 10^(SNRdB(iSNRdB)/10); % Convert to non-dB value
    
        %%%%%%%%%% RANDOM SOURCE
        x_txbits = (1-sign(randn(N_bits,1)))/2; % ~Bernoulli source(p=1/2)
        
        %--Generate random symbols and normalize it to unit symbol energy
        tx_re = randsrc(1,N_bits,a_MQAM);
        tx_im = randsrc(1,N_bits,a_MQAM);
        q_txsyms = (tx_re + j*tx_im);
        q_txsyms_norm = q_txsyms/sqrt(Es);
        %%%RANDOM SOURCE-ED
        
        %++++++++++
        %%%%%%%%%% THE CHANNEL TO CONVOLVE
        re = randn(1,N_bits);
        im = randn(1,N_bits);
        noise = 1/sqrt(2)*(re + j*im); % WGN with 0dB variance
        q_rxsyms_noisy = q_txsyms_norm + (1/sqrt(SNR))*noise; % AWGN
        %%%CHANNEL CONVOLUT-ED 
        %++++++++++
        
        %%%%%%%%%% QAM DEMODULATOR (rectangular)
        %--Segregate Re-Im & De-Norm Rx symbols
        q_rxsyms_re = real(q_rxsyms_noisy)* sqrt(Es);
        q_rxsyms_im = imag(q_rxsyms_noisy)* sqrt(Es);
        %--Mapping to the nearest alphabet
        rx_re = floor(q_rxsyms_re/2)*2+1;
        rx_im = floor(q_rxsyms_im/2)*2+1;
        %--mapping boundary conditions for RE and IM
        rx_re(find(max(a_MQAM)<rx_re)) = max(a_MQAM);
        rx_re(find(min(a_MQAM)>rx_re)) = min(a_MQAM);
        rx_im(find(max(a_MQAM)<rx_im)) = max(a_MQAM);
        rx_im(find(min(a_MQAM)>rx_im)) = min(a_MQAM);
        %%%QAM DEMODULATED -ED
        
        q_rxsyms = rx_re+j*rx_im;
        N_symerrs = nnz(q_rxsyms(:)-q_txsyms(:));
        SER(iSIMU,iSNRdB) = N_symerrs/N_bits;
        
        %--Calculating theoretical SER
        EbN0 = SNR/k; % EbN0 = EsN0/k
        Qfn = q(sqrt(d_minsq * EbN0));
        SER_theo(iSIMU, iSNRdB) = (4*(1-(1/sqrt(M_ary)))*Qfn) - ...
                                  (4*(1-(1/sqrt(M_ary)))^2*(Qfn^2));
                
        %-- Scatter plots
        figure(iSIMU*10);
        plot(q_rxsyms,'.k');
        title(SIMU{iSIMU}.legend);
        axis square;
        axis equal;
        grid on;
    end
end
toc;
%%%%%%%%%% SER PLOTS
figure(1)
for iSIMU = 1:length(SIMU)     % iterating through Simulations
    semilogy(SNRdB, SER_theo(iSIMU, :),'b:','LineWidth',2);
    hold on
    semilogy(SNRdB, SER(iSIMU, :), SIMU{iSIMU}.linestyle, 'Linewidth',1);
    hold on
end
axis([0 30 10^-6 1]);
grid on
legend('theoretical', 'Simulated');
xlabel('Es/No, dB');
ylabel('SER(P_s)');
title('Ps versus E_b/N_0');