How to solve error?? Too many output arguments (water fill code)

%% Main function to calculate MIMO channel capacity and associated figures
%% of merit like complementary CDF & Outage probability
% Objective: (1)This function compares Mean channel capacity for different
% MIMO realizations (SISO, SIMO,MISO, MIMO)as a function of SNR
% (2) compares complementary CDF for different MIMO realizations
% (SISO, SIMO,MISO, MIMO)as a function of capacity in bps/Hz
% (3) compares Outage probability for different MIMO realizations
% (SISO, SIMO,MISO, MIMO)as a function of SNR
% Author: Raza Umar as part of EE 575 Information Theory Assignment
% Date: May 08, 2010
%% >>>>>>>>>>>>>>>>>>>>>>>> ... CLEANING ... <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
close all;
clear all;
clc;
%% >>>>>>>>>>>>>>>>> ... SIMULTAION PARAMS ... <<<<<<<<<<<<<<<<<<<<<<<<<<<<
SNR_dB=-10:30;
SNR=10.^(SNR_dB./10);
ch_realizations=10000; % Monte Carlo sim. of 10,000 channel realizations
c_outage=4;
epsilon=1e-6;
%% >>>>>>>>>>>>>>>>>>>>>>> ... SISO ... <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% initialization and param setting
c11_i=zeros(length(SNR),ch_realizations);
c11=zeros(length(SNR),1);
% matrix initializations to avoid matrix growing inside loop
% Mean channel capacity calculations
for i1=1:length(SNR)
SNR_i= SNR(i1); % for one specific value of SNR
h11=1/sqrt(2).*complex(randn(1,ch_realizations),randn(1,ch_realizations));
% complex normal r.v. with var=1/2 per dim.
h11_mag_sq=abs(h11).^2;
c11_i(i1,1:ch_realizations)=log2(1+SNR_i.*h11_mag_sq); % instantaneous cap.
c11(i1)=mean(c11_i(i1,:)); % mean capacity
end
%complementary CDF at SNR=10dB
c11_all=c11_i(21,:); % channel realizations at 10dB SNR
range_11=0:0.1:max(c11_all);
count11=histc(c11_all,range_11);
count11_norm=cumsum(count11)/max(cumsum(count11));
comp_CDF_11=1-count11_norm;
%outage probability
j1= 1;
for i1=find(SNR_dB==2):find(SNR_dB==20)
c11_all=abs(c11_i(i1,:));
range11=0:0.1:max(c11_all);
count11=histc(c11_all,range11);
count11_norm=cumsum(count11)/max(cumsum(count11));
if(isempty(find(abs(range11-c_outage)<epsilon)))
outage11(j1)=1;
else
outage11(j1)=count11_norm(find(abs(range11-c_outage)<epsilon));
end
j1=j1+1;
end
%% >>>>>>>>>>>>>>>>>>>>> ... MIMO (4x4) ... <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% initialization and param setting
Pt=1; % power budget for MIMO water filling
Nr=4;
Nt=4;
I_Nr=eye(Nr);
c44_i=zeros(length(SNR),ch_realizations);
c44=zeros(length(SNR),1);
% matrix initializations to avoid matrix growing inside loop
% Mean channel capacity calculations
for i1=1:length(SNR)
SNR_i= SNR(i1); % for one specific value of SNR
H_vec=1/sqrt(2).*complex(randn(Nr*Nt,ch_realizations),randn(Nr*Nt,ch_realizations));
% Channel matrix with elements as complex normal r.v. having var=1/2 perdim.
H=reshape(H_vec,Nr,Nt,ch_realizations);
for i2=1:ch_realizations
H_i=H(1:Nr,1:Nt,i2);
H_H_hermt=H_i*H_i';
arg=I_Nr+ (SNR_i/Nt).*H_H_hermt;
arg2= det(arg);
c44_i(i1,i2)=log2(arg2); % instantaneous cap.
%% alternate implementation to find inst. cap. (using eigen values)
% lamda=eig(H_H_hermt);
%arg3=(1+(SNR_i/Nt)*lamda(1))*(1+(SNR_i/Nt)*lamda(2))*(1+(SNR_i/Nt)*lamda(3))*(1+(SNR_i/Nt)*lamda(4));
% c442i(i1,i2)=log2(arg3); % instantaneous cap.
%% MIMO water filling
    lamda=eig(H_H_hermt);
    N=1./(SNR_i.*lamda);
    N_wf=N';
    pt=1;
    [v P_wf]=wfill(Pt,4,N_wf);
arg_wf=(1+P_wf(1)/N_wf(1)).*(1+P_wf(2)/N_wf(2)).*(1+P_wf(3)/N_wf(3)).*(1+P_wf(4)/N_wf(4));
c44_i_wf(i1,i2)=log2(arg_wf); % instantaneous cap. using wf
%%
end
c44(i1)=mean(c44_i(i1,:)); % mean capacity
c44_wf(i1)=mean(c44_i_wf(i1,:)); % mean capacity
end
%%
%%complementary CDF at SNR=10dB
c44_all=abs(c44_i(21,:)); % channel realizations at 10dB SNR
range_44=0:0.1:max(c44_all);
count44=histc(c44_all,range_44);
count44_norm=cumsum(count44)/max(cumsum(count44));
comp_CDF_44=1-count44_norm;
%%MIMO WF
c44_all_wf=abs(c44_i_wf(21,:)); % channel realizations at 10dB SNR
range_44_wf=0:0.1:max(c44_all_wf);
count44_wf=histc(c44_all_wf,range_44_wf);
count44_norm_wf=cumsum(count44_wf)/max(cumsum(count44_wf));
comp_CDF_44_wf=1-count44_norm_wf;
%%
%outage prob
j1= 1;
for i1=find(SNR_dB==2):find(SNR_dB==20)
c44_all=abs(c44_i(i1,:));
range44=0:0.1:max(c44_all);
count44=histc(c44_all,range44);
count44_norm=cumsum(count44)/max(cumsum(count44));
if(isempty(find(abs(range44-c_outage)<epsilon)))
outage44(j1)=1;
else
outage44(j1)=count44_norm(find(abs(range44-c_outage)<epsilon));
end
j1=j1+1;
end
%%MIMO WF
j1= 1;
for i1=find(SNR_dB==2):find(SNR_dB==20)
c44_all_wf=abs(c44_i_wf(i1,:));
range44_wf=0:0.1:max(c44_all_wf);
count44_wf=histc(c44_all_wf,range44_wf);
count44_norm_wf=cumsum(count44_wf)/max(cumsum(count44_wf));
if(isempty(find(abs(range44_wf-c_outage)<epsilon)))
outage44_wf(j1)=1;
else
outage44_wf(j1)=count44_norm_wf(find(abs(range44_wfc_outage)<epsilon));
end
j1=j1+1;
end
%% >>>>>>>>>>>>>>>>>>>>> ... SIMO (1x4) ... <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% initialization and param setting
Nr=4;
Nt=1;
c14_i=zeros(length(SNR),ch_realizations);
c14=zeros(length(SNR),1);
% matrix initializations to avoid matrix growing inside loop
% Mean channel capacity calculations
for i1=1:length(SNR)
    SNR_i= SNR(i1); % for one specific value of SNR
h14=1/sqrt(2).*complex(randn(Nr*Nt,ch_realizations),randn(Nr*Nt,ch_realizations));
% complex normal r.v. having var=1/2 per dim.
h14_mag_sq=abs(h14).^2;
h14_mag_sq_sum=sum(h14_mag_sq);
c14_i(i1,1:ch_realizations)=log2(1+SNR_i.*h14_mag_sq_sum); %
instantaneous cap.
c14(i1)=mean(c14_i(i1,:)); % mean capacity
end
%complementary CDF at SNR=10dB
c14_all=c14_i(21,:); % channel realizations at 10dB SNR
range_14=0:0.1:max(c14_all);
count14=histc(c14_all,range_14);
count14_norm=cumsum(count14)/max(cumsum(count14));
comp_CDF_14=1-count14_norm;
%outage probability
j1= 1;
for i1=find(SNR_dB==2):find(SNR_dB==20)
c14_all=abs(c14_i(i1,:));
range14=0:0.1:max(c14_all);
count14=histc(c14_all,range14);
count14_norm=cumsum(count14)/max(cumsum(count14));
if(isempty(find(abs(range14-c_outage)<epsilon)))
outage14(j1)=1;
else
outage14(j1)=count14_norm(find(abs(range14-c_outage)<epsilon));
end
j1=j1+1;
end
%% >>>>>>>>>>>>>>>>>>>>> ... MISO (4x1) ... <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% initialization and param setting
Nr=1;
Nt=4;
c41_i=zeros(length(SNR),ch_realizations);
c41=zeros(length(SNR),1);
% matrix initializations to avoid matrix growing inside loop
% Mean channel capacity calculations
for i1=1:length(SNR)
SNR_i= SNR(i1); % for one specific value of SNR
h41=1/sqrt(2).*complex(randn(Nr*Nt,ch_realizations),randn(Nr*Nt,ch_realizations));
% complex normal r.v. having var=1/2 per dim.
h41_mag_sq=abs(h41).^2;
h41_mag_sq_sum=sum(h41_mag_sq);
c41_i(i1,1:ch_realizations)=log2(1+(SNR_i./Nt).*h41_mag_sq_sum); %
instantaneous cap.
c41(i1)=mean(c41_i(i1,:)); % mean capacity
end
%complementary CDF at SNR=10dB
c41_all=c41_i(21,:); % channel realizations at 10dB SNR
range_41=0:0.1:max(c41_all);
count41=histc(c41_all,range_41);
count41_norm=cumsum(count41)/max(cumsum(count41));
comp_CDF_41=1-count41_norm;
%outage probability
j1= 1;
for i1=find(SNR_dB==2):find(SNR_dB==20)
c41_all=abs(c41_i(i1,:));
range41=0:0.1:max(c41_all);
count41=histc(c41_all,range41);
count41_norm=cumsum(count41)/max(cumsum(count41));
if(isempty(find(abs(range41-c_outage)<epsilon)))
outage41(j1)=1;
else
outage41(j1)=count41_norm(find(abs(range41-c_outage)<epsilon));
end
j1=j1+1;
end
%% >>>>>>>>>>>>>>>>>>>> ... Plotting ... <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% Mean channel capacity
figure('Name','Mean Capacity comparisons (vs SNR) for Flat Fading MIMO channels');
plot(SNR_dB,abs(c44_wf),'bx-',SNR_dB,abs(c44),'ro',SNR_dB,c41,'c+-',SNR_dB,c14,'gs-',SNR_dB,c11,'b:');
legend('4x4 MIMO WF','4x4 MIMO','4x1 MISO','1x4 SIMO','1X1SISO','Location','NorthWest');
title('Mean Capacity vs SNR');
xlabel('SNR (dB) --->');
ylabel('Mean Capacity bps/Hz --->');
axis([-10 30 0 40]);
%complementary CDF at SNR=10dB
figure('Name','Complementary CDF at SNR = 10dB');
plot(range_44_wf,comp_CDF_44_wf,'bx-',range_44,comp_CDF_44,'ro',range_41,comp_CDF_41,'c+-',range_14,comp_CDF_14,'gs-',range_11,comp_CDF_11,'b:');
legend('4x4 MIMO WF','4x4 MIMO','4x1 MISO','1x4 SIMO','1X1SISO','Location','NorthEast');
title('Complementary CDF comparisons (vs capacity) at SNR=10dB');
xlabel('Mean Capacity bps/Hz --->');
ylabel('1 - Outage Probability --->');
axis([0 15 0.9 1]);
%outage probability vs SNR for 4 bps/Hz
figure('Name','Outage probability comparisons (vs SNR) for Flat Fading Channels');
%plot(range_44,comp_CDF_44,'bx-',range_44,comp_CDF_44,'ro',range_41,comp_CDF_41,'c+-',range_14,comp_CDF_14,'gs-',range_11,comp_CDF_11,'b:');
semilogy(SNR_dB(find(SNR_dB==2):find(SNR_dB==20)),outage44_wf,'bx-',SNR_dB(find(SNR_dB==2):find(SNR_dB==20)),outage44,'ro',SNR_dB(find(SNR_dB==2):find(SNR_dB==20)),outage41,'c+-',SNR_dB(find(SNR_dB==2):find(SNR_dB==20)),outage14,'gs-',SNR_dB(find(SNR_dB==2):find(SNR_dB==20)),outage11,'b:');
legend('4x4 MIMO WF','4x4 MIMO','4x1 MISO','1x4 SIMO','1X1SISO','Location','SouthEast');
title('Outage probability vs SNR for 4 bps/Hz');
xlabel('SNR in dB --->');
ylabel('Outage Probability --->');
axis([2 20 1e-6 1]);
% MIMO water-filling capacity gain
figure('Name','capacity gain of water-filling');
plot(SNR_dB,abs(c44_wf)'-abs(c44));
grid on;
title('WF gain in capacity');
xlabel('SNR in dB --->');
ylabel('WF gain in bps/Hz --->');