# Why does my simulation of the binary symmetric channel disagree with theory?

Could someone please answer me why the simulated performance (using matlab) of a Hamming code over a binary symmetric channel does not fit exacly with the theoretical? Here is my code:

clear all

n=15;
k=11;
I=eye(k);

%define the parity matrix for the hamming code
P=[1 1 0 0;1 0 1 0;0 1 1 0;1 0 0 1;0 1 0 1;0 0 1 1;1 1 1 0;1 1 0 1;1 0 1 1;0 1 1 1;1 1 1 1];

%The generator matrix combination of parity and identity matrix
G=[I,P];

%Define parity check matrix
HI=eye(n-k);
H=[P',HI];

%verify that rows of H represent the dual space for the rows of G
%(verification = 0)
verification=rem(G*H',2);

UoB=1043162;

%convert UoB number to binary
ci=dec2bin(UoB);
ci=ci-('0');

%take the 11 rightmost digits
cw=ci(:,10:20);
cw=rem(cw*G,2);

%monte carlo simulation index
mc=10000;
index=1;
Pe_matrix=zeros(21,mc);
Pe_vector=zeros(1,21);


In the procedure demonstrated bellow random error patterns error=rand(1,15)<p; are generated for 21 different crossover probabilities.Consequently these patterns are subtracted from the original codeword c=xor(cw,error); in order to find the errored codeword (in any place where the error pattern has an 1 the codeword will flip its bit. After that, the probability of error for the current error patern is calculated Pe=(sum(xor(cw,c)))./k; by summing up the ones returned (the places where the errored codeword and the original differ) from the subtraction of the errored codeword and the original one, and by dividing them by the total length of the codeword k=numel(c);.The Pe_matrix stores all the Pe values for each of the cross over probabilities and each of the mc iterations.

for i=1:mc

for p=0:0.01:.2
error=rand(1,15)<p;
c=xor(cw,error);
k=numel(c);
Pe=(sum(xor(cw,c)))./k;
Pe_matrix(index,i)=Pe;
index=index+1;
end
index=1;

end


In order to obtain smoother results Monte Carlo simulations were performed. In each of the Monte Carlo iterations,each row of the Pe_matrix has the Pe values for the different crossover probabilities.Now by summing up the elements of each row and by dividing them with the total number of the Monte Carlo simulations we acquire a strong estimate of the BER for the Binary symmetric channel.

for j=1:21
Pe_vector(1,j)=sum(Pe_matrix(j,:))/mc;
end

%In the code demonstrated bellow the theoretical performance of the channel
%is calculated in order to compare it with the simulated.

for pe=0:0.01:.2
index=index+1;
Pblkerror=1-(1-pe)^15-15*pe*(1-pe)^14;
Pd(1,index)=(3/15)*Pblkerror;
end

%finally by creating a legend one can compare the theoretical with the
%simulated performance by the depicted BER graphs.
semilogy(Pe_vector,'r-o')
hold on;
semilogy(Pd,'k-*')
hold on;
legend('simulation','theoretical performance');
ylabel('Probability of error,Pe');
xlabel('cross over probability');
title('BsC simulation');
grid on;

• Perhaps if you could say a little about how your simulated results don't match the theoretical ones, someone might be able to help you instead of trying to debug your program. Is the simulated error probability larger than theoretical result? smaller? larger or smaller depending on SNR or channel error probability? – Dilip Sarwate Feb 11 '12 at 22:33
• The plot in y axis has the error probability and in x the crossover of the bsc.The theoretica curve starts from approximately 10^-3 and reaches just above 10^-1.The simulated curve instead of starting from 10^-3 starts from 10^-2 and follows the same direction with the theoretical while it is a bit shhifted upwards. – Rizias Feb 12 '12 at 16:16
• I read MATLAB very poorly, but are you comparing word error probability Pblkerror=1-(1-pe)^15-15*pe*(1-pe)^14 to average bit error probability Pe=(sum(xor(cw,c)))./k or using Pd(1,index)=(3/15)*Pblkerror as the theoretical bit error probability? This last seems to effectively assume that when a block error occurs, exactly $3$ of the $15$ channel bits are in error, but what you really need is the probability of data bit error, not channel bit error. – Dilip Sarwate Feb 12 '12 at 16:30
• I edited your code to make it a little more readable and took longer comments outside. If you disagree with this edit, by all means please roll it back. – Phonon Feb 12 '12 at 19:41
• pls i want to plot BSC, BEC and entropy function using matlab. Can someone help me with the codes? – Yusuf Ahmed Jan 8 '17 at 19:09