%% Fig.6 in Paper: Mean Value of Minimum Transmit Power (TDD with Spherical Uncertainty Regions)
close all;
clear all;

% Parameters
SNR0 = 10; % target SNR
varG = 1; % variance of AWGN at receiver
varH = 2; % variance of channel matrix
estSNRdB = 8:14; % estimation SNR (dB)
estSNR = 10.^(estSNRdB/10); % estimation SNR
nT = 2; % number of transmit antennas
nR = 2; % number of receive antennas

% Algorithms
P0 = zeros(4,7); % minimum required transmit power
N = 2; % number of iterations
for a = 1:N
    H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
    E = randn(nR,nT) + i*randn(nR,nT); % estimation error matrix without specifying variance
    
    qosP = 0.85; % probability of providing required QoS
    r = sqrt(chi2inv(qosP, 2*nR*nT)*varH./(1+estSNR)); % radius of spherical uncertainty region
    for b = 1:7
        varE = varH/estSNR(b); % variance of error matrix
        estH = (H+E*sqrt(varE/2))*varH/(varH+varE); % estimated CSI matrix
        [U, ~] = eig(estH'*estH); % normalized eigenvectors and eigenvalues (ascending order)
        
        % Robust with all eigenmodes
        cvx_begin
            variable t;
            variable delta(nR,nT) complex;     
            dual variables gamma{nT}; % dual variables
            minimize (t);
            norm(delta(:)) <= r(b);
            for c = 1:nT
                gamma{c} : t >= quad_form((estH+delta)*U(:,c), diag(ones(nR,1)));
            end
        cvx_end
        P0(1,b) = P0(1,b)+SNR0*varG/real(cell2mat(gamma)'*diag(U'*(estH+delta)'*(estH+delta)*U));
    
        % Non-robust with maximum eigenmode
        cvx_begin
            variable delta(nR,nT) complex;
            minimize (norm((estH+delta)*U(:,nT), 'fro'));
            norm(delta(:)) <= r(b);
        cvx_end
        P0(2,b) = P0(2,b)+SNR0*varG/real(cvx_optval^2);
    
        % OSTBC with equal power over all eigenmodes
        cvx_begin
            variable delta(nR,nT) complex;
            f = 0;
            for c = 1:nT
                f = f+quad_form((estH+delta)*U(:,c), diag(ones(nR,1)))/nT;
            end
            minimize (f);
            norm(delta(:)) <= r(b);
        cvx_end
        P0(3,b) = P0(3,b)+SNR0*varG/real(cvx_optval);
    end
    
    % Robust for qosP = 0.6
    qosP = 0.6; % probability of providing required QoS
    r = sqrt(chi2inv(qosP, 2*nR*nT)*varH./(1+estSNR)); % radius of spherical uncertainty region
    for b = 1:7
        varE = varH/estSNR(b); % variance of error matrix
        estH = (H+E*sqrt(varE/2))*varH/(varH+varE); % estimated CSI matrix
        [U, ~] = eig(estH'*estH); % normalized eigenvectors and eigenvalues (ascending order)
        cvx_begin
            variable t;
            variable delta(nR,nT) complex;     
            dual variables gamma{nT}; % dual variables
            minimize (t);
            norm(delta(:)) <= r(b);
            for c = 1:nT
                gamma{c} : t >= quad_form((estH+delta)*U(:,c), diag(ones(nR,1)));
            end
        cvx_end
        P0(4,b) = P0(4,b)+SNR0*varG/real(cell2mat(gamma)'*diag(U'*(estH+delta)'*(estH+delta)*U));
    end
end
P0 = 10*log10(P0/N); % minimum required transmit power (dB)
save('data/minPowerTDD.mat', 'estSNRdB', 'P0'); % save data

% Figures
output_fig = 'figure/minPowerTDD';
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
set(gca, 'Fontsize', 16);
hold on;
grid on;
plot(estSNRdB, P0(1,:), '-o', estSNRdB, P0(2,:), '--o', estSNRdB, P0(3,:), '-.o', estSNRdB, P0(4,:), '-*', 'linewidth', 2);
legend('Robust: P_{in}=0.85, n_T=2, n_R=2', 'Non-Robust: P_{in}=0.85, n_T=2, n_R=2', ...
    'OSTBC: P_{in}=0.85, n_T=2, n_R=2', 'Robust: P_{in}=0.6, n_T=2, n_R=2');
xlabel('Estimation SNR (dB)');
ylabel('Mean Value of Minimum Transmit Power (dB)');
eval(['print -depsc2 ' output_fig]);