map.m

卷积码的最大后验概率译码算法的一个验证matlab算法

clc;clear
t=cputime;
disp('islem yapiliyor. Lutfen bekleyiniz...')
son=4;
for snr=0:son
    hata(snr+1)=0;
    snr
    carpan=[100 100 1000 1000 10000 100000];
    for j=1:carpan(snr+1)        
        n=100;
        u=randint(1,n-2);   %最后两位为0
        u(n)=0;
        kafes=poly2trellis(3,[5 7]);%此生成多项式与书上是反的，因此后面有些是不同的。
        v=convenc(u,kafes);
        for i=1:2*n
            s(i)=2*v(i)-1;
        end
        r = awgn(s,snr,'measured');
        %Asagidaki x degerleri BPSK mod黮asyona g鰎e bulunmustur. 
        %每个状态输出两个符号（以BPSK调制方式）,状态s1:00, s2:10, s3:11， s4:01
        x111=-1;x112=-1;x121=1;x122=1; %从1状态到1状态的输出为00，从1状态到2状态的输出为11   
        x231=-1;x232=1;x241=1;x242=-1; %从2状态到3状态的输出为01，从2状态到4状态的输出为10 
        x311=1;x312=1;x321=-1;x322=-1;
        x431=1;x432=-1;x441=-1;x442=1;
        %gamalarin hesabi
        k=0;
       Eb=2;N0=Eb/(10^(snr/10));sigma=N0/2;
       for i=1:2:2*n
            k=k+1;
            gama11(k)=exp((-(r(i)-x111)^2-(r(i+1)-x112)^2)/2/sigma);  %计算输入r(i),r(i+1),可能是某个输出（x,x)的似然函数
            gama12(k)=exp((-(r(i)-x121)^2-(r(i+1)-x122)^2)/2/sigma);
            gama23(k)=exp((-(r(i)-x231)^2-(r(i+1)-x232)^2)/2/sigma);
            gama24(k)=exp((-(r(i)-x241)^2-(r(i+1)-x242)^2)/2/sigma);
            gama31(k)=exp((-(r(i)-x311)^2-(r(i+1)-x312)^2)/2/sigma);
            gama32(k)=exp((-(r(i)-x321)^2-(r(i+1)-x322)^2)/2/sigma);
            gama43(k)=exp((-(r(i)-x431)^2-(r(i+1)-x432)^2)/2/sigma);
            gama44(k)=exp((-(r(i)-x441)^2-(r(i+1)-x442)^2)/2/sigma);
        end
        %alfalarin hesabi:  %
        %初始状态为s00，所以只有alfa1，alfa2有值,alfa1(0)=1,其它alfa（0）=0
        alfa1(1)=gama11(1);alfa2(1)=gama12(1);alfa3(1)=0;alfa4(1)=0;
        for i=2:n
            alfa1(i)=alfa1(i-1)*gama11(i)+alfa3(i-1)*gama31(i); %后续状态转移图，从状态1出发
            alfa2(i)=alfa1(i-1)*gama12(i)+alfa3(i-1)*gama32(i);
            alfa3(i)=alfa2(i-1)*gama23(i)+alfa4(i-1)*gama43(i);
            alfa4(i)=alfa2(i-1)*gama24(i)+alfa4(i-1)*gama44(i);
        end
        %Betalarin hesabi %最后的状态为s00，只有状态1，3有值,beta1(n)=1,betax(n)=0
        beta1(n-1)=gama11(n);beta3(n-1)=gama31(n);beta2(n-1)=0;beta4(n-1)=0;
        for i=n-2:-1:1
            beta1(i)=beta1(i+1)*gama11(i+1)+beta2(i+1)*gama12(i+1);%前续状态转移图，进入状态1
            beta2(i)=beta3(i+1)*gama23(i+1)+beta4(i+1)*gama24(i+1);
            beta3(i)=beta1(i+1)*gama31(i+1)+beta2(i+1)*gama32(i+1);
            beta4(i)=beta3(i+1)*gama43(i+1)+beta4(i+1)*gama44(i+1);
        end
        %L(ak)larin hesabi  %输入的似然比，分子为输入为1的状态转移，分母为输入为0的状态转移
        payda=gama12(1)*beta2(1); %初始状态为s00，只有下面式子的第一项
                L_a(1)=log((gama11(1)*beta1(1))/(payda));
        for i=2:n-2
            pay=(alfa1(i-1)*gama11(i)*beta1(i)+alfa2(i-1)*gama23(i)*beta3(i)+...
                 alfa3(i-1)*gama31(i)*beta1(i)+alfa4(i-1)*gama43(i)*beta3(i));
            payda=(alfa1(i-1)*gama12(i)*beta2(i)+alfa2(i-1)*gama24(i)*beta4(i)+...
                   alfa3(i-1)*gama32(i)*beta2(i)+alfa4(i-1)*gama44(i)*beta4(i));
            L_a(i)=log(pay/payda);
        end 
        for i=1:n-2
            if L_a(i)>0
                a(i)=0;
            else
                a(i)=1;
            end
        end
        a(n)=0; %怎么最后一个不计算，缺省为0,a(n-1)=?
        hata(snr+1)=hata(snr+1)+sum(abs(u-a));
    end    
    hata(snr+1);
    bit_hata_orani(snr+1)=hata(snr+1)/100/carpan(snr+1)
end
snr=0:son;
semilogy(snr,bit_hata_orani),grid
title('MAP');xlabel('Eb/N0(dB)');ylabel('Bit Error Rate');
zaman=cputime-t

%For further infornation: [e-mail]