hammingcodes3.m

A simple m-file to compute the experiemntal BER of coded BPSK

% file demonstrating basic block code encoding/decoding

%===============================================
% Hamning code specifications
%===============================================
% create generator and parity check matrix
P=[1 1 0
    0 1 1 
    1 0 1
    1 1 1];
[A,B]=size(P);
Ik=eye(A);

G=[Ik P]; % Generator Matrix (systematic form)

I=eye(B);
H=[P' I]; % Parity Check Matrix (systematic form)

[k,N]=size(G);
coderate=k./N;

% construct the standard array (decoding table)
m=N-k;
table = zeros(2^m,N);
emptyRows = [2:2^N]';
% Each row contains a single-error pattern

errorpattern=zeros(2.^m,N);

for i=0:2.^m-1,
    if i==0,
        e=zeros(1,N);
    elseif i==1
        e=zeros(1,N);
        e(1)=1;
    else
        elast=e(end);
        e(2:end)=e(1:end-1);
        e(1)=elast;
    end
    % e is one of the 2^m possible single error pattern
    errorpattern(i+1,:)=e;
end
% errorpattern   % uncomment to see contents of errorpattern

%compute syndrome table for each single error pattern
syndrome_table=rem(errorpattern*H',2);    
%===============================================

NoofCodewords=15000; % increse this to extend the BER curve to lower snrs
EbNodBVals=2:1:18;

BER=zeros(size(EsNodBVals));


for SNR=1:length(EbNodBVals),
    NoErrorTotal=0;
    NoBitsTotal=0;
    for Noword=1:NoofCodewords,
        %===============================================
        % encoding/decoding
        %===============================================
        % create random input message block (1*k)
        m=floor(rand(1,k).*2);
        % encoding c=mG (mod 2)
        c=rem(m*G,2); % codeword

        % simulate BPSK modulator
        signalsent=2.*c-1;

        % noise
        EbNo=10.^(EbNodBVals(SNR)./10);
        EsNo=EbNo.*coderate;
        Es=1;
        No=Es./EsNo;
        sigma=sqrt(No)./2;
        n=sigma.*randn(size(signalsent));
        sr=signalsent+n; %awgn, sr= received signal
        % hard quantize demodulator
        r=zeros(size(sr));
        r(sr>=0)=1;

        %compute syndrome
        s=rem(r*H',2);
        % find syndrome in syndrome table
        for i=1:2.^m,
            s_table=syndrome_table(i,:);
            if s_table==s, % syndrome found
                e=errorpattern(i);
            end
        end
        % compute corrected codeword
        c_estimate=rem(r+e,2);

        % computes No of bit errors
        No_errors=length(find(c~=c_estimate));
        No_bits=length(c_estimate);
        NoErrorTotal=NoErrorTotal+No_errors;
        NoBitsTotal= NoBitsTotal+No_bits;
    end
    BER(SNR)=NoErrorTotal./NoBitsTotal;
end
% plot results
if ishold, hold,end
semilogy(EbNodBVals, BER,'r')
figure(gcf)
hold
EbNo=10.^(EbNodBVals./10);
%EbNo=EsNo;

p=1/2*erfc(sqrt(2.*coderate.*EbNo)/sqrt(2));

%p=q(sqrt(2.*coderate.*EbNo));
Pb=p-p.*(1-p).^(N-1);
semilogy(EbNodBVals,Pb,'g:')
axis([2 12 10.^(-10) 1])
grid