0
$\begingroup$

I simply want to visualize 【Capacity bounds and mutual information】 versus 【SNR $\gamma$(dB) or probability $p$ (The input to the channel is a binomial distribution with parameter $p$)】

But I cannot reproduce the results in a paper (On the Capacity of FSO Links over Gamma-Gamma Atmospheric Turbulence Channels Using OOK Signaling) for some reason.

Is there some mistake in my Mathematica code?

\[Gamma] = 10;(*\[Gamma]=P_{opt}/\sqrt{N_o W} is the SNR \
definition*)Clear["Global`*"];
p = 1/2;
\[Xi] = 1;
(*\[Xi] represents the square of the increment in Euclidean distance \
due to the use of a pulse shape of high PAOPR(peak-to-average optical \
power ratio),alternative to the classical rectangular pulse.*)
\[Theta] = 0.9028;
\[Kappa] = 2;(*with 0<\[Theta]<1,representing the fact that the \
channel under study is constrained to \[Kappa]=2 W T_b degrees of \
freedom*)\[Alpha] = 4;
\[CapitalTheta] = \[Alpha]^(3 \[Alpha]/
       2)*(1 - \[Alpha])^((1 - \[Alpha])/
      2)*(1 - \[Alpha]/2)^(1 - \[Alpha]/2);
\[Beta] = 4;
(*\[Alpha] and \[Beta] are related to the atmospheric conditions*)
SI = 1/\[Alpha] + 1/\[Beta] + 
  1/(\[Alpha]*\[Beta]);(*scintillation index,a parameter of interest \
used to describe the strength of atmospheric fading*)Capacity = 
 MeijerG[{{1, 
      1, (1 - \[Alpha])/2, (2 - \[Alpha])/2, (1 - \[Beta])/
       2, (2 - \[Beta])/2}, {}}, {{1}, {0}}, 
    16*(1/p - 
       1)*\[Xi]*\[Theta]*\[Kappa]*\[Gamma]^2/(\[Alpha]^2*\[Beta]^2)]*(\
2^(\[Alpha] + \[Beta])/
      Log[E, 2]/(2*Pi*Gamma[\[Alpha]]*Gamma[\[Beta]])) /. {\[Gamma] ->
     10*10^\[Gamma]};
Subscript[C, "nt"] = 
  1/2*Log[2, 
     1 + (1/p - 1)*\[Kappa]*\[Xi]*\[Theta]*\[Gamma]^2] /. {\[Gamma] ->
      10*10^\[Gamma]};
(*nt means:no atmospheric turbulence*)
Subscript[C, "AWGN"] = 
  1/2*Log[2, 1 + 2*\[Xi]*\[Theta]*\[Kappa]*\[Gamma]^2] /. {\[Gamma] ->
      10*10^\[Gamma]};
Subscript[C, "H"] = 
  Log[2, (Sqrt[\[Xi]*\[Theta]*\[Kappa]]*\[Gamma] + 2)*
     Sqrt[E/(2 Pi)]] /. {\[Gamma] -> 10*10^\[Gamma]};
Subscript[C, "F"] = 
  Log[2, Sqrt[
       E^2/(4 Pi)]^\[Alpha]*(Sqrt[\[Xi]*\[Theta]*\[Kappa]]*\[Gamma])^\
\[Alpha]*1/\[CapitalTheta]] /. {\[Gamma] -> 10*10^\[Gamma]};
Subsuperscript[C, "AWGN", "turb"] = 
 MeijerG[{{1, 
      1, (1 - \[Alpha])/2, (2 - \[Alpha])/2, (1 - \[Beta])/
       2, (2 - \[Beta])/2}, {}}, {{1}, {0}}, 
    32*\[Xi]*\[Theta]*\[Kappa]*\[Gamma]^2/(\[Alpha]^2*\[Beta]^2)]*(2^(\
\[Alpha] + \[Beta] - 2)/
      Log[E, 2]/(2*Pi*Gamma[\[Alpha]]*Gamma[\[Beta]])) /. {\[Gamma] ->
     10*10^\[Gamma]}; plot1 = 
 Plot[{Capacity, Subscript[C, "nt"], Subscript[C, "AWGN"], 
   Subscript[C, "H"], Subscript[C, "F"], 
   Subsuperscript[C, "AWGN", "turb"]}, {\[Gamma], 0.01, 50}, 
  PlotLegends -> "Expressions", PlotRange -> Full, 
  AxesLabel -> {"SNR \[Gamma](in dB)", "Capacity(bits/channel use)"}];
Subscript[f, "I"] = 
  2*(\[Alpha]*\[Beta])^((\[Alpha] + \[Beta])/2)/(Gamma[\[Alpha]]*
      Gamma[\[Beta]])*i^((\[Alpha] + \[Beta])/2 - 1)*
   BesselK[\[Alpha] - \[Beta], 2 Sqrt[\[Alpha]*\[Beta]*i]];
CGammap[\[Gamma]_] := 
  NIntegrate[
   1/2*Log[2, (1 + (1/p - 1)*\[Kappa]*\[Xi]*\[Theta]*\[Gamma]^2*i^2)]*
     Subscript[f, "I"] /. {\[Gamma] -> 10*10^\[Gamma]}, {i, 0, 
    Infinity}];
plot2 = ListPlot[{#, CGammap@#} & /@ Range[5, 50, 0.1], 
   PlotLegends -> "anotherBound"];
Show[plot1, plot2]

enter image description here

Clear["Global`*"];
\[Xi] = 1;
\[Theta] = 0.9028;
\[Kappa] = 2;

fYwhenxequals1[y_, p_, i_, \[Gamma]_] := 
  1/Sqrt[2*Pi]*
   Exp[-(y - 1/p*\[Gamma]*Sqrt[\[Xi]*\[Theta]*\[Kappa]]*i^2)^2/2];
fYwhenxequals0[y_] := 1/Sqrt[2*Pi]*Exp[-y^2/2];

MutualInformation[i_, p_, \[Gamma]_] := 
  Module[{denominator}, 
   denominator[y_] := (1 - p)*fYwhenxequals0[y] + 
     p*fYwhenxequals1[y, p, i, \[Gamma]];
   (1 - p)*
     NIntegrate[
      fYwhenxequals0[y]*
       Log[2, fYwhenxequals0[y]/denominator[y]], {y, -65, 109}] + 
    p*NIntegrate[
      fYwhenxequals1[y, p, i, \[Gamma]]*
       Log[2, fYwhenxequals1[y, p, i, \[Gamma]]/
         denominator[y]], {y, -65, 109}]];
listplotWhenSNRequalsMinus5dB = 
  ListPlot[{#, MutualInformation[1, #, 10^(-5/10)]} & /@ 
    Range[0, 1, 0.00625]];
listplotWhenSNRequalsMinus10dB = 
  ListPlot[{#, MutualInformation[1, #, 10^(-10/10)]} & /@ 
    Range[0, 1, 0.00625]];
listplotWhenSNRequalsMinus1dB = 
  ListPlot[{#, MutualInformation[1, #, 10^(-1/10)]} & /@ 
    Range[0, 1, 0.00625]];
Show[listplotWhenSNRequalsMinus5dB, listplotWhenSNRequalsMinus10dB, \
listplotWhenSNRequalsMinus15dB, 
 PlotLegends -> {"SNR is -5 dB", "SNR is -10 dB", "SNR is -1 dB"}]

enter image description here

Clear["Global`*"];
GammaGammaDistributionPDF[\[Alpha]_, \[Beta]_] := 
  2*(\[Alpha]*\[Beta])^((\[Alpha] + \[Beta])/2)/(Gamma[\[Alpha]]*
      Gamma[\[Beta]])*i^((\[Alpha] + \[Beta])/2 - 1)*
   BesselK[\[Alpha] - \[Beta], 2 Sqrt[\[Alpha]*\[Beta]*i]];
Plot[{GammaGammaDistributionPDF[4, 4], 
  GammaGammaDistributionPDF[4, 1], 
  GammaGammaDistributionPDF[1, 1]}, {i, 0.001, 10}, PlotRange -> Full,
  PlotLabel -> 
  Row[{"Gamma-Gamma distribution for different \[Alpha] and \[Beta]", 
    ""}], AxesLabel -> {"i", Row[{Subscript["f", "I"], "(i)"}]}, 
 PlotLegends -> "Expressions"]


(*For plane wave,here sigmaR^2 is the variance of laser \
scintillation (Rytov variance?),and its magnitude determines the strength of \
atmospheric turbulence.The larger its value,the stronger the \
turbulence intensity.*)
GammaGammaDistributionPDFForPlaneWave[sigmaR_] := 
 GammaGammaDistributionPDF[
  1/(Exp[(0.49 sigmaR^2)/((1 + 1.11 sigmaR^(12/5))^(7/6))] - 1), 
  1/(Exp[(0.51 sigmaR^2)/((1 + 0.69 sigmaR^(12/5))^(5/6))] - 1)]
Plot[{GammaGammaDistributionPDFForPlaneWave[0.4], 
  GammaGammaDistributionPDFForPlaneWave[0.8], 
  GammaGammaDistributionPDFForPlaneWave[1.4]}, {i, 0.001, 6}, 
 PlotRange -> Full, 
 PlotLabel -> 
  Row[{"Gamma-Gamma distribution for different ", 
    Subsuperscript[\[Sigma], "R", ""], ""}], 
 AxesLabel -> {"i", Row[{Subscript["f", "I"], "(i)"}]}, 
 PlotLegends -> "Expressions"]

enter image description here

enter image description here


By the way,

LogNormal model for weak turbulence

LogNormalDistributionForWeakTurbulence[MeanOfLogorithmLightIntensity_,
   StandardDeviationOfLogorithmLightIntensity_] :=
 Module[{MeanOfChi, VarianceOfChi, StandardDeviationOfChi}, 
  MeanOfChi = MeanOfLogorithmLightIntensity/2; 
  VarianceOfChi = StandardDeviationOfLogorithmLightIntensity^2/4;
  StandardDeviationOfChi = Sqrt[VarianceOfChi];
  PDF[LogNormalDistribution[MeanOfChi, StandardDeviationOfChi], i]]

Plot[Table[
   LogNormalDistributionForWeakTurbulence[0.05, 
    StandardDeviationOfLogorithmLightIntensity], \
{StandardDeviationOfLogorithmLightIntensity, {0.2, 0.4, 0.6, 0.8} // 
     Sqrt}] // Evaluate, {i, 0, 4}, Filling -> Axis, 
 PlotRange -> Full, PlotLegends -> {0.2, 0.4, 0.6, 0.8}, 
 PlotLabel -> 
  Row[{"LogNormal distribution for different ", 
    Subsuperscript[\[Sigma], "l", "2"], ""}]]

enter image description here

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.