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); %% 產(chǎn)生kasami序列
% M=2^(length(fbconnection)/2); %%系統(tǒng)的最大用戶數(shù)
M=6;%%系統(tǒng)的最大用戶數(shù)
kasamiseq1=kasamiseq(1:M,:);
clear kasamiseq
kasamiseq=kasamiseq1;
ind_ka=find(kasamiseq==0);
kasamiseq(ind_ka)=-1;% 極性變換
N=2^length(fbconnection)-1; % N為kasami序列長度，即為擴頻倍數(shù)（擴頻增益）
% kasamiseq=(1/sqrt(N))*kasamiseq; %% 歸一化

SNR=3:9;
n1=length(SNR);
a1=SNR(1);
a2=SNR(n1);
sgma=1;
R=kasamiseq*kasamiseq';
iR=inv(R);
alfa=[0.1 0.3 1];      %%用戶活躍概率
total_T=1000;
tic

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


% for alfa_i=1:length(alfa)
alfa_i=3;
        for TT=1:total_T
                identity=rand(1,M);
                identity(1)=1; %%% reference user
                for i=2:M
                    if (identity(i)>alfa(alfa_i))
                        identity(i)=0;
                    else
                        identity(i)=1;
                    end
                end
                identity_yuan(TT,:)=identity;
        end
%%%%%%%%%%%%%%%%經(jīng)過驗證是對的%%%%%%%%%%%%%%%%%%%



    
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);
        for i=2:M
            A(i,i)=sqrt(Ec);
        end
        clear i
        ss=[];
        data_yuan=[];
       for TT_i=1:total_T
              ind_iden=find(identity_yuan(TT_i,:)==1);  %%the identity of user
%               identity(ind_iden)=ind_iden;
              num_fact=length(ind_iden);  %%實際用戶數(shù)
              num_data=1; %%每個離散時間間隔內(nèi)傳輸?shù)亩M制數(shù)
              T=1;   %%觀測時間
              data_number=num_data*T;
              data=randsrc(num_fact,data_number);
              data_yuan=[data_yuan data(1)];%%%傳送的數(shù)據(jù)都是隨機的，但是用戶數(shù)已經(jīng)定下來了%%%%
    
              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;  %%%經(jīng)匹配濾波
    %%%  ML檢測
         for j=1:2^M
             base2=dec2base(j-1,2,M);
            temp=zeros(1,length(base2));
             for jj=1:length(base2)
                 temp(jj)=base2(jj)-48;
             end
            ind_temp=find(temp==0);
            temp(ind_temp)=-1;
            b(:,j)=temp';
        end
        clear base2 j temp ind_temp jj;
        %%%%%這里是認為每個用戶都存在的經(jīng)典ML估計法%%%%%%%%%
        %%%%%可以寫成這樣的原因是用戶存在與否每個可能性都考慮到了，且用戶傳輸?shù)臄?shù)據(jù)是+1還是-1也考慮進去了

         for i=1:total_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(alfa_i,cyc_snr)=ber;%%%是一個三維的矩陣%%%%%
%         BER(cyc_snr)=ber;
%         data_ml(:,:,cyc_snr,alfa_i)=est_b1;
   
        clear i j f ind_f 

  %%%   FRST和ML、FRST和MAL兩者結(jié)合檢測 
  
        for cyc_T=1:total_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));%%%y中的列上有些元素是空的%%%%%
                preP=2^(-M+length(ind_Xb))*alfa(alfa_i)^(M-length(ind_Xb))*(1-alfa(alfa_i))^(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(alfa_i,cyc_snr)=ber1;
        BER2(alfa_i,cyc_snr)=ber2;
        BER1(cyc_snr)=ber1;
        BER2(cyc_snr)=ber2;
        data_frstml(:,:,cyc_snr,alfa_i)=final_Xfmlb;
        data_frstmal(:,:,cyc_snr,alfa_i)=final_Xfmalb;
        
       cyc_snr=cyc_snr+1
end
% clear identity_yuan
% end
toc 



figure(1)
semilogy(SNR,BER(1,:),'-kd',SNR,BER(2,:),'-k*',SNR,BER(3,:),'-ks',SNR,BER1(1,:),'-rx',SNR,BER1(2,:),'-ro',SNR,BER1(3,:),'-r+',SNR,BER2(1,:),'-g<',SNR,BER2(2,:),'-g>',SNR,BER2(3,:),'-g^')
% title('ML、FRST-ML和FRST-MAP的比較');
% semilogy(SNR,BER,'-kd',SNR,BER1,'-k*',SNR,BER2,'-ks')
xlabel('Eb/No(dB)');
ylabel('BER');
legend('ML alpha=0.1','ML alpha=0.3','ML alpha=1.0','RST-ML alpha=0.1','RST-ML alpha=0.3','RST-ML alpha=1.0','RST-MAP alpha=0.1','RST-MAP alpha=0.3','RST-MAP alpha=1.0')
% legend('ML alfa=0.1','RST-ML alfa=0.1','RST-MAP alfa=0.1')
 grid on