block_turbo_decode.m

这是一个很好的关于turbo乘积码的matlab代码

function  ber = block_turbo_decode;
clear all;
% (64,57,4)*(64,57,4) BCH code  
n1 = 64;
k1 = 57;
delta1 = 4;

n2 = 64;
k2 = 57;
delta2 = 4;

x = round(rand(1, k1*k2));    % info. bits
% place k1*k2 information bits in a matrix 
for i = 1:k1
    for j = 1:k2
        infor(i,j) = x(k2*(i-1)+j);
    end
end
% coding the k1 row using code epslon2
% g2(x)=x^6 + x + 1
% define D2=[D0 D1 D2 D3 D4 D5]
D2 = zeros(1,6);
for i = 1:k1
    for j = 1:k2
        D2(1,1:6) = [infor(i,j) xor(infor(i,j),D2(1,1)) D2(1,2:5)];
    end
    for k = 1:6
        r2(i,k) = D2(1,7-k);
    end
    
    c_temp2(i,:) =[infor(i,:) r2(i,:)];
    temp2 = sum(c_temp2(i,:));
    temp2 = mod(temp2,2);
    c2(i,:) = [c_temp2(i,:) temp2];
end
% coding the n2 columns using code epslon1
% g1(x)=x^6 + x + 1
D1 = zeros(1,6);
for j = 1:n2
    for i = 1:k1
        D1(1,1:6) = [c2(i,j) xor(c2(i,j),D1(1,1)) D1(1,2:5)];
    end
    
    for k = 1:6
        r1(k,j) = D1(1,7-k);
    end
    c_temp1(:,j) =[c2(:,j); r1(:,j)];
    temp1 = sum(c_temp1(:,j));
    temp1 = mod(temp1,2);
    c(:,j) = [c_temp1(:,j); temp1];
end
% bpsk 
x = 2*c - ones(n1,n2);
% awgn channel
EbN0db = 0.5;
rate = 0.793;
en = 10^(EbN0db/10);
sigma = 1/sqrt(2*rate*en);
r = x + sigma*randn(n1,n2);

%for k =1:8
    w(:,:) = zeros(n1,n2);
    wrot(:,:) = zeros(n2,n1);
%end
    
for i=1:n1
    R(i,:) = r(i,:);
    if i > k1
        R(i,:) = zeros(1,n2);
    end
end

belta = [0 0.5 0.7 0.9 1 1 1 1];
sclar = [0.2 0.3 0.5 0.7 0.9 1 1 1];
for m =1:8
    % 行译码
    w(1:k1,:) = decode(k1,n2,belta(m),R(1:k1,:),m);
    R(:,:) = r(:,:) + w(:,:)*sclar(m);
    
    % 列译码
    wrot(:,:) = decode(n2,n1,belta(m),R(:,:)',m);
    w(:,:) = wrot(:,:)';
    R(:,:) = r(:,:) + w(:,:)*sclar(m);
end

for i = 1:k1
    for j = 1:k2
        if R(i,j) > 0
            xhat(i,j) = 1;
        else
            xhat(i,j) = 0;
        end
    end
end

error = 0;
for i = 1:k1
    for j = 1:k2
        if xhat(i,j) ~= infor(i,j)
            error = error + 1;
        end
    end
end
ber = error/(k1*k2);