% No direct Link
% written by Peilan Wang 10/28/2019
% test beamforming gain and nmse of OMP,BPDN, GAMP for reflected link only
%% 
clear all, clc,close all,
test_mode = 0; % 0--> nmse vs. SNR ; 1---> nmse vs T, SNR = 20; 2--> Suc_rate vs T
global M
global Nr
global Nt
global T




%% fixed
M_h = 8; 
M_v = 8;
M = M_h*M_v;
Nt = 16;% 32; 
grid_rh = 32;
grid_rv = 32;
grid_r = grid_rh*grid_rv;
grid_t = 64;% 32;
Nr = M;

numiter = 2e2;
snr = 10;
noise = 1;

T_L = 5:25:135


save_Channel_T;

Lr =length(T_L);

NMSE_bp = zeros(length(T_L),numiter);
NMSE_omp = zeros(length(T_L),numiter);
NMSE_gamp = zeros(length(T_L),numiter);
NMSE_ls = zeros(length(T_L),numiter);
NMSE_ls_psur = zeros(length(T_L),numiter);


Suc_rec_bp = zeros(length(T_L),numiter);
Suc_rec_omp = zeros(length(T_L),numiter);
Suc_rec_gamp = zeros(length(T_L),numiter);
Suc_rec_ls= zeros(length(T_L),numiter);
Suc_rec_ls_psur= zeros(length(T_L),numiter);

Beam_Gain_0 = zeros(length(T_L),numiter);
Beam_Gain_omp = zeros(length(T_L),numiter);
Beam_Gain_bp = zeros(length(T_L),numiter);
Beam_Gain_gamp = zeros(length(T_L),numiter);
Beam_Gain_ls = zeros(length(T_L),numiter);
Beam_Gain_ls_psur = zeros(length(T_L),numiter);
% VecH0 = zeros(length(T_L),numiter,Nr*Nt);
% VecH_omp = zeros(length(T_L),numiter,Nr*Nt);
% VecH_bp = zeros(length(T_L),numiter,Nr*Nt);
% VecH_gamp = zeros(length(T_L),numiter,Nr*Nt);

for it = 1: Lr
      
       T = T_L(it);
%         K_S = T_L(it)
  for im = 1: numiter

       


%             [bGain_0,bGain_omp,bGain_bp,bGain_gamp,bGain_ls,nmse_omp,nmse_bp,nmse_gamp,nmse_ls] = test_noDirec_fen(T,SNR,noise)
            
%             [bGain_0,bGain_omp,bGain_bp,bGain_gamp,bGain_ls,bGain_ls_psur,nmse_omp,nmse_bp,nmse_gamp,nmse_ls,nmse_ls_psur] = test_UPA_L0323(T,snr,noise)
 [bGain_0,bGain_omp,bGain_bp,bGain_gamp,bGain_ls,bGain_ls_psur,...
   nmse_omp,nmse_bp,nmse_gamp,...
    nmse_ls,nmse_ls_psur] = test_UPA_L0331_T(snr,noise,ATRR_ls,H_cas,Atrr,t0,vecH_0,al_Sig_vec);

            Beam_Gain_0(it,im) = bGain_0;
            Beam_Gain_omp(it,im) = bGain_omp;
            Beam_Gain_bp(it,im) = bGain_bp;
            Beam_Gain_gamp(it,im) = bGain_gamp;
            Beam_Gain_ls(it,im) = bGain_ls;
            Beam_Gain_ls_psur(it,im) = bGain_ls_psur;
            
            
            NMSE_ls_psur(it,im) = nmse_ls_psur;
            NMSE_ls(it,im) = nmse_ls;
            NMSE_bp(it,im) = nmse_bp;
            NMSE_omp(it,im) = nmse_omp;
            NMSE_gamp(it,im) = nmse_gamp;

            if nmse_bp<1e-2
                Suc_rec_bp(it,im) = 1;
            end
            if nmse_omp<1e-2
                Suc_rec_omp(it,im) = 1;
            end
            if nmse_gamp<1e-2
                Suc_rec_gamp(it,im) = 1;
            end
            if nmse_ls<1e-2
                Suc_rec_ls(it,im) = 1;
            end
            if nmse_ls_psur<1e-2
                Suc_rec_ls_psur(it,im) = 1;
            end
            
    end
end

NMSE_bp_mean = mean(NMSE_bp,2);
NMSE_omp_mean = mean(NMSE_omp,2);
NMSE_gamp_mean = mean(NMSE_gamp,2);
NMSE_ls_mean = mean(NMSE_ls,2);
NMSE_ls_psur_mean = mean(NMSE_ls_psur,2);

Beam_Gain_0_mean = (mean(Beam_Gain_0,2));
Beam_Gain_omp_mean = (mean(Beam_Gain_omp,2));
Beam_Gain_bp_mean = (mean(Beam_Gain_bp,2));
Beam_Gain_gamp_mean =( mean(Beam_Gain_gamp,2));
Beam_Gain_ls_mean =( mean(Beam_Gain_ls,2));
Beam_Gain_ls_psur_mean = ( mean(Beam_Gain_ls_psur,2));

% save   %% 2019/11/12
save qw%newT_0331_path2
figure(1)
% 
% semilogy(T_L,NMSE_bp_mean,'o--',T_L,NMSE_omp_mean,'+--',T_L,NMSE_gamp_mean,'^--',T_L,NMSE_ls_mean,'*--','linewidth',2)
% legend('BPDN','OMP','GAMP','LS')
% xlabel('T')
% ylabel('NMSE')
% grid on


semilogy(T_L,NMSE_bp_mean,'o--',T_L,NMSE_omp_mean,'+--',T_L,NMSE_gamp_mean,'^--',T_L,NMSE_ls_mean,'*--',T_L,NMSE_ls_psur_mean,'x--','linewidth',2)
legend('BPDN','OMP','GAMP','LS','LS-perfect')
xlabel('T')
ylabel('NMSE')
grid on


figure(2)
semilogy(T_L,Beam_Gain_0_mean,'*--',T_L,Beam_Gain_omp_mean,'o--',T_L,Beam_Gain_bp_mean,'+--',T_L,Beam_Gain_gamp_mean,'^--',T_L,Beam_Gain_ls_mean,'*--',T_L,Beam_Gain_ls_psur_mean,'x--','linewidth',2)
legend('Perfect Channel','OMP','BPDN','GAMP','LS','LS-perfect')
xlabel('T')
ylabel('Receive Signal Power')
grid on
