%這個程序由rayleigh 的方法產生由Nakagami 分布的隨機變量；
function aa=newnakagamirnd(m,n,fdt)

%m為nakagami 衰落系數
%n為nakagami樣本個數
%fdt為多普勒頻移fd和時間t的乘積
if m==0.65
    a1=-0.0828;
    a2=-0.45634;
    a3=-15.8891;
    b1=63.1955;
    b2=23.2981;
elseif m==0.75
    a1=-0.0547;
    a2=-0.3679;
    a3=-1.0336;
    b1=6.2107;
    b2=1.8533;
elseif m==0.85
    a1=-0.0336;
    a2=-0.1543;
    a3=-0.4733;
    b1=4.9250;
    b2=1.2082;
elseif m==1.0
    a1=0;
    a2=0;
    a3=0;
    b1=0;
    b2=0;
elseif m==1.5
    a1=0.0993;
    a2=0.0560;
    a3=0.2565;
    b1=0.5276;
    b2=0.0770;
elseif m==2.0
    a1=0.1890;
    a2=-0.0128;
    a3=0.2808;
    b1=-0.0809;
    b2=0.0638;
elseif m==3.0
    a1=0.3472;
    a2=-0.2145;
    a3=0.2626;
    b1=-0.6779;
    b2=0.1690;
elseif m==4.0
    a1=0.4846;
    a2=-0.4231;
    a3=0.2642;
    b1=-0.9729;
    b2=0.2727;
elseif m==5.0
    a1=0.6023;
    a2=-0.6238;
    a3=0.2789;
    b1=-1.1798;
    b2=0.3732;
elseif m==6.0
    a1=0.7139;
    a2=-0.8305;
    a3=0.3223;
    b1=-1.3232;
    b2=0.4558;
elseif m==7.0
    a1=0.8167;
    a2=-1.0244;
    a3=0.3761;
    b1=-1.4233;
    b2=0.5192;
elseif m==8.0
    a1=0.9260;
    a2=-1.2350;
    a3=0.4557;
    b1=-1.4872;
    b2=0.5628;
elseif m==10.0
    a1=1.1088;
    a2=-1.6095;
    a3=0.6015;
    b1=-1.6046;
    b2=0.6488;
end

fil=besselj(0,(0:1:2*pi*fdt));
x1=randn(1,n);
y1=randn(1,n);


x22=conv(x1,fil);
x2=x22/sqrt(mean(x22.^2));
y22=conv(y1,fil);
y2=y22/sqrt(mean(y22.^2));
sita=atan(y2./x2);
r=sqrt(x2.^2+y2.^2);
u=1-exp((-r.^2)./2);
z=(sqrt(log(1./(1-u)))).^(1/m);
bb=z+(a1.*z+a2.*z.^2+a3.*z.^3)./(1+b1.*z+b2.*z.^2);
w=mean(bb.^2);
aa0=bb/sqrt(w);
aa1=aa0.*cos(sita);
aa2=aa0.*sin(sita);
aa3=aa1+sqrt(-1)*aa2;

for i=1:n
    aa(i)=aa3(randint(1,1,[1,n+length(fil)-1]));
end