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!