opt_mutual_mimo.m

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

%----------------------------------------------------------------------
% Optimization program to find the transmit covariance matrix that achieves
% the ergodic capacity of a MIMO channel given the channel mean and transmit
% correlation matrices. The program then calculates the mutual information
% using the covariance matrix that maximizes Jensen's upper bound on
% the average mutual information, and compares that with the capacity itself.
%
%----------------------------------------------------------------------
% Author: Mai Vu
% Date: Mar 02, 2004
% Modified: Apr 19, 2004
%
%----------------------------------------------------------------------
% INPUT: The channel mean and transmit correlation matrices
%        The program either takes inputs from channel_file, or generates the
%        matrices arbitrarily.
%
% OUTPUT: It saves the results in save_file, which includes: 
%         - The channel parameters: mean H0, transmit correlation Rt, SNR
%         - The optimal transmit covariance matrix: Q_opt, x_opt (the
%           vectorized form of Q_opt)
%         - The channel ergodic capacity: mi_opt
%         - Optimality gap (a bound on the gap between the numerical
%           solution and the true capacity): fgap_opt
%         - The mutual information obtained using the input covariance matrix
%           which is optimal for the Jensen bound: mutual_jen_fw
%         - The mutual information with iid equal power input: mutual_equal
%
%        The program also produces two graph files:
%         - graph_mi_file: plots the capacity and the other two mutual
%           information curves
%         - graph_eig_file: plots the input power distribution (i.e. the
%           eigenvalue of the transmit covariance matrix) for the three
%           cases: optimal, Jensen bound and equal power.
%
% PARAMETERS:
%         N: Number of transmit antennas
%         M: Number of receive antennas
%         NSAMPLE: number of channel samples used in each Monte-Carlo simulation
%         SNR_dB: input signal to noise ratio in dB
%         K: Rician K factor
%         random_data: 0 means generate random channel mean and correlation matrices
%                      1 means take input from channel_file
%         spec_case: 0 means channel has both mean and tx corr
%                    1 means channel has zero mean
%                    2 means channel has no transmit correlation
%         file_version: file version
%         channel_file: input mat file for channel parameters (.mat)
%         save_file: output mat file (.mat)
%         graph_mi_file: output graph file (.eps)
%         graph_eig_file: output graph file (.eps)
%----------------------------------------------------------------------  

clear all;
close all;
clc;

disp('Initializing...');

%--------------------PARAMETERS--------------------
% change these parameters to run different simulations
N = 4;  % number of Tx antennas
M = 4;  % number of Rx antennas
NSAMPLE = 10000;

SNR_dB = [-10:5:15]; % SNR in dB
K = 5;  % channel K factor

random_data = 0; % generate random channel mean and correlation matrices
spec_case = 0; % 1 is no mean, 2 is no tx corr

% setting up various file names
if spec_case == 1
  suffix = 'corr';
elseif spec_case == 2
  suffix = 'mean';
else
  suffix = 'both';
end

file_version = 'a';
channel_file = ['mutual_mimo_channel_4_4_', file_version];
save_file = ['opt_mutual_mimo_', suffix, '_', file_version];
graph_mi_file = ['mutual_mimo_compare_', suffix, '_', file_version, '.eps'];
graph_eig_file = ['eigval_cov_compare_', suffix, '_', file_version, '.eps'];

%--------------------PROGRAM--------------------
I = eye(N);
SNR = 10.^(SNR_dB/10);

% use the test data or create random channel statistics
if random_data
  Rt_half = randn(N,N) + j*randn(N,N);
  Rt = Rt_half*Rt_half';
  Rt = Rt/(K+1);
  
  H0 = randn(M,N) + j*randn(M,N);
  H0 = H0/norm(H0, 'fro')*sqrt(K*M*N/(K+1));  
  
else
  load(channel_file);
  [M N] = size(H0);
end

if (spec_case == 2)
  % check for mean only case
  Rt = I*real(trace(Rt))/N; % this is to retain the same K factor
elseif (spec_case == 1)
  % check for correlation only case
  Rt = Rt/trace(Rt)*N; % full correlation
  H0 = zeros(M,N);
end

Rt_sqrt = sqrtm(Rt);

% Jensen upper bound waterfill matrix
S0 = (H0'*H0 + M*Rt);
[U lambda] = eig(S0);
lambda = real(diag(lambda));

% mutual information at the initial input covariance matrix
% for k=1:length(SNR_db)
%   fval_init(k) = mutual_mimo_func(M,N,H0,Rt_sqrt,gamma(k),Q,NSAMPLE,I);
% end

disp('Optimizing the average mutual information...');
[x_opt, Q_opt, fval_opt, fgap_opt] = mutual_mimo_opt(H0, Rt, M, N, SNR_dB);

lambda_equal = ones(N,1)/N;  

disp('Calculating Jensen''s upper bound on the average mutual information...');
for k = 1:length(SNR_dB)
  disp(['  SNR = ', num2str(SNR_dB(k)), 'dB']);  
  lambda_jen_wf(:,k) = waterfill(SNR(k), lambda);
  mutual_jen_wf(k) = ...
      - mutual_mimo_func_diag(M,N,H0*U,U'*Rt_sqrt*U,...
                              SNR(k),lambda_jen_wf(:,k),NSAMPLE,I);
  mutual_equal(k) = ...
      - mutual_mimo_func_diag(M,N,H0,Rt_sqrt,...
                              SNR(k),lambda_equal,NSAMPLE,I);
  if (size(x_opt,1) > N)
    lambda_opt(:,k) = eig(Q_opt(:,:,k));
  else
    lambda_opt(:,k) = x_opt(:,k);
  end
end

mi_opt = -fval_opt;

%--------------------SAVE AND DISPLAY RESULTS--------------------
% save data and plot graphs
save(save_file, 'H0', 'Rt', 'SNR_dB', 'x_opt', 'Q_opt', 'mi_opt', 'fgap_opt', ...
     'lambda_jen_wf', 'mutual_jen_wf', 'mutual_equal');

disp('________________________');
disp('Optimization successful.');
disp('  The optimum covariance matrix is in Q_opt');
disp('  The optimum power allocation is in lambda_opt');
disp('  The ergodic capacity is in mi_opt');

figure;
plot(SNR_dB, mi_opt/log(2), 'r');
hold on;
plot(SNR_dB, mutual_jen_wf/log(2), '-.');
plot(SNR_dB, mutual_equal/log(2), 'g:');
xlabel('SNR in dB');
ylabel('Mutual information (bps/Hz)');
legend('Ergodic capacity', 'MI w/Jensen cov', 'MI w/equal power',2);
saveas(gcf, graph_mi_file, 'psc2');

figure;
fig_leg(:, 1) = plot(SNR_dB, lambda_opt, 'r');
hold on;
fig_leg(:, 2) = plot(SNR_dB, lambda_jen_wf, 'b--');
fig_leg(:, 3) = plot(SNR_dB, lambda_equal*ones(1, length(SNR_dB)), 'g:');
xlabel('SNR in dB');
ylabel('Eigenvalues of the transmit covariance matrix');
legend(fig_leg(1, :), 'Optimum power', 'Jensen power', 'Equal power',1);
saveas(gcf, graph_eig_file, 'psc2');