Utils#

compute_beta2(lamb: float, cd_coefficient: float, speed_of_light: float) → float[source]#

Compute the Chromatic Dispersion coefficient β₂ in ps²/km

Parameters:

lamb (float) – Wavelength (nm)
cd_coefficient (float) – Chromatic Dispersion coefficient (ps/nm/km)
speed_of_light (float) – Speed of light (m/s)
by (The formula is given)
:: (.. math) – beta_2 = -frac{10^3 cdot D cdot lambda^2}{2pi c}
WARNING:: (..) –
- All input values must be in the units specified above.
- Output β₂ is given in picosecond squared per kilometer (ps²/km).

Returns:

beta2 (float) – Group velocity dispersion β₂ (in ps²/km)

Example 1#

>>> lamb = 1550         # in nm
>>> cd_coefficient = 17 # in ps/nm/km
>>> c = 299792458       # in m/s
>>> beta2 = compute_beta2(lamb, cd_coefficient, c)
>>> print(beta2)
-21.682619391414896

Notes

The formula is derived from Eq. (4) and (5) of [Savory, 2008].

References

[1] Savory, Seb J. “Digital filters for coherent optical receivers.” Optics express 16.2 (2008): 804-817.
[2] Eghbali, Amir, et al. “Optimal least-squares FIR digital filters for compensation of chromatic dispersion in digital coherent optical receivers.” Journal of lightwave technology 32.8 (2014): 1449-1456.

apply_chromatic_dispersion(x: ndarray, z: float, beta2: float, alpha_dB: float | None = None, fs: float = 1, direction: int = 1) → ndarray[source]#

Apply chromatic dispersion effects in optical fiber communications.

Apply chromatic dispersion effects in the frequency domain for fiber-optic communication systems. It applies a dispersion-induced phase shift to the input signal in the frequency domain and considers signal attenuation [1].

Parameters:

x (numpy array) – Complex signal
z (float) – Step length in km.
beta2 (float) – coefficient in ps**2/km
alpha_dB (float, optional) – attenuation in dB/km (default: None)
fs (float, optional) – Sampling frequency in hertz (Hz).
direction (int, optional) – Propagation direction, 1 for forward and -1 for backward. Defaults to 1.

Returns:

y (numpy array) – Complex signal after chromatic dispersion

References

[1] Savory, Seb J. “Digital filters for coherent optical receivers.” Optics express 16.2 (2008): 804-817.
[2] Eghbali, Amir, et al. “Optimal least-squares FIR digital filters for compensation of chromatic dispersion in digital coherent optical receivers.” Journal of lightwave technology 32.8 (2014): 1449-1456.

apply_kerr_nonlinearity(x: ndarray, z: float, gamma: float, gain: float = 1, direction: int = 1, intensity: ndarray | None = None) → ndarray[source]#

Apply Kerr nonlinearity phase rotation to a signal.

The rotation is \(\exp(j \gamma z I)\) where \(I\) is the intensity driving it. For a single field that is \(|x|^2\); for the Manakov model of two polarizations it is the total intensity \(|E_x|^2 + |E_y|^2\) shared by both, which is what intensity is for.

Parameters:

x (np.ndarray) – Input complex signal.
z (float) – Fiber length in km.
gamma (float) – Kerr coefficient in rad/W/km. The Manakov model passes \(8\gamma/9\) here: the factor belongs to the model, not to the fibre.
gain (float, optional) – Gain factor. Default is 1.
direction (int, optional) – Propagation direction (1=forward, -1=backward). Default is 1.
intensity (np.ndarray, optional) – Intensity driving the rotation, broadcastable against x. Default (None) uses |x|**2, the scalar NLSE.

Returns:

np.ndarray – Signal after Kerr nonlinear phase rotation.

compute_erbium_doped_fiber_amplifier_gain(alpha_dB: float, L_span: float) → float[source]#

Compute the amplitude gain of an EDFA that compensates for fiber loss.

Parameters:

alpha_dB (float) – Fiber loss in dB/km.
L_span (float) – Span length in km.

Returns:

float – Amplitude gain factor.

compute_erbium_doped_fiber_N_ase(alpha_dB: float, L_span: float, NF_dB: float, h: float = 6.62607015e-34, nu: float = 193414489032258.06) → float[source]#

Compute ASENoise params

Parameters:

alpha_dB (float) – Fiber loss (dB/km)
L_span (float) – Length of the link (in km)
NF_dB (float) – Noise Figure (dB)

The formula is given by :

\[N_{ASE} = (e^{\alpha L}-1) h \nu n_{sp}\]

where

\[\alpha = \frac{\alpha_{dB}}{10 \log_{10} e} = \frac{\ln 10}{10}\,\alpha_{dB}\]

References

[1] Essiambre, René-Jean, Gerhard Kramer, Peter J. Winzer, Gerard J. Foschini, and Bernhard Goebel. “Capacity limits of optical fiber networks.” Journal of Lightwave technology 28, no. 4 (2010): 662-701.

get_linear_step_size(L_span: float, StPS: int) → ndarray[source]#

Compute uniformly spaced step sizes for the split-step Fourier method.

Parameters:

L_span (float) – Span length in km.
StPS (int) – Number of steps per span.

Returns:

np.ndarray – Array of step sizes, each equal to L_span / StPS.

get_logarithmic_step_size(L_span: float, StPS: int, alpha_dB: float = 0, step_log_factor: float = 0.4) → ndarray[source]#

Compute logarithmically spaced step sizes for the split-step Fourier method.

Parameters:

L_span (float) – Span length in km.
StPS (int) – Number of steps per span.
alpha_dB (float, optional) – Fiber loss in dB/km. Default is 0.
step_log_factor (float, optional) – Logarithmic step factor. Default is 0.4.

Returns:

np.ndarray – Array of logarithmically spaced step sizes.

References

[1] O. V. Sinkin et al., “Optimization of the split-step Fourier method in modeling optical-fiber communications systems,” J. Lightwave Technol., vol. 21, no. 1, pp. 61-68, Jan. 2003.

itu_grid_frequency(n: int, m: int = 1) → tuple[float, float][source]#

Center frequency and slot width of an ITU-T G.694.1 flexible-grid channel.

A WDM channel is described by the integer couple \((n, m)\) (decision D19): the center frequency is \(193.1 + n \times 0.00625\) THz and the slot width is \(12.5 \times m\) GHz. The fixed grids are particular cases (e.g. the 50 GHz grid uses even \(n\) multiples of 8 and \(m = 4\)).

Parameters:

n (int) – Signed channel index on the 6.25 GHz granularity.
m (int, optional) – Slot width multiplier (slot width = 12.5 * m GHz). Default is 1.

Returns:

center_Hz (float) – Channel center frequency in Hz.
width_Hz (float) – Slot width in Hz.

References

ITU-T Recommendation G.694.1 (2020), section 7 (flexible DWDM grid).

Examples

>>> center, width = itu_grid_frequency(0, m=4)
>>> print(f"{center/1e12:.4f} THz, {width/1e9:.1f} GHz")
193.1000 THz, 50.0 GHz
>>> center, width = itu_grid_frequency(-32, m=4)
>>> print(f"{center/1e12:.4f} THz")
192.9000 THz

is_polarization_pair(x: ndarray, block: str) → bool[source]#

Read the propagation model off the shape of the field.

A (N,) or (..., 1, N) field is one polarization and obeys the scalar NLSE; a (..., 2, N) field is a polarization pair and obeys the Manakov equation. Anything else on the antenna axis is refused, because a pointwise Kerr step applied row by row would describe parallel fibres, not one fibre carrying several signals.

The leading axes are a batch. Only -2 and -1 are read, so (B, P, N) propagates B independent realizations in one call and every block downstream broadcasts over them – the same rule the MIMO blocks follow on their antenna axis (D2). That is also why the polarization axis must be written even when there is one of them: (B, N) is indistinguishable from (P, N), so a batch of single-polarization fields is (B, 1, N).

Parameters:

x (np.ndarray) – Field the block is about to propagate.
block (str) – Name of the calling block, for the message.

Returns:

bool – True when the Manakov model applies.

Raises:

ShapeError – If the polarization axis has any size other than 1 or 2.

manakov_kerr(x: ndarray, gamma: float, manakov: bool) → tuple[float, ndarray | None][source]#

Kerr coefficient and driving intensity of one split step.

The Manakov equation shares the total intensity \(|E_x|^2 + |E_y|^2\) between the two polarizations and scales the coefficient by \(8/9\), the average of the nonlinearity over the fast random birefringence of a real fibre. That factor is a property of the model, not of the glass, which is why it lives here and not in FiberSpec.

Parameters:

x (np.ndarray) – Field at the current step.
gamma (float) – Kerr coefficient of the fibre, in rad/W/km.
manakov (bool) – Whether the field is a polarization pair.

Returns:

gamma (float) – Coefficient to use in the rotation.
intensity (np.ndarray or None) – Intensity driving it, or None for the scalar NLSE.

dbm_to_watt(power_dBm: ndarray | float) → ndarray | float[source]#

Optical power from dBm to watts.

Signal Model#

\[P[\mathrm{W}] = 10^{-3} \times 10^{P[\mathrm{dBm}] / 10}\]

Optical data sheets, standards and papers quote power in dBm while every equation in this package takes watts, so the conversion sits at the boundary of nearly every optical script. Written out by hand it is easy to get subtly wrong – the factor is \(10^{-3}\) and the divisor is 10, not 20, because this is a power and not an amplitude.

Axes: element-wise.

param power_dBm:: Power in dBm, i.e. decibels relative to one milliwatt.
type power_dBm:: float or np.ndarray
returns:: float or np.ndarray – Power in watts.

References

G. P. Agrawal, Fiber-Optic Communication Systems, 4th ed., Wiley, 2010, Section 1.2.

Examples

>>> print(dbm_to_watt(0.0))
0.001
>>> print(round(dbm_to_watt(-3.0), 7))
0.0005012
>>> print(np.round(dbm_to_watt(np.array([-10.0, 0.0, 10.0])), 6))
[0.0001 0.001  0.01  ]

watt_to_dbm(power_W: ndarray | float) → ndarray | float[source]#

Optical power from watts to dBm, the inverse of dbm_to_watt().

Signal Model#

\[P[\mathrm{dBm}] = 10 \log_{10}\!\left( \frac{P[\mathrm{W}]}{10^{-3}}\right)\]

Axes: element-wise.

param power_W:: Power in watts. Must be positive: zero power is \(-\infty\) dBm and negative power is not a power.
type power_W:: float or np.ndarray
returns:: float or np.ndarray – Power in dBm.
raises ValueError:: If any value is not strictly positive. A silent -inf or nan propagating into a link budget is worse than a stop.

References

G. P. Agrawal, Fiber-Optic Communication Systems, 4th ed., Wiley, 2010, Section 1.2.

Examples

>>> print(watt_to_dbm(1e-3))
0.0
>>> print(round(watt_to_dbm(5.0118723e-4), 4))
-3.0
>>> print(np.round(watt_to_dbm(np.array([1e-4, 1e-3, 1e-2])), 6))
[-10.   0.  10.]

launch_amplitude(power_W: ndarray | float, *, polarizations: int = 1) → ndarray | float[source]#

Field amplitude that launches a given optical power.

Signal Model#

An Amplifier multiplies the field, and a power is the squared modulus of it, so a launch power is set through a square root:

\[a = \sqrt{\frac{P}{P_{\mathrm{pol}}}}\]

The division is the point. A dual-polarization signal carries the channel power split between its two polarizations, so each one is launched at \(\sqrt{P/2}\) – and that factor, written by hand at each call site, is the same one that makes an ASE budget or a nonlinear coefficient disagree by 3 dB (see budget()).

Axes: element-wise – an array of powers gives an array of amplitudes, which is what a launch-power sweep needs.

param power_W:: Channel power in watts, summed over the polarizations.
type power_W:: float or np.ndarray
param polarizations:: 1 (default) for a single-polarization field, 2 for a polarization pair.
type polarizations:: int, optional, keyword-only
returns:: float or np.ndarray – Amplitude gain to give an Amplifier.
raises ValueError:: If a power is negative, or if polarizations is neither 1 nor 2.

Examples

>>> print(f"{launch_amplitude(dbm_to_watt(0.0)):.6f}")
0.031623

Two polarizations share the channel power, so each carries half of it – 3 dB down, which is \(\sqrt{2}\) in amplitude:

>>> single = launch_amplitude(1e-3)
>>> pair = launch_amplitude(1e-3, polarizations=2)
>>> print(f"{single / pair:.6f}")
1.414214

linear_step_transfer(n_samples: int, z: float, beta2: float, *, fs: float = 1, alpha_dB: float | None = None, direction: int = 1, like: ndarray | None = None) → ndarray[source]#

Transfer function of one linear step of the split-step method.

Signal Model#

The linear operator of the NLSE is diagonal in frequency, and it is the whole of it – attenuation and dispersion together:

\[H[k] = e^{-\frac{\alpha}{2} z d}\, e^{j \frac{\beta_2}{2} z \omega_k^2 d}\]

with \(\alpha = \frac{\ln 10}{10}\alpha_{dB}\), \(\omega_k = 2\pi f_s f[k]\) over the FFT grid, and \(d = \pm 1\) the propagation direction. The two factors are the same operator \(D\), one real and constant across the band, the other complex and frequency-dependent; separating them would make a caller carry a scalar alongside the array for nothing.

Why this is a function of its own. It depends only on \((n, z, \beta_2, f_s, \alpha_{dB}, d)\) – never on the signal – while a split-step loop applies it thousands of times. Building it costs 0.42 ms at 12288 samples against 0.77 ms for the FFT pair it serves, of which the complex exponential alone is 88 %, so recomputing it per step is most of the cost of propagating. FiberLink builds one per distinct step length in prepare and reuses them.

Axes: element-wise – returns a 1-D array of n_samples bins, which broadcasts against any field shaped (..., n_samples).

param n_samples:: Length of the FFT, i.e. of the field’s last axis.
type n_samples:: int
param z:: Step length in km.
type z:: float
param beta2:: Group-velocity dispersion in ps^2/km.
type beta2:: float
param fs:: Sampling frequency in Hz. Default 1.
type fs:: float, optional, keyword-only
param alpha_dB:: Attenuation in dB/km. Default None, a lossless step.
type alpha_dB:: float, optional, keyword-only
param direction:: 1 forward, -1 backward (back-propagation). Default 1.
type direction:: int, optional, keyword-only
param like:: Array whose type the FFT grid should follow, so that a cupy field gets a cupy transfer function (D3).
type like:: np.ndarray, optional, keyword-only
returns:: np.ndarray – Transfer function on the FFT grid, in FFT order.

Examples

A lossless step is unit modulus – dispersion moves phase, not energy:

>>> H = linear_step_transfer(8, 10.0, -21.7, fs=100e9)
>>> print(np.round(np.abs(H), 12))
[1. 1. 1. 1. 1. 1. 1. 1.]

Backward is the conjugate of forward, which is what makes back-propagation exact in the absence of noise:

>>> back = linear_step_transfer(8, 10.0, -21.7, fs=100e9, direction=-1)
>>> print(np.round(np.max(np.abs(back - np.conj(H))), 12))
0.0

apply_frequency_response(x: ndarray, H: ndarray, *, axis: int = -1) → ndarray[source]#

Filter a signal by a transfer function, in the frequency domain.

Signal Model#

\[y = \mathcal{F}^{-1}\left\{ H \odot \mathcal{F}\{x\} \right\}\]

The multiplication is periodic: no zero padding, so the result is a circular convolution and y has the length of x. That is what a split-step wants, and what a channel mask or a Raman tilt wants; a filter whose impulse response must not wrap around – pulse shaping – pads instead, which is a different operation and lives in comnumpy.core.filters.

Axes: declared axis – filters along axis and broadcasts over the rest, so a (B, P, N) field is filtered by a length-N response without a loop.

param x:: Signal to filter.
type x:: np.ndarray
param H:: Transfer function on the FFT grid, broadcastable against x.
type H:: np.ndarray
param axis:: Axis to filter along. Default -1.
type axis:: int, optional, keyword-only
returns:: np.ndarray – Filtered signal, same shape as x.

Examples

>>> x = np.array([1.0, 0.0, 0.0, 0.0])
>>> print(np.round(np.real(apply_frequency_response(x, np.ones(4))), 12))
[1. 0. 0. 0.]

step_transfers(n_samples: int, lengths: ndarray, *, beta2: float, fs: float, alpha_dB: float | None = None, direction: int = 1, previous: Tuple[ndarray, ndarray, Tuple[int, float, float, float | None, int, Tuple[float, ...]]] | None = None) → Tuple[ndarray, ndarray, Tuple[int, float, float, float | None, int, Tuple[float, ...]]][source]#

Transfer functions of a split-step schedule, one per distinct length.

Signal Model#

A split step applies linear_step_transfer() once per step, and the operator depends only on the step length. A linear schedule makes every step the same length, so one transfer function serves all of them whatever StPS is; a logarithmic schedule makes them all different and needs StPS. Both cases are found the same way – by asking which lengths are distinct – so the propagation loop indexes a table and never branches on the schedule.

Rebuilding costs a complex exponential per distinct length, which is 88 % of the cost of a transfer function and 60 % of the FFT pair it serves, so the table is keyed on everything it depends on and returned unchanged when the key matches. A Monte-Carlo runs the same block at the same length over and over, and then it is built once for the whole sweep instead of once per pass.

Axes: element-wise – the table is (n_distinct, n_samples) and each row broadcasts against a field shaped (..., n_samples).

param n_samples:: Length of the FFT, i.e. of the field’s last axis.
type n_samples:: int
param lengths:: Length of every step, in km, in the order they are applied.
type lengths:: np.ndarray
param beta2:: Group-velocity dispersion in ps^2/km.
type beta2:: float, keyword-only
param fs:: Sampling frequency in Hz.
type fs:: float, keyword-only
param alpha_dB:: Attenuation in dB/km. Default None, a lossless step.
type alpha_dB:: float, optional, keyword-only
param direction:: 1 forward, -1 for back-propagation. Default 1.
type direction:: int, optional, keyword-only
param previous:: The (table, index, key) returned by the last call, returned as is when nothing it depends on has changed.
type previous:: tuple, optional, keyword-only
returns:: tuple – (table, index, key): the (n_distinct, n_samples) transfer functions, the row each step reads, and the key the three were built from.

Examples

A linear schedule collapses to one row, however many steps:

>>> table, index, _ = step_transfers(64, np.full(50, 1.6), beta2=-21.7,
...                                  fs=100e9)
>>> print(table.shape, index[:4], index[-1])
(1, 64) [0 0 0 0] 0

A logarithmic one keeps them apart:

>>> table, index, _ = step_transfers(64, np.array([0.5, 1.5, 4.0]),
...                                  beta2=-21.7, fs=100e9)
>>> print(table.shape, index)
(3, 64) [0 1 2]

Asking again with the same schedule returns the same object, which is what keeps a sweep from rebuilding it:

>>> again = step_transfers(64, np.array([0.5, 1.5, 4.0]), beta2=-21.7,
...                        fs=100e9, previous=(table, index, _))
>>> print(again[0] is table)
True