非线性BD-GMD均分功率用户排序误码率

clc;
clear all;
close all;
warning off;
%------------Attention 没有必要做功率归一化--------
%------------计算误码率时 要除以M,而不是log2(M)-----
%------------GMD分解自带不与 论文中一致 ------------
%------------h最好每星座点变一次，不然Frame_num过大的话造成BER不稳定---------
factor=0;
%设置仿真参数
Nt =16; % Tx 天线
Nr =16; % Rx 天线
%UserAntennaNum的累加和为H的行数
UserNum=4;
UserAntennaNum_old=[ 4 4 4 4];%对比不同天线数目下，BER性能
    stream_num_old=[ 4 4 4 4];
    %stream_num=[3 3 3;2 2 5;2 3 4;3 4 2;2 4 3; 5 2 2];
    %stream_num=[3 3 4;2 2 6;2 3 4;4 4 2;2 4 4; 6 2 2];
[row_test col_test]=size(UserAntennaNum_old);

SymPerFrame = 1; 
Frame_num = 50; % 仿真中帧的数目
M =16;% QAM调制 
P_all=5;
% the sum power of all stream is P_all
%PowerScale= sqrt(2*(M-1)/3);
ModPeriod  = 2*sqrt(M) ;% mod操作数，用于THP预编码                
% 设置信噪比范围
    SNR = [-5,-3:3:24];
%初始化误比特率
Len = length(SNR);
ber_bdgmd_noorder=zeros(row_test,Len);
ber_bdgmd_optimalorder=zeros(row_test,Len);

        
            %PowerScale = 2*M/3;% 功率归一化

            PowerScale_linear_GMD=2*M/3;
            

            %由于星座点映射之后的功率时2(M-1)/3,再经过mod之后放大M/(M-1),所以mod后功率为PowerScale，而每个流上的最终功率为
            %SigPowerDbw_of_each_stream，所以beta为功率控制系数，可为系数也可为对角矩阵

%-----------------循环过程开始-----------------
for index = 1:Len;%SNR取值变化循环
%--------------------- to send the program's progress-------------------
     fprintf('%s %d %s \n', 'The progress of this program is ',SNR(index)/max(SNR)*100,'%');
  
  snr = 10.^(SNR(index)/10);% 将SNR转为线性域的变化范围    
 %   noise_power = Nt/snr;% 噪声功率

%-----------------帧的变化范围，在每一个帧内信道保持不变----------------- 
    for cnt = 1:Frame_num
        
         H1 = sqrt(1/2)*(randn(UserAntennaNum_old(1),Nt)+j*randn(UserAntennaNum_old(1),Nt));%随机生成信道矩阵
         H2 = sqrt(1/2)*(randn(UserAntennaNum_old(2),Nt)+j*randn(UserAntennaNum_old(2),Nt));%随机生成信道矩阵
         H3 = sqrt(1/2)*(randn(UserAntennaNum_old(3),Nt)+j*randn(UserAntennaNum_old(3),Nt));%随机生成信道矩阵
         H4 = sqrt(1/2)*(randn(UserAntennaNum_old(3),Nt)+j*randn(UserAntennaNum_old(3),Nt));%随机生成信道矩阵
   %%  no user order 均分功率
         %用户排序
%        
         [H UserAntennaNum stream_num]=random_order(H1,H2,H3,H4,UserAntennaNum_old,stream_num_old);
         
%-----------------每一帧中包含2*100个符号，循环从1到100-----------------
        for sym_index = 1:SymPerFrame;
           
             H_est=H; %完美信道信息
         for n=1:row_test
             %n代表同时配置的几种天线接收方案
             
             %--------------------产生信号源------------------------              
            Source = randint(sum(stream_num(n,:)),1,M);%生成2*1的符号，符号取值范围是0~15
            Sym = qammod(Source,M,0,'gray'); %进行QAM调制和功率归一化
            SigPower_of_each_stream=P_all/sum(stream_num(n,:));%按流均分功率时，将要发射时stream功率       
            
            SigPowerDbw_of_each_stream=10*log10(SigPower_of_each_stream); 
   
            %----------------非线性BD-GMD预编码----------------------------
            % 功率分配
            beta_nonlinear_GMD=sqrt(SigPower_of_each_stream/PowerScale_linear_GMD);
            % 预编码
            [d_thp,L,Q,G,PPP]= GMD_encoder(Sym,H_est,ModPeriod,UserNum,UserAntennaNum(n,:),stream_num(n,:),beta_nonlinear_GMD);
            % @para Sym ,just is one MQAM symbol,it is a complex number.
            % @para H_est ,the radnom channel which has been added the
            %  noise,here it is a 4*4 complex matrix.
            
            % beta_water_BD_GMD功率分配矩阵求倒数
            %[ beta_water_BDGMD_mine_inverse nonzero_stream_index]=water_inverse(diag(beta_water_BDGMD_mine));
            % MIMO channel
             HH=H*d_thp;
             y_thp = awgn(HH,SNR(index),SigPowerDbw_of_each_stream);
             % 注水定理只体现在噪声上
                   
            % Caculate channel capacity
            

            
            G=diag(1./diag(G));
            
            % receiver end 
            r_thp = G*Q'*y_thp/beta_nonlinear_GMD;         
            Rec_thp = mod_thp(r_thp,ModPeriod);% mod
            % QAM demodulator 
           % Rec_thp = Rec_thp_1;% Scale the symbol to QAM symbol            
            Rec_Data_thp = qamdemod(Rec_thp,M,0,'gray');
            % BER 只统计功率分配不为0 的流的误码率
            [err ratio] = biterr(Rec_Data_thp,Source,log2(M));% ber           
            ber_bdgmd_noorder(n,index) = ber_bdgmd_noorder(n,index) + err;
           
            % 计算总流数
            %stream_num_BD_GMD_mine_total(n,index)=stream_num_BD_GMD_mine_total(n,index)+length(nonzero_stream_index);
            
         

          end

        end; % loop for num   
       
                  
    %% 功率分配 遍历排序 
                                   
       [H,UserAntennaNum,stream_num]=optimal_order(H1,H2,H3,H4,UserNum,UserAntennaNum_old,stream_num_old);
       
       for sym_index = 1:SymPerFrame;
           
             H_est=H; %完美信道信息
         for n=1:row_test
             %n代表同时配置的几种天线接收方案
             
             %--------------------产生信号源------------------------              
            Source = randint(sum(stream_num(n,:)),1,M);%生成2*1的符号，符号取值范围是0~15
            Sym = qammod(Source,M,0,'gray'); %进行QAM调制和功率归一化
            SigPower_of_each_stream=P_all/sum(stream_num(n,:));%按流均分功率时，将要发射时stream功率       
            
            SigPowerDbw_of_each_stream=10*log10(SigPower_of_each_stream); 
            % 功率分配
            beta_nonlinear_GMD=sqrt(SigPower_of_each_stream/PowerScale_linear_GMD);
            % 预编码
            [d_thp,L,Q,G,PPP]=GMD_encoder(Sym,H_est,ModPeriod,UserNum,UserAntennaNum(n,:),stream_num(n,:),beta_nonlinear_GMD);
            % @para Sym ,just is one MQAM symbol,it is a complex number.
            % @para H_est ,the radnom channel which has been added the
            %  noise,here it is a 4*4 complex matrix.
            
            
            % beta_water_BD_GMD功率分配矩阵求倒数
            %[ beta_water_BDGMD_mine_inverse nonzero_stream_index]=water_inverse(diag(beta_water_BDGMD_mine));
            % MIMO channel
             HH=H*d_thp;
             y_thp = awgn(HH,SNR(index),SigPowerDbw_of_each_stream);
             % 注水定理只体现在噪声上
             %noise=y_thp-HH;
             %y_thp_wf=HH+diag(beta_water_BDGMD_mine_inverse)*noise;
            
            % Caculate channel capacity
  
            
            G=diag(1./diag(G));
            
            % receiver end 
            r_thp = G*Q'*y_thp/beta_nonlinear_GMD;         
            Rec_thp = mod_thp(r_thp,ModPeriod);% mod
            % QAM demodulator 
           % Rec_thp = Rec_thp_1;% Scale the symbol to QAM symbol            
            Rec_Data_thp = qamdemod(Rec_thp,M,0,'gray');
            % BER 只统计功率分配不为0 的流的误码率
            [err ratio] = biterr(Rec_Data_thp,Source,log2(M));% ber           
            ber_bdgmd_optimalorder(n,index) = ber_bdgmd_optimalorder(n,index) + err;
            %err_rem((cnt-1)*SymPerFrame+sym_index)=err;
            % 计算总流数
            %stream_num_BD_GMD_mine_total(n,index)=stream_num_BD_GMD_mine_total(n,index)+length(nonzero_stream_index);
            
         end
       end
    end; % loop for iteration          
end % loop for snr
%----------------------------循环结束----------------------------
stream_num_of_per_loop=sum(stream_num);
ber_bdgmd_noorder=ber_bdgmd_noorder./(Frame_num*SymPerFrame*stream_num_of_per_loop*sqrt(M));%计算误比特率

ber_bdgmd_optimalorder=ber_bdgmd_optimalorder./(Frame_num*SymPerFrame*stream_num_of_per_loop*sqrt(M));

%[3 3 4;2 2 6;2 3 5;4 4 2;2 4 4; 6 2 2]

%----------------保存当前workspace中所有变量值到filename.mat文件中---------

filename='Linear-BD-GMD-SVD-THP-16QAM-4user-uesrodrder-streamreduced-40000';%-----需要改，根据程序不同改名
nowtime=clock;
filename_time=sprintf('%s--%d.%d.%d.%d.%d.mat',filename,nowtime(1),nowtime(2),nowtime(3),nowtime(4),nowtime(5));
save(filename_time);
%-------------------------------画图-------------------------------------
figure(1)
semilogy(SNR,ber_bdgmd_noorder(1,:),'LineWidth',2,'Color',[1 0 0],'Marker','*');
hold on;

semilogy(SNR,ber_bdgmd_optimalorder(1,:),'LineWidth',2,'Color',[1 0 1],'Marker','+');

grid on;
xlabel('SNR(dB)');ylabel('BER');title('非线性BD-GMD均分功率用户排序误码率图');
%stream_num=[ 3 3 3;3 2 3;2 2 5;2 3 4;3 4 2;2 4 3;5 2 2];
legend('BD-GMD-noorder','BD-GMD-THP-optimalorder');

 

function  [Q,R,P] =BD_GMD_5_stream(H,UserNum,UserAntennaNum,stream_num,flag) 

% H 为信道 Nr*Nt
% H 按行分成Num个字h对应每个用户，UserAntennaNum为向量，对应每用户的天线数
% flag标示是 真二用户 1 假二用户 0 主要为了减少计算量

%将H分成各子信道矩阵
[m,n]=size(H);
I=eye(n);
%Q_end=[];
%R_end=zeros(m);
%P_end=[];
%

if UserNum~=2
    
    [Q1,R1,P1] = BD_GMD_5_stream(H,2,[UserAntennaNum(1) sum(UserAntennaNum(2:end))],[stream_num(1) sum(stream_num(2:end))],0);%假二用户
    H2=H(UserAntennaNum(1)+1:end,:);
%    rank_of_rightmatrix=rank((I-P(:,1:stream_num(1))*P(:,1:stream_num(1))'));
%    rank_of_old_H=rank(H(UserAntennaNum(1)+1:end,:));   
%    rank_of_new_H=rank(H2);
    [Q2,R2,P2] = BD_GMD_5_stream(H2*(I-P1*P1'),UserNum-1,UserAntennaNum(2:end),stream_num(2:end),1);
 
%    合成最后分解矩阵
    Q=blkdiag(Q1,Q2);
    R=blkdiag(R1,R2);
    R(stream_num(1)+1:end,1:stream_num(1))=Q2'*H2*P1;
    P=[P1 P2];
%      Q(UserAntennaNum(1)+1:end,stream_num(1)+1:end)=Q2;
%      R(stream_num(1)+1:end,stream_num(1)+1:end)=R2;
%      R(stream_num(1)+1:end,1:stream_num(1))=Q2'*H(UserAntennaNum(1)+1:end,:)*P(:,1:stream_num(1));
%      P(:,stream_num(1)+1:end)=P2;
      
else
if flag==0% 假二用户，不必对H2进行分解，只对H1进行分解即可 
    H1=H(1:UserAntennaNum(1),:);
    %求BD-GMD
   [Q,R,P]=gmd_zcy_streamreduce(H1,stream_num(1));
   % Q:Nr1*r1, R:r1*r1,P:Nt*r1
    
    
else %真二用户 需要对H1和H2都要分解
    
%%求二分矩阵的BD-GMD
H1=H(1:UserAntennaNum(1),:);
H2=H(UserAntennaNum(1)+1:UserAntennaNum(1)+UserAntennaNum(2),:);

%求BD-GMD
[Q1,R1,P1]=gmd_zcy_streamreduce(H1,stream_num(1));

HH=H2*(I-P1*P1');
[Q2,R2,P2]=gmd_zcy_streamreduce(HH,stream_num(2));
r=(Q2)'*H2*P1;
R=blkdiag(R1,R2);
[row col]=size(R1);
R(stream_num(1)+1:stream_num(1)+stream_num(2),1:col)=r;

P=[P1  P2];

Q=blkdiag(Q1,Q2);
end

%%

end