clc
clear
radium=600;                              %小區(qū)半徑
fc=5e9;                                  %中心載波頻率
n1=2;                                    %衰落系數(shù)(LOS)
n2=4;                                    %衰落系數(shù)(NLOS)
sample_number=1;
B=120000;                                %每個(gè)子信道的帶寬大小
subcarrier_number=128;                    %子載波數(shù)目
N0=3.9811*1e-21;                          %功率譜密度-174dbm
noise=N0*B;                              %每個(gè)子信道的噪聲
ber=1.0e-3;
gamma=-log2(5*ber)/(1.6*log2(exp(1)));
number=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%基站和中繼端的坐標(biāo)
BS_x=0;
BS_y=0;
Relay_x=[-300,0];%各個(gè)中繼的橫坐標(biāo);2個(gè)中繼在一個(gè)扇區(qū)。小區(qū)分為3個(gè)扇區(qū) 
Relay_y=[100*sqrt(3),200*sqrt(3)];%各個(gè)中繼的縱坐標(biāo)
SNR=156;
p=10^(SNR/10)*noise;
user_number=[5,10,15,20,25];

total_number=100;


% total_number=10000;
R1_ave=zeros(total_number,length(user_number));
R2_ave=zeros(total_number,length(user_number));
R3_ave=zeros(total_number,length(user_number));
R4_ave=zeros(total_number,length(user_number));
R5_ave=zeros(total_number,length(user_number));
for cycle_number=1:1:total_number
    cycle_number
    for for_user_number=1:1:length(user_number) 
        user_x=zeros(1,user_number(for_user_number));
        user_y=zeros(1,user_number(for_user_number));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %確定每個(gè)用戶的坐標(biāo)
        for i=1:1:user_number(1,for_user_number)
            user_x(i)=-600+900*rand(1);
            if abs(user_x(i))<=300
                if user_x(i)<=0
                    user_y(i)=600*rand(1); %應(yīng)該為300*sqrt(3)*rand(1)
                else
                    user_y(i)=rand(1)*(sqrt(3)*300-sqrt(3)*user_x(i))+sqrt(3)*user_x(i);;
                end
            else
                ceil1=sqrt(3)*(600-abs(user_x(i)));
                user_y(i)=ceil1*rand(1);
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        distance_BS_RN=400*ones(1,2);
        distance_RN_UT1=zeros(1,user_number(for_user_number));
        distance_BS_UT=zeros(1,user_number(for_user_number));
        for k=1:1:number %%%%  k表征撒點(diǎn)數(shù)
            [I1,Q1]=chtry(2,sample_number,20);
            [I2,Q2]=chtry(user_number(for_user_number),sample_number,20);
            [I3,Q3]=chtry(user_number(for_user_number),sample_number,20);
            [I4,Q4]=chtry(user_number(for_user_number),sample_number,20);
            for i=1:1:2  %%%%2個(gè)中繼，BS-RS的信道矩陣
                relay=I1(i,:)+sqrt(-1)*Q1(i,:);
                middle1=fft(relay,8*subcarrier_number);
                for jj=1:1:subcarrier_number
                    H_BS_RN(k,i,jj)=sum(middle1(8*(jj-1)+1:8*jj))/8;
                end
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           中繼1到用戶的信道矩陣
            for i=1:1:user_number(for_user_number)
                distance_RN_UT1(i)=sqrt((user_x(i)-Relay_x(1))^2+(user_y(i)-Relay_y(1))^2);
                middle2_1=fft(I2(i,:)+sqrt(-1)*Q2(i,:),8*subcarrier_number);
                for jj=1:1:subcarrier_number
                    H_RN_UT1(k,i,jj)=sum(middle2_1(8*(jj-1)+1:8*jj))/8;
                end
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           中繼2到用戶的信道矩陣
            for i=1:1:user_number(for_user_number)
                distance_RN_UT1(i)=sqrt((user_x(i)-Relay_x(2))^2+(user_y(i)-Relay_y(2))^2);
                middle2_1=fft(I2(i,:)+sqrt(-1)*Q2(i,:),8*subcarrier_number);
                for jj=1:1:subcarrier_number
                    H_RN_UT2(k,i,jj)=sum(middle2_1(8*(jj-1)+1:8*jj))/8;
                end
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             基站到用戶的信道矩陣
            for i=1:1:user_number(for_user_number)
                distance_BS_UT(i)=sqrt((user_x(i))^2+(user_y(i))^2);
                middle3=fft(I3(i,:)+sqrt(-1)*Q3(i,:),8*subcarrier_number);
                for jj=1:1:subcarrier_number
                    H_BS_UT(k,i,jj)=sum(middle3(8*(jj-1)+1:8*jj))/8;
                end
            end
        end
        R1=zeros(number,user_number(for_user_number));    
        R2=zeros(number,user_number(for_user_number));   
        R3=zeros(number,user_number(for_user_number));   
        R4=zeros(number,user_number(for_user_number));
        R5=zeros(number,user_number(for_user_number)); 
        for t=1:1:number
            for i=1:1:2         
                for j=1:1:subcarrier_number        
                    H_BS_RN(t,i,j)=sqrt(10^(-4.6)*1/(200*sqrt(3))^2)*H_BS_RN(t,i,j);  %加入大尺度       
                end
            end
            for i=1:1:user_number(for_user_number)      
                for j=1:1:subcarrier_number      
                    H_RN_UT1(t,i,j)=sqrt(10^(-4.6)*1/(distance_RN_UT1(i))^4)*H_RN_UT1(t,i,j);    
                end
            end
            for i=1:1:user_number(for_user_number)      
                for j=1:1:subcarrier_number      
                    H_RN_UT2(t,i,j)=sqrt(10^(-4.6)*1/(distance_RN_UT1(i))^4)*H_RN_UT2(t,i,j); %%error!   
                end
            end
            for i=1:1:user_number(for_user_number)     
                for j=1:1:subcarrier_number          
                    H_BS_UT(t,i,j)=sqrt(10^(-4.6)*1/(distance_BS_UT(i))^4)*H_BS_UT(t,i,j);                       
                end
            end 
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %proposed algorithm
            R_least=12800*ones(1,user_number(for_user_number));
            occupy=zeros(1,subcarrier_number);
            C=zeros(user_number(for_user_number),subcarrier_number,3);
            for i=1:1:user_number(for_user_number)
                for j=1:1:subcarrier_number
                    for k=1:1:3
                        if k==3
                            C(i,j,k)=B*log2(1+p*(abs(H_BS_UT(t,i,j)))^2/(noise*gamma));
                        elseif k==1
                            C(i,j,k)=1/2*min(B*log2(1+p*(abs(H_BS_RN(t,1,j)))^2/(noise*gamma)),B*log2(1+p*(abs(H_RN_UT1(t,i,j)))^2/(noise*gamma)));
                        else
                            C(i,j,k)=1/2*min(B*log2(1+p*(abs(H_BS_RN(t,1,j)))^2/(noise*gamma)),B*log2(1+p*(abs(H_RN_UT2(t,i,j)))^2/(noise*gamma)));%%error!!
                        end
                    end
                end
            end
            C_all=C;
            C_max=zeros(user_number(for_user_number),subcarrier_number);
            for i=1:1:user_number(for_user_number)
                for j=1:1:subcarrier_number
                    C_max(i,j)=max(C(i,j,:));
                end
            end
            user_index=ones(1,user_number(for_user_number));
            subcarrier_index=ones(1,subcarrier_number);
            C_max1=C_max;
            while sum(subcarrier_index)>0
                if sum(user_index)>0
                    [max1,index]=max(max(C_max1'));
                    [max2,index2]=max(max(C_max1));
                    subcarrier_index(index2)=0;
                    R1(t,index)=R1(t,index)+max1;
                    occupy(1,index2)=index;
                    if R1(t,index)>=R_least(1,index)
                        C_max1(index,:)=0;
                        user_index(1,index)=0;
                    end
                    C_max1(:,index2)=0;
                else
                    for j=1:1:subcarrier_number
                        if subcarrier_index(1,j)==1
                            [max1,index]=max(C_max(:,j));
                            R1(t,index)=R1(t,index)+max1;
                            subcarrier_index(1,j)=0;
                            occupy(1,j)=index;
                        end
                    end
                end
            end
            for i=1:1:subcarrier_number