gofdm_4qam_rayli.m

4 16 64 QAM plots for simulated and theoretical results for OFDM system under a standard Rayleigh ch

%  Title       : 4 QAM Constellation - OFDM - Channels
%  Author      : goutam.ramamurthy.06<AT>student.lth.se
%  History     : 2009-02-12      created
%  References  : lots of material
%  Description :
%
%        1         2         3         4         5         6         7         
%23456789012345678901234567890123456789012345678901234567890123456789012345
%

close all;
clear all;
clc;

%%%%%%%%%% GLOBAL SYSTEM PARAMETERS

%--General Channel Parameters 
ON_noise        = 1;         % 1 = ON; 0 = OFF
ON_esti         = 0;          % 1 = Pilots estimate; 0 = Known CH
channel         = 'PURE_AWGN'; % TID_DISP,RAYL,PURE_AWGN

%--OFDM Parameters 
N_OSyms  = 500;      % Choose total # of symbols in 1 frame to simulateTXRX
N_AvgItr = 1;      % Choose Number of average iterations considered

N_FFT       = 16;      % Hardcoded 
N_c         = N_FFT;  % loading all carriers with data

% keep L_CP less than the maximum dispersion
L_CP = 4;             % Normal CP first symbol. L_CP=160, L_CP_2 = 144;(LTE)

%--Case-Specific Simulation Parameters
idx = 1;
% SIM{idx}.name       = '4QAM';
% SIM{idx}.equ        = 'XX';   % XX = dont care
% SIM{idx}.M          = 4;
% idx = idx+1;

SIM{idx}.name       = '4QAM';
SIM{idx}.equ        = 'ZF';   % XX = dont care
SIM{idx}.M          = N_FFT;

%--Other Common Parameters 
SNRb_dB = (-2:2:16).';     % Eb/N0 range for simulation
N_Snrs = length(SNRb_dB);
EbN0 = 10.^(SNRb_dB/10);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
tic;
for iSIM = 1:idx        % iterating through Simulations
    fprintf('SIMULATION No.%d\n', idx);
    
    M = SIM{iSIM}.M;    % No. of constellation points
    k = log2(M);
    Es = (M-1)*2/3;     % Avg. Ener per QAM sym, norm = sqrt(Es)
    Eb = Es/k;
    d_minsq = 3*log2(M)/(M-1); % for M-QAM BER theo calc
        
    N_QSyms = N_c * N_OSyms ;
    N_bits  = N_QSyms * k;        % an integral multiple of k channels
    
    s = (1:sqrt(M)/2);        % In rect QAMs, for M=2^n & n ~ {even I}
    a_MQAM = [-(2*s-1) (2*s-1)];  % alphabets

    for iSNR = 1:N_Snrs
        fprintf('Eb/N0: %d dB\n', SNRb_dB(iSNR));        
        N0 =Eb/EbN0(iSNR);
            
        N_totsymerr = 0;
        N_totbiterr = 0;
        
        %%%%RND SRC:Generate a column of random binary bits
        txbits = (1-sign(randn(N_bits, 1)))/2; % ~Bernoulli src(p=1/2)
        
        for iN_AvgItr = 1:N_AvgItr
            
            %%%%QAM MAPPER
            %--using matlab's modem object generate QAM syms and normalize
            qam_mod = modem.qammod(M);    % create a modulator object
            qam_mod.symbolorder = 'gray' ;  % Enable/disable gray coding
            qam_mod.inputtype = 'bit';
            %qam_mod.disp;                   % Print the properties of modulator
            q_txsyms = modulate(qam_mod, txbits);  % QAM modulate
            q_txsyms_norm = q_txsyms/sqrt(Es);
            
            %%%%OFDM MOD
            O_txseq= [];
            for i = 1:N_OSyms                 % iterate over N OFDMsyms
                offset = (i-1)*N_c;
                txsyms_ext = q_txsyms_norm(offset+1:offset+N_c);
                O_sym = (sqrt(N_FFT))*ifft(txsyms_ext, N_FFT);
                %--ADD CP  
                O_txseq = [O_txseq;O_sym(N_FFT-L_CP+1:end);O_sym];
            end

            %%%%THE CHANNEL TO CONVOLVE:
            noise_len = size(O_txseq);
            re = randn(noise_len);
            im = randn(noise_len);
            awgn = 1/sqrt(2)*(re + j*im); % WGN with 0dB variance
            if strcmp(channel, 'RAYL') ~= 0
                h = (randn(noise_len,1)+j*randn(noise_len,1))/sqrt(2);
                rxdata = h.*O_txseq + ON_noise*sqrt(N0/2)*awgn;
                %%%% needs lot of realizations to give good results so,
                %%%% have atleast N_AvgItr = 100
%                 h = (randn+j*randn)/sqrt(2);
%                 rxdata = filter(h,1,q_txsyms_norm) + ON_noise*sqrt(N0/2)*awgn;
            else % if strcmp(SIM{idx}.ch, 'PURE_AWGN')~=0
                rxdata = O_txseq + ON_noise*sqrt(N0/2)*awgn;
            end


            
            %%%%RECEIVER =============>
            %--ESTI and EQZ first
            if (ON_esti)
                %%% do pilot estimation and stuff here
            else
                if strcmp(channel, 'RAYL')~=0
                    H = fft(h, N_FFT);
                else % if strcmp(SIM{idx}.ch, 'PURE_AWGN')~=0
                    H = ones(N_FFT,1);
                end
            end
            
            if strcmp(SIM{idx}.equ, 'ZF')~=0
                EQZ = 1./H;
                %EQZ =(1./H)/sqrt(sum(abs(H).^2)); 
                EQZ = EQZ(:);
            end
            
            %%%%OFDM DEMOD
            q_rxsyms = [];            % pre-declaration
            for i = 1:N_OSyms     % iterate over N OFDMsyms
                offset = (i-1)*N_FFT + (i*L_CP);
                %--REM CP
                O_rxsym = (1/sqrt(N_FFT))*fft(rxdata(offset+1:offset+N_FFT),N_FFT);
                if strcmp(SIM{idx}.equ, 'ZF')~=0
                    O_rxsym = EQZ.*O_rxsym;                      
                end  
                q_rxsyms = [q_rxsyms; O_rxsym(1:N_FFT)];
            end % N_OFDMsyms 
            
            if strcmp(channel, 'PURE_AWGN')~=0 %denormalize->amplifies s+n
                q_rxsyms_denorm = q_rxsyms*sqrt(Es); 
            else
                q_rxsyms_denorm = q_rxsyms;
            end
            
            %%%%RECT QAM MAPPING:
            %--Segregate Re-Im
            q_rxsyms_re = real(q_rxsyms_denorm);
            q_rxsyms_im = imag(q_rxsyms_denorm);
            %--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);
            q_rxsyms = rx_re+j*rx_im;
            
            qam_demod = modem.qamdemod(qam_mod);  % copy props from mod object
            %qam_demod.disp;                     % Prints properties of DEMOD
            rxbits = demodulate(qam_demod, q_rxsyms_denorm); % QAM modulate
            
            N_totsymerr = N_totsymerr + nnz(q_rxsyms(:)-q_txsyms(:));
            N_totbiterr = N_totbiterr + nnz(rxbits(:)-txbits(:));
        end % for iN_AvgItr 
        
        %--Simulated SER
        N_symerrs = N_totsymerr/N_AvgItr;
        SER(iSIM,iSNR) = N_symerrs/N_QSyms;
        
        %--Simulated BER
        N_biterrs = N_totbiterr/N_AvgItr;
        BER(iSIM,iSNR) = N_biterrs/N_bits;
    end  % for iSNR
    
    %%%%THEORETICAL BER CALC:
    if strcmp(channel, 'RAYL') ~= 0
        %--BER(TH) for a rayleigh fading channel
        ber = berfading(SNRb_dB,'qam',M,1);
        BER_theo(iSIM, :)= ber(:);
    else% if strcmp(SIM{iSIM}.ch, 'PURE_AWGN') ~= 0
        %--SER (TH)
        Qfn = qfunc(sqrt(d_minsq * EbN0));
        for iSER = 1 : length(Qfn)
            SER_theo(iSIM,iSER) = (4*(1-(1/sqrt(M)))*Qfn(iSER)) - ...
                                  (4*(1-(1/sqrt(M)))^2*(Qfn(iSER)^2));
        end
        BER_theo(iSIM, :) = (1/k)*SER_theo(iSIM, :); %--BER (TH)
    end
end  % for iSIM
toc;


%%%%Performance Plots:
for iSIM = 1:length(SIM)     % iterating through Simulations
    figure(iSIM) % BER v/s EbN0
    hdl = semilogy(SNR_dB, SER_theo(iSIM,:),'b:','LineWidth',2);
    hold on
    hdl = semilogy(SNRb_dB, SER(iSIM,:), 'rs-', 'Linewidth',1);
    grid on
    xlabel('Eb/No, dB');
    ylabel('BER(P_b)');
    title('BER Performance');
end
axis([-2 26 10^-7 1]);