tic
clc;
clear;
T_Antenna=7;    
R_Antenna=8;
code_g=[1 0 1 1;1 1 1 1]   
code_rate=1/2;
[c1,constrain_length]=size(code_g);
data_size=5600;     
source_size=data_size-constrain_length+1;  %有一部分比特用于歸零.
Eb_NO_dB=[-5:0];
a=length(Eb_NO_dB);
source_num=[10,10,10,10,10,10];
puncture_mode=2;
code_rate=1/puncture_mode;
joint_iteration_num=4;
mapping=[1,j,-j,-1] % QPSK 星座映射圖
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mapping=([-3,-3,-3,-3,-1,-1,-1,-1,+3,+3,+3,+3,+1,+1,+1,+1]...
%           +j*[-3,-1,+3,+1,-3,-1,+3,+1,-3,-1,+3,+1,-3,-1,+3,+1])/sqrt(10); % 16 QAM 星座映射圖
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 64 QAM 星座映射圖
% mapping1=repmat([-7,-5,-1,-3,+7,+5,+1,+3]',1,8).*ones(8,8);
% mapping2=repmat([-7,-5,-1,-3,+7,+5,+1,+3]',1,8).*ones(8,8);
% mapping1=mapping1';
% mapping1=mapping1(:)';
% mapping2=mapping2(:)';
% mapping=(mapping1+j*mapping2)/sqrt(42);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mod_num=length(mapping);
MMSE_PIC_BER=zeros(a,joint_iteration_num);
ZF_SIC_BER=zeros(a,joint_iteration_num);
MMSE_SIC_BER=zeros(a,joint_iteration_num);
for i1=1:length(Eb_NO_dB)   
    fprintf('@@@@@@@@@@@@@@@@@@@@@@@@@Now the SNR is %fdb.@@@@@@@@@@@@@@@@@@@@@@@@@@\n',Eb_NO_dB(i1));  
    bit_er_snr_rate_MMSE_PIC=zeros(source_num(i1),joint_iteration_num);
    bit_er_snr_rate_ZF_SIC=zeros(source_num(i1),joint_iteration_num);
    bit_er_snr_rate_MMSE_SIC=zeros(source_num(i1),joint_iteration_num);
    for i2=1:source_num(i1)
        i2
        inf_source=randint(1,source_size); 
        [coder] = myencoder(inf_source,code_g); 
        channel_Alpha=interleaver(2*data_size);                              % 加信道交織
        coder=coder(1,channel_Alpha);
        coder_lay=reshape(coder,T_Antenna,length(coder)/T_Antenna);
        tran_sig=modulate(coder_lay,mapping);   
        [rec_sig,channel_matrix,noise_var]=channel(tran_sig,Eb_NO_dB(i1),R_Antenna,mod_num,code_rate);
        aprior_det=zeros(1,2*data_size); 
        extrin_det=zeros(1,2*data_size);
        for i3=1:joint_iteration_num
            err_num=0;
            aprior_det_MMSE_PIC=extrin_det_MMSE_PIC(1,channel_Alpha);
            aprior_det_ZF_SIC=extrin_det_ZF_SIC(1,channel_Alpha);
            aprior_det_MMSE_SIC=extrin_det_MMSE_SIC(1,channel_Alpha);
            if i3==1
               LLR_MMSE_PIC=MMSE(rec_sig,channel_matrix,noise_var,mapping);
           else
               LLR_MMSE_PIC=MMSE_PIC(rec_sig,channel_matrix,noise_var,mapping,aprior_det);
           end;
           LLR_ZF_SIC=ZF_SIC(rec_sig,channel_matrix,mapping,noise_var,aprior_det);
           LLR_ZF_SIC=LLR_ZF_SIC-aprior_det;
           LLR_MMSE_SIC=MMSE_SIC(rec_sig,channel_matrix,mapping,noise_var,aprior_det);
           LLR_MMSE_SIC=LLR_MMSE_SIC-aprior_det;
           LLR_MMSE_PIC(1,channel_Alpha)=LLR_MMSE_PIC;
           LLR_ZF_SIC(1,channel_Alpha)=LLR_ZF_SIC;
           LLR_MMSE_SIC(1,channel_Alpha)=LLR_MMSE_SIC;      
           [LLR_MMSE_PIC_Con,LLR_MMSE_PIC_I]=decode(LLR_MMSE_PIC,code_g);
           [LLR_ZF_SIC_Con,LLR_ZF_SIC_I]=decode(LLR_ZF_SIC,code_g);
           [LLR_MMSE_SIC_Con,LLR_MMSE_SIC_I]=decode(LLR_MMSE_SIC,code_g);
           extrin_det_MMSE_PIC=LLR_MMSE_PIC_Con;
           extrin_det_ZF_SIC=LLR_ZF_SIC_Con-LLR_ZF_SIC;
           extrin_det_MMSE_SIC=LLR_MMSE_SIC_Con-LLR_MMSE_SIC;
           %%%%%%%%%%%%%%%%%% 算各次跌代的誤碼率 %%%%%%%%%%%%%%%%%%%%%%%%
           %%%%%%% MMSE_PIC %%%%%%%%
           LLR_MMSE_PIC_I=LLR_MMSE_PIC_I(1:source_size);
           inf_dec_MMSE_PIC=(sign(LLR_MMSE_PIC_I)+1)/2;
           err_num_MMSE_PIC=sum(inf_dec_MMSE_PIC~=inf_source);
           bit_er_snr_rate_MMSE_PIC(i2,i3)=err_num_MMSE_PIC/(source_size);
           %%%%%%% ZF_SIC %%%%%%%%%
           LLR_ZF_SIC_I=LLR_ZF_SIC_I(1:source_size);
           inf_dec_ZF_SIC=(sign(LLR_ZF_SIC_I)+1)/2;
           err_num_ZF_SIC=sum(inf_dec_ZF_SIC~=inf_source);
           bit_er_snr_rate_ZF_SIC(i2,i3)=err_num_ZF_SIC/(source_size);
           %%%%%%% MMSE_SIC %%%%%%%%
           LLR_MMSE_SIC_I=LLR_MMSE_SIC_I(1:source_size);
           inf_dec_MMSE_SIC=(sign(LLR_MMSE_SIC_I)+1)/2;
           err_num_MMSE_SIC=sum(inf_dec_MMSE_SIC~=inf_source);
           bit_er_snr_rate_MMSE_SIC(i2,i3)=err_num_MMSE_SIC/(source_size);
       end;
    end;
    MMSE_PIC_BER(i1,:)=sum(bit_er_snr_rate_MMSE_PIC,1)/source_num(i1);
    ZF_SIC_BER(i1,:)=sum(bit_er_snr_rate_ZF_SIC,1)/source_num(i1);
    MMSE_SIC_BER(i1,:)=sum(bit_er_snr_rate_MMSE_SIC,1)/source_num(i1);
end;
semilogy(Eb_NO_dB,MMSE_PIC_BER(:,1),'b-*');
hold on; 
semilogy(Eb_NO_dB,ZF_SIC_BER(:,1),'b.-');
semilogy(Eb_NO_dB,MMSE_SIC_BER(:,1),'b-*');
semilogy(Eb_NO_dB,MMSE_PIC_BER(:,2),'r-+');
semilogy(Eb_NO_dB,ZF_SIC_BER(:,2),'r-+');
semilogy(Eb_NO_dB,MMSE_SIC_BER(:,2),'r-+');
semilogy(Eb_NO_dB,MMSE_PIC_BER(:,3),'g-o');
semilogy(Eb_NO_dB,ZF_SIC_BER(:,3),'g-o');
semilogy(Eb_NO_dB,MMSE_SIC_BER(:,3),'g-o');
semilogy(Eb_NO_dB,MMSE_PIC_BER(:,3),'m-p');
semilogy(Eb_NO_dB,ZF_SIC_BER(:,3),'m-p');
semilogy(Eb_NO_dB,MMSE_SIC_BER(:,3),'m-p');
legend('MMSE_PIC 1 iteration','ZF_SIC 1 iteration','MMSE_SIC 1 iteration'...
       'MMSE_PIC 2 iteration','ZF_SIC 2 iteration','MMSE_SIC 2 iteration'...
       'MMSE_PIC 3 iteration','ZF_SIC 3 iteration','MMSE_SIC 3 iteration'...
       'MMSE_PIC 4 iteration','ZF_SIC 4 iteration','MMSE_SIC 4 iteration')
toc