diff --git a/docs/json.rst b/docs/json.rst index ddc0bdd06..e313d00b0 100644 --- a/docs/json.rst +++ b/docs/json.rst @@ -523,6 +523,10 @@ See ``delta_power_range_db`` for more explaination. | | | ``ggn_spectrally_separated`` (see eq. 21 | | | | from `arXiv:1710.02225 | | | | `_). | +| | | ``ggn_approx`` (see eq. 24-25 | +| | | from `jlt:9741324 | +| | | `_). | +---------------------------------------------+-----------+---------------------------------------------+ | ``dispersion_tolerance`` | (number) | Optional. Pure number. Tuning parameter for | | | | ggn model solution. Default value is 1. | diff --git a/gnpy/core/science_utils.py b/gnpy/core/science_utils.py index 612eb47b2..41d1d1d29 100644 --- a/gnpy/core/science_utils.py +++ b/gnpy/core/science_utils.py @@ -12,7 +12,7 @@ from numpy import interp, pi, zeros, cos, array, append, ones, exp, arange, sqrt, trapz, arcsinh, clip, abs, sum, \ concatenate, flip, outer, inner, transpose, max, format_float_scientific, diag, sort, unique, argsort, cumprod, \ - polyfit + polyfit, log, reshape, swapaxes, full, nan from logging import getLogger from scipy.constants import k, h from scipy.interpolate import interp1d @@ -276,6 +276,7 @@ class NliSolver: List of implemented methods: 'gn_model_analytic': eq. 120 from arXiv:1209.0394 'ggn_spectrally_separated': eq. 21 from arXiv: 1710.02225 + 'ggn_approx': eq. 24-25 jlt:9741324 """ SPM_WEIGHT = (16.0 / 27.0) @@ -324,6 +325,28 @@ def compute_nli(spectral_info: SpectralInformation, srs: StimulatedRamanScatteri g_nli = sum(g_nli, 1) g_nli = interp(spectral_info.frequency, cut_frequency, g_nli) nli = spectral_info.baud_rate * g_nli # Local white noise + elif 'ggn_approx' in sim_params.nli_params.method: + if sim_params.nli_params.computed_channels is not None: + cut_indices = array(sim_params.nli_params.computed_channels) - 1 + elif sim_params.nli_params.computed_number_of_channels is not None: + nb_ch_computed = sim_params.nli_params.computed_number_of_channels + nb_ch = len(spectral_info.channel_number) + cut_indices = array([round(i * (nb_ch - 1) / (nb_ch_computed - 1)) for i in range(0, nb_ch_computed)]) + else: + cut_indices = array(spectral_info.channel_number) - 1 + + eta = NliSolver._ggn_approx(cut_indices, spectral_info, fiber, srs) + + # Interpolation over the channels not indicated as computed channels in simulation parameters + cut_power = outer(spectral_info.signal[cut_indices], ones(spectral_info.number_of_channels)) + cut_frequency = spectral_info.frequency[cut_indices] + pump_power = outer(ones(cut_indices.size), spectral_info.signal) + cut_baud_rate = outer(spectral_info.baud_rate[cut_indices], ones(spectral_info.number_of_channels)) + + g_nli = eta * cut_power * pump_power ** 2 / cut_baud_rate + g_nli = sum(g_nli, 1) + g_nli = interp(spectral_info.frequency, cut_frequency, g_nli) + nli = spectral_info.baud_rate * g_nli # Local white noise else: raise ValueError(f'Method {sim_params.nli_params.method} not implemented.') @@ -526,6 +549,89 @@ def _frequency_offset_threshold(beta2, symbol_rate): freq_offset_th = ((k_ref * delta_f_ref) * rs_ref * beta2_ref) / (beta2 * symbol_rate) return freq_offset_th + @staticmethod + def _ggn_approx(cut_indices, spectral_info: SpectralInformation, fiber, srs, spm_weight=SPM_WEIGHT, + xpm_weight=XPM_WEIGHT): + """Computes the nonlinear interference power evaluated at the fiber input. + The method uses eq. 24-25 of https://ieeexplore.ieee.org/document/9741324 + """ + # Spectral Features + nch = spectral_info.number_of_channels + frequency = spectral_info.frequency + baud_rate = spectral_info.baud_rate + slot_width = spectral_info.slot_width + roll_off = spectral_info.roll_off + df = spectral_info.df + diag(full(nch, nan)) + + # Physical fiber parameters + alpha = fiber.alpha(frequency) + beta2 = fiber.beta2(frequency) + gamma = outer(fiber.gamma(frequency[cut_indices]), ones(nch)) + + identity = diag(ones(nch)) + weight = spm_weight * identity + xpm_weight * (ones([nch, nch]) - identity) + weight = weight[cut_indices, :] + + dispersion_tolerance = sim_params.nli_params.dispersion_tolerance + phase_shift_tolerance = sim_params.nli_params.phase_shift_tolerance + max_slot_width = max(slot_width) + max_beta2 = max(abs(beta2)) + delta_z = sim_params.raman_params.result_spatial_resolution + + # Approximation psi + loss_profile = srs.loss_profile[:nch] + z = srs.z + psi = NliSolver._approx_psi(df=df, frequency=frequency, beta2=beta2, baud_rate=baud_rate, + loss_profile=loss_profile, z=z) + + # GGN for SPM + for cut_index in cut_indices: + dn = 0 + cut_frequency = frequency[cut_index] + cut_baud_rate = baud_rate[cut_index] + cut_roll_off = roll_off[cut_index] + cut_beta2 = beta2[cut_index] + cut_alpha = alpha[cut_index] + k_tol = dispersion_tolerance * abs(cut_alpha) + phi_tol = phase_shift_tolerance / delta_z + f_cut_resolution = min(k_tol, phi_tol) / abs(max_beta2) / (4 * pi ** 2 * (1 + dn) * max_slot_width) + f_pump_resolution = min(k_tol, phi_tol) / abs(max_beta2) / (4 * pi ** 2 * max_slot_width) + psi[cut_index, cut_index] = NliSolver._generalized_psi(cut_frequency, cut_frequency, cut_baud_rate, + cut_roll_off, cut_frequency, cut_baud_rate, + cut_roll_off, f_cut_resolution, f_pump_resolution, + srs, cut_alpha, cut_beta2, 0, cut_frequency) + psi = psi[cut_indices, :] + cut_baud_rate = outer(baud_rate[cut_indices], ones(nch)) + pump_baud_rate = outer(ones(cut_indices.size), baud_rate) + + eta_cut_central_frequency = \ + gamma ** 2 * weight * psi / (cut_baud_rate * pump_baud_rate ** 2) + eta = cut_baud_rate * eta_cut_central_frequency # Local white noise + + return eta + + @staticmethod + def _approx_psi(df, frequency, baud_rate, beta2, loss_profile, z): + """Computes the approximated psi function similarly to the one used in the GN model. + The method uses eq. 25 of https://ieeexplore.ieee.org/document/9741324""" + pump_baud_rate = outer(ones(frequency.size), baud_rate) + cut_beta = outer(beta2, ones(frequency.size)) + pump_beta = outer(ones(frequency.size), beta2) + delta_z = abs(z[:-1] - z[1:]) + + loss_lin = log(loss_profile) + pump_alpha = (loss_lin[:, 1:] - loss_lin[:, :-1]) / delta_z + leff = abs((loss_profile[:, 1:] - loss_profile[:, :-1]) / sqrt(abs(pump_alpha))) * pump_alpha / abs(pump_alpha) + leff = reshape(outer(leff, ones(z.size - 1)), newshape=[leff.shape[0], leff.shape[1], leff.shape[1]]) + leff2 = leff * swapaxes(leff, 2, 1) + leff2 = sum(leff2, axis=(1, 2)) + z_int = outer(ones(frequency.size), leff2) + + delta_beta = (cut_beta + pump_beta) / 2 + psi = z_int * pump_baud_rate / (4 * pi * abs(delta_beta * df)) + return psi + + def estimate_nf_model(type_variety, gain_min, gain_max, nf_min, nf_max): if nf_min < -10: diff --git a/tests/test_science_utils.py b/tests/test_science_utils.py index fdf0a80af..355cfa8be 100644 --- a/tests/test_science_utils.py +++ b/tests/test_science_utils.py @@ -10,6 +10,7 @@ from pandas import read_csv from numpy.testing import assert_allclose from numpy import array +from copy import deepcopy import pytest from gnpy.core.info import create_input_spectral_information, create_arbitrary_spectral_information @@ -132,3 +133,26 @@ def test_fiber_lumped_losses_srs(set_sim_params): stimulated_raman_scattering = RamanSolver.calculate_attenuation_profile(spectral_info_input, fiber) power_profile = stimulated_raman_scattering.power_profile assert_allclose(power_profile, expected_power_profile, rtol=1e-3) + + +@pytest.mark.usefixtures('set_sim_params') +def test_nli_solver(): + """Test the accuracy of NLI solver.""" + fiber = Fiber(**load_json(TEST_DIR / 'data' / 'test_science_utils_fiber_config.json')) + fiber.ref_pch_in_dbm = 0.0 + # fix grid spectral information generation + spectral_info_input = create_input_spectral_information(f_min=191.3e12, f_max=196.1e12, roll_off=0.15, + baud_rate=32e9, spacing=50e9, tx_osnr=40.0, + tx_power=1e-3) + + SimParams.set_params(load_json(TEST_DIR / 'data' / 'sim_params.json')) + sim_params = SimParams() + + spectral_info_output = fiber(deepcopy(spectral_info_input)) + + sim_params.nli_params.method = 'ggn_approx' + sim_params.raman_params.result_spatial_resolution = 10e3 + spectral_info_output_approx = fiber(deepcopy(spectral_info_input)) + + assert_allclose(spectral_info_output.nli, spectral_info_output_approx.nli, rtol=1e-1) +