%%%%-----此段程式碼完成PNC下行鏈路LDPC譯碼，中繼端資料假設是理想接收無差錯---%%%%
clc
clear
code_length=24;  %碼長
code_rate=1/2;   %位元速率
max_iter=20;     %譯碼最大迭代次數
SNR=3;
%%----LDPC檢驗矩陣和生成矩陣的構造----%%%
H=genH(code_length,code_rate);  %生成校驗矩陣
G=genG(H);       %獲得生成矩陣
G_gf1=gf(G,2);   %轉換到有限域
G_gf2=gf(G,2);   %轉換到有限域
dim=size(H);
d=dim(2)-dim(1);
%%---隨機產生A端和B端原始資料，並進行編碼
%------A端-----%
msg1=randi([0 1],1,d);
s_gf1=gf(msg1,2);          %轉換到有限域
code_bit1=s_gf1*G_gf1;    %%編碼
x1=code_bit1;
code_bit1=double(code_bit1.x);
%---- B端-----%
msg2=randi([0 3],1,d);
s_gf2=gf(msg2,2);          %轉換到有限域
code_bit2=s_gf2*G_gf2;    %%編碼
x2=code_bit2;
code_bit2=double(code_bit2.x);
%-----假設在中繼端無差錯接收-----%
yr=x1+x2;
yr=double(yr.x);

%%下行鏈路分兩路進行BPSK調製傳輸
yr=de2bi(yr);
yr1=yr(:,1);
yr2=yr(:,2);

y_bp1=-2*yr1+1;
y_bp2=-2*yr2+1;
%--加噪聲--%
yaw1=awgn(y_bp1,SNR,'measured');
yaw2=awgn(y_bp2,SNR,'measured');
%--接收端軟解碼，獲取碼字的概率資訊--%
sigma=sqrt(10^(-0.1*SNR));
%第一路
f10=1./(1+exp(-2*yaw1/sigma^2));
f11=1-f10;
%第二路
f20=1./(1+exp(-2*yaw2/sigma^2));
f21=1-f20;
%---合併兩路，計算GF（4）上碼字的概率
p0=f10.*f20;
p1=f11.*f20;
p2=f10.*f21;
p3=f11.*f21;
%%
%%
%B端求解A端的資訊
len=length(yr);
Pa0=zeros(1,len);
for loop=1:len
    if code_bit2(loop)==0
        Pa0(loop)=p0(loop);
    elseif code_bit2(loop)==1
        Pa0(loop)=p1(loop);
    elseif code_bit2(loop)==2
        Pa0(loop)=p2(loop);
    elseif code_bit2(loop)==3
        Pa0(loop)=p3(loop);
    end
end

Pa1=1-Pa0;
[yhat_a] = decode(H,max_iter,Pa1,Pa0);
ya=yhat_a';
err1=nnz(code_bit1-ya);
%%
%%
%A端解B端的資訊
Pb0=zeros(1,len);
Pb1=zeros(1,len);
Pb2=zeros(1,len);
Pb3=zeros(1,len);
for loop=1:len
    if code_bit1(loop)==0;
        Pb0(loop)=p0(loop);
        Pb1(loop)=p1(loop);
        Pb2(loop)=p2(loop);
        Pb3(loop)=p3(loop);
    elseif code_bit1(loop)==1
        Pb0(loop)=p1(loop);
        Pb1(loop)=p0(loop);
        Pb2(loop)=p3(loop);
        Pb3(loop)=p2(loop);
    end
end

Pb10=Pb0+Pb1;
Pb11=Pb2+Pb3;

Pb20=Pb0+Pb2;
Pb21=Pb1+Pb3;

yhat1= decode(H,max_iter,Pb11,Pb10);
yhat2= decode(H,max_iter,Pb21,Pb20);

yhat_b=[yhat2,yhat1];  %這裡位元發生顛倒，需注意
yb=(bi2de(yhat_b))';    %%二進位制位元流轉換為4元
err2=nnz(code_bit2-yb);