function Pb=awgn_theo(M,K,rho_in_db)
rho=10^(rho_in_db/10);
part=0;
t=2^(K-1)/(2^K-1);
for i=2:M
    comb=factorial(M)/(factorial(i)*factorial(M-i));
    part=part+(-1)^i*comb*exp(-K*rho*(1-1/i));
end
Ps=part/M;
Pb=t*Ps;