% This script is used to evaluate Fig. 1b) in the paper "Adaptive Removal of Multipath Interference in Short Reach 112 GBd PAM-4 IM/DD Systems"

%% Parameters
df       = 1e6;        % Laser linewidth [Hz]
SIR_dB   = 20;                % Interference attenuation [dB]
alpha    = 10^(-SIR_dB/20);   % Interference attenuation [linear]
n_fiber  = 1.467;             % Refractive index
c        = physconst('lightspeed'); % [m/s]

L        = linspace(0,250,50); % Interference delay [m]  
tau      = n_fiber./c.*L;       % Interference time (= tau) [s]

tau_c = 1/(pi*df);              % laser coherence time [s]
L_c   = (c/n_fiber)*tau_c;      % laser coherence length [m]

var_sat = 2*alpha^2;            % Analytical saturation of variance

%% Monte–Carlo Simulation
fs               = 100e9;            % sampling rate [Hz]
Tsim             = 50e-6;            % sim duration [s]
N                = round(Tsim*fs);   % number of samples for each realization
max_delay_samples  = round(max(tau)*fs); % largest delay that is evaluated (based on max. Interference delay)
phase_noise_std   = sqrt(2*pi*df/fs); % standard dev. phase noise

num_realizations = 50;             % number of parallel runs
monte_carlo_variance = zeros(num_realizations, length(L));
parfor r = 1:num_realizations

    % generate a realization of phase noise random walk
    dphi = phase_noise_std * randn(1, N + max_delay_samples); % matlab randn process has std = 1
    phi  = cumsum(dphi);
    phi_direct  = phi(max_delay_samples+1 : max_delay_samples+N);
    var_k = zeros(1, length(L));
    for t = 1:length(tau)

        nd = round( tau(t)*fs ); % delay in samples for current interference time
        phi_delayed = phi(max_delay_samples+1-nd : max_delay_samples+N-nd); %cut out interfering signal part (was earlier)

        E = exp(1j*phi_direct) + alpha*exp(1j*phi_delayed); % E-fields combined
        I = abs(E).^2; % photo current as magnitude square of E-field
        var_k(t) = var(I);

    end
    monte_carlo_variance(r, :) = var_k;
end

avg_of_mc_variances = mean(monte_carlo_variance, 1);
std_of_mc_variances = std(monte_carlo_variance, 0, 1);

%% Analytic variance
L_        = linspace(0,250,500); % Interference delay [m]  
tau_      = n_fiber./c.*L_;
analytic_variance = 2*alpha^2 * (1 - exp(-2*pi*df.*tau_)).^2;

%% Plot 
cols =     [0.3467    0.5360    0.6907
            0.9153    0.2816    0.2878
            0.4416    0.7490    0.4322];

coherence_length_multiples = 0.5:0.5:ceil(L(end)/L_c);

figure();
hold on;
plot(L, avg_of_mc_variances, 'LineWidth',2, 'DisplayName','Simulation','Color',cols(1,:),'LineStyle','-');
errorbar(L, avg_of_mc_variances,std_of_mc_variances, 'LineWidth',0.7,'LineStyle','none', 'DisplayName','Simulation','Color',cols(1,:),'HandleVisibility','off');

plot(L_, analytic_variance, 'LineWidth',2, 'DisplayName','Analytic','Color',cols(2,:),'LineStyle','-');
xticks(coherence_length_multiples.*L_c);
xticklabels(round(coherence_length_multiples.*L_c,1));

norm_to_coherence_len = 1;
if norm_to_coherence_len
    xticklabels(coherence_length_multiples);
    xlabel('$n \cdot L_c$', 'FontSize',12);
else
    xlabel('Interference Delay [m]', 'FontSize',12);
end

xline(L_c.*coherence_length_multiples, 'LineWidth',1.5, 'DisplayName','Coh. Length','HandleVisibility','off','Color',[0.7,0.7,0.7],'LineStyle','-');
xlim([0,L(end)]);
yline(var_sat, '-.k','LineWidth',1.5, 'DisplayName','Saturation: 2$\alpha ^2$');
grid on;
ylabel('Intensity Variance', 'FontSize',12);
title(sprintf('MPI Variance; %d MHz; SIR: %d dB',df.*1e-6,SIR_dB), 'FontSize',14);
legend('Location','southeast');