MATLAB Code
%% ==============================
% FSO Simulation: Weak / Moderate / Strong Turbulence
% Toolbox-Free, Memory-Safe
% Nbits up to 1e6
%% ==============================
clear; clc; close all;
%% ------------------------------
% Parameters
%% ------------------------------
Nbits = 1e6; % number of transmitted bits
Pt = 1; % transmitted power (normalized)
sigmaR2_weak = 0.3; % Rytov variance for weak turbulence
sigmaR2_mod = 1; % moderate
sigmaR2_str = 10; % strong
noiseVar = 0.01; % AWGN variance
% FSO link parameters
L = 1000; % link distance in meters
lambda = 1550e-9; % laser wavelength (m)
theta = 1e-3; % beam divergence (rad)
Ar = 1e-4; % receiver aperture (m^2)
alpha_att = 0.2/1000; % atmospheric attenuation (1/m)
% Path loss (geometric + atmospheric)
PL = (Ar/(pi*(theta*L)^2)) * exp(-alpha_att*L);
fprintf('Path loss factor = %.4e\n', PL);
%% ------------------------------
% Transmit signal (OOK)
%% ------------------------------
bits = randi([0 1], Nbits, 1);
tx = Pt * bits;
%% ------------------------------
% Weak Turbulence (Log-Normal)
%% ------------------------------
% Log-amplitude variance
sigma2 = sigmaR2_weak;
% Generate 1D log-normal intensity
X = sqrt(sigma2) * randn(Nbits,1) - sigma2/2;
I_weak = exp(X);
% Received signal
rx_weak = tx .* I_weak + sqrt(noiseVar)*randn(Nbits,1);
% BER (threshold = 0.5*Pt)
rx_bits = rx_weak > 0.5*Pt;
BER_weak = sum(bits ~= rx_bits)/Nbits;
fprintf('Weak turbulence BER = %.4e\n', BER_weak);
%% ------------------------------
% Moderate Turbulence (Gamma-Gamma)
%% ------------------------------
sigmaR2 = sigmaR2_mod;
% Approx alpha and beta (small integers for memory)
alpha = max(1, round((exp((0.49*sigmaR2)/(1+1.11*sigmaR2^(12/5))^(7/6))-1)^(-1)));
beta = max(1, round((exp((0.51*sigmaR2)/(1+0.69*sigmaR2^(12/5))^(5/6))-1)^(-1)));
% Generate Gamma RVs (toolbox-free, 1D)
IL = mean(-log(rand(alpha,Nbits)),1); % Large-scale
IS = mean(-log(rand(beta,Nbits)),1); % Small-scale
I_mod = IL .* IS; % Gamma-Gamma intensity
% Received signal
rx_mod = tx .* I_mod' + sqrt(noiseVar)*randn(Nbits,1); % I_mod' to match column vector
% BER
rx_bits = rx_mod > 0.5*Pt;
BER_mod = sum(bits ~= rx_bits)/Nbits;
fprintf('Moderate turbulence BER = %.4e\n', BER_mod);
%% ------------------------------
% Strong Turbulence (Negative Exponential)
%% ------------------------------
% Exponential intensity
I_strong = exprnd(1,Nbits,1);
% Received signal
rx_strong = tx .* I_strong * PL + sqrt(noiseVar)*randn(Nbits,1);
% BER
rx_bits = rx_strong > 0.5*Pt;
BER_strong = sum(bits ~= rx_bits)/Nbits;
fprintf('Strong turbulence BER = %.4e\n', BER_strong);
%% ------------------------------
% PDF Comparison Plot
%% ------------------------------
edges = linspace(0,5,300);
figure;
hold on;
histogram(I_weak, edges,'Normalization','pdf','DisplayStyle','stairs');
histogram(I_mod, edges,'Normalization','pdf','DisplayStyle','stairs');
histogram(I_strong,edges,'Normalization','pdf','DisplayStyle','stairs');
hold off;
legend('Weak (Log-Normal)','Moderate (Gamma-Gamma)','Strong (Exponential)');
xlabel('Intensity'); ylabel('PDF');
title('FSO Turbulence Regimes (Toolbox-Free, 1D)');
grid on;
Output
