clear all; clc;

tic
warning off;

%% 基本参数
f=2.4*10^9;         % 频率
lambda=3*10^8/f;
kappa0=2*pi/lambda;

Z0=376.73;                              %真空阻抗
sigma_2=5.6*10^(-3)/kappa0^2/Z0^2;      %噪声功率 V2/m2
P_T=10^(-4);                            %电流密度功率P_T

K = 8;                                  %用户数量

user_position=zeros(3,K);               %用户坐标r_k

L=30;                                   %用户中心的Z轴坐标
user_position(:,1)=[1 1 L];
user_position(:,2)=[1 -1 L];
user_position(:,3)=[-1 -1 L];
user_position(:,4)=[-1 1 L];
user_position(:,5)=[5 5 L];
user_position(:,6)=[5 -5 L];
user_position(:,7)=[-5 -5 L];
user_position(:,8)=[-5 5 L];

%% 发射机面积 坐标
A_T = [0:0.1:1];                            %发射机面积
A_T(1) = 0.001;

length_L=length(A_T);                       %不同设置的个数

Iteration=10;                                %重复试验次数

sum_rate_PDM=zeros(length_L,Iteration);         %所提方法和速率
sum_rate_MF=zeros(length_L,Iteration);         %所提方法和速率
sum_rate_ditial=zeros(length_L,Iteration);      %离散MIMO和速率
sum_rate_single=zeros(length_L,Iteration);      %单极化离散MIMO和速率

L_x=zeros(length_L,1);                      %X轴方向长度设置
L_y=zeros(length_L,1);                      %Y轴方向长度设置

%% 主函数
for q=1:Iteration
    parfor m=1:length_L

        sprintf('面积，第%d轮',m)

        %% 空间域采样设置
        L_x=sqrt(A_T(m));                    %发射机X轴方向长度
        L_y=sqrt(A_T(m));

        delta_x=33;                             %X轴量化次数
        delta_y=33;                             %y轴量化次数

        S_T=zeros(3,delta_x*delta_y);           %发射机表面坐标
        for a=1:delta_x
            for b=1:delta_y
                S_T(1,a+delta_x*(b-1))=L_x/delta_x*a-L_x/2;
                S_T(2,a+delta_x*(b-1))=L_y/delta_y*b-L_y/2;
                S_T(3,a+delta_x*(b-1))=0;
            end
        end

        dS=L_x*L_y/delta_x/delta_y;             %发射机量化单元面积

        NumdS=delta_x*delta_y;                  %发射机量化单元数
        %% 波数域采样设置
        N_x=9;                                  %波数域x轴基展开量化
        N_y=9;                                  %波数域y轴基展开量化
        N_F=N_x*N_y;                            %波束域展开总个数
        %% 傅里叶基函数生成
        f_base=zeros(N_F,NumdS);                %傅里叶基函数
        for n_x=1:N_x
            for n_y=1:N_y
                for n=1:NumdS
                    f_base(n_x+N_x*(n_y-1),n)=1/sqrt(L_x*L_y)          *exp(1j*2*pi*(n_x-(N_x/2))*(S_T(1,n)-L_x/2)/L_x);
                    f_base(n_x+N_x*(n_y-1),n)=f_base(n_x+N_x*(n_y-1),n)*exp(1j*2*pi*(n_y-(N_y/2))*(S_T(2,n)-L_y/2)/L_y);
                end
            end
        end
        %% 发射机表面电流初始化
        theta_k = zeros(3,K,NumdS);
        temp=0;
        for n=1:NumdS
            for k=1:K
                theta_k(:,k,n)=(2*rand(3,1)-1+2j*(rand(3,1)-0.5))/sqrt(6);
                temp=temp+norm(theta_k(:,k,n))^2*dS;
            end
        end
        theta_k=theta_k*sqrt(P_T)/sqrt(temp);
        %% 接收机Combiner初始化
        Psi_k=zeros(3,K);   %combiner初始化
        for k=1:K
            Psi_k(:,k)=(2*rand(3,1)-1+2j*(rand(3,1)-0.5))/sqrt(6);
            Psi_k(:,k)=Psi_k(:,k)/norm(Psi_k(:,k));
        end
        %% 离散MIMO参数
        Mx = max(floor(L_x*2/lambda),1);
        My = max(floor(L_y*2/lambda),1);
        
        M = Mx*My;                          %离散天线总数
        
        S_m = zeros(3,M);
        for mx = 1:Mx
            for my = 1:My
                S_m(1,mx + Mx*(my - 1)) = mx*lambda/2-L_x/2;
                S_m(2,mx + Mx*(my - 1)) = my*lambda/2-L_y/2;
            end
        end
        
        H_k = MIMO_channel_generate(lambda,user_position,S_m,K,M);
        
        v = 2*rand(3*M*K,1)-1+1j*2*(rand(3*M*K,1)-0.5);           %离散MIMO预编码向量初始化
        v = sqrt(P_T/(lambda/2)^2)*v/norm(v);
        
        sum_rate_ditial(m,q) = Digital_MIMO_precoding(K,M,H_k,Psi_k,v,P_T,sigma_2,lambda);
        
        %% 容量优化
        G_k = Freespace_G_k_generate(lambda,user_position,S_T,K,NumdS,delta_x,delta_y);
        [~,~,sum_rate_MF(m,q)] = MF_MIMO_precoding(K,NumdS,dS,G_k,Psi_k,theta_k,P_T,sigma_2);

    end
end
sum_rate_single = mean(sum_rate_single,2);
sum_rate_MF = mean(sum_rate_MF,2);
sum_rate_ditial = mean(sum_rate_ditial,2);
sum_rate_PDM = mean(sum_rate_PDM,2);

figure;
hold on;
grid on;
box on;
plot(A_T,sum_rate_PDM,'-r','LineWidth',1.5);
plot(A_T,sum_rate_MF,'-g','LineWidth',1.5);
plot(A_T,sum_rate_ditial,'-b','LineWidth',1.5);

xlabel('Aperture area $A_{\rm T}$ (${\rm m}^2$)','Interpreter','latex');
ylabel('Capacity (bps/Hz)','Interpreter','latex');
legend('CAP-MIMO','MF','Fully-digital MIMO','Interpreter','latex','FontSize',12);
set(gca,'FontName','Times','FontSize',14);

save('Digital_MIMO_area.mat','A_T','sum_rate_MF','sum_rate_ditial');

toc


