MATLAB.rar_cdama_nakagami_nakagami m

%*******************************************************************************% % % % ------- Discontinuous Nakagami-m simulator for arbitrary m ------ % % % % By Giuseppe Thadeu Freitas de Abreu, 2006 (All rights reserved) % % Center for Wireless Communnications, University of Oulu % % % %*******************************************************************************% function [Naka,f1,f2] = DiscontinuousNakagami(m,Omega,Nsp,Nfr) % This program generates piecewise-continuous Nakagami-m variates specified by: % OBS: If f1 and f2 are called, envelope and phase distributions are also generated % ------------------------------------------------------- % % m : Vector of m-parameters % Omega : Vector of Omega-parameters % Nsp : Number of samples % Tch : Coherence time of the channel % Tsb : Symbol duration (sampling period) % Nrd : Number of randoms used in the sum of sinusoids % % ------------------------------------------------------- % Examples: % 1) To generate 100 frames of 50 samples of a Nakagami-m chanel % with m = 1.35 and unitary power, type: % < Naka = PiecewiseContinuousNakagami(1.35,1,50,100); > % 2) To see also the distributions of the channel described, type, e.g.: % < [Naka,f1,f2] = PiecewiseContinuousNakagami(1.35,1,500,10000); > % OBS: Accurate distributions plots require a substantial amount of samples! M = length(m); mL = floor(2*m)/2; mU = mL + 0.5; p = 2*mL.*(mU-m)./m; testm = (floor(mL) == mL); for m_index = 1:M, if testm(m_index), mI(m_index) = mL(m_index); mH(m_index) = mU(m_index); else mI(m_index) = mU(m_index); mH(m_index) = mI(m_index); end end % Standard deviation of variates with interger m sigmaI = sqrt(Omega./(2*mI)); Naka = []; % Standard deviation of variates with half-interger m for m_index = 1:M, if mH(m_index) == 0.5, sigmaH(m_index) = sqrt(Omega(m_index)./(2*mH(m_index))) else sigmaH(m_index) = sqrt(Omega(m_index)./(floor(mH(m_index))+ceil(mH(m_index)))); end Naka_m = []; for frame = 1:Nfr, testp(m_index) = (rand <= p(m_index)); if (testm(m_index) & testp(m_index))|(~testm(m_index) & ~testp(m_index)), % Selection of variates with integer m % Variates with integer m X = sigmaI(m_index)*randn(mI(m_index),Nsp); Y = sigmaI(m_index)*randn(mI(m_index),Nsp); SignX = (2*randint-1); SignY = (2*randint-1); NakaFr = (SignX*sqrt(sum(X.^2,1)) + j*SignY*sqrt(sum(Y.^2,1))).*exp(j*randint(1,1,[0 3])*pi/2); else % Variates with half-integer m if mH == 0.5, % Case wehere m = 0.5 (Real gaussian variates) NakaFr = sigmaH(m_index)*randn(1,Nsp).*exp(j*randint(1,1,[0 3])*pi/2); else % Case wehere m > 0.5 (Gamma variates) X = sigmaH(m_index)*randn(ceil(mH(m_index)),Nsp); Y = sigmaH(m_index)*randn(ceil(mH(m_index))-1,Nsp); SignX = (2*randint-1); SignY = (2*randint-1); NakaFr = (SignX*sqrt(sum(X.^2,1)) + j*SignY*sqrt(sum(Y.^2,1))).*exp(j*randint(1,1,[0 3])*pi/2); end end Naka_m = [Naka_m NakaFr]; end Naka = [Naka; Naka_m]; end % --- Envelope Distributions --- if nargout >= 2, close all MarkerString = {'x','+','o','v','s','d','^','p'}; MarkerSize = [10 ,8 ,8 ,8 ,8 ,8 ,8 ,10 ]; Nbins = 49; figure(1) hold on grid on box on % Ploting auxiliary data for Legend LegendLabels = {'Exact'}; plot(-10,-10,'k','LineWidth',1.5); for m_index = 1:M, plot(-10,-10,[char(MarkerString(m_index)) 'k-'],'MarkerSize',MarkerSize(m_index),'MarkerFaceColor','w','LineWidth',1); LegendLabels(m_index+1) = {['m = ' num2str(m(m_index)) '; \Omega = ' num2str(Omega(m_index))]}; end for m_index = 1:M, % Histogram of Nakagami envelope distribution NakaEnv = abs(Naka(m_index,:)); [Henv,Xenv,Lenv,Uenv(m_index)] = Histogram(NakaEnv,Nbins); [XnormEnv,YnormEnv(m_index)] = NakaEnvelopePeak(m(m_index),Omega(m_index)); Henv = Henv*YnormEnv(m_index)/max(Henv); % Plotting simulated and theoretical distributions HtheoEnv = NakagamiEnvelopeDistribution(m(m_index),Omega(m_index),Xenv); plot(Xenv,Henv,[char(MarkerString(m_index)) 'k-'],'MarkerSize',MarkerSize(m_index),'MarkerFaceColor','w','LineWidth',1); plot(Xenv,HtheoEnv,'k','LineWidth',1.5); end top = ceil(10*max(YnormEnv))/10; axis([0 max(Uenv) 0 top]); h=get(gcf,'CurrentAxes'); set(h,'FontSize',12,'XTick',[0:0.5:max(Uenv)],'YTick',[0:0.1:top]); xlabel('Envelope: |Z|','FontSize',12); ylabel('p_{|Z|}(z)','FontSize',12); title('Nakagami-{\it m} Envelope Distribution','FontSize',14,'pos',[3,1.01*top]); h = legend(LegendLabels,1); set(h,'FontSize',12); FileName = ['NakagamiEnvDistribution']; print('-deps2',[FileName '.eps']); f1=figure(1); saveas(f1,[FileName '.fig']); end % --- Envelope Distributions --- if nargout == 3, TtheoPhase = (0:0.0005:0.5)*pi*4; for m_index = 1:M, RtheoPhase(m_index,:) = NakagamiPhaseDistribution(m(m_index),TtheoPhase); RtheoPhaseL = NakagamiPhaseDistribution(mL(m_index),TtheoPhase); RtheoPhaseU = NakagamiPhaseDistribution(mU(m_index),TtheoPhase); RmixPhase(m_index,:) = RtheoPhaseL*p(m_index) + (1-p(m_index))*RtheoPhaseU; end XnormPhase = pi/4; % Peak for m > 1, valley for m <= 1 Nbins = 79; for m_index= 1:M, % Histogram of Nakagami phase distribution NakaPhase = angle(Naka(m_index,:)); [T,R] = rose(NakaPhase,Nbins); Tphase = T(2:4:end); Rphase(m_index,:) = R(2:4:end); r1 = floor((1-0.9)*pi*Nbins/2) + 1; r2 = floor((1-0.6)*pi*Nbins/2) + 1; YnormPhase = gamma(m(m_index))*(abs(sin(2*XnormPhase)).^(m(m_index)-1))/((2^m(m_index))*gamma(m(m_index)/2)^2); if m(m_index) >=1, Rphase(m_index,:) = Rphase(m_index,:)*YnormPhase/max(Rphase(m_index,r1:r2)); % Peak at pi/4 else if min(Rphase(m_index,r1:r2)) > 0, Rphase(m_index,:) = Rphase(m_index,:)*YnormPhase/min(Rphase(m_index,r1:r2)); % Valley at pi/4 end end end figure(2) mmpolar(TtheoPhase,RtheoPhase,'-k',TtheoPhase,RmixPhase,':k',Tphase,Rphase,'--k','RLimit',[0 0.4],'RTickAngle',37.5,'RTickOffset',0.07); annotation('arrow',[0.52,0.8],[0.52,0.895]); text(0.85,0.8,'m > 1','FontSize',14,'FontName','courier'); h=get(gcf,'CurrentAxes'); set(h,'FontSize',12,'FontName','courier'); xlabel(''); ylabel(''); title('Nakagami-{\itm} Phase pdf (Mixture Method)','FontSize',14,'pos',[0,1.11],'FontName','courier'); FileName = ['NakagamiPhaseDistribution']; print('-deps2',[FileName '.eps']); f2=figure(2); saveas(f2,[FileName '.fig']); end % -------------------------------------------------------------------------------- % % Below are functions used to construct the envelope distribution plots % % -------------------------------------------------------------------------------- function [Xpeak,Ypeak] = NakaEnvelopePeak(m,Omega) Xpeak = sqrt(Omega*(1-0.5/m)); if Xpeak == 0, Ypeak = (2/gamma(m))*(m/Omega)^m; else Ypeak = NakagamiEnvelopeDistribution(m,Omega,Xpeak); end % -------------------------------------------------------------------