% 2017.4.14
% 算法比较: Proposed WSR[17] Sum-MMSE w/o DL ZF-CI[13] MMSE-CI-based[13]
% 4用户对,全部用户功率限制条件(Ps=Pr=1,信号功率归一化)
% 仿真平均误码率(Average BER)性能曲线:
clc;
clear all;
Nb = 8;
Nr = 8;
Nk = 2;
K = 4;
a = 1;
Ps = a;
Pr = Ps; % 信号功率归一化
SNR = 0:2:14; % 信噪比相同设置
frame_length = 10000; % 传输信号的帧结构
Imax = 20; % 最大迭代次数
Block = 100; % 最大蒙特卡罗试验次数设置
for i_p = 1:length(SNR)
snr = 10^(SNR(i_p)/10);
b = a/snr;
sigma = sqrt(b/2); % 噪声功率随着SNR的改变
%sigma = 1;
loop = 0;
err(i_p) = 0;
ber(i_p) = 0;
err2(i_p) = 0;
ber2(i_p) = 0;
err5(i_p) = 0;
ber5(i_p) = 0;
err_zf(i_p) = 0;
ber_zf(i_p) = 0;
err_mmse(i_p) = 0;
ber_mmse(i_p) = 0;
while loop < Block
loop = loop + 1;
% 发射端产生信号调制QPSK
data0 = randint(Nk, frame_length);% 随机产生发送数据
data = data0(:)';
I = zeros(Nk,frame_length/2);
Q = zeros(Nk,frame_length/2);
I = (2*data0(:,1:2:frame_length)-1)/sqrt(2);
Q = (2*data0(:,2:2:frame_length)-1)/sqrt(2);
symbols = I+sqrt(-1)*Q; % sqrt(-1)将实数变成复数
s = symbols(:,:); % QPSK调制后的输入信号
s1 = s;
s2 = s;
s3 = s;
s4 = s;
Hrb = sqrt(1/2) * (randn(Nr,Nb) + sqrt(-1) * randn(Nr,Nb));
for k = 1:K
Hdb{k} = sqrt(1/2) * (randn(Nk,Nb) + sqrt(-1) * randn(Nk,Nb));
Hdr{k} = sqrt(1/2) * (randn(Nk,Nr) + sqrt(-1) * randn(Nk,Nr));
end
% getting noise power
n_relay = sigma *(randn(Nr,frame_length/2) + sqrt(-1) * randn(Nr,frame_length/2)); % 中继处噪声向量
for k = 1:K
n_user{k} = sigma *(randn(Nk,frame_length/2) + sqrt(-1) * randn(Nk,frame_length/2));% 用户处噪声向量
end
%% Proposed algorithm
i1 = 1;
B1(:,:,i1) = sqrt(Ps/Nb).* eye(Nb);
B1_1(:,:,i1) = B1(:,2*1-1:2*1); % 发射端用户预编码矩阵初始化
B1_2(:,:,i1) = B1(:,2*2-1:2*2);
B1_3(:,:,i1) = B1(:,2*3-1:2*3);
B1_4(:,:,i1) = B1(:,2*4-1:2*4);
F1(:,:,i1) = sqrt(Pr/trace(Hrb * B1(:,:,i1) * B1(:,:,i1)'* Hrb'+ b .* eye(Nr))).* eye(Nr);
for k1 = 1:K
H1{k1}(:,:,i1) = [Hdb{k1};Hdr{k1} * F1(:,:,i1) * Hrb];
Z_z1{k1}(:,:,i1) = [b .* eye(Nk) zeros(Nk) ;zeros(Nk) b .* Hdr{k} * F1(:,:,i1) * F1(:,:,i1)'* Hdr{k}'+ b .* eye(Nk)];
W1{k1}(:,:,i1) = B1(:,2*k1-1:2*k1)'* H1{k1}(:,:,i1)'* inv(H1{k1}(:,:,i1) * B1(:,:,i1) * B1(:,:,i1)'* H1{k1}(:,:,i1)'+ Z_z1{k1}(:,:,i1));
W1_1{k1}(:,:,i1) = W1{k1}(:,1:Nk);
W1_2{k1}(:,:,i1) = W1{k1}(:,Nk+1:2*Nk);
E1{k1}(:,:,i1) = eye(Nk) - W1{k1}(:,:,i1) * H1{k1}(:,:,i1) * B1(:,2*k1-1:2*k1) - B1(:,2*k1-1:2*k1)'* H1{k1}(:,:,i1)'* W1{k1}(:,:,i1)'+ W1{k1}(:,:,i1) * H1{k1}(:,:,i1) * B1(:,:,i1) * B1(:,:,i1)'* H1{k1}(:,:,i1)'* W1{k1}(:,:,i1)'+ W1{k1}(:,:,i1) * Z_z1{k1}(:,:,i1) * W1{k1}(:,:,i1)';
R1{k1}(i1) = -log2(det(E1{k1}(:,:,i1)));
end
while(i1<=Imax)
% 拉格朗日函数法—KKT条件(二分法)
% 引入中间变量
V0_1(:,:,i1) = zeros(Nb);
for k1 = 1:K
V1(:,:,i1) = V0_1(:,:,i1) + H1{k1}(:,:,i1)'* W1{k1}(:,:,i1)'* W1{k1}(:,:,i1)* H1{k1}(:,:,i1);
V0_1(:,:,i1) = V1(:,:,i1);
C1{k1}(:,:,i1) = H1{k1}(:,:,i1)'* W1{k1}(:,:,i1)';
end
%
e1 = 0.00001;
min1 = 0.001;
max1 = 1e5;
mid_b1 = (max1 + min1)/2;
bb1_1(i1) = trace((inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{1}(:,:,i1)) * (inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{1}(:,:,i1))');
bb2_1(i1) = trace((inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{2}(:,:,i1)) * (inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{2}(:,:,i1))');
bb3_1(i1) = trace((inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{3}(:,:,i1)) * (inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{3}(:,:,i1))');
bb4_1(i1) = trace((inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{4}(:,:,i1)) * (inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{4}(:,:,i1))');
while (abs(bb1_1(i1) + bb2_1(i1) + bb3_1(i1) + bb4_1(i1) - Ps)>e1)
if bb1_1(i1) + bb2_1(i1) + bb3_1(i1) + bb4_1(i1) > Ps
min1 = mid_b1;
else
max1 = mid_b1;
end
mid_b1 = (max1 + min1)/2;
bb1_1(i1) = trace((inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{1}(:,:,i1)) * (inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{1}(:,:,i1))');
bb2_1(i1) = trace((inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{2}(:,:,i1)) * (inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{2}(:,:,i1))');
bb3_1(i1) = trace((inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{3}(:,:,i1)) * (inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{3}(:,:,i1))');
bb4_1(i1) = trace((inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{4}(:,:,i1)) * (inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{4}(:,:,i1))');
end
B1_1(:,:,i1+1) = inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{1}(:,:,i1);
B1_2(:,:,i1+1) = inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{2}(:,:,i1);
B1_3(:,:,i1+1) = inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{3}(:,:,i1);
B1_4(:,:,i1+1) = inv(V0_1(:,:,i1) + mid_b1.*eye(Nb)) * C1{4}(:,:,i1);
B1(:,:,i1+1) = [B1_1(:,:,i1+1) B1_2(:,:,i1+1) B1_3(:,:,i1+1) B1_4(:,:,i1+1)];
% 第四步:更新中继预编码矩阵F
% 方法一:(QCQP问题)—CVX优化工具箱
% M1(:,:,i1) = W1_2{1}(:,:,i1) * Hdr{1};
% M2(:,:,i1) = W1_2{2}(:,:,i1) * Hdr{2};
% M3(:,:,i1) = W1_2{3}(:,:,i1) * Hdr{3};
% M4(:,:,i1) = W1_2{4}(:,:,i1) * Hdr{4};
% U1(:,:,i1) = Hrb * B1_1(:,:,i1+1);
% U2(:,:,i1) = Hrb * B1_2(:,:,i1+1);
% U3(:,:,i1) = Hrb * B1_3(:,:,i1+1);
% U4(:,:,i1) = Hrb * B1_4(:,:,i1+1);
% S1(:,:,i1) = Hrb * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hdb{1}'* W1_1{1}(:,:,i1)';
% S2(:,:,i1) = Hrb * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hdb{2}'* W1_1{2}(:,:,i1)';
% S3(:,:,i1) = Hrb * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hdb{3}'* W1_1{3}(:,:,i1)';
% S4(:,:,i1) = Hrb * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hdb{4}'* W1_1{4}(:,:,i1)';
% A1(:,:,i1) = kron(conj((Hrb * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hrb')')+ b .* eye(Nr),M1(:,:,i1)'* M1(:,:,i1));
% A2(:,:,i1) = kron(conj((Hrb * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hrb')')+ b .* eye(Nr),M2(:,:,i1)'* M2(:,:,i1));
% A3(:,:,i1) = kron(conj((Hrb * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hrb')')+ b .* eye(Nr),M3(:,:,i1)'* M3(:,:,i1));
% A4(:,:,i1) = kron(conj((Hrb * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hrb')')+ b .* eye(Nr),M4(:,:,i1)'* M4(:,:,i1));
% c1(:,i1) = vec(M1(:,:,i1)'* U1(:,:,i1)');
% c2(:,i1) = vec(M2(:,:,i1)'* U2(:,:,i1)');
% c3(:,i1) = vec(M3(:,:,i1)'* U3(:,:,i1)');
% c4(:,i1) = vec(M4(:,:,i1)'* U4(:,:,i1)');
% d1(:,i1) = vec(M1(:,:,i1)'* S1(:,:,i1)');
% d2(:,i1) = vec(M2(:,:,i1)'* S2(:,:,i1)');
% d3(:,i1) = vec(M3(:,:,i1)'* S3(:,:,i1)');
% d4(:,i1) = vec(M4(:,:,i1)'* S4(:,:,i1)');
% t1(i1) = trace(eye(Nk) - W1_1{1}(:,:,i1) * Hdb{1} * B1_1(:,:,i1+1) - B1_1(:,:,i1+1)'* Hdb{1}'* W1_1{1}(:,:,i1)'+ W1_1{1}(:,:,i1) * Hdb{1} * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hdb{1}'*W1_1{1}(:,:,i1)'+ W1_1{1}(:,:,i1) * W1_1{1}(:,:,i1)'+ W1_2{1}(:,:,i1) * W1_2{1}(:,:,i1)');
% t2(i1) = trace(eye(Nk) - W1_1{2}(:,:,i1) * Hdb{2} * B1_2(:,:,i1+1) - B1_2(:,:,i1+1)'* Hdb{2}'* W1_1{2}(:,:,i1)'+ W1_1{2}(:,:,i1) * Hdb{2} * B1(:,:,i1+1) * B1(:,:,i1+1)'* Hdb{2}'*W1_1{2}(:,:,i1)'+ W1_1{2}(:,:,i1) * W1_1{2}(:,:,i1)'+ W1_2{2}(:,:,i1) * W1_2{2}(:,:,i1)');
% t3(i1