From 555bf0130c25ec9e801dd66c8a534291f3d60b0f Mon Sep 17 00:00:00 2001 From: FenghaoZhu <61408487+FenghaoZhu@users.noreply.github.com> Date: Thu, 4 Jan 2024 16:01:04 +0800 Subject: [PATCH] Add files via upload --- R_WMMSE.m | 93 ++++++++++++++++++++++++++++ Test.m | 164 +++++++++++++++++++++++++++++++++++++++++++++++++ Test_R_WMMSE.m | 78 +++++++++++++++++++++++ Test_WMMSE.m | 54 ++++++++++++++++ WMMSE.m | 71 +++++++++++++++++++++ find_U.m | 14 +++++ find_U1.m | 17 +++++ find_V1.m | 47 ++++++++++++++ find_W.m | 6 ++ find_W1.m | 9 +++ find_X.m | 68 ++++++++++++++++++++ sum_rate.m | 14 +++++ sum_rate1.m | 17 +++++ 13 files changed, 652 insertions(+) create mode 100644 R_WMMSE.m create mode 100644 Test.m create mode 100644 Test_R_WMMSE.m create mode 100644 Test_WMMSE.m create mode 100644 WMMSE.m create mode 100644 find_U.m create mode 100644 find_U1.m create mode 100644 find_V1.m create mode 100644 find_W.m create mode 100644 find_W1.m create mode 100644 find_X.m create mode 100644 sum_rate.m create mode 100644 sum_rate1.m diff --git a/R_WMMSE.m b/R_WMMSE.m new file mode 100644 index 0000000..7331f8d --- /dev/null +++ b/R_WMMSE.m @@ -0,0 +1,93 @@ +clc;clear; +K = 1; % 基站个数，此版本固定为1 +T = 512; % 发射天线个数 +R = 4; % 接收天线个数 +epsilon = 1e-3; % 收敛条件 +sigma2 = 1; % 噪声功率 +snr = 10; % 信噪比 +P = db2pow(snr)*sigma2; % 发射功率 +I = 16; % 用户个数 +alpha1 = ones(I,K); % 权重系数，都假设相同 +d = 4; % 假设每个用户都有d条路独立的数据流 +max_iter = 100; +rate = []; % 初始化一个空速率记录rate +time = []; % 初始化一个空运行时间记录rate +tic; %开始计时 +begin_time = toc; % 标记开始时间，结束时间减去开始时间就是使用的时间 +time = [time 0]; +% 初始化信道向量 +H = zeros(R,T,I); % 信道系数 用户数量*每个用户天线数量*基站天线数量 +% 把所有用户的信道矩阵拼接以第一个维度拼接起来 +H_full = zeros(I*R,T); +for i=1:I + H(: , :, i)=sqrt(1/2)*(randn(R,T)+1i*randn(R,T)); % 圆复高斯信道 + H_full((i-1)*R+1:i*R,:)=H(:,:, i); +end + +% 计算H_hat = H_full*H_full' +H_hat = H_full*H_full'; + + + + +% 初始化W和U矩阵 +U =randn(R,d,I) + 1j*randn(R,d,I); +W = zeros(d,d,I); +for i=1:I + W(:,:,i)=eye(d,d); +end + +% 用低维度的坍塌矩阵X代替波束赋形矩阵，V = H'*X, 随机初始化X +X = zeros(R*I,d,I); % 每个基站的坍塌波束赋形矩阵 +V = zeros(T,d,I); % 每个基站的波束赋形矩阵 +for i=1:I + x = sqrt(1/2)*(randn(R*I,d)+1i*randn(R*I,d)); % 随机初始化低维度的坍塌矩阵X + v = H_full'*x; + X(:,:, i) = sqrt(P/(I*trace(H_hat*x*x')))*x; + V(:,:, i) = H_full'*X(:,:,i); +end + +% 求初始化发射波束V后求系统和速率 +rate_old = sum_rate(H,V,sigma2,R,I,alpha1); +rate = [rate rate_old]; + +iter1 = 1; +while(1) + U = find_U(H_hat,X,sigma2, P, R,I,d); % R-WMMSE公式24 + W = find_W(U,H_hat,X, R , I,d); % R-WMMSE公式25 + X = find_X(alpha1,H_hat,sigma2,U,W,T, R, I,d ,P); % R-WMMSE公式26 + + % 计算放缩系数beta + beta = 0; + + for i = 1:I + beta = beta + trace(H_hat*X(:,:,i)*(X(:,:,i)')); + end + + beta = P / beta; + + % 计算真正的波束赋形向量 + for i = 1:I + V(:,:,i) = sqrt(beta) * H_full'*X(:,:,i); + end + + + rate_new = sum_rate(H,V,sigma2,R,I,alpha1); % 计算和速率 + rate = [rate rate_new]; + elapsed_time = toc; + elapsed_time = elapsed_time - begin_time; + time = [time elapsed_time]; + iter1 = iter1 + 1; + if abs(rate_new-rate_old) / rate_old < epsilon || iter1 > max_iter + break; + end + rate_old = rate_new; +end + +% plot(0:iter1-1,rate,'r-o') +% grid on +% xlabel('Iterations') +% ylabel('Sum rate (bits per channel use)') +% set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1) +% title('R-WMMSE, K=1, T=128, R=4, D=4, \epsilon=1e-3','Interpreter','tex') +% title('SISO-IFC, K=3, T=1, R=1, \epsilon=1e-3','Interpreter','tex') \ No newline at end of file diff --git a/Test.m b/Test.m new file mode 100644 index 0000000..34ca48c --- /dev/null +++ b/Test.m @@ -0,0 +1,164 @@ +% 此脚本用来运行WMMSE和R-WMMSE这两个算法的性能差距，指标包括运行时间，迭代次数以及最终速率 +% 暂时此脚本只支持单基站的仿真情景 + +% 定义运行参数 +clc;clear; +rng(1); % 设置随机数种子 +K=1; % 基站个数，目前由于R-WMMSE程序只支持一个基站，故只能固定为1 +T=256; % 基站天线数量 +R=4; % 每个用户天线数量 +epsilon=0.001; % 收敛设定的限制 +sigma2=1; % 噪声功率 +snr=10; % 用户的信噪比 +I=16; % 用户数量 +alpha1=ones(I,K); % 用户权重 +d=4; % 每个用户流数 +max_iter=100; % 最大的迭代次数 +num_sample = 100; % 信道样本数量 + + +% 生成信道，以便于在相同信道的条件下测试R-WMMSE和WMMSE的性能差别 +H_WMMSE = cell(I,K,K); % 生成原版WMMSE所需的信道 +for i=1:I + for k = 1:K + for j=1:K + H_WMMSE{i,k,j}=sqrt(1/2)*(randn(R,T)+1i*randn(R,T)); % 圆复高斯信道 + end + end +end + +H_R_WMMSE = zeros(R,T,I); % 生成R-WMMSE所需的信道 +for i=1:I + H_R_WMMSE(:,:,i) = H_WMMSE{i}; +end + +% 运行对比测试,生成前两张图 +[iter_R_WMMSE, time_R_WMMSE, rate_R_WMMSE] = Test_R_WMMSE(K,T,R,epsilon,sigma2,snr,I,alpha1,d,max_iter); +[iter_WMMSE, time_WMMSE, rate_WMMSE] = Test_WMMSE(K,T,R,epsilon,sigma2,snr,I,alpha1,d,max_iter); + +figure(1); +plot(0:iter_R_WMMSE-1,rate_R_WMMSE, '-sb') +hold on +plot(0:iter_WMMSE-1,rate_WMMSE, '-*r') +grid on +xlabel('Iterations') +ylabel('Sum rate (bits per channel use)') +set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1) +title(['Compare WMMSE with R-WMMSE, K=',num2str(K), ',', 'T=', num2str(T), ',', 'R=', num2str(R), ',','\epsilon=', num2str(epsilon)]) +legend('R-WMMSE','WMMSE') +savefig(['./figs/T=',num2str(T),', I=',num2str(I),', d=',num2str(d),',',num2str(snr),'dB, sumrate.fig']) + +figure(2); +plot(0:iter_R_WMMSE-1,time_R_WMMSE, '-bs') +hold on +plot(0:iter_WMMSE-1,time_WMMSE, 'd-r', 'MarkerFaceColor', 'r') +grid on +xlabel('Iterations') +ylabel('Time (s)') +set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1) +title(['Compare WMMSE with R-WMMSE, K=',num2str(K), ',','T=', num2str(T), ',','R=', num2str(R), ',','\epsilon=', num2str(epsilon)]) +legend('R-WMMSE','WMMSE') +savefig(['./figs/T=',num2str(T),', I=',num2str(I),', d=',num2str(d),',',num2str(snr),'dB, runtime.fig']) + +% 运行对比测试，生成运行时间关于用户数量变化的图 +% begin_user=5; +% end_user=45; +% len=end_user-begin_user+1; +% time1 = zeros(len,1); % WMMSE的不同用户数量的运行时间 +% time2 = zeros(len,1); % R-WMMSE的不同用户数量的运行时间 +% +% bar = waitbar(0,'开始测试'); % waitbar显示进度条 +% for num_user=begin_user:end_user +% for f = 1:num_sample +% [iter_R_WMMSE, time_R_WMMSE, rate_R_WMMSE] = Test_R_WMMSE(K,T,R,epsilon,sigma2,snr,num_user,ones(num_user,1),d,max_iter); +% [iter_WMMSE, time_WMMSE, rate_WMMSE] = Test_WMMSE(K,T,R,epsilon,sigma2,snr,num_user,ones(num_user,1),d,max_iter); +% time1(num_user-begin_user+1)=time1(num_user-begin_user+1)+time_WMMSE(iter_WMMSE); +% time2(num_user-begin_user+1)=time2(num_user-begin_user+1)+time_R_WMMSE(iter_R_WMMSE); +% str=['计算中...',num2str(100*f/100),'%'];% 百分比形式显示处理进程 +% waitbar((num_user-begin_user+1)/len,bar,str) % 更新进度条bar +% end +% time1(num_user-begin_user+1) = time1(num_user-begin_user+1) / (num_sample); +% time2(num_user-begin_user+1) = time2(num_user-begin_user+1) / (num_sample); +% end +% close(bar); % 循环结束关闭进度条 +% +% figure(1); +% plot(begin_user:end_user,time2, '-sb') % R-WMMSE的不同用户数量的运行时间 +% hold on +% plot(begin_user:end_user,time1, '-*r') % WMMSE的不同用户数量的运行时间 +% grid on +% xlabel('Number of users') +% ylabel('Average CPU Time (s)') +% set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1) +% title(['Compare WMMSE with R-WMMSE, K=',num2str(K), ',', 'T=', num2str(T), ',', 'R=', num2str(R), ',','\epsilon=', num2str(epsilon)]) +% legend('R-WMMSE','WMMSE') +% savefig('./figs/T=128, I=12, d=2, 10dB, runtime.fig') + +% 运行对比测试，生成运行时间关于基站天线数量变化的图 +% antenna_num_pool = [32, 64, 128, 256, 512]; +% len = length(antenna_num_pool); +% time1 = zeros(len,1); % WMMSE的不同用户数量的运行时间 +% time2 = zeros(len,1); % R-WMMSE的不同用户数量的运行时间 +% bar = waitbar(0,'开始测试'); % waitbar显示进度条 +% for num_antenna_index=1:len +% for f = 1:num_sample +% [iter_R_WMMSE, time_R_WMMSE, rate_R_WMMSE] = Test_R_WMMSE(K,antenna_num_pool(num_antenna_index),R,epsilon,sigma2,snr,I,ones(I,1),d,max_iter); +% [iter_WMMSE, time_WMMSE, rate_WMMSE] = Test_WMMSE(K,antenna_num_pool(num_antenna_index),R,epsilon,sigma2,snr,I,ones(I,1),d,max_iter); +% time1(num_antenna_index)=time1(num_antenna_index)+time_WMMSE(iter_WMMSE); +% time2(num_antenna_index)=time2(num_antenna_index)+time_R_WMMSE(iter_R_WMMSE); +% str=['计算中...',num2str(100*f/100),'%'];% 百分比形式显示处理进程 +% waitbar((num_antenna_index)/len,bar,str) % 更新进度条bar +% end +% time1(num_antenna_index) = time1(num_antenna_index) / (num_sample); +% time2(num_antenna_index) = time2(num_antenna_index) / (num_sample); +% end +% close(bar); % 循环结束关闭进度条 +% +% figure(1); +% semilogy(antenna_num_pool,time2, '-sb') % R-WMMSE的不同用户数量的运行时间 +% hold on +% semilogy(antenna_num_pool,time1, '-*r') % WMMSE的不同用户数量的运行时间 +% grid on +% xlabel('Number of BS Antennas') +% ylabel('Average CPU Time (s)') +% set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1) +% title(['Compare WMMSE with R-WMMSE, K=',num2str(K), ',', 'T=', num2str(T), ',', 'R=', num2str(R), ',','\epsilon=', num2str(epsilon)]) +% legend('R-WMMSE','WMMSE') +% savefig('./figs/I=16, d=4, 10dB, runtime1.fig') + +% 运行遍历测试，测试所有的可能组合并且保存数据以及图 +% T_pool = [32, 64, 128, 256, 512, 1024]; % 天线数量范围 +% I_num_pool = 5:5:45; % 用户数量范围 +% d_num_pool = 1:8; % 用户数据流数数量范围 +% snr_num_pool = -10:5:30; % 信噪比范围 +% [iter_R_WMMSE, time_R_WMMSE, rate_R_WMMSE] = Test_R_WMMSE(K,T,R,epsilon,sigma2,snr,I,alpha1,d,max_iter); +% [iter_WMMSE, time_WMMSE, rate_WMMSE] = Test_WMMSE(K,T,R,epsilon,sigma2,snr,I,alpha1,d,max_iter); +% +% figure(1); +% plot(0:iter_R_WMMSE-1,rate_R_WMMSE, '-sb') +% hold on +% plot(0:iter_WMMSE-1,rate_WMMSE, '-*r') +% grid on +% xlabel('Iterations') +% ylabel('Sum rate (bits per channel use)') +% set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1) +% title(['Compare WMMSE with R-WMMSE, K=',num2str(K), ',', 'T=', num2str(T), ',', 'R=', num2str(R), ',','\epsilon=', num2str(epsilon)]) +% legend('R-WMMSE','WMMSE') +% savefig(['./figs/T=',num2str(T),', I=',num2str(I),', d=',num2str(d),',',num2str(snr),'dB, sumrate.fig']) +% +% figure(2); +% plot(0:iter_R_WMMSE-1,time_R_WMMSE, '-bs') +% hold on +% plot(0:iter_WMMSE-1,time_WMMSE, 'd-r', 'MarkerFaceColor', 'r') +% grid on +% xlabel('Iterations') +% ylabel('Time (s)') +% set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1) +% title(['Compare WMMSE with R-WMMSE, K=',num2str(K), ',','T=', num2str(T), ',','R=', num2str(R), ',','\epsilon=', num2str(epsilon)]) +% legend('R-WMMSE','WMMSE') +% savefig(['./figs/T=',num2str(T),', I=',num2str(I),', d=',num2str(d),',',num2str(snr),'dB, runtime.fig']) + + + + + diff --git a/Test_R_WMMSE.m b/Test_R_WMMSE.m new file mode 100644 index 0000000..4d2561b --- /dev/null +++ b/Test_R_WMMSE.m @@ -0,0 +1,78 @@ +function [iter1, time, rate] = Test_R_WMMSE(K, T, R, epsilon, sigma2, snr, I, alpha1, d, max_iter) +%TEST_R_WMMSE 是用来测试R-WMMSE性能的函数 +% 输入函数运行参数，返回迭代次数，运行时间与速率信息 +P = db2pow(snr)*sigma2; % 发射功率 +rate = []; % 初始化一个空速率记录rate +time = []; % 初始化一个空运行时间记录rate +tic; %开始计时 +begin_time = toc; % 标记开始时间，结束时间减去开始时间就是使用的时间 +time = [time 0]; +% 初始化信道向量 +H = zeros(R,T,I); % 信道系数 用户数量*每个用户天线数量*基站天线数量 +% 把所有用户的信道矩阵拼接以第一个维度拼接起来 +H_full = zeros(I*R,T); +for i=1:I + H(: , :, i)=sqrt(1/2)*(randn(R,T)+1i*randn(R,T)); % 圆复高斯信道 + H_full((i-1)*R+1:i*R,:)=H(:,:, i); +end + +% 计算H_hat = H_full*H_full' +H_hat = H_full*H_full'; + +% 初始化W和U矩阵 +U =randn(R,d,I) + 1j*randn(R,d,I); +W = zeros(d,d,I); +for i=1:I + W(:,:,i)=eye(d,d); +end + +% 用低维度的坍塌矩阵X代替波束赋形矩阵，V = H'*X, 随机初始化X +X = zeros(R*I,d,I); % 每个基站的坍塌波束赋形矩阵 +V = zeros(T,d,I); % 每个基站的波束赋形矩阵 +for i=1:I + x = sqrt(1/2)*(randn(R*I,d)+1i*randn(R*I,d)); % 随机初始化低维度的坍塌矩阵X + v = H_full'*x; + X(:,:, i) = sqrt(P/(I*trace(H_hat*x*x')))*x; + V(:,:, i) = H_full'*X(:,:,i); +end + +% 求初始化发射波束V后求系统和速率 +rate_old = sum_rate(H,V,sigma2,R,I,alpha1); +rate = [rate rate_old]; + +iter1 = 1; +while(1) + U = find_U(H_hat,X,sigma2, P, R,I,d); % R-WMMSE公式24 + W = find_W(U,H_hat,X, R , I,d); % R-WMMSE公式25 + X = find_X(alpha1,H_hat,sigma2,U,W,T, R, I,d ,P); % R-WMMSE公式26 + + % 计算放缩系数beta + beta = 0; + + for i = 1:I + beta = beta + trace(H_hat*X(:,:,i)*(X(:,:,i)')); + end + + beta = P / beta; + + % 计算真正的波束赋形向量 + for i = 1:I + V(:,:,i) = sqrt(beta) * H_full'*X(:,:,i); + end + + + rate_new = sum_rate(H,V,sigma2,R,I,alpha1); % 计算和速率 + rate = [rate rate_new]; + elapsed_time = toc; + elapsed_time = elapsed_time - begin_time; + time = [time elapsed_time]; + iter1 = iter1 + 1; + if abs(rate_new-rate_old) / rate_old < epsilon || iter1 > max_iter + break; + end + rate_old = rate_new; +end + + +end + diff --git a/Test_WMMSE.m b/Test_WMMSE.m new file mode 100644 index 0000000..032bf02 --- /dev/null +++ b/Test_WMMSE.m @@ -0,0 +1,54 @@ +function [iter1, time, rate] = Test_WMMSE(K, T, R, epsilon, sigma2, snr, I, alpha1, d, max_iter) +%TEST_WMMSE 是用来测试WMMSE性能的函数 +% 输入函数运行参数，返回迭代次数，运行时间与速率信息 +P = db2pow(snr)*sigma2; % 发射功率 +rate = []; % 初始化一个空向量记录rate +time = []; % 初始化一个空运行时间记录rate +tic; %开始计时 +begin_time = toc; % 标记开始时间，结束时间减去开始时间就是使用的时间 +time = [time 0]; +% 初始化信道向量 +H = cell(I,K,K); % 信道系数 +for i=1:I + for k = 1:K + for j=1:K + H{i,k,j}=sqrt(1/2)*(randn(R,T)+1i*randn(R,T)); % 圆复高斯信道 + end + end +end + +U = cell(I,K); +U(:)={zeros(R,d)}; + +% 随机初始化发射波束向量 +V = cell(I,K); % 算法第一行 +for i=1:I + for k=1:K + v = randn(T,d)+1i*randn(T,d); % 随机初始化 + V{i,k}=sqrt(P/I)*v/norm(v,"fro"); + end +end + +% 求初始化发射波束V后求系统和速率 +rate_old = sum_rate1(H,V,sigma2,R,I,K,alpha1); +rate = [rate rate_old]; + + +iter1 = 1; +while(1) + U = find_U1(H,V,sigma2,R,I,K,d); % Tbale I line 4 in p.4435 + W = find_W1(U,H,V,I,K,d); % Tbale I line 5 in p.4435 + V = find_V1(alpha1,H,U,W,T,I,K,P); % Tbale I line 6 in p.4435 + rate_new = sum_rate1(H,V,sigma2,R,I,K,alpha1); + rate = [rate rate_new]; + elapsed_time = toc; + elapsed_time = elapsed_time - begin_time; + time = [time elapsed_time]; + iter1 = iter1 + 1; + if abs(rate_new-rate_old) / rate_old < epsilon || iter1 > max_iter + break; + end + rate_old = rate_new; +end +end + diff --git a/WMMSE.m b/WMMSE.m new file mode 100644 index 0000000..f1e1654 --- /dev/null +++ b/WMMSE.m @@ -0,0 +1,71 @@ +clc;clear; +K = 1; % 基站个数 +T = 128; % 发射天线个数 +R = 4; % 接收天线个数 +epsilon = 1e-3; % 收敛条件 +sigma2 = 1; % 噪声功率 +snr = 30; % 信噪比 +P = db2pow(snr)*sigma2; % 发射功率 + +I = 4; % 每个基站服务的用户个数 +alpha1 = ones(I,K); % 权重系数，都假设相同 + +d = 1; % 假设每个用户有d条路独立的数据流 +max_iter = 100; +tic; +rate = []; % 初始化一个空向量记录rate +time = []; % 初始化一个空运行时间记录rate +tic; %开始计时 +begin_time = toc; % 标记开始时间，结束时间减去开始时间就是使用的时间 +time = [time 0]; +% 初始化信道向量 +H = cell(I,K,K); % 信道系数 +for i=1:I + for k = 1:K + for j=1:K + H{i,k,j}=sqrt(1/2)*(randn(R,T)+1i*randn(R,T)); % 圆复高斯信道 + end + end +end + +U = cell(I,K); +U(:)={zeros(R,d)}; + +% 随机初始化发射波束向量 +V = cell(I,K); % 算法第一行 +for i=1:I + for k=1:K + v = randn(T,d)+1i*randn(T,d); % 随机初始化 + V{i,k}=sqrt(P/I)*v/norm(v,"fro"); + end +end + +% 求初始化发射波束V后求系统和速率 +rate_old = sum_rate1(H,V,sigma2,R,I,K,alpha1); +rate = [rate rate_old]; + + +iter1 = 1; +while(1) + U = find_U1(H,V,sigma2,R,I,K,d); % Tbale I line 4 in p.4435 + W = find_W1(U,H,V,I,K,d); % Tbale I line 5 in p.4435 + V = find_V1(alpha1,H,U,W,T,I,K,P); % Tbale I line 6 in p.4435 + rate_new = sum_rate1(H,V,sigma2,R,I,K,alpha1); + rate = [rate rate_new]; + elapsed_time = toc; + elapsed_time = elapsed_time - begin_time; + time = [time elapsed_time]; + iter1 = iter1 + 1; + if abs(rate_new-rate_old) / rate_old < epsilon || iter1 > max_iter + break; + end + rate_old = rate_new; +end +toc +% plot(0:iter1-1,rate,'r-o') +% grid on +% xlabel('Iterations') +% ylabel('Sum rate (bits per channel use)') +% set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1) +% title('MIMO-IFC, K=4, T=3, R=2, \epsilon=1e-3','Interpreter','tex') +%title('SISO-IFC, K=3, T=1, R=1, \epsilon=1e-3','Interpreter','tex') \ No newline at end of file diff --git a/find_U.m b/find_U.m new file mode 100644 index 0000000..74abc26 --- /dev/null +++ b/find_U.m @@ -0,0 +1,14 @@ +function U = find_U(H_hat,X,sigma2, P, R,I,d) + + J = zeros(R,R,I); %计算不含噪声项的矩阵 + IN = zeros(R,R,I); %计算含噪声项的矩阵 + U = zeros(R,d,I); + + for i=1:I + for l=1:I + J(:,:,i) = J(:,:,i) + H_hat(R*(i-1)+1:R*i, :)*X(:,:,l)*(X(:,:,l)')*(H_hat(R*(i-1)+1:R*i, :)'); % R-WMMSE公式24括号里面求信号项 + IN(:,:,i) = IN(:,:,i) + sigma2 / P * trace(H_hat*X(:,:,l)*(X(:,:,l)'))* eye(R); % R-WMMSE公式24括号里面求干扰项 + end + U(:,:,i) = (J(:,:,i) + IN(:,:,i)) \ (H_hat(R*(i-1)+1:R*i, :)*X(:,:,i)); + end +end \ No newline at end of file diff --git a/find_U1.m b/find_U1.m new file mode 100644 index 0000000..30ab09f --- /dev/null +++ b/find_U1.m @@ -0,0 +1,17 @@ +function U = find_U1(H,V,sigma2,R,I,K,d) + J = cell(I,K); + J(:)={zeros(R,R)}; + U = cell(I,K); + U(:)={zeros(R,d)}; + for i=1:I + for k=1:K + for j=1:K + for l=1:I + J{i,k} = J{i,k} + H{i,k,j}*V{l,j}*(V{l,j}')*(H{i,k,j}'); % 算法Table I, 第四行括号求和的部分 + end + end + J{i,k} = J{i,k} + sigma2*eye(R); % 算法Table I, 第四行括号求和加上噪声项的部分 + U{i,k} = J{i,k}\H{i,k,k}*V{i,k}; % 算法Table I, 第四行括号求逆，然后乘以H乘以V + end + end +end \ No newline at end of file diff --git a/find_V1.m b/find_V1.m new file mode 100644 index 0000000..18bb4ea --- /dev/null +++ b/find_V1.m @@ -0,0 +1,47 @@ +function V = find_V1(alpha1,H,U,W,T,I,K,P) + V = cell(I,K); + A = cell(K); + A(:) = {zeros(T,T)}; + + + for k=1:K % 公式15括号内求和部分 + for j=1:K + for l=1:I + A{k} = A{k} + alpha1(l,j)*H{l,j,k}'*U{l,j}*W{l,j}*(U{l,j}')*H{l,j,k}; + end + end + end + + max_iter = 100; % 二分法查找最优对偶变量\mu + mu = zeros(K,1); + for k=1:K % 对每个基站迭代寻找最优\mu + mu_min = 0; + mu_max = 10; + iter = 0; + while(1) + mu1 = (mu_max+mu_min) / 2; + P_tem = 0; + for i=1:I % 计算功率和 + V_tem = inv((A{k}+mu1*eye(T)))*alpha1(i,k)*((H{i,k,k}')*U{i,k}*W{i,k}); % 公式15 + P_tem = P_tem + real(trace(V_tem*V_tem')); + end + if P_tem > P + mu_min = mu1; + else + mu_max = mu1; + end + iter = iter + 1; + + if abs(mu_max - mu_min) < 1e-5 || iter > max_iter + break + end + end + mu(k) = mu1; + end + + for i=1:I + for k=1:K + V{i,k} = inv((A{k}+mu(k)*eye(T)))*alpha1(i,k)*((H{i,k,k}')*U{i,k}*W{i,k}); % 公式15 + end + end +end \ No newline at end of file diff --git a/find_W.m b/find_W.m new file mode 100644 index 0000000..9facb79 --- /dev/null +++ b/find_W.m @@ -0,0 +1,6 @@ +function W = find_W(U,H_hat,X, R, I,d) + W = zeros(d,d,I); + for i=1:I + W(:,:,i) = inv(eye(d)-U(:,:,i)'*H_hat(R*(i-1)+1:R*i, :)*X(:,:,i)); + end +end \ No newline at end of file diff --git a/find_W1.m b/find_W1.m new file mode 100644 index 0000000..6f856ce --- /dev/null +++ b/find_W1.m @@ -0,0 +1,9 @@ +function W = find_W1(U,H,V,I,K,d) + W = cell(I,K); + W(:) = {zeros(d,d)}; + for i=1:I + for k=1:K + W{i,k} = inv(eye(d)-U{i,k}'*H{i,k,k}*V{i,k}); % 算法Table I, 第五行 + end + end +end \ No newline at end of file diff --git a/find_X.m b/find_X.m new file mode 100644 index 0000000..ee5a64d --- /dev/null +++ b/find_X.m @@ -0,0 +1,68 @@ +function X = find_X(alpha1, H_hat,sigma2, U, W, T , R ,I ,d ,P ) + + eta = alpha1(1,1) * trace(U(:,:,1)*W(:,:,1)*(U(:,:,1)')); + W_hat = alpha1(1,1) * W(:,:,1); + U_hat = U(:,:,1); + % 计算分块矩阵，为了进一步降低复杂度 + for i = 2:I + W_hat = blkdiag(W_hat, alpha1(i,1) * W(:,:,i)); + U_hat = blkdiag(U_hat, U(:,:,i)); + eta = eta + alpha1(i,1) * trace(U(:,:,i)*W(:,:,i)*(U(:,:,i)')); + end + + eta = eta * sigma2 / P; + % 计算低维度的坍塌矩阵X_hat + X_hat = U_hat / (eta * inv(W_hat) + U_hat' * H_hat * U_hat); + + % 将低维度的坍塌矩阵X_hat由二维矩阵形式转换为三维矩阵形式以兼容代码数据格式 + X = zeros(R*I,d,I); % 每个基站的坍塌波束赋形矩阵 + for i = 1:I + X(:,:,i) = X_hat(:, d*(i-1) + 1:d*i); + end + + % % 另外一种稍慢的表达方式 + % J=zeros(I*R, I*R, I); + % IN=zeros(I,1); + % X=zeros(R*I,d, I); + % for i=1:I + % + % for l=1:I + % J(:,:,i) = J(:,:,i) + alpha1(l, 1) * H_hat(R*(l-1)+1:R*l, :)'*U(:,:,l)*W(:,:,l)*(U(:,:,l)')*(H_hat(R*(l-1)+1:R*l, :)); % R-WMMSE公式24括号里面求信号项 + % IN(i) = IN(i) + alpha1(l, 1) * trace(U(:,:,l)*W(:,:,l)*(U(:,:,l)')); % R-WMMSE公式26括号里面求干扰项 + % end + % + % X(:,:,i) = (J(:,:,i) + IN(i) * H_hat * sigma2 / P) \ (alpha1(i, 1) * (H_hat(R*(i-1)+1:R*i, :)'*U(:,:,i)*W(:,:,i))); + % + % end +end + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/sum_rate.m b/sum_rate.m new file mode 100644 index 0000000..03da9b5 --- /dev/null +++ b/sum_rate.m @@ -0,0 +1,14 @@ +function system_rate = sum_rate(H,V,sigma2,R,I,alpha1) + rate = zeros(I,1); + for i=1:I + denominator = zeros(R,R); + for l=1:I + denominator = denominator + H(:,:,i)*V(:,:,l)*V(:,:,l)'*H(:,:,i)'; + end + numerator = H(:,:,i)*V(:,:,i)*V(:,:,i)'*H(:,:,i)'; + denominator = denominator - numerator + sigma2*eye(R); + + rate(i) = log2(det(eye(R)+numerator / denominator)); + end + system_rate = real(sum(rate.*alpha1,'all')); +end \ No newline at end of file diff --git a/sum_rate1.m b/sum_rate1.m new file mode 100644 index 0000000..aa32083 --- /dev/null +++ b/sum_rate1.m @@ -0,0 +1,17 @@ +function system_rate = sum_rate1(H,V,sigma2,R,I,K,alpha1) + rate = zeros(I,K); + for i=1:I + for k = 1:K + temp = zeros(R,R); + for l=1:I + for j=1:K + if l~=i || j~=k + temp = temp + H{i,k,j}*V{l,j}*(V{l,j}')*(H{i,k,j}'); + end + end + end + rate(i,k) = log2(det(eye(R)+H{i,k,k}*V{i,k}*(V{i,k}')*(H{i,k,k}')*inv(temp + sigma2*eye(R)))); % 公式2 + end + end + system_rate = real(sum(rate.*alpha1,'all')); +end \ No newline at end of file