why this codes showing Undefined function 'wf' for input arguments of type 'double'. [v P_wf]=wf(Pt,4,N_wf); and sometimes faulting the value. the code is not running?

clc;

clear all;

close all;

%% >>>>>>>>>>>>>>>>> ... 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, 1)))

outage11(j1)=1;

else

outage11(j1)=count11_norm(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=sqrt(SNR_i)*((randn(1,1000)+1j*randn(1,1000))/sqrt(2));

% Channel matrix with elements as complex normal r.v. having var=1/2 per dim.

H=reshape(H_vec,4, []);

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';

[v, P_wf]=wf(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, 1)))

outage44(j1)=1;

else

outage44(j1)=count44_norm(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, 1)))

outage44_wf(j1)=1;

else

outage44_wf(j1)=count44_norm_wf(abs(range44_wf-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, 1)))

outage14(j1)=1;

else

outage14(j1)=count14_norm(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, 1)))

outage41(j1)=1;

else

outage41(j1)=count41_norm(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','1X1 SISO','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','1X1 SISO','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','1X1 SISO','Location','SouthEast');

title('Outage probability vs SNR for 4 bps/Hz');

xlabel('SNR in dB');

ylabel('1-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');