%% Fig.7 in Paper: Mean Value of Minimum Transmit Power (FDD with Cubic Uncertainty Regions)
close all;
clear all;

% Parameters
SNR0 = 10; % target SNR
varG = 1; % variance of AWGN at receiver
varH = 2; % variance of channel matrix
qSNRdB = 5:14; % quantization SNR (dB)
qSNR = 10.^(qSNRdB/10); % quantization SNR
qSTEP = sqrt(6*varH./qSNR); % quantization step

% Algorithms
P0 = zeros(4,10); % minimum required transmit power
N = 300; % number of iterations
for a = 1:N
    nT = 4; % number of transmit antennas
    nR = 4; % number of receive antennas
    H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
    for b = 1:10
        estH = qSTEP(b).*round(H./qSTEP(b)); % 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);
            max(abs(real(delta(:)))) <= qSTEP(b)/2;
            max(abs(imag(delta(:)))) <= qSTEP(b)/2;
            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'));
            max(abs(real(delta(:)))) <= qSTEP(b)/2;
            max(abs(imag(delta(:)))) <= qSTEP(b)/2;
        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);
            max(abs(real(delta(:)))) <= qSTEP(b)/2;
            max(abs(imag(delta(:)))) <= qSTEP(b)/2;
        cvx_end
        P0(3,b) = P0(3,b)+SNR0*varG/real(cvx_optval);
    end
    
    % Robust for nT = 6 and nR = 6
    nT = 6; % number of transmit antennas
    nR = 6; % number of receive antennas
    H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
    for b = 1:10
        estH = qSTEP(b).*round(H./qSTEP(b)); % 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);
            max(abs(real(delta(:)))) <= qSTEP(b)/2;
            max(abs(imag(delta(:)))) <= qSTEP(b)/2;
            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/minPowerFDD.mat', 'qSNRdB', 'P0'); % save data

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