mutual_mimo_opt.m

实现mimo预编码环境中容量的所有计算,包括不同的准则

function [x_opt, Q_opt, fval_opt, fgap_opt] = mutual_mimo_opt(H0, Rt, M, N, SNR_db)

%----------------------------------------------------------------------
% Optimization function to find the ergodic capacity of a MIMO channel given
% the channel mean and transmit correlation matrices. The program returns
% the optimum variable values (which can be the eigenvalues of the transmit
% covariance matrix or the vectorized version of this matrix), the optimum
% transmit covariance matrix itself, the mutual information value and the
% function gap to optimality.
%
% Author: Mai Vu
% Date: Dec 04, 2003
% Modified: Mar 02, 2004; Apr 19, 2004.
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% INPUTS:
%
%   N: Number of transmit antennas
%   M: Number of receive antennas
%   H0: channel mean (MxN)
%   Rt: transmit correlation
%   SNR_db: SNR in dB
%
% THE PROBLEM
%
%   max f = E[logdet(I + SQ)]
%   s.t  trace(Q) = 1
%
%   where  S = H'*H
%          H = H0 + Hw*Rt_sqrt
%
% OUTPUTS:
%   
%   x_opt: In symmetric Q case, the vector contains re-order elements of Q
%   Q_opt: Optimum input covariance matrix for each SNR input
%   fval_opt: Mutual information value at each SNR input
%   fgap_opt: Gap to optimality of the MI value at each SNR input
%----------------------------------------------------------------------
  
% optimization parameters
NSAMPLE = 10000;
MAXITER = 10;
MAXBARR = 10;
MAXLINES = 50;
epsilon = 5*10^-6;
mu = 100;
barr_t = 100;

% signal to noise ratio
gamma = 10.^(SNR_db/10);

I = eye(N,N);
if (H0 == 0)
  % special case correlation only
  [U Rt] = eig(Rt);
  diag_flag = 1;
elseif (Rt/trace(Rt)*N == I)
  % special case mean only
  [V H0 U] = svd(H0);
  diag_flag = 1;
else
  diag_flag = 0;
end
Rt_sqrt = sqrtm(Rt);

% initial transmit covariance matrix
S0 = H0'*H0 + M*Rt;
Q = S0/trace(S0);
lambda = real(eig(Q));

% indices for calculating the gradient/hessian and reorganize the unknowns
% into a matrix
ind1 = [];
ind2 = [];
for n=1:N
  ind1 = [ind1, n:1:N];
  ind2 = [ind2, 1:1:N-n+1];
end

ind3 = [];
for n=1:N
  if n>1
    for k=1:n-1
      ind3 = [ind3 k+(n-k)*N-(n-k-1)*(n-k)/2 + N*(N-1)/2];
    end
  end
  for k=n:N
    ind3 = [ind3 n+(k-n)*N-(k-n-1)*(k-n)/2];
  end
end

% number of unknowns
Nr = (N+1)*N/2;
Nx = N^2;

% initial unknown values
Q12 = diag(Q(ind1,ind2));
x = [real(Q12(1:N))/2; real(Q12(N+1:Nr)); imag(Q12(N+1:Nr))];

% initial line search counter array
ls_cnt = [];
% maximum number of barrier outter steps
barr_cnt = ceil(log(1/epsilon/barr_t)/log(mu)) + 1;

% solve for diagonal Q
if (diag_flag)
  
  for k=1:length(gamma)
    disp(['  SNR = ', num2str(SNR_db(k)), 'dB']);
    [lambda_new, fval2, fgap2, cnt2(k)] = ...
        mutual_mimo_newton_diag(M, N, H0, Rt_sqrt, gamma(k), lambda, ...
                                NSAMPLE, MAXITER, MAXLINES, epsilon);
    x_opt(:,k) = lambda_new(:,cnt2(k));
    Q_opt(:,:,k) = U*diag(x_opt(:,k))*U';
    fval_opt(k) = fval2(cnt2(k));
    fgap_opt(k) = fgap2(cnt2(k));
  end

else
% solve for symmetric Q
  for k=1:length(gamma)
    disp(['  SNR = ', num2str(SNR_db(k)), 'dB']);
    [x_new, fval, fgap, cnt_temp, ls_cnt_new] = ...
        mutual_mimo_newton_symm_choose(M, N, H0, Rt_sqrt, gamma(k), Q, NSAMPLE, ...
                                       MAXITER, MAXLINES, epsilon, ind1, ind2, ...
                                       ind3, barr_t, mu, MAXBARR);

    x_opt(:,k) = x_new(:,cnt_temp(length(cnt_temp)));
    Q_opt(:,:,k) = form_Q_symm(x_opt(:,k), N, Nr, Nx, ind3);
    fval_opt(k) = fval(length(fval));
    fgap_opt(k) = fgap(length(fgap));
  end

end