clear all; clc;

tic

f=2.4*10^9;
lambda=3*10^8/f;
kappa=2*pi/lambda;

Z0=376.73;              %真空阻抗

K=8;                    %用户数量

%% 发射机面积 坐标
Ls_x=sqrt(0.25);                           %发射机板面宽度
Ls_y=sqrt(0.25);
delta_x=33;                         %x轴量化
delta_y=33;                         %y轴量化
s_center= [0 0 0]';                 %发射机中心

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

dS=Ls_x*Ls_y/delta_x/delta_y;       %发射机量化单元面积
NumdS=delta_x*delta_y;              %发射机量化单元数

%% 离散MIMO
Mx = max(floor(Ls_x*2/lambda),1);
My = max(floor(Ls_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-Ls_x/2;
        S_m(2,mx + Mx*(my - 1)) = my*lambda/2-Ls_y/2;
    end
end

%% 傅里叶基展开
delta_u=3;                         %u轴基展开量化
delta_v=3;                         %v轴基展开量化
NumF=delta_u*delta_v;

f_base=zeros(NumF,NumdS);           %傅里叶基生成
for u=1:delta_u
    for v=1:delta_v
        for n=1:NumdS
            f_base(u+delta_u*(v-1),n)=1/sqrt(Ls_x*Ls_y)*exp(1j*2*pi*(u-(delta_u/2))*(Vs(1,n)-Ls_x/2)/Ls_x)*exp(1j*2*pi*(v-(delta_v/2))*(Vs(2,n)-Ls_y/2)/Ls_y);
        end
    end
end

%% 接收机面积 坐标
Vr=lambda^2/4/pi;
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

%% 信道（格林函数）生成
clear k u v a b m n
%% 容量计算
sigma_2_dB = [-40:2.5:-10];
P_T=10^(-4)*ones(length(sigma_2_dB),1);
sigma_2=5.6*10.^(sigma_2_dB/10)/kappa^2/Z0^2;    %V2/m2

length_L=length(sigma_2);

Iteration=20;

sum_rate_H_MIMO_Upper=zeros(length_L,Iteration);
sum_rate_H_MIMO_freespace=zeros(length_L,Iteration);
sum_rate_H_MIMO_wavenumber=zeros(length_L,Iteration);
sum_rate_ditial = zeros(length_L,Iteration);
sum_rate_MF = zeros(length_L,Iteration);

SNR=zeros(length_L,1);

for a=1:Iteration
    parfor m=1:length_L
        
        sprintf('功率，第%d轮',m)

        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(m))/sqrt(temp);

        user_position=zeros(3,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];

        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(m)/(lambda/2)^2)*v/norm(v);

        sum_rate_ditial(m,a) = Digital_MIMO_precoding(K,M,H_k,Psi_k,v,P_T(m),sigma_2(m),lambda);

        G_k = Freespace_G_k_generate(lambda,user_position,Vs,K,NumdS,delta_x,delta_y);

        [~,~,sum_rate_MF(m,a)] = MF_MIMO_precoding(K,NumdS,dS,G_k,Psi_k,theta_k,P_T(m),sigma_2(m));

        SNR(m)=P_T(m)/sigma_2(m);
    end

end
toc
%%
sum_rate_H_MIMO_Upper=mean(sum_rate_H_MIMO_Upper,2);
sum_rate_H_MIMO_freespace=mean(sum_rate_H_MIMO_freespace,2);
sum_rate_H_MIMO_wavenumber=mean(sum_rate_H_MIMO_wavenumber,2);
sum_rate_ditial=mean(sum_rate_ditial,2);
sum_rate_MF = mean(sum_rate_MF,2);

P_T=SNR*5.6*10.^(-3)/kappa^2/Z0^2;  %等效增大功率为

P_T=P_T*10^6;


figure;
semilogx(P_T,sum_rate_ditial,'-rp',P_T,sum_rate_MF,'-bs','LineWidth',1.5);
hold on;
grid on;
box on;
xlabel('Maximum transmit power $P_{\rm T}$ (${\rm mA}^2$)','Interpreter','latex');
ylabel('Sum-rate (bps/Hz)','Interpreter','latex');
legend('Fully-digital MIMO','MF','Interpreter','latex','FontSize',12);
set(gca,'FontName','Times','FontSize',14);
ylim([0 20])

save('Digital_MIMO_power.mat','P_T','sum_rate_MF','sum_rate_ditial');


