computebersoftdecisionquantized.m

Computes BER v EbNo curve for convolutional encoding / soft decision Viterbi decoding scheme assum

function Pb=ComputeBerSoftDecisionQuantized(EbNodBvals, Q, Bd, coderate)

% Pb=ComputeBerSoftDecisionQuantized(EbNodBvals, Q, Bd, coderate)
% Q denotes the number of quantization levels at the input of the
% demodulator
% Choosing Q=2 yelds hard decision VA decoding
% computes BER for quantized VA decoding
% based on the information-bit WEF
% This function also plots the BER v Eb/No curve
% Uniform quantization is assumed to derive the DMC channel model
% 
% See also EstimateBitWEF ComputeBerSoftDecisionQuantized

% (c) Dr B. Gremont, 2007

I=find(Bd>0);
dfree=I(1);
% Q=4; % number fo quantization levels
% EbNodBvals=0:0.5:10;
Pb=zeros(size(EbNodBvals));

for k=1:length(EbNodBvals),
    EbNodB=EbNodBvals(k);
    M=2; %BPSK
    EbNo=10.^(EbNodB./10);
    EsNo=EbNo.*coderate.*log2(M);
    Es=1;
    No=0.5.*Es./EsNo;
    sigma=sqrt(No./2);
    
    if Q>2,
        Thresholds=zeros(1,Q-1);
        MeanValues=zeros(1,Q);
        for j=1:Q-1, % find the thresholds
            Thresholds(j)=1./(2.*(Q-1))+(j-1)./(Q-1);
        end
        MeanValues=(0:Q-1).*1./(Q-1);
        % evalaute transition prob. matrix of DMC
        Pij=zeros(2,Q);
        for j=1:Q, % q-ary output
            for i=1:2, % binary input
                if i==1, % bit zero
                    Mean=-1;
                elseif i==2 % bit one
                    Mean=+1;
                end
                if j==1, % bottom threshold
                    Pij(i,j)=1-q((Thresholds(1)-Mean)./sigma);
                elseif j==Q,
                    Pij(i,j)=q((Thresholds(end)-Mean)./sigma);
                else
                    Pij(i,j)=q((Thresholds(j-1)-Mean)./sigma)-q((Thresholds(j)-Mean)./sigma);
                end
            end
        end
        % Compute Bhattacharyya Parameter D0
        Do=0;
        for j=1:Q,
            Do=Do+sqrt(Pij(1,j).*Pij(2,j));
        end
        DoPoly=zeros(size(Bd));
        for j=1:length(Bd),
            DoPoly(j)=(Do).^(j-1);
        end
        Pb(k)=sum(DoPoly.*Bd);
    elseif Q==2,
        I=find(Bd>0);
        P2=zeros(size(Bd));
        p=q(sqrt(2.*coderate.*EbNo));
        for count=1:length(I),
            d=I(count);
            if rem(d,2)==1, % if d odd
                for k2=(d+1)./2:1:d,
                    P2(I(count))=P2(I(count))+nchoosek(d,k2).*p.^k2.*(1-p).^(d-k2); % see Proakis p 490
                end
            else % if d even
                for k2=(d./2+1):d,
                    P2(I(count))=P2(I(count))+nchoosek(d,k2).*p.^k.*(1-p).^(d-k2);
                end
                P2(I(count))=P2(I(count))+0.5.*nchoosek(d,d./2).*p.^(d./2).*(1-p).^(d./2);
            end
        end
        Pb(k)=sum(Bd.*P2);
    end
end
% semilogy(EbNodBvals,Pb,'r')