% Figure 3.9: Achievable rate of MIMO receivers for a fast fading channel. 
% Elisabeth de Carvalho, April 2012

clear all


% number of channel realization
It = 10000;

% Number of transmit and receive antennas
M=4;

% SNR range in dB
% SNR = ratio between total transmit power and noise variance
SNRdBvalues =  [-10:2:30];

% initialize variables: maximal achievable rate for one SNR value and one channel realization
Cmf     = zeros(1,It);
Czf     = zeros(1,It);
Cmmse   = zeros(1,It);
Csic    = zeros(1,It);


SNRidx = 0;
for SNRdB=SNRdBvalues
    SNRdB
    SNRidx = SNRidx + 1;
    SNR=10^(SNRdB/10);
    
    % compute maximal achievable rate 
    for kk=1:It
        
        % generate channel realization
        H = (randn(M)+j*randn(M))/sqrt(2);
        
        % Matched filter
        for ii=1:M
            % output SNR for each stream
            hi  = H(:,ii);
            Hib = H; Hib(:,ii)=zeros(M,1);
            SNRi(ii)  = real( ( abs(hi'*hi)^2) / (hi'*Hib*Hib'*hi + hi'*hi/SNR*M) );
        end
        Cmf(kk) = sum ( log2( 1 + SNRi) );                  % maximal achievable rate
        
        
        %% Zero-Forcing
        Dzf     = inv(H'*H)/SNR*M;
        SNRzf   = real( diag(1./diag(Dzf)) );               % output SNR (diagonal matrix)
        Czf(kk) = log2( det(eye(M) + SNRzf) );              % maximal achievable rate
        
        
        %% MMSE
        Dmmse   = inv(H'*H + eye(M)/SNR*M);
        SNRmmse = real( diag(1./diag(Dmmse))*SNR/M -eye(M) );       % output SNR (diagonal matrix)
        Cmmse(kk) = real(log2(det(eye(M) + SNRmmse)));              % maximal achievable rate
        
        
        %% SIC
        Csic(kk) = real(log2(det(eye(M) + SNR/M*H*H')));
        
    end
    
    % average over all channel realizations for a given SNR value
    Mmf(SNRidx)     = mean(Cmf);
    Mzf(SNRidx)     = mean(Czf);
    Mmmse(SNRidx)   = mean(Cmmse);
    Msic(SNRidx)    = mean(Csic);
    
end

Mmf
Mzf
Mmmse
Msic

% plot
figure(1)
plot(SNRdBvalues, Mmf,'k')
hold on
plot(SNRdBvalues, Mzf,'b')
hold on
plot(SNRdBvalues, Mmmse,'m')
plot(SNRdBvalues, Msic,'r')
title('Achievable rate of MIMO receivers for i.i.d. fast fading channel')
xlabel('Average SNR \rho (dB)')
ylabel('Achievable rate (bits/transmission)')
legend('Matched Filter','Zero-Forcing', 'MMSE','SIC',2)
grid
hold off
