diff --git a/Functions/Theory/calcFWM/FWM_products.m b/Functions/Theory/calcFWM/FWM_products.m new file mode 100644 index 0000000..5429d0e --- /dev/null +++ b/Functions/Theory/calcFWM/FWM_products.m @@ -0,0 +1,69 @@ + + + +w0 = [1290:2:1290+15*2]'; +w0 = [1290:2:1290+15*2]'; +w0 = [ 1302 1304 1306 1308]'; +%w0 = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]'; +m = 0.5*(numel(w0)^3 - numel(w0)^2); +a = [1,1,1]'; + +w = w0; + + + + +for o = 2:3 + + p = nchoosek(w,3); + + q = []; + parfor i = 1:size(p,1) + + q_ = perms(p(i,:)); + q = [q;q_]; + + end + p = q; + %p = unique(q,"rows"); + w_ = p(:,1) + p(:,2) - p(:,3); + a_ = ones(size(w_)).* 1/o; + + w = [w ; w_]; + a = [a ; a_]; +% w = unique(w); +% a = unique(a); + + m(end+1) = 0.5*(numel(w)^3 - numel(w)^2); + +end + +figure(11) +hold on + +lambda = min(w):max(w); + +gen = sum(lambda == w,1); +gen(gen==0) = NaN; +stem(lambda,gen,"filled",'LineWidth',1,'Marker','o','MarkerSize',2,'LineStyle',':') + +initial = sum(lambda == w0,1); +initial(initial==0) = NaN; +stem(lambda,initial,"filled",'LineWidth',1.5,'MarkerSize',5,'Marker','^'); + +xlabel('Wavelength'); +ylabel('number of FWM products'); + +grid minor +legend('Generated Products', 'Initial Channel Position') +AxesMain = gca; +fig = gcf; +fontsize(AxesMain,8,"points") + +fig.Units = "centimeters"; +fig.Position = [2 2 8.5 7]; + + + + + diff --git a/Functions/Theory/calcFWM/JLT_statistics_laser_and_zdw.m b/Functions/Theory/calcFWM/JLT_statistics_laser_and_zdw.m new file mode 100644 index 0000000..79d052f --- /dev/null +++ b/Functions/Theory/calcFWM/JLT_statistics_laser_and_zdw.m @@ -0,0 +1,43 @@ +%% Laser Offset Statistics + + figure +for i = 1 + res = 1.7e6; + n_chann = 16; + df_T_exact = (-n_chann/2+0.5:n_chann/2).* 200e9; + + for key = 1:100 + laser_frequency_imperfection(key,:) = res .* round(randn(1,n_chann)*i*100); + df_T(key,:) = df_T_exact + laser_frequency_imperfection(key,:); + end + + hold on + histogram(laser_frequency_imperfection.*1e-6,100,"Normalization","probability","EdgeColor","none","FaceAlpha",0.3,'DisplayName',['Std. Dev.: ',num2str(mean(std(laser_frequency_imperfection))*1e-6),' MHz']); + xlabel('Laser Frequency Offset in MHz'); + ylabel('Probability'); + title(['Laser deviations from exact grid.']) + +end + +%% ZDW Statistics + +for k = 1:1000 + + % Set parameters + meanUniformMin = 1309; + meanUniformMax = 1315; + meanValue = 1310; + sigma = 2; + + % Seed the random number generator (assuming Mersenne Twister) + rng(k); + + % Generate normally distributed random numbers + randomNumbers(k,:) = normrnd(meanValue, sigma, [n_chann, 1]).'; + +end + +figure; +histogram(randomNumbers,100,"Normalization","probability","EdgeColor","none"); +xlabel('ZDW in nm'); +ylabel('Probability'); \ No newline at end of file diff --git a/Functions/Theory/calcFWM/analytical_calculation_paper.m b/Functions/Theory/calcFWM/analytical_calculation_paper.m new file mode 100644 index 0000000..d088424 --- /dev/null +++ b/Functions/Theory/calcFWM/analytical_calculation_paper.m @@ -0,0 +1,65 @@ +%%FWM analysis from "Analytical Calculation of the Number of +%%Four-Wave-Mixing Products in Optical Multichannel Communication Systems" + +N_ = [4,8,16]; + +df_hz = 400e9; +center_nm = 1310; + +figure() +for i = 1:length(N_) + + vec = (2*N_(i)-1:-1:2-N_(i)) -(N_(i)/2+0.5); + channelplan_hz = nm2hz(center_nm) + (vec * df_hz) ; + channelplan_nm = hz2nm(channelplan_hz); + + [Mndg,Mdg] = getProducts(N_(i)); + total(i) = sum(Mndg) + sum(Mdg); + subplot(1,length(N_),i) + xline(calcWavelengthPlan(N_(i), df_hz, center_nm)); + hold on + stem(channelplan_nm,(Mndg+Mdg),'filled','LineWidth',1,'Marker','o','MarkerSize',2) + stem(channelplan_nm,(Mdg),'filled','LineWidth',1,'Marker','none'); + ylim([0,100]) + xlim([1260, 1365]); + grid off + xlabel('O-band wavelength region in nm'); + ylabel('Number of FWM products'); + title([num2str(N_(i)),' ch.']) + +end + + +function [Mndg,Mdg] = getProducts(N) + + s = abs(2-N-1) ; + + for n = 2-N:2*N-1 + + if n<-N + Mdg(n+s) = 0; + elseif (-N <= n)&&(n <= 0) + Mdg(n+s) = N - ceil((N-n)/2); + elseif (1 <= n)&&(n <= N) + Mdg(n+s) = N - 1 - floor(n/2) - ceil((N-n)/2); + elseif (N < n)&&(n <= 2*N) + Mdg(n+s) = N - floor(n/2); + elseif n > 2*N + Mdg(n+s) = 0; + end + + if n<-N + Mndg(n+s) = 0; + elseif (-N <= n)&&(n < 1) + Mndg(n+s) = ceil((N^2 + n^2 - 2*N - 2*n + 2*N*n)/4); + elseif (1 <= n)&&(n <= N) + Mndg(n+s) = ceil(((N^2 - 6*N - 2*n^2 + 2*n + 4)/4) + floor((N*n)/2)); + elseif (N < n)&&(n <= 2*N) + Mndg(n+s) = floor(N^2 + n^2 /4 - N*n); + elseif n > 2*N + Mndg(n+s) = 0; + end + + end +end + diff --git a/Functions/Theory/calcFWM/calcFwmEfficiency.m b/Functions/Theory/calcFWM/calcFwmEfficiency.m new file mode 100644 index 0000000..b1a3560 --- /dev/null +++ b/Functions/Theory/calcFWM/calcFwmEfficiency.m @@ -0,0 +1,23 @@ +function [eta, deltaBeta] = calcFwmEfficiency(f_i, f_j, f_k, f_0, Ds, alphaDbPerKm, Lkm) +% FWM efficiency including phase mismatch and attenuation. +% alphaDbPerKm is the power attenuation in dB/km, Lkm is the fiber length in km. + +deltaBeta = calcPhaseMatching(f_i, f_j, f_k, f_0, Ds); + +alphaNpPerM = alphaDbPerKm .* log(10) ./ 10 ./ 1e3; +Lm = Lkm .* 1e3; + +denominator = alphaNpPerM.^2 + deltaBeta.^2; +term1 = alphaNpPerM.^2 ./ denominator; + +loss_term = 1 - exp(-alphaNpPerM .* Lm); + +if abs(alphaNpPerM) < eps + term2 = 4 .* sin(deltaBeta .* Lm ./ 2).^2 ./ max((alphaNpPerM .* Lm).^2, eps); +else + term2 = 1 + 4 .* exp(-alphaNpPerM .* Lm) .* sin(deltaBeta .* Lm ./ 2).^2 ./ (loss_term.^2); +end + +eta = term1 .* term2; + +end diff --git a/Functions/Theory/calcFWM/calcFwmPower.m b/Functions/Theory/calcFWM/calcFwmPower.m new file mode 100644 index 0000000..a6b65b7 --- /dev/null +++ b/Functions/Theory/calcFWM/calcFwmPower.m @@ -0,0 +1,48 @@ +function [P_fwm, eta, deltaBeta, Leff] = calcFwmPower( ... + f_i, f_j, f_k, f_0, Ds, alphaDbPerKm, Lkm, ... + P_i, P_j, P_k, gammaWInvKmInv, degeneracyFactor) +% Calculate FWM power for a fiber with attenuation and phase mismatch. +% +% Inputs: +% f_i, f_j, f_k, f_0 : frequencies in Hz +% Ds : dispersion slope in ps / (nm^2 km) +% alphaDbPerKm : attenuation in dB/km +% Lkm : fiber length in km +% P_i, P_j, P_k : launch powers in W +% gammaWInvKmInv : nonlinear coefficient in 1/(W km) +% degeneracyFactor : typically 3 for degenerate FWM, 6 for non-degenerate + +if nargin < 8 || isempty(P_i) + P_i = 1; +end +if nargin < 9 || isempty(P_j) + P_j = P_i; +end +if nargin < 10 || isempty(P_k) + P_k = 1; +end +if nargin < 11 || isempty(gammaWInvKmInv) + gammaWInvKmInv = 1; +end +if nargin < 12 || isempty(degeneracyFactor) + degeneracyFactor = 1; +end + +[eta, deltaBeta] = calcFwmEfficiency(f_i, f_j, f_k, f_0, Ds, alphaDbPerKm, Lkm); + +alphaNpPerM = alphaDbPerKm .* log(10) ./ 10 ./ 1e3; +Lm = Lkm .* 1e3; +gammaWInvMInv = gammaWInvKmInv ./ 1e3; + +if abs(alphaNpPerM) < eps + Leff = Lm; +else + Leff = (1 - exp(-alphaNpPerM .* Lm)) ./ alphaNpPerM; +end + +P_fwm = degeneracyFactor .* eta .* ... + (gammaWInvMInv .* Leff).^2 .* ... + P_i .* P_j .* P_k .* ... + exp(-alphaNpPerM .* Lm); + +end diff --git a/Functions/Theory/calcFWM/calcPhaseMatching.m b/Functions/Theory/calcFWM/calcPhaseMatching.m new file mode 100644 index 0000000..4c2fd4f --- /dev/null +++ b/Functions/Theory/calcFWM/calcPhaseMatching.m @@ -0,0 +1,20 @@ +function deltaBeta = calcPhaseMatching(f_i, f_j, f_k, f_0, Ds) +% Approximate phase mismatch for degenerate FWM close to the ZDW. +% Inputs are frequencies in Hz. +% f_i, f_j : pump frequencies (equal in the degenerate case) +% f_k : signal frequency +% f_0 : zero-dispersion frequency +% Ds : dispersion slope in ps / (nm^2 km) + +c = physconst('LightSpeed'); + +pump_frequency = 0.5 .* (f_i + f_j); +lambda_zdw_m = c ./ f_0; +dispersion_slope_si = Ds .* 1e3; + +deltaBeta = -(2 .* pi .* lambda_zdw_m.^4 ./ c.^2) .* ... + dispersion_slope_si .* ... + (pump_frequency - f_0) .* ... + (pump_frequency - f_k).^2; + +end diff --git a/Functions/Theory/calcFWM/scriptFWM.m b/Functions/Theory/calcFWM/scriptFWM.m new file mode 100644 index 0000000..abe3685 --- /dev/null +++ b/Functions/Theory/calcFWM/scriptFWM.m @@ -0,0 +1,104 @@ +clear; +clc; + +% Sweep the degenerate pump frequency around its nominal wavelength. +pump_detuning_hz = (-800:0.01:800) .* 1e9; + +% Degenerate FWM setup: two pump photons at f_p and one signal at f_s +% generate an idler at f_i = 2*f_p - f_s. +pump_wavelength_nm = 1310; +signal_wavelength_nm = 1308; +zdw_wavelength_nm = 1310; + +f_pump_nominal = wavelength2frequency(pump_wavelength_nm, 'nm'); +f_signal_scalar = wavelength2frequency(signal_wavelength_nm, 'nm'); +f_zdw_scalar = wavelength2frequency(zdw_wavelength_nm, 'nm'); + +f_pump = f_pump_nominal + pump_detuning_hz; +f_signal = f_signal_scalar .* ones(size(f_pump)); +f_zdw = f_zdw_scalar .* ones(size(f_pump)); +f_idler = 2 .* f_pump - f_signal; + +% Fiber parameters +dispersion_slope_ps_nm2_km = 0.07; +attenuation_db_per_km = 0.21; +fiber_length_km = 10; + +% Launch powers and nonlinear coefficient +pump_power_dbm = 10; +signal_power_dbm = 10; +pump_power_w = dbm2watt(pump_power_dbm); +signal_power_w = dbm2watt(signal_power_dbm); +gamma_w_inv_km_inv = 1.3; +degeneracy_factor = 3; + +[P_fwm, eta, delta_beta, L_eff_m] = calcFwmPower( ... + f_pump, f_pump, f_signal, f_zdw, ... + dispersion_slope_ps_nm2_km, attenuation_db_per_km, fiber_length_km, ... + pump_power_w, pump_power_w, signal_power_w, ... + gamma_w_inv_km_inv, degeneracy_factor); + +f_pump_thz = f_pump .* 1e-12; +f_zdw_thz = f_zdw_scalar .* 1e-12; +f_signal_thz = f_signal_scalar .* 1e-12; +idler_power_dbm = 10 .* log10(max(P_fwm, realmin) ./ 1e-3); + +figure; +tiledlayout(2,1); + +ax1 = nexttile; +plot(ax1, f_pump_thz, eta, 'LineWidth', 2); +hold(ax1, 'on'); +xline(ax1, f_zdw_thz, '--r', 'ZDW', 'LineWidth', 1.2, ... + 'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'bottom'); +xline(ax1, f_signal_thz, '--k', 'Signal', 'LineWidth', 1.2, ... + 'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'middle'); +ylabel(ax1, 'FWM efficiency'); +grid(ax1, 'on'); +title(ax1, 'FWM Efficiency and Idler Power versus Pump Frequency'); + +ax2 = nexttile; +plot(ax2, f_pump_thz, idler_power_dbm, 'LineWidth', 2); +hold(ax2, 'on'); +xline(ax2, f_zdw_thz, '--r', 'ZDW', 'LineWidth', 1.2, ... + 'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'bottom'); +xline(ax2, f_signal_thz, '--k', 'Signal', 'LineWidth', 1.2, ... + 'LabelOrientation', 'horizontal', 'LabelVerticalAlignment', 'middle'); +xlabel(ax2, 'Pump frequency (THz)'); +ylabel(ax2, 'FWM idler power (dBm)'); +grid(ax2, 'on'); + +fprintf('Pump wavelength : %.3f nm -> %.6f THz\n', ... + pump_wavelength_nm, f_pump_nominal .* 1e-12); +fprintf('Signal wavelength : %.3f nm -> %.6f THz\n', ... + signal_wavelength_nm, f_signal_scalar .* 1e-12); +fprintf('ZDW wavelength : %.3f nm -> %.6f THz\n', ... + zdw_wavelength_nm, f_zdw_scalar .* 1e-12); +fprintf('Pump launch power : %.2f dBm -> %.4g W\n', ... + pump_power_dbm, pump_power_w); +fprintf('Signal launch power : %.2f dBm -> %.4g W\n', ... + signal_power_dbm, signal_power_w); +fprintf('Peak FWM efficiency : %.4g\n', max(eta)); +fprintf('Peak FWM idler power : %.4g W (%.2f dBm)\n', ... + max(P_fwm), 10 .* log10(max(P_fwm) ./ 1e-3)); +fprintf('Effective fiber length : %.4f km\n', L_eff_m ./ 1e3); +fprintf('Idler wavelength range : %.3f nm to %.3f nm\n', ... + min(frequency2wavelength(f_idler, 'nm')), max(frequency2wavelength(f_idler, 'nm'))); +fprintf('Max |delta beta| : %.4g 1/m\n', max(abs(delta_beta))); + +function wavelength = frequency2wavelength(frequency, outputUnit) + c = physconst('LightSpeed'); + wavelength = c ./ frequency; + + switch lower(outputUnit) + case 'm' + case 'nm' + wavelength = wavelength .* 1e9; + otherwise + error('Unsupported output unit "%s". Use "m" or "nm".', outputUnit); + end +end + +function power_w = dbm2watt(power_dbm) + power_w = 1e-3 .* 10.^(power_dbm ./ 10); +end diff --git a/Functions/Theory/calcFWM/validate_fwm.m b/Functions/Theory/calcFWM/validate_fwm.m new file mode 100644 index 0000000..928a232 --- /dev/null +++ b/Functions/Theory/calcFWM/validate_fwm.m @@ -0,0 +1,124 @@ +%% Validate the analytical FWM product count against a brute-force reference +% The paper counts channel combinations, not only unique output frequencies: +% non-degenerate: i < j, k ~= i, k ~= j, n = i + j - k +% degenerate: i == j, k ~= i, n = 2*i - k + +clear; +clc; + +N_values = [4, 8, 16]; +plot_N = 8; + +fprintf('Validating analytical FWM product count from Goebel and Hanik (2008)\n'); + +for idxN = 1:numel(N_values) + N = N_values(idxN); + fprintf('\nN = %d\n', N); + + [Mndg_ana, Mdg_ana] = getProducts(N); + [Mndg_brute, Mdg_brute, n_values] = getProductsBruteForce(N); + + diff_ndg = Mndg_ana - Mndg_brute; + diff_dg = Mdg_ana - Mdg_brute; + + fprintf(' Analytical total : %d\n', sum(Mndg_ana) + sum(Mdg_ana)); + fprintf(' Brute-force total : %d\n', sum(Mndg_brute) + sum(Mdg_brute)); + + if all(diff_ndg == 0) && all(diff_dg == 0) + fprintf(' Match : yes\n'); + else + fprintf(' Match : no\n'); + fprintf(' Non-degenerate diff: %s\n', mat2str(diff_ndg)); + fprintf(' Degenerate diff : %s\n', mat2str(diff_dg)); + end + + if N == plot_N + plotComparison(n_values, Mndg_ana, Mdg_ana, Mndg_brute, Mdg_brute, N); + end +end + +function [Mndg, Mdg, n_values] = getProductsBruteForce(N) + n_values = (2 - N):(2*N - 1); + Mndg = zeros(size(n_values)); + Mdg = zeros(size(n_values)); + + for idx = 1:numel(n_values) + n = n_values(idx); + + % Degenerate products: two identical pumps and one different channel. + for i = 1:N + k = 2*i - n; + if isValidChannel(k, N) && (k ~= i) + Mdg(idx) = Mdg(idx) + 1; + end + end + + % Non-degenerate products: unordered pump pair plus one third channel. + for i = 1:N + for j = (i + 1):N + k = i + j - n; + if isValidChannel(k, N) && (k ~= i) && (k ~= j) + Mndg(idx) = Mndg(idx) + 1; + end + end + end + end +end + +function tf = isValidChannel(channel_idx, N) + tf = (channel_idx >= 1) && (channel_idx <= N) && (channel_idx == round(channel_idx)); +end + +function plotComparison(n_values, Mndg_ana, Mdg_ana, Mndg_brute, Mdg_brute, N) + figure; + + subplot(1,2,1); + stem(n_values, Mndg_ana + Mdg_ana, 'filled', 'LineWidth', 1, 'Marker', 'o', 'MarkerSize', 2); + hold on; + stem(n_values, Mdg_ana, 'filled', 'LineWidth', 1, 'Marker', 'none'); + title(['Analytical (N=', num2str(N), ')']); + xlabel('Product index n'); + ylabel('Number of FWM products'); + legend('Total', 'Degenerate'); + grid on; + + subplot(1,2,2); + stem(n_values, Mndg_brute + Mdg_brute, 'filled', 'LineWidth', 1, 'Marker', 'o', 'MarkerSize', 2); + hold on; + stem(n_values, Mdg_brute, 'filled', 'LineWidth', 1, 'Marker', 'none'); + title(['Brute force (N=', num2str(N), ')']); + xlabel('Product index n'); + ylabel('Number of FWM products'); + legend('Total', 'Degenerate'); + grid on; +end + +function [Mndg, Mdg] = getProducts(N) + s = abs(2 - N - 1); + + for n = 2 - N:2*N - 1 + if n < -N + Mdg(n + s) = 0; + elseif (-N <= n) && (n <= 0) + Mdg(n + s) = N - ceil((N - n)/2); + elseif (1 <= n) && (n <= N) + Mdg(n + s) = N - 1 - floor(n/2) - ceil((N - n)/2); + elseif (N < n) && (n <= 2*N) + Mdg(n + s) = N - floor(n/2); + elseif n > 2*N + Mdg(n + s) = 0; + end + + if n < -N + Mndg(n + s) = 0; + elseif (-N <= n) && (n < 1) + Mndg(n + s) = ceil((N^2 + n^2 - 2*N - 2*n + 2*N*n)/4); + elseif (1 <= n) && (n <= N) + Mndg(n + s) = ceil(((N^2 - 6*N - 2*n^2 + 2*n + 4)/4) + floor((N*n)/2)); + elseif (N < n) && (n <= 2*N) + Mndg(n + s) = floor(N^2 + n^2/4 - N*n); + elseif n > 2*N + Mndg(n + s) = 0; + end + end +end diff --git a/Functions/Theory/calcFWM/wavelength2frequency.m b/Functions/Theory/calcFWM/wavelength2frequency.m new file mode 100644 index 0000000..fe0f1eb --- /dev/null +++ b/Functions/Theory/calcFWM/wavelength2frequency.m @@ -0,0 +1,18 @@ +function frequency = wavelength2frequency(wavelength, inputUnit) +% Convert wavelength to optical frequency. +% Supported units: m, nm. + +c = physconst('LightSpeed'); + +switch lower(inputUnit) + case 'm' + wavelength_m = wavelength; + case 'nm' + wavelength_m = wavelength .* 1e-9; + otherwise + error('Unsupported input unit "%s". Use "m" or "nm".', inputUnit); +end + +frequency = c ./ wavelength_m; + +end diff --git a/Libs/mutual information rate TUM/LICENSE.txt b/Functions/Theory/mutual information rate TUM/LICENSE.txt similarity index 100% rename from Libs/mutual information rate TUM/LICENSE.txt rename to Functions/Theory/mutual information rate TUM/LICENSE.txt diff --git a/Libs/mutual information rate TUM/README.md b/Functions/Theory/mutual information rate TUM/README.md similarity index 100% rename from Libs/mutual information rate TUM/README.md rename to Functions/Theory/mutual information rate TUM/README.md diff --git a/Libs/mutual information rate TUM/air.m b/Functions/Theory/mutual information rate TUM/air.m similarity index 100% rename from Libs/mutual information rate TUM/air.m rename to Functions/Theory/mutual information rate TUM/air.m diff --git a/Libs/mutual information rate TUM/example_mi_cg_complex_16qam.m b/Functions/Theory/mutual information rate TUM/example_mi_cg_complex_16qam.m similarity index 100% rename from Libs/mutual information rate TUM/example_mi_cg_complex_16qam.m rename to Functions/Theory/mutual information rate TUM/example_mi_cg_complex_16qam.m diff --git a/Libs/mutual information rate TUM/example_mi_cg_real_2D_awgn_vs_phasenoise.m b/Functions/Theory/mutual information rate TUM/example_mi_cg_real_2D_awgn_vs_phasenoise.m similarity index 100% rename from Libs/mutual information rate TUM/example_mi_cg_real_2D_awgn_vs_phasenoise.m rename to Functions/Theory/mutual information rate TUM/example_mi_cg_real_2D_awgn_vs_phasenoise.m diff --git a/Libs/mutual information rate TUM/example_silas.m b/Functions/Theory/mutual information rate TUM/example_silas.m similarity index 100% rename from Libs/mutual information rate TUM/example_silas.m rename to Functions/Theory/mutual information rate TUM/example_silas.m diff --git a/Libs/mutual information rate TUM/mi_cg.m b/Functions/Theory/mutual information rate TUM/mi_cg.m similarity index 100% rename from Libs/mutual information rate TUM/mi_cg.m rename to Functions/Theory/mutual information rate TUM/mi_cg.m diff --git a/Libs/mutual information rate TUM/silas_example_mi_cg_pam.m b/Functions/Theory/mutual information rate TUM/silas_example_mi_cg_pam.m similarity index 100% rename from Libs/mutual information rate TUM/silas_example_mi_cg_pam.m rename to Functions/Theory/mutual information rate TUM/silas_example_mi_cg_pam.m