clear all, close all, clc
%% parameter
SNR_dB = [-20:5:20];
SNR_linear = 10.^(SNR_dB/10.);
N_iter = 100;
N = 64; % number of beams (transmit antennas)
N_RF = 4; % number of RF chians
K = 6; % number of users
L = 3; % number of paths per user
lamada = 1; % wavelength
d = lamada/2;
P = 30;%K; % total transmitted power
Pconv = 0.9; % energy conversion efficiency
Pmin = 0.1; % energy harvesting threshold
Imax = 20; % iteration times of power allocation
bit = 4; % quantification bit of phase shifter
Rtemp0 = zeros(K,1);
temp_it0 = [];
temp_it1 = [];
temp_it2 = [];
temp_it3 = [];
temp_it4 = [];
temp_it5 = [];
temp_it6 = [];
ite = 1:Imax;
%%
for i_snr = 1:length(SNR_dB)
SNR = SNR_linear(i_snr);
sigma2 = P/SNR;
noise = sigma2;
temp0 = zeros(1,Imax);
temp1 = zeros(1,Imax);
temp2 = zeros(1,Imax);
temp3 = zeros(1,Imax);
temp4 = zeros(1,Imax);
temp5 = zeros(1,Imax);
temp6 = zeros(1,Imax);
iter = 1;
while iter <= N_iter
H = generate_channel(N,K,L,lamada,d); % generate the channel
%%% for initialization
F = H*inv(H'*H);
beta = sqrt(P/trace(F'*F));
F = F*beta;
H_eq = H'*F;
for k=1:K
sum_power = sum(abs(H_eq(k,:)).^2);
psl = 1-Pmin/Pconv/(sum_power+sigma2);
sum_inf = sum(abs(H_eq(k,:)).^2)-abs(H_eq(k,k))^2;
Rtemp0(k) = log2(1+abs(H_eq(k,k))^2/(sum_inf+sigma2+noise/psl));
end
Rmin = min(Rtemp0)/10;
%%% full digital - ZF
p0 = sum(abs(F).^2,1);
p0 = p0'*(P/sum(p0));
F = F./repmat(sqrt(sum(abs(F).^2,1)),N,1);
[SE0,power0,feasbile_point0] = PA0(H',F,p0,K,sigma2,noise,P,Pconv,Pmin,Rmin,Imax);
%%% Two-stage hybrid precoding: rf_num = K
[AP1,DP1,p0] = two_stage(H,K,bit,P);
[SE1,power1,feasbile_point1] = PA0(H'*AP1,DP1,p0,K,sigma2,noise,P,Pconv,Pmin,Rmin,Imax);
[AP1,DP1,p0] = two_stage_sub(H,K,bit,P);
[SE6,power6,feasbile_point6] = PA0(H'*AP1,DP1,p0,K,sigma2,noise,P,Pconv,Pmin,Rmin,Imax);
%%% HP-NOMA
if K > N_RF
[HA_full,AP_full,HA_sub,AP_sub,setu] = A_precoder(H,N,N_RF,K,bit); % analog precoding
setf_full = user_grouping(HA_full,N_RF,K,setu); % user grouping
DP_full = D_precoder(HA_full,AP_full,K,N_RF,setf_full); % digital precoding
[SE2,power2,feasbile_point2] = PA(HA_full,DP_full,setf_full,K,N_RF,sigma2,noise,P,Pconv,Pmin,Rmin,Imax);
setf_sub = user_grouping(HA_sub,N_RF,K,setu); % user grouping
DP_sub = D_precoder(HA_sub,AP_sub,K,N_RF,setf_sub); % digital precoding
[SE4,power4,feasbile_point4] = PA(HA_sub,DP_sub,setf_sub,K,N_RF,sigma2,noise,P,Pconv,Pmin,Rmin,Imax);
else
feasbile_point2 = 1;
feasbile_point4 = 1;
end
%%% HP-OMA
if K > N_RF
[SE3,power3,feasbile_point3] = PA00(HA_full,DP_full,setf_full,K,N_RF,sigma2,noise,P,Pconv,Pmin,Rmin,Imax);
[SE5,power5,feasbile_point5] = PA00(HA_sub,DP_sub,setf_sub,K,N_RF,sigma2,noise,P,Pconv,Pmin,Rmin,Imax);
else
feasbile_point3 = 1;
feasbile_point5 = 1;
end
if (feasbile_point0 == 0) || (feasbile_point1 == 0) || (feasbile_point2 == 0) || (feasbile_point3 == 0) || (feasbile_point4 == 0) || (feasbile_point5 == 0) || (feasbile_point6 == 0)
continue;
end
temp0 = temp0 + SE0;
temp1 = temp1 + SE1;
temp6 = temp6 + SE6;
if K <= N_RF
temp2 = temp1;
temp3 = temp1;
temp4 = temp6;
temp5 = temp6;
else
temp2 = temp2 + SE2;
temp3 = temp3 + SE3;
temp4 = temp4 + SE4;
temp5 = temp5 + SE5;
end
iter = iter+1;
end
C0(i_snr) = temp0(end)/N_iter;
C1(i_snr) = temp1(end)/N_iter;
C2(i_snr) = temp2(end)/N_iter;
C3(i_snr) = temp3(end)/N_iter;
C4(i_snr) = temp4(end)/N_iter;
C5(i_snr) = temp5(end)/N_iter;
C6(i_snr) = temp6(end)/N_iter;
C00(i_snr) = C0(i_snr)/(sum(power0)+N*300+200)*10^3;
C11(i_snr) = C1(i_snr)/(sum(power1)+K*300+N*K*40+200)*10^3;
C22(i_snr) = C2(i_snr)/(sum(power2)+N_RF*300+N*N_RF*40+200)*10^3;
C33(i_snr) = C3(i_snr)/(sum(power3)+N_RF*300+N*N_RF*40+200)*10^3;
C44(i_snr) = C4(i_snr)/(sum(power4)+N_RF*300+N*40+200)*10^3;
C55(i_snr) = C5(i_snr)/(sum(power5)+N_RF*300+N*40+200)*10^3;
C66(i_snr) = C6(i_snr)/(sum(power6)+K*300+N*40+200)*10^3;
temp_it0 = [temp_it0;temp0];
temp_it1 = [temp_it1;temp1];
temp_it2 = [temp_it2;temp2];
temp_it3 = [temp_it3;temp3];
temp_it4 = [temp_it4;temp4];
temp_it5 = [temp_it5;temp5];
temp_it6 = [temp_it6;temp6];
end
CI0 = temp_it0(1,:)/N_iter;
CI1 = temp_it1(1,:)/N_iter;
CI2 = temp_it2(1,:)/N_iter;
CI3 = temp_it3(1,:)/N_iter;
CI4 = temp_it4(1,:)/N_iter;
CI5 = temp_it5(1,:)/N_iter;
CI6 = temp_it6(1,:)/N_iter;
%%
figure;
plot(SNR_dB,C0,'k-*','Linewidth',1.5);
hold on
plot(SNR_dB,C2,'b-*','Linewidth',1.5);
hold on
plot(SNR_dB,C3,'g-*','Linewidth',1.5);
hold on
plot(SNR_dB,C4,'y-*','Linewidth',1.5);
hold on
plot(SNR_dB,C5,'m-*','Linewidth',1.5);
grid on
xlabel('SNR (dB)');
ylabel('Spectral efficiency (bps/Hz)');
legend('SWIPT-Fully Digital ZF Precoding', 'SWIPT-Fully-Connected HP-NOMA', 'SWIPT-Fully-Connected HP-OMA', 'SWIPT-Sub-Connected HP-NOMA', 'SWIPT-Sub-Connected HP-OMA');
figure;
plot(SNR_dB,C00,'k-','Linewidth',1.5);
hold on
plot(SNR_dB,C22,'b-','Linewidth',1.5);
hold on
plot(SNR_dB,C33,'g-*','Linewidth',1.5);
hold on
plot(SNR_dB,C44,'y-','Linewidth',1.5);
hold on
plot(SNR_dB,C55,'m-*','Linewidth',1.5);
grid on
xlabel('SNR (dB)');
ylabel('Energy efficiency (bps/Hz/W)');
legend('SWIPT-Fully Digital ZF Precoding', 'SWIPT-Fully-Connected HP-NOMA', 'SWIPT-Fully-Connected HP-OMA', 'SWIPT-Sub-Connected HP-NOMA', 'SWIPT-Sub-Connected HP-OMA');
% figure;
% plot(ite,CI2,'b-','Linewidth',1.5);
% hold on
% plot(ite,CI4,'y-','Linewidth',1.5);
% grid on
% xlabel('Iteration');
% ylabel('Spectral efficiency (bps/Hz)');
% legend('SWIPT-Fully-Connected HP-NOMA', 'SWIPT-Sub-Connected HP-NOMA');
% filename = ['SR_', num2str(N), '_', num2str(N_RF), '_', num2str(K), '.mat'];
% save(filename, 'C0', 'C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'temp_it2', 'temp_it4', 'SNR_dB', 'ite');
