fading channel simulation with Monte-Carlo method

I'm working on a fading simulation in Matlab based on the work of Nikolay Kostov :

http://www.radioeng.cz/fulltexts/2003/03_04_12_16.pdf

but I find different results in the experimental BER from the theoretical one while performing simulations for M-PAM and M-QAM modulations. the first function generates the fade coefficients to simulate the fading channel ( Rayleigh / Rician ), the second one is the calculation of the BER through the Monte-Carlo method for a M-PAM modulation.

I think the problem stays in the fade coefficients, I tried to change 'sigma' to values less than 1 and I get some results, but since this
parameter is strictly bond to 'K' ( the Rician coefficient ) these settings aren't good for every 'K' I put in the "fading" function. I don't want to use the
'rayleighchan' and 'ricianchan' functions in matlab since my goal is to transpose these code in Python to work with the GNURadio deployment kit.
I hope you can help me. Thx in advance :D.

I'm working with MATLAB R2009b

#################################################################################################
% Generate the fade coefficients

function r = rice_fading_sigma(Kdb, N, Mi,sigma)
K = 10^(Kdb/10);
const = 1/(2*(K+1));
x = sigma*randn(1,N);
y = sigma*randn(1,N);
r = sqrt(const*((x + sqrt(2*K)).^2 + y.^2));
rt = zeros(1,Mi*length(r));
ki = 1;
for i=1:length(r)
rt(ki:i*Mi) = r(i);
ki = ki+Mi;
end
r = rt;

##################################################################################################
% The Monte-Carlo simulation

function BER = pam_rice(logEbNo,M,Kdb,sigma)

k = log2(M);
% M is the number of symbols
EbNolin = 10.^(logEbNo/10);
Nit = 100000;
Ns 0000;
Tstop = 100;
k_lin = 10^(Kdb/10);
ber_th = berfading(logEbNo,'pam',M,1,k_lin)
% Plot theoretical results.
figure;
semilogy(logEbNo,ber_th,'b');
xlabel('E_b/N_0 (dB)');
ylabel('Bit Error Rate');
grid on;
drawnow;

snr = logEbNo+10*log10(k);

for it_snr = 1:length(EbNolin)
bit_err = 0;
for it = 1:Nit
a = randint(Ns,1,M);
s = pammod(a,M,0,'gray');
r = rice_fading_sigma(Kdb,Ns,1,sigma);
v_r = s.*r';
v = awgn(v_r,snr(it_snr),'measured');
z = pamdemod(v,M,0,'gray');
errors = biterr(a,z);
bit_err = bit_err + errors;
if bit_err >= Tstop
break
end
end
BER(it_snr) = bit_err/(it*Ns*k);
end
hold on
semilogy(logEbNo, BER,'r-','LineWidth',1);
legend('Theoretical BER','Empirical BER');
title('Comparing Theoretical and Empirical Error Rates');
hold off;
%grid;

#################################################################################################### Start a New Thread