BER graph MATLAB calibration

Question

I'm doing a 16QAM system (transmitter, channel and receiver), and BER and PER curves of the results. However, I'm having some problems with noise at the receiver. I'm running the system inside two loops: for all the Eb/No values and for all the packets and I sent 200 symbols and 1000 packets but this still happens. I would like to check whether the result from this code is correct or not:

clear all
clc
numPkts=1000;

N = 200; % number of symbols
M = 16;   % constellation size
k = log2(M); % bits per symbol
pv=4; %prefix length


% defining the real and imaginary PAM constellation
% for 16-QAM
alphaRe = [-(2*sqrt(M)/2-1):2:-1 1:2:2*sqrt(M)/2-1];
alphaIm = [-(2*sqrt(M)/2-1):2:-1 1:2:2*sqrt(M)/2-1];
k_16QAM = 1/sqrt(10);

Eb_N0_dB  = [0:15]; % multiple Es/N0 values
Es_N0_dB  = Eb_N0_dB + 10*log10(k);
erTot=zeros(1,length(Eb_N0_dB));

% Mapping for binary <--> Gray code conversion
ref = [0:k-1];
map = bitxor(ref,floor(ref/2));
[tt ind] = sort(map);                                

for ii = 1:length(Eb_N0_dB)
for pktX=1:numPkts    
% symbol generation
% ------------------
ipBit = rand(1,N*k,1)>0.5; % random 1's and 0's
ipBitReshape = reshape(ipBit,k,N).';
bin2DecMatrix = ones(N,1)*(2.^[(k/2-1):-1:0]) ; % conversion from binary to decimal
% real
ipBitRe =  ipBitReshape(:,[1:k/2]);
ipDecRe = sum(ipBitRe.*bin2DecMatrix,2);
ipGrayDecRe = bitxor(ipDecRe,floor(ipDecRe/2));
% imaginary
ipBitIm =  ipBitReshape(:,[k/2+1:k]);
ipDecIm = sum(ipBitIm.*bin2DecMatrix,2);
ipGrayDecIm = bitxor(ipDecIm,floor(ipDecIm/2)); 
% mapping the Gray coded symbols into constellation
modRe = alphaRe(ipGrayDecRe+1);
modIm = alphaIm(ipGrayDecIm+1);
% complex constellation
mod = modRe + j*modIm;
s1 = k_16QAM*mod; % normalization of transmit power to one 

s=[s1(length(s1)-pv+1:end) s1]; %add prefix


% noise
% -----
EsNo=10^(Es_N0_dB(ii)/10);
stanDevNoise=sqrt((1)/(2*EsNo));

n =stanDevNoise *[randn(1,length(s)) + j*randn(1,length(s))]; % white guassian noise, 0dB variance 



h=(1/sqrt(2))*(randn+j*randn);
y1= conv(s,h) + n; % additive white gaussian noise



%removes prefix
        y1(1:pv) = [];   

y=y1/h;
% demodulation
% ------------
y_re = real(y)/k_16QAM; % real part
y_im = imag(y)/k_16QAM; % imaginary part

% rounding to the nearest alphabet
ipHatRe = 2*floor(y_re/2)+1;
ipHatRe(find(ipHatRe>max(alphaRe))) = max(alphaRe);
ipHatRe(find(ipHatRe<min(alphaRe))) = min(alphaRe);
ipHatIm = 2*floor(y_im/2)+1;
ipHatIm(find(ipHatIm>max(alphaIm))) = max(alphaIm);
ipHatIm(find(ipHatIm<min(alphaIm))) = min(alphaIm);

% Constellation to Decimal conversion
ipDecHatRe = ind(floor((ipHatRe+4)/2+1))-1; % LUT based
ipDecHatIm = ind(floor((ipHatIm+4)/2+1))-1; % LUT based

% converting to binary string
ipBinHatRe = dec2bin(ipDecHatRe,k/2);
ipBinHatIm = dec2bin(ipDecHatIm,k/2);

% converting binary string to number
ipBinHatRe = ipBinHatRe.';
ipBinHatRe = ipBinHatRe(1:end).';
ipBinHatRe = reshape(str2num(ipBinHatRe).',k/2,N).' ;

ipBinHatIm = ipBinHatIm.';
ipBinHatIm = ipBinHatIm(1:end).';
ipBinHatIm = reshape(str2num(ipBinHatIm).',k/2,N).' ;

% counting errors for real and imaginary
nBitErr(pktX) = size(find([ipBitRe- ipBinHatRe]),1) + size(find([ipBitIm - ipBinHatIm]),1) ;

end
erTot(ii)=erTot(ii)+sum(nBitErr); %total errors in all packets

simBer(ii)=(erTot(ii)/(N*k*numPkts)); %bit error rate

totPktErRate(ii)=(erTot(ii)/(numPkts)); 
end

theoryBer = (1/k)*3/2*erfc(sqrt(k*0.1*(10.^(Eb_N0_dB/10))));

close all; figure
semilogy(Eb_N0_dB,theoryBer,'bs-','LineWidth',2);
hold on
semilogy(Eb_N0_dB,simBer,'mx-','LineWidth',2);
axis([0 15 10^-5 1])
grid on
legend('theory', 'simulation');
xlabel('Eb/No, dB')
ylabel('Bit Error Rate')
title('Bit error probability curve for 16-QAM modulation')

Thanks!

Answer 1

The code provided makes the following assumptions:

• 16-QAM modulation using Gray-coding bit mapping
• a flat slow/block Rayleigh fading channel model.
• coherent decoding under perfect channel state information estimation

Due to it's similarity with the Additive-White-Gaussian-Noise (AWGN) channel, a logical first step in understanding and calibrating the system performance under the assumptions stated above is to evaluate its performance without fading (i.e. substituting the channel model with an AWGN channel by setting h=1 in the provided code).

AWGN channel

You may want to verify the calibration of Symbol-Error-Rate (SER) performance as this can have a large impact on the (BER) performance, and SER curves are readily available for coherent decoding of uncoded 16-QAM constellation (see e.g dsplog, these lecture slides, this book, etc.) Those references also include the following approximation to the SER of 16-QAM:

1.5*erfc(sqrt(EsN0/10))

where EsN0 = 10.^(0.1*EsN0_dB).

Note that results may be equivalently provided in terms of either Es/N0 (the average energy per symbol) or Eb/N0 (the average energy per bit). For a k-bits signal constellation (constellation size of 2^k), the relationship between Es/N0 and Eb/N0 is given as

Es/N0 = k*Eb/N0

Thus for 16-QAM, Es/N0 = 4Eb/N0 (or Es/N0_dB = Eb/N0_dB + 6dB).

For a Gray coded scheme, a BER approximation for sufficiently high Eb/N0 can then be obtained from the fact that a symbol error translates into 1 bit in error (out of the k-bits in the symbol) most of the time, thus BER ~ SER/k (or again for 16-QAM: BER ~ SER/4).

Es/N0 (dB)   Eb/N0 (dB)   SER      BER approx
15           9            1.8e-2   4.5e-3
16           10           7.0e-3   1.8e-3
18           12           5.5e-4   1.4e-4
20           14           1.2e-5   3.0e-6
25           19           2.7e-15  6.7e-16

As a side note, the confidence interval of simulation results using 2,000,000 symbols at SERs below approximately 10^-5 can start to be quite significant. As an illustration, the following graph shows the SER of 16-QAM in blue, with the expected 95% confidence interval of a 2,000,000 symbols simulation in red:

Rayleigh block fading channel

Once performance calibration has been established for the AWGN channel, we can get back to the Rayleigh block fading channel used in the posted code.

Assuming perfect channel state information estimation at the receiver and if there were no noise, it is possible to scale the received signal back exactly onto the original transmitted symbols using the transformation:

y = y1/h;

When noise is present, this transformation unfortunately also scales the noise. Fortunately, the noise remains white and Gaussian such that the basic derivation of AWGN channel equations can be reused with some work. Over independent packets, the statistical distribution of the scaling abs(h) follows a Rayleigh distribution (with parameter sigma²=1/2). Thus to get the average effect of this scaling on SER, it is possible to compute the weighted sum (where the weight is the probability density function of the Rayleigh distribution) of the effects over the range of possible scaling values using the integral:

This can be done numerically with MATLAB using:

function SER = AwgnSer(EsN0)
  SER = 1.5*erfc(sqrt(0.1*EsN0));
end
function f = WeightedAwgnSer(x)
  weight = 2*x.*exp(-x.*conj(x));
  f = weight*AwgnSer(EsN0*x.*conj(x));
end
function SER = BlockRayleighFadingSer(EsN0)
  for ii=1:length(EsN0)
    SER(ii) = quad(inline('WeightedAwgnSer(EsN0(ii),s)','s'), 0, inf);
  end
end

A similar derivation can be obtained for the BER:

function BER = AwgnBer(EsN0)
  x = sqrt(0.1*EsN0);
  q1 = 0.5*erfc(x);
  q3 = 0.5*erfc(3*x);
  q5 = 0.5*erfc(5*x);
  BER = (12*q1+8*q3-4*q5 - q1*(q1+q3-2*q5)+(q3-q5)*q5)/16;
end
function f = WeightedAwgnBer(x)
  weight = 2*x.*exp(-x.*conj(x));
  f = weight*AwgnBer(EsN0*x.*conj(x));
end
function SER = BlockRayleighFadingBer(EsN0)
  for ii=1:length(EsN0)
    SER(ii) = quad(inline('WeightedAwgnBer(EsN0(ii),s)','s'), 0, inf);
  end
end

Note that I've use an exact formula for the BER since the weighted average tends to be affected by low signal-to-noise ratio where the approximation is not very good. It does not make a huge difference on the curve (~0.3dB at Eb/N0=10dB) but it is not something I want to worry about when calibrating performance curves.

This yields the following performance curves:

Other considerations

Decoding performance can be affected by a number of other factors which are beyond the scope of this answer. The following thus only briefly touches on a few common ones and links to external references which may be used for additional information.

The decoder in the posted code uses explicit knowledge of the fading effect (as evidenced from the line y=y1/h;). It is generally not the case, and fading must first be estimated. The estimation of the channel effect at the receiver is beyond the scope of this answer, but generally imperfect estimation result in some performance loss. Performance curves of perfect knowledge are often used as a practical benchmark against which to compare performance under imperfect channel estimation.

Channel coding is often done to improve system performance. Common benchmarks used for coded modulation over the AWGN channel are:

• Shannon's channel capacity
• Union upper bound (eg. these lectures notes), and multiple other bounds found in research literature
• Uncoded modulation performance (which we derived here)

Similarly for coded modulation over flat block Rayleigh fading channel, the following benchmarks are commonly used:

• Outage probability (see section 5.4.1 of this book)
• Uncoded modulation performance (which we derived here)