clc
clear all
close all
format long;
NMC = 1e6;
Om2 = 1; TH = 10.^(0./10);
M2 = 3;N = 3; Po = 0.8^1; %Po = sqrt(0.7);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 产生随机数%%%%%%%%%%
% ym = randraw('nakagami', [M2, Om2], N, NMC ); ym = ym';
% yc = randraw('nakagami', [M2, Om2], N, NMC ); yc = yc';
% ym = sqrt(Om2./2./M2).*chi2rnd(2*M2, NMC,N);
% C1 = randn(NMC,N); C2 = randn(NMC,N);
% % yc =sqrt(VAR(1))./sqrt(2).*complex(C1,C2);
% yc =sqrt(Om2)./sqrt(2).*complex(C1,C2);
% var(ym)
% var((yc))
% ym1= Po.*ym + sqrt(1-Po.^2).*yc;
% A=corrcoef(ym1,ym)
% ym = abs(ym); ym1 = abs(ym1);
% Ym = ym.*ym;
% Ym1 = ym1.*ym1;
% A=corrcoef(ym1,ym)
Ym= (gamrnd(M2,Om2./M2,NMC,N)); VAR = var(Ym)./1;
% Yc=gamrnd(M2,Om2./M2,NMC,N);
% Yc = randraw('nakagami', [M2, Om2], N, NMC ); Yc = Yc';
C1 = randn(NMC,N); C2 = randn(NMC,N);
% Yc = sqrt(Om2)./sqrt(2.*1).*complex(C1,C2);
Yc =sqrt(VAR(1))./sqrt(2).*complex(C1,C2);
Ym1=sqrt(Po).*(Ym).^(0.5) + sqrt(1-Po).*(Yc).^(0.5);
mean(Ym)
mean(Yc)
var(Ym)
var(Yc)
var(Yc.*Yc)
% Ym1= Po.*Ym + sqrt(1-Po.^2).*Yc;
A=corrcoef(Ym1,Ym)
% Ym = abs(Ym); Ym1 = abs(Ym1);
A=corrcoef(Ym1,Ym)
Ym1 = Ym1.^2;
for j = 1:NMC
[YM(j),nYM] = max((Ym(j,:))); % 得到最大信道的序号nYM
YM1(j) = Ym1(j,nYM); % 根据被选中的过时信道,得到实时信道
end
% for j = 1:NMC
% [YM1(j),nYM] = max((Ym1(j,:))); % 得到最大信道的序号nYM
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 随机数PDF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
step = 0.05;XXX = 0:step:5; %I=0.01:0.01:5;
hDD = hist(YM1,XXX);
NCorrGamma_pdf = hDD/(step*sum(hDD)); % 仿真得到统计意义上的 Gamama PDF
NCorrGammaRNCDF = cumsum(hDD)/NMC;
% NGammaRNPDF = N.*(GammaRNCDF).^(N-1).*Gamma_pdf; % N个数中最大的数的PDF
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 理论 PDF %%%%%%%%%%%%%%
syms xx yy zz;
for iii = 1:1:length(XXX)
xxx = XXX(iii);
PDF_XY = @(yy) (M2./Om2).^(M2+1).*(xxx.*yy./Po).^((M2-1)./2)./((1-Po).*gamma(M2)).*exp(-(xxx+yy)./(1-Po).*M2./Om2)...
.*mfun('BesselI',M2-1, 2.*M2./Om2./(1-Po).*sqrt(Po.*xxx.*yy));
PDF_G = @(yy) (M2./Om2).^M2./gamma(M2).*yy.^(M2-1).*exp(-M2./Om2.*yy); % 得到数学表达式上的 Gamama PDF
CDF_G = @(yy) 1 - (mfun('GAMMA',M2, M2./Om2.*yy)./gamma(M2)); % 得到数学表达式上的 Gamama CDF
PDF_NG = @(yy) N.*PDF_G(yy).*(CDF_G(yy)).^(N-1);
PDF_XdY = @(yy) (M2./Om2)./(1-Po).*(xxx./yy./Po).^((M2-1)./2).*exp(-(xxx+Po.*yy)./(1-Po)./1.*M2./Om2)...
.*mfun('BesselI',M2-1, 2.*M2./Om2./(1-Po).*sqrt(Po.*xxx.*yy));
NCorrGammaPDF(iii) = quadgk(@(yy) PDF_XY(yy)./PDF_G(yy).*PDF_NG(yy),0,1e2);
NCorrGammaPDF1(iii) = quadgk(@(yy) PDF_XdY(yy).*PDF_NG(yy), 0,1e2);
end
figure
plot(XXX,NCorrGamma_pdf,'ro-');
hold on
plot(XXX,NCorrGammaPDF,' bp')
plot(XXX,NCorrGammaPDF1,' gd')
xlabel('x')
ylabel('Gamma pdf, p(I)')
grid;
评论0