clear all;
clc;       
fbconnection=[1 0 0 0 0 1];
% fbconnection=[0 0 1 0 0 0 0 0 0 1]; %%本原多項式
kasamiseq=kasami_seq(fbconnection); %% 產生kasami序列
% M=2^(length(fbconnection)/2); %%系統的最大用戶數
M=6;%%系統的最大用戶數
kasamiseq1=kasamiseq(1:M,:);
clear kasamiseq
kasamiseq=kasamiseq1;
ind_ka=find(kasamiseq==0);
kasamiseq(ind_ka)=-1;% 極性變換
N=2^length(fbconnection)-1; % N為kasami序列長度，即為擴頻倍數（擴頻增益）
% kasamiseq=(1/sqrt(N))*kasamiseq; %% 歸一化
sgma=1;
snr_exp=6;
snr_exp=10^(snr_exp/10);
Eb_exp=2*sgma^2*snr_exp;
Ec_exp=Eb_exp/N;

R=kasamiseq*kasamiseq';
iR=inv(R);
    
T=8000;   %%觀測時間
tic
for T_i=1:T
    identity=randsrc(1,M);
    
   identity_yuan(T_i,:)=identity;
end

b=confi_ml(M);    
[X,Xb]=gen_X(M);
[row_Xb,col_Xb]=size(Xb);

SNR=0:12;
n1=length(SNR);
a1=SNR(1);
a2=SNR(n1);   
    cyc_snr=1;
    for snrindb=a1:a2
        snr=10^(snrindb/10);
        Eb=2*sgma^2*snr;
        Ec=Eb/N;
        A=zeros(M,M);
        A(1,1)=sqrt(Ec_exp);
        for i=2:M
            A(i,i)=sqrt(Ec);
        end
        clear i
        ss=[];
        data_yuan=[];
       for TT_i=1:T
               %%the identity of user
%               identity(ind_iden)=ind_iden;
              num_fact=M;  %%實際用戶數
              num_data=1; %%每個離散時間間隔內傳輸的二進制數
%               T=1;   %%觀測時間
              data_number=num_data;
              data=randsrc(num_fact,data_number);
              data_yuan=[data_yuan data(1)];
    
              s=0;
              for i=1:num_fact
                     temp=A(ind_iden(i),ind_iden(i))*kasamiseq(ind_iden(i),:)'*data(i,:);
                     s=s+temp;
              end
             ss=[ss,s];
       end
        [ss_row,ss_col]=size(ss);
        r=ss+sgma*randn(ss_row,ss_col);
        clear temp i
        y=kasamiseq*r;  %%%經匹配濾波
        
    %%%  ML檢測

         for i=1:T
             for j=1:2^M
                 f(j)=(y(:,i)-R*A*b(:,j))'*iR*(y(:,i)-R*A*b(:,j));
             end
             [value_f,ind_f]=min(f);
             est_b1(:,i)=b(:,ind_f);
        end
        [number,ber]=symerr(est_b1(1,:),data_yuan);   
        BER(cyc_snr)=ber;
        data_ml(:,:,cyc_snr)=est_b1;
   
        clear i j f ind_f est_b1

  %%%   FRST和ML、FRST和MAL兩者結合檢測 
  
        for cyc_T=1:T
            for i=1:col_Xb
                ind_Xb=find(Xb(:,i)==0);
                fml(i)=(y(:,cyc_T)-R*A*Xb(:,i))'*iR*(y(:,cyc_T)-R*A*Xb(:,i));
                preP=2^(-M+length(ind_Xb))*alfa^(M-length(ind_Xb))*(1-alfa)^(length(ind_Xb));
                fmal(i)=preP*exp((-1/(2*sgma^2))*((y(:,cyc_T)-R*A*Xb(:,i))'*iR*(y(:,cyc_T)-R*A*Xb(:,i))));
            end
            [value_fml,ind_fml]=min(fml);
            [value_fmal,ind_fmal]=max(fmal);
            final_Xfmlb(:,cyc_T)=Xb(:,ind_fml); 
            final_Xfmalb(:,cyc_T)=Xb(:,ind_fmal);         
            clear i fml fmal  
        end
        [number,ber1]=symerr(final_Xfmlb(1,:),data_yuan); 
        [number,ber2]=symerr(final_Xfmalb(1,:),data_yuan); 
        BER1(cyc_snr)=ber1;
        BER2(cyc_snr)=ber2;
        data_frstml(:,:,cyc_snr)=final_Xfmlb;
        data_frstmal(:,:,cyc_snr)=final_Xfmalb;
        cyc_snr=cyc_snr+1;
    end

toc 



figure(1)
semilogy(SNR-6,BER,'-d',SNR-6,BER1,'-*',SNR-6,BER2,'-s')
% title('ML、FRST-ML和FRST-MAP的比較');
axis([-6 6 10^(-4) 1])
xlabel('interference user minus the reference user of Eb/No(dB)');
ylabel('BER');
legend('ML','RST-ML ','RST-MAP')
 grid on