function zd1690249673022886508
% 关于nakagami-m分布matlab仿真的问题
% 举报违规检举侵权投诉|今天 00:03 提问者悬赏:100分 | wuyueyi625 | 分类:其他编程语言 | 浏览18次
%
% 我想出这个图 我需要做的就是做一个关于m值对网络性能影响的仿真 现在我在网上找到了代码 但是我对matlab实在不是很会用
% http://wenku.baidu.com/link?url=lGx6IkPPyeyYWtVKleA3hGnT0nPR090q4Iva_FsTnMZ0PlgiXFi-VxWn9S6hu_bjKcaX2dO5sSIjYZwD8CjrnBdgL1Tf9tmaJMNfkuBAnA_
% 这个网站就是我找到的代码 但是直接放进matlab里无法运行出图 求大神指点 我的毕设就靠各位的力量了
% 1.用monte carlo仿真nakagami分布:
figure(1)
m = 1;
omaga = 1;
r = [0:0.1:2];
for n = 1:6000
[realamp] = nakagami(m, omaga);
amp(n) = realamp;
end
for n4 = 1:21
cnt1(n4) = 0;
end
for n1 = 1:21
for n2 = 1:6000
if round(amp(n2)*10) == round(r(n1)*10)
cnt1(n1) = cnt1(n1) +1;
end
end
end
for n3 = 1:21
cnt1(n3) = cnt1(n3)/600;
end
plot(r,cnt1,'b');
hold on;
m = 2;
omaga = 1;
r = [0:0.1:2];
for n = 1:6000
[realamp] = nakagami(m, omaga);
amp(n) = realamp;
end
for n4 = 1:21
cnt2(n4) = 0;
end
for n1 = 1:21
for n2 = 1:6000
if round(amp(n2)*10) == round(r(n1)*10)
cnt2(n1) = cnt2(n1) +1;
end
end
end
for n3 = 1:21
cnt2(n3) = cnt2(n3)/600;
end
plot(r,cnt2,'r');
xlabel('r');
ylabel('pdf');
legend('m=1','m=2');
hold off
% 2.理想的仿真: 先将m=1,2,带入PDF的公式,简化概率密度函数,然后画出理想图形:
figure(2)
x=0:0.01:2.0;
y1=2.*x.*exp(-x.*x);
plot(x,y1);
hold on;
y2=8.*(x.^3).*exp(-2.*x.*x);
plot(x,y2,'r');
xlabel('r');
ylabel('pdf');
legend('m=1','m=2');
% 对于模拟的图形和仿真的图形相比较,可以看出用monte carlo方法通过取随机数近似的
% nakagami概率密度函数大致形状是一样的,m=1时比m=2时有更大的拖尾。
function [realamp] = nakagami(m, omaga)
% Nakagami的主函数;
r = [0:0.1:2];
dt = 0.01;
gam = 0;
for n = 0:0.01:20
gam = gam + (n.^(m-1)).*exp(-n).*dt;
end
tosam = 0;
tocnt = 0;
for rr = 1:21
pdf(rr) = ((2*m^m)*((r(rr))^(2*m-1))*(exp(-(m/omaga)*(r(rr))^2)))/(gam*omaga^m);
sam(rr) = round(pdf(rr)*10);
tosam = tosam + sam(rr);
cnt(rr) = 0;
end
toamp = 0;
for n1 = 1:400
amp(n1) = round(rand*20);
for n2 = 1:21
if amp(n1)==round(r(n2)*10)
cnt(n2) = cnt(n2) + 1;
break;
end
end
if cnt(n2) > sam(n2)
cnt(n2) = cnt(n2) - 1;
amp(n1) = -1;
end
end
for n3 = 1:21
tocnt = tocnt + cnt(n3);
end
for n4 = 1:200
number = round(rand*400);
if number == 0
number = round(rand*400);
end
realamp = amp(number)/10;
if realamp ~= -0.1
break;
end
end
return