diff --git a/Constellations/ReshapedPSKConstellation.py b/Constellations/ReshapedPSKConstellation.py new file mode 100644 index 0000000..12cf566 --- /dev/null +++ b/Constellations/ReshapedPSKConstellation.py @@ -0,0 +1,81 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Sander Wahls (TU Delft) 2018 +# Shrinivas Chimmalgi (TU Delft) 2018 + +import numpy as np +import scipy.integrate as integrate +import warnings + +from Constellations import MPSKConstellation + +class ReshapedPSKConstellation(MPSKConstellation): + """Reshaped phase shift keying (PSK) constellation (implements + BaseConstellation). + + The goal of the reshaping procedure is to fit the average (normalized) + energy of the generated pulses to a desired value Ed. For details, see + Gui et al., Opt. Express 26(21), 2018. + """ + + def __init__(self, m, b0_fun, Ed, bnds): + """Constructor for a reshaped m x n PSK constellation. + + Parameters + ---------- + m : int + + b0_fun : function + Carrier waveform. The function should maps any input vector of the + type numpy.array(float), which represents a vector of nonlinear + frequencies xi, to a output vector of the type numpy.array(complex), + which represents the values of the carrier waveform at these xi. + Ed : float + Desired average energy of the generated pulses (with respect to + normalized units). Should be positive. + bnds : numpy.array(float) + Vector with two entries [a,b], which are used as initial bounds in + the bisection produdure based on which the constellation is + reshaped. It should be 0`_ Optical Networking and Communication Conference & Exhibition (OFC), paper M3Z.13, Mar. 2019. +Please consider citing this paper if NFDMLab contributes to your research. + First Steps ----------- @@ -37,6 +39,7 @@ Acknowledgements ---------------- - This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 716669). +- This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 766115. - Fiber-topic transmissions are simulated using a slightly modified port of `SSPROP `_. - Alexander Geisler has contributed a Raman gain profile. diff --git a/Documentation/modulators.rst b/Documentation/modulators.rst index 8939f31..1d93ed9 100644 --- a/Documentation/modulators.rst +++ b/Documentation/modulators.rst @@ -12,5 +12,9 @@ Modulators :special-members: .. autoclass:: Modulators.DiscSpecModulator + :members: + :special-members: + +.. autoclass:: Modulators.TimeDomainPulseShapingModulator :members: :special-members: \ No newline at end of file diff --git a/Example_b_Modulation_SMF_DDF.py b/Example_b_Modulation_SMF_DDF.py new file mode 100644 index 0000000..f5840f6 --- /dev/null +++ b/Example_b_Modulation_SMF_DDF.py @@ -0,0 +1,58 @@ +import Examples +import sys +import numpy as np +import scipy.io as scio + +Fiber_type = 'DDF' # 'DDF' or 'SSMF' + +if Fiber_type == 'DDF': + example = Examples.BajajEtAl2020_b_modulation() + example.n_steps_per_span = 500 + example.path_average = False + example.beta2 = -25.491e-27 + example.gamma = 1.3e-3 +elif Fiber_type == 'SSMF': + example = Examples.GuiEtAl2018() + example.n_steps_per_span = 200 + Factor_gamma_eff_80km = 0.2646 + example.beta2 = -25.49e-27*Factor_gamma_eff_80km #-8.414e-27 + example.gamma = 1.3e-3 + example.path_average = True +else: + raise Exception('Specify the correct fiber type') +example.n_spans = 1 +#example.gamma = 1.3e-3 #2.2e-3 +#example.alpha = 0 +example.constellation_type = 'PSK' +example.constellation_level = 8 +example.n_symbols_per_block = 9 +example.fiber_span_length = 80e3 +example.noise = True +example.noise_figure = 6 +example.tx_bandwidth = 40e9 +example.rx_bandwidth = 40e9 +example.Ed = 16 +example.reconfigure() +runs = 10 +tx_data, rx_data = example.run(runs) + + +#Rotation compensation +########################################################################################### +#Compensate for the rotation mismatch +RotationCompensation = 1 +if RotationCompensation == 1: + #Find the average rotation error + Rotation_Error = (np.angle(rx_data["symbols"])-np.angle(tx_data["symbols"])-np.pi)%(2*np.pi)+np.pi + MeanRotation = np.mean(Rotation_Error) + #Compensate for the average rotation + rx_data["symbols"] = rx_data["symbols"] * np.exp(-1j*MeanRotation) +########################################################################################### + +#Plot the result + +example.evaluate_results(tx_data, rx_data) + + +print('Energy per carrier',example.Ed) +print('Fiber type : ', Fiber_type) diff --git a/Examples/BajajEtAl2020_b_modulation.py b/Examples/BajajEtAl2020_b_modulation.py new file mode 100644 index 0000000..9bb0730 --- /dev/null +++ b/Examples/BajajEtAl2020_b_modulation.py @@ -0,0 +1,178 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Vinod Bajaj (TU Delft) 2019, 2021 +# Sander Wahls (TU Delft) 2018-2019 +# Shrinivas Chimmalgi (TU Delft) 2018 +# Sander Wahls (KIT), 2024 + +import numpy as np +import math + +from Examples import BaseExample + +class BajajEtAl2020_b_modulation(BaseExample): + '''This example recreates the b-modulation transmission in dispersion-decreasing + fiber from the paper + + "Exact NFDM Transmission in the Presence of Fiber-Loss" by V. Bajaj, + S. Chimmalgi, V. Aref, and S. Wahls, + + published in the Journal of Lightwave Technology 38(11), 2020.''' + + def __init__(self): + # Fiber parameters + + self.beta2 = -25.491e-27 + """Dispersion parameter in s**2/m.""" + + self.gamma = 1.3e-3 + """Nonlinearity parameter in (W m)**(-1).""" + + self.fiber_type = "DDF" + """Fiber type: "DDF" for dispersion decreasing fiber or "SSMF" for standard single mode fiber.""" + + self.Tscale = 1.25e-9 # s + """Time scale used during normalization in s.""" + + self.alpha = 0.2e-3 + """Loss coefficient in 1/m.""" + + self.post_boost = True + """Boost at the end of each span (lumped amplification). True or + False.""" + + self.path_average = False + """Use path-averaged fiber parameters during normalization. True or + False.""" + + self.noise = True + """Add ASE noise during amplification (only if post_boost == True). + True or False.""" + + self.noise_figure = 3 + """Noise figure in dB.""" + + self.n_spans = 8 + """Number of fiber spans.""" + + self.fiber_span_length = 80e3 + """Length of a fiber span in m.""" + + self.n_steps_per_span = 500 + """Number of spatial steps per span during the numerical simulation of + the fiber transmission.""" + + # Modulator parameters + self.constellation_type = 'PSK' + """We support 'QAM' and 'PSK'.""" + + self.constellation_level = 16 + """Level of the QAM constellation used as input for the reshaping + process (4, 16, 256, ...).""" + + self.n_symbols_per_block = 9 + """Number of carriers.""" + + self.Ed = 4 + """Desired average pulse energy in normalized units. Used during the + reshaping of the constellation. (See the paper for details.)""" + + # Filters + + self.tx_bandwidth = 33e9 + """Bandwidth of the (ideal) low-pass filter at the transmitter in Hz.""" + + self.rx_bandwidth = 33e9 # in GHz + """Bandwidth of the (ideal) low-pass filter at the receiver in Hz.""" + + self.reconfigure() + + def reconfigure(self): + distance = self.n_spans*self.fiber_span_length #m + from Modulators.CarrierWaveforms import flat_top + carrier_waveform = lambda xi : flat_top(xi, T0) + self._carrier_spacing = 15.0 + T0 = 4.5 + T = np.array([-2.25, 2.25]) # changed from -1 1 + + # Normalization + + if self.path_average == True: + from Normalization import Lumped + self._normalization = Lumped(self.beta2, + self.gamma, + self.Tscale, + alpha=np.mean(self.alpha)*np.log(10)*0.1, + zamp=self.fiber_span_length) + else: + from Normalization import Lossless + self._normalization = Lossless(self.beta2, self.gamma, self.Tscale) + + # Constellation + if self.constellation_type == 'QAM': + from Constellations import ReshapedQAMConstellation + m = int(math.sqrt(self.constellation_level)) + assert self.constellation_level == m*m + bnds = np.array([0, 4/np.abs(carrier_waveform(0.0))]) + self._constellation = ReshapedQAMConstellation(m, m,carrier_waveform,self.Ed, bnds) + elif self.constellation_type == 'PSK': + from Constellations import ReshapedPSKConstellation + bnds = np.array([0,4/np.abs(carrier_waveform(0.0))]) + self._constellation = ReshapedPSKConstellation(self.constellation_level,carrier_waveform,self.Ed,bnds) + else: + raise Exception('Constellation not supported') + + + from Links.DDF_profile import Get_Beta_Gamma_Profile + dz = self.fiber_span_length/self.n_steps_per_span + profile = Get_Beta_Gamma_Profile(self.alpha * np.log(10) * 0.1, self.beta2, self.gamma, dz, + self.n_steps_per_span) + + + + # Modulator + + from Modulators import ContSpecModulator + normalized_distance = self.normalization.norm_dist(distance) + normalized_duration = T[1] - T[0] + required_normalized_dt = normalized_duration/self.n_symbols_per_block/128 # originally it was 8 instead of 128 in denominator + required_dxi = self._carrier_spacing / 8 + self._modulator = ContSpecModulator(carrier_waveform, + self._carrier_spacing, + self.n_symbols_per_block, + normalized_distance*profile['avg_D_z'], + T, + required_normalized_dt, + required_dxi, + "b", + matched_filter_percentage=70) + + # Link + + from Links import SMFSplitStep + dt = self._normalization.denorm_time(self.modulator.normalized_dt) + self._link = SMFSplitStep(dt, dz, self.n_steps_per_span, self.fiber_type, + self.alpha, self.beta2, self.gamma, False, self.n_spans, + self.post_boost, self.noise, self.noise_figure) + + # Filters + + from Filters import PassThrough + self._tx_filter = PassThrough() + from Filters import FFTLowPass + self._rx_filter = FFTLowPass(self.rx_bandwidth/2, dt) diff --git a/Examples/BajajEtAl2020_multi_soliton.py b/Examples/BajajEtAl2020_multi_soliton.py new file mode 100644 index 0000000..ab2b992 --- /dev/null +++ b/Examples/BajajEtAl2020_multi_soliton.py @@ -0,0 +1,159 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Vinod Bajaj (TU Delft) 2019, 2021 +# Sander Wahls (TU Delft) 2018-2019 +# Shrinivas Chimmalgi (TU Delft) 2018 +# Sander Wahls (KIT), 2024 + +import numpy as np +import math + +from Examples import BaseExample + +class BajajEtAl2020_multi_soliton(BaseExample): + '''This example recreates the multi-soliton transmission in dispersion-decreasing + fiber from the paper + + "Exact NFDM Transmission in the Presence of Fiber-Loss" by V. Bajaj, + S. Chimmalgi, V. Aref, and S. Wahls, + + published in the Journal of Lightwave Technology 38(11), 2020.''' + + + def __init__(self): + # Fiber parameters + + self.beta2 =-25.491e-27 + """Dispersion coefficient in s**2/m.""" + + self.gamma = 1.36e-3 + """Nonlinearity coefficient in (W m)**(-1).""" + + self.fiber_type = "DDF" + """Fiber type: "DDF" dispersion for decreasing fiber or "SSMF" for standard single mode fiber.""" + + self.Tscale = 4.1022e-11 + + """Time scale used during normalization in s.""" + + self.alpha = 0.2e-3 + """Loss coefficient in 1/m.""" + + self.n_spans = 8 + """Number of fiber spans.""" + + self.n_steps_per_span = 500 + """Number of spatial steps per fiber span during simulations. Use more than 2 steps per km length for the transmission over the DDF for better accuracy.""" + + self.fiber_span_length = 80e3 #If you change it, then change the average dispersion (D_z) in the demodulator + """Length of a fiber span in m.""" + + self.post_boost = True + """Boost at end of each span (lumped amplification). True or False.""" + + self.path_average = False + """Use path-average fiber parameters during normalization. True or False. False for the transmission over DDF.""" + + self.noise = False + """Add ASE noise (lumped amplification only). True or False.""" + + self.noise_figure = 6 + """Noise figure in dB.""" + + # Receiver bandwidth + + self.tx_bandwidth = 33*1e9 + """Bandwidth of the ideal low-pass at the transmitter in Hz.""" + + self.rx_bandwidth = 33*1e9 # Hz + """Bandwidth of the ideal low-pass at the receiver in Hz.""" + + # Modulator parameters + + self.constellation_level = 4 + """Level of the QAM constellation (4, 16, 256, ...).""" + + self.eigenvalues = np.array([0.45j-0.6, 0.3j-0.4, 0.45j-0.2, 0.3j, 0.45j+0.2, 0.3j+0.4, 0.45j+0.6]) + """Eigenvalue pattern.""" + + self.residues_amplitude = np.exp(np.array([11.85, 7.06, 7.69, 3.81, 1.93, -0.62, -5.43])) + """Spectral amplitudes for each of the eigenvalues. These values get + multiplied with symbols drawn from the constellation before pulse + generation). """ + + + self.reconfigure() + + def reconfigure(self): + distance = self.n_spans*self.fiber_span_length #m + assert(np.size(self.eigenvalues) == np.size(self.residues_amplitude)) + T = np.array([-9*np.pi, 9*np.pi]) + + + # Normalization + + if self.path_average == True: + from Normalization import Lumped + self._normalization = Lumped(self.beta2, self.gamma, + self.Tscale, + alpha=np.mean(self.alpha)*np.log(10)*0.1, + zamp=self.fiber_span_length) + else: + from Normalization import Lossless + self._normalization = Lossless(self.beta2, self.gamma, self.Tscale) + + + # Constellation + + from Constellations import QAMConstellation + m = int(math.sqrt(self.constellation_level)) + assert self.constellation_level == m*m + self._constellation = QAMConstellation(m, m) + + from Links.DDF_profile import Get_Beta_Gamma_Profile + dz = self.fiber_span_length/self.n_steps_per_span + nz = self.n_spans*self.n_steps_per_span + profile = Get_Beta_Gamma_Profile(self.alpha * np.log(10) * 0.1, self.beta2, self.gamma, dz, + self.n_steps_per_span) + + # Modulator + + from Modulators import DiscSpecModulator + normalized_distance = self.normalization.norm_dist(distance) + required_normalized_dt = (T[1] - T[0])/512 + self._modulator = DiscSpecModulator(self.eigenvalues, + self.residues_amplitude, + normalized_distance*profile['avg_D_z'], + T, required_normalized_dt) + + # Link + + from Links import SMFSplitStep # Propagation in Dispersion Decreasing Fiber + dt = self.normalization.denorm_time(self.modulator.normalized_dt) + + self._link = SMFSplitStep(dt, dz, self.n_steps_per_span, self.fiber_type, + self.alpha, self.beta2, self.gamma, + False, self.n_spans, self.post_boost, + self.noise, self.noise_figure) + + # Filters + + from Filters import PassThrough + from Filters import FFTLowPass + self._tx_filter = PassThrough() + self._rx_filter = FFTLowPass(self.rx_bandwidth/2, dt) diff --git a/Examples/BaseExample.py b/Examples/BaseExample.py index b00e41b..431784a 100644 --- a/Examples/BaseExample.py +++ b/Examples/BaseExample.py @@ -270,6 +270,8 @@ def evaluate_results(self, tx_data, rx_data): rx_data["nfspecs"][0].show(legend=["tx", "rx"]) constellation_diagram = ConstellationDiagram(self.constellation) + # scaling the Rx symbols, to the rms power of Tx symbols, This will isolate the effect of scaling from EVM computation + rx_data["symbols"] = np.sqrt(np.mean(np.abs(tx_data["symbols"])**2))*rx_data["symbols"]/np.sqrt(np.mean(np.abs(rx_data["symbols"])**2)) constellation_diagram.plot(rx_data["symbols"], new_fig=True) error_vector_magnitude = ErrorVectorMagnitude() @@ -282,7 +284,7 @@ def evaluate_results(self, tx_data, rx_data): eff, br, bw_rx = modulation_efficiency.compute(rx_data["t"], rx_data["q"], rx_data["q"], nerr, nbits) _, _, bw_tx = modulation_efficiency.compute(tx_data["t"], tx_data["q"], tx_data["q"], 0, 1) - tx_power_level = np.mean(np.abs(rx_data["q"])**2) + tx_power_level = np.mean(np.abs(tx_data["q"])**2) tx_power_level_in_dBm = 10*np.log10(tx_power_level / 0.001) link_length = self.link.span_length * self.link.n_spans diff --git a/Examples/BuelowArefIdler2016.py b/Examples/BuelowArefIdler2016.py index a8d5b3d..0917c9e 100644 --- a/Examples/BuelowArefIdler2016.py +++ b/Examples/BuelowArefIdler2016.py @@ -17,6 +17,8 @@ # Contributors: # Sander Wahls (TU Delft) 2018-2019 # Shrinivas Chimmalgi (TU Delft) 2018 +# Vinod Bajaj (TU Delft) 2021 +# Sander Wahls (KIT) 2024 import numpy as np import math @@ -41,10 +43,13 @@ def __init__(self): self.gamma = 1.6e-3 """Nonlinearity coefficient in (W m)**(-1).""" + self.fiber_type = "SSMF" + """Fiber type: "DDF" dispersion decreasing fiber or "SSMF" standard single mode fiber.""" + self.Tscale = 4.5473e-11 """Time scale used during normalization in s.""" - self.alpha = np.array([0.2e-3]) + self.alpha = 0.2e-3 """Loss coefficient in 1/m.""" self.n_spans = 20 @@ -132,7 +137,7 @@ def reconfigure(self): dt = self.normalization.denorm_time(self.modulator.normalized_dt) dz = self.fiber_span_length/self.n_steps_per_span nz = self.n_spans*self.n_steps_per_span - self._link = SMFSplitStep(dt, dz, self.n_steps_per_span, + self._link = SMFSplitStep(dt, dz, self.n_steps_per_span, self.fiber_type, self.alpha, self.beta2, self.gamma, False, self.n_spans, self.post_boost, self.noise, self.noise_figure) diff --git a/Examples/GuiEtAl2018.py b/Examples/GuiEtAl2018.py index f99cb6f..03fb31c 100644 --- a/Examples/GuiEtAl2018.py +++ b/Examples/GuiEtAl2018.py @@ -17,6 +17,8 @@ # Contributors: # Sander Wahls (TU Delft) 2018-2019 # Shrinivas Chimmalgi (TU Delft) 2018 +# Vinod Bajaj (TU Delft) 2019, 2021 +# Sander Wahls (KIT) 2024 import numpy as np import math @@ -40,10 +42,13 @@ def __init__(self): self.gamma = 1.2e-3 """Nonlinearity parameter in (W m)**(-1).""" + self.fiber_type = "SSMF" + """Fiber type: "DDF" for dispersion decreasing fiber or for "SSMF" standard single mode fiber.""" + self.Tscale = 1.25e-9 # s """Time scale used during normalization in s.""" - self.alpha = np.array([0.2e-3]) + self.alpha = 0.2e-3 """Loss coefficient in 1/m.""" self.post_boost = True @@ -71,6 +76,8 @@ def __init__(self): the fiber transmission.""" # Modulator parameters + self.constellation_type = 'QAM' + """We support 'QAM' and 'PSK'.""" self.constellation_level = 16 """Level of the QAM constellation used as input for the reshaping @@ -91,6 +98,11 @@ def __init__(self): self.rx_bandwidth = 33e9 # in GHz """Bandwidth of the (ideal) low-pass filter at the receiver in Hz.""" + self.dX_factor = 8 + """Factor used when determining the time domain step size dt and the + nonlinear frequency domain step size dxi. Both are proportional to + 1/dX_factor.""" + self.reconfigure() def reconfigure(self): @@ -115,22 +127,26 @@ def reconfigure(self): self._normalization = Lossless(self.beta2, self.gamma, self.Tscale) # Constellation - - from Constellations import ReshapedQAMConstellation - m = int(math.sqrt(self.constellation_level)) - assert self.constellation_level == m*m - bnds = np.array([0, 4/np.abs(carrier_waveform(0.0))]) - self._constellation = ReshapedQAMConstellation(m, m, - carrier_waveform, - self.Ed, bnds) + if self.constellation_type == 'QAM': + from Constellations import ReshapedQAMConstellation + m = int(math.sqrt(self.constellation_level)) + assert self.constellation_level == m*m + bnds = np.array([0, 4/np.abs(carrier_waveform(0.0))]) + self._constellation = ReshapedQAMConstellation(m, m,carrier_waveform,self.Ed, bnds) + elif self.constellation_type == 'PSK': + from Constellations import ReshapedPSKConstellation + bnds = np.array([0,4/np.abs(carrier_waveform(0.0))]) + self._constellation = ReshapedPSKConstellation(self.constellation_level,carrier_waveform,self.Ed,bnds) + else: + raise Exception('Constellation not supported') # Modulator from Modulators import ContSpecModulator normalized_distance = self.normalization.norm_dist(distance) normalized_duration = T[1] - T[0] - required_normalized_dt = normalized_duration/self.n_symbols_per_block/8 - required_dxi = self._carrier_spacing / 8 + required_normalized_dt = normalized_duration/self.n_symbols_per_block/self.dX_factor + required_dxi = self._carrier_spacing / self.dX_factor self._modulator = ContSpecModulator(carrier_waveform, self._carrier_spacing, self.n_symbols_per_block, @@ -145,8 +161,8 @@ def reconfigure(self): from Links import SMFSplitStep dt = self._normalization.denorm_time(self.modulator.normalized_dt) dz = self.fiber_span_length/self.n_steps_per_span - self._link = SMFSplitStep(dt, dz, self.n_steps_per_span, self.alpha, - self.beta2, self.gamma, False, self.n_spans, + self._link = SMFSplitStep(dt, dz, self.n_steps_per_span, self.fiber_type, + self.alpha, self.beta2, self.gamma, False, self.n_spans, self.post_boost, self.noise, self.noise_figure) # Filters diff --git a/Examples/LeArefBuelow2017.py b/Examples/LeArefBuelow2017.py index 607b344..a676937 100644 --- a/Examples/LeArefBuelow2017.py +++ b/Examples/LeArefBuelow2017.py @@ -18,6 +18,7 @@ # Sander Wahls (TU Delft) 2018-2019 # Christoph Mahnke 2018 # Marius Brehler (TU Dortmund) 2019 +# Sander Wahls (KIT) 2024 import numpy as np import matplotlib.pyplot as plt @@ -47,7 +48,7 @@ def __init__(self): # Link parameters - self.alpha = np.array([0.2e-3]) + self.alpha = 0.2e-3 """Loss coefficient in 1/m.""" self.beta2 = -5.75e-27 @@ -56,6 +57,9 @@ def __init__(self): self.gamma = 1.6e-3 """Nonlinearity coefficient in (W m)**(-1).""" + self.fiber_type = "SSMF" + """Fiber type: "DDF" for dispersion decreasing fiber or "SSMF" for standard single mode fiber.""" + self.Tscale = 4e-10 """Time scale used during normalization in s.""" @@ -215,6 +219,7 @@ def reconfigure(self): self._link = SMFSplitStep(dt, dz, self.n_steps_per_span, + self.fiber_type, self.alpha, self.beta2, self.gamma, diff --git a/Examples/TimeDomainPulseShaping.py b/Examples/TimeDomainPulseShaping.py new file mode 100644 index 0000000..088d893 --- /dev/null +++ b/Examples/TimeDomainPulseShaping.py @@ -0,0 +1,156 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Sander Wahls (TU Delft) 2019, 2021 + +import numpy as np +import math + +from Examples import BaseExample +from Constellations import QAMConstellation +from Links import SMFSplitStep +from Normalization import Lumped +from Filters import FFTLowPass, DispersionCompensation, Concatenate +from Modulators import TimeDomainPulseShapingModulator, CarrierWaveforms + +class TimeDomainPulseShaping(BaseExample): + """Conventional time domain pulse shaping with raised cosines and a linear + dispersion compensation filter.""" + + def __init__(self): + + self.roll_off_factor = 0.5 + """Roll-off factor for the raised cosine pulses.""" + + self.T0 = 2.5e-11 + """Symbol duration in s.""" + + self.n_symbols_per_block = 128 + """Number of symbols per block.""" + + self.n_guard_symbols = 32 + """The length of the guard interval is given by this number times T0.""" + + self.beta2 = -5e-27 + """Dispersion parameter in s**2/m.""" + + self.gamma = 1.2e-3 + """Nonlinearity parameter in (W m)**(-1).""" + + self.Tscale = 1.25e-9 # s + """Time scale used during normalization in s.""" + + self.alpha = np.array([0.2e-3]) + """Loss coefficient in 1/m.""" + + self.post_boost = True + """Boost at the end of each span (lumped amplification). True or + False.""" + + self.path_average = True + """Use path-averaged fiber parameters during normalization. True or + False.""" + + self.noise = True + """Add ASE noise during amplification (only if post_boost == True). + True or False.""" + + self.noise_figure = 6 + """Noise figure in dB.""" + + self.n_spans = 8 + """Number of fiber spans.""" + + self.fiber_span_length = 80e3 + """Length of a fiber span in m.""" + + self.n_steps_per_span = 40 + """Number of spatial steps per span during the numerical simulation of + the fiber transmission.""" + + self.tx_bandwidth = 40e9 + """Bandwidth of the (ideal) low-pass filter at the transmitter in Hz.""" + + self.rx_bandwidth = 40e9 # in GHz + """Bandwidth of the (ideal) low-pass filter at the receiver in Hz.""" + + self.constellation_level = 4 + """Level of the QAM constellation (4, 16, 256, ...).""" + + self.dispersion_compensation = True + """Compensate for chromatic dispersion at the end of the link. True or + False.""" + + self.tx_gain = 0.065 + """Gain (in linear units) that is applied at the transmitter before + transmission. Use for power control.""" + + self.dt_factor = 16 + """Factor that determines the time domain step size, which is + proportional to 1/dt_factor.""" + + self.make_n_samples_pow2 = False + """If true, it is ensured that the number of samples is a power of two.""" + + self.reconfigure() + + def reconfigure(self): + + self._normalization = Lumped(self.beta2, + self.gamma, + self.T0, + self.alpha, + self.fiber_span_length) + + normalized_T0 = self._normalization.norm_time(self.T0) + self.pulse_fun = lambda t : self.tx_gain*CarrierWaveforms.raised_cosine(t, self.roll_off_factor, normalized_T0) + self.pulse_spacing = normalized_T0 + requested_normalized_dt = normalized_T0 / self.dt_factor + + m = int(math.sqrt(self.constellation_level)) + assert self.constellation_level == m*m + self._constellation = QAMConstellation(m, m) + + self._modulator = TimeDomainPulseShapingModulator(self.pulse_fun, + self.pulse_spacing, + requested_normalized_dt, + self.n_symbols_per_block, + self.n_guard_symbols, + self.make_n_samples_pow2) + self.dt = self._normalization.denorm_time(self._modulator.normalized_dt) + + dz = self.fiber_span_length/self.n_steps_per_span + self._link = SMFSplitStep(self.dt, + dz, + self.n_steps_per_span, + "SSMF", + self.alpha, + self.beta2, + self.gamma, + False, + self.n_spans, + self.post_boost, + self.noise, + self.noise_figure) + + self._tx_filter = FFTLowPass(self.tx_bandwidth/2, self.dt) + rx_lowpass = FFTLowPass(self.rx_bandwidth/2, self.dt) + if self.dispersion_compensation: + dcf = DispersionCompensation(self.beta2, self.fiber_span_length*self.n_spans, self.dt) + self._rx_filter = Concatenate(rx_lowpass, dcf) + else: + self._rx_filter = rx_lowpass diff --git a/Examples/__init__.py b/Examples/__init__.py index 2f8ec2a..0253b96 100644 --- a/Examples/__init__.py +++ b/Examples/__init__.py @@ -16,8 +16,13 @@ # # Contributors: # Sander Wahls (TU Delft) 2019 +# Vinod Bajaj (TU Delft) 2019 +# Sander Wahls (KIT) 2024 from Examples.BaseExample import BaseExample from Examples.BuelowArefIdler2016 import BuelowArefIdler2016 from Examples.LeArefBuelow2017 import LeArefBuelow2017 from Examples.GuiEtAl2018 import GuiEtAl2018 +from Examples.BajajEtAl2020_multi_soliton import BajajEtAl2020_multi_soliton +from Examples.BajajEtAl2020_b_modulation import BajajEtAl2020_b_modulation +from Examples.TimeDomainPulseShaping import TimeDomainPulseShaping diff --git a/Filters/Concatenate.py b/Filters/Concatenate.py new file mode 100644 index 0000000..ef2061e --- /dev/null +++ b/Filters/Concatenate.py @@ -0,0 +1,12 @@ +from Filters import BaseFilter + +class Concatenate(BaseFilter): + '''Concatenates two filters.''' + + def __init__(self, filter1, filter2): + """Constructor.""" + self._filter1 = filter1 + self._filter2 = filter2 + + def filter(self, input): + return self._filter2.filter(self._filter1.filter(input)) diff --git a/Filters/DispersionCompensation.py b/Filters/DispersionCompensation.py new file mode 100644 index 0000000..c714ddc --- /dev/null +++ b/Filters/DispersionCompensation.py @@ -0,0 +1,58 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Sander Wahls (TU Delft) 2018-2019 + +import numpy as np +from Filters import BaseFilter + +class DispersionCompensation(BaseFilter): + '''Dispersion compensation fiter implemented using the FFT.''' + + def __init__(self, beta2, link_length, h_sec): + """Constructor. + + Parameters + ---------- + beta2 : float + link_length : float + h_sec : float + """ + self._beta2 = beta2 + self._link_length = link_length + self._h_sec = h_sec + + def filter(self, input): + """Removes chromatic disperion effects from a signal. + + Parameters + ---------- + input : numpy.array(complex) + Vector of equi-distant time domain samples q[n]=q(t0+n*h_sec), + n=0,1,2,...,N-1, where h_sec is the sampling interval that was + provided to __init__(...). + + Returns + ------- + numpy.array(complex) + Vector of the same properties as the input. + """ + ft = np.fft.fft(input) + N = np.size(ft) + freq_Hz = np.fft.fftfreq(N, d=self._h_sec) + ft *= np.exp(-0.5j*(2*np.pi*freq_Hz)**2*self._link_length*self._beta2) + return np.fft.ifft(ft) diff --git a/Filters/__init__.py b/Filters/__init__.py index c7c7d0b..7c80e95 100644 --- a/Filters/__init__.py +++ b/Filters/__init__.py @@ -20,3 +20,5 @@ from Filters.BaseFilter import BaseFilter from Filters.PassThrough import PassThrough from Filters.FFTLowPass import FFTLowPass +from Filters.DispersionCompensation import DispersionCompensation +from Filters.Concatenate import Concatenate diff --git a/Helpers/SumOfShiftedWaveforms.py b/Helpers/SumOfShiftedWaveforms.py new file mode 100644 index 0000000..38864f7 --- /dev/null +++ b/Helpers/SumOfShiftedWaveforms.py @@ -0,0 +1,115 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Sander Wahls (TU Delft) 2018-2019 + +import numpy as np + +class SumOfShiftedWaveforms: + """Generates a waveform by adding shifted and scaled copies of a base + waveform. The generated waveform is of the form + + s(X) = sum_n symbol[n]*waveform_fun(X - n*waveform_spacing + offset). + """ + + def __init__(self, waveform_fun, waveform_spacing, n_waveforms, grid): + """Constructor. + + Parameters + ---------- + + waveform_fun : function that maps np.array(float) to np.array(complex) + The base waveform. + waveform_spacing : float + The spacing between the copies of the base waveform. + n_waveforms : int + The number of base waveforms that enter each generated waveform. + grid : np.array(float) + Sort array of locations at which the waveform is generated. + """ + + self._waveform_fun = waveform_fun + self._waveform_spacing = waveform_spacing + self._n_waveforms = n_waveforms + self._grid = grid + + self._centers = (np.arange(0, n_waveforms)-(n_waveforms-1)/2.0)*waveform_spacing + self._center_idx = np.zeros(n_waveforms, dtype=int) + for n in range(0, n_waveforms): + self._center_idx[n] = np.argmin(abs(grid - self._centers[n])) + + # Determine interval [left, right] outside of which the waveform is + # effectively zero. This is later used to speed up the modulation + # process when the number of carriers is getting larger. + + vals = waveform_fun(grid) + tol = 10 * np.finfo(vals[0]).eps * np.max(np.abs(vals)) + idx = np.argwhere(np.abs(vals)>tol) + i1 = idx[0][0] + i2 = idx[-1][0] + self._left = grid[i1] + self._right = grid[i2] + + @property + def waveform_fun(self): + """Returns the waveform_fun that was provided to the constructor + (read-only).""" + return self._waveform_fun + + @property + def waveform_spacing(self): + """Returns the waveform_spacing that was provided to the constructor + (read-only).""" + return self._waveform_spacing + + @property + def grid(self): + """Returns the grid that was provided to the constructor (read-only).""" + return self._grid + + @property + def centers(self): + """Returns a vector of grid points that contains the locations at which + the centers of base waveforms are placed when they are summed up. Equal + to the terms n*waveform_spacing+offset in the description of this + class (read-only).""" + return self._centers + + @property + def center_idx(self): + """Returns a vector of indices i[n] such that centers[n] = grid[i[n]] + (read-only).""" + return self._center_idx + + def generate_waveform(self, symbols): + """Generates a waveform as described in the description of this class.""" + nc = np.size(symbols) + assert nc == self._n_waveforms + vals = np.zeros(np.size(self._grid), dtype=complex) + for n in range(0, nc): + shifted_grid = self.grid - self.centers[n] + idx = np.logical_and(shifted_grid>=self._left, shifted_grid<=self._right) + vals[idx] = vals[idx] + symbols[n]*self._waveform_fun(shifted_grid[idx]) + return vals + + def extract_symbols(self, waveform): + """Extracts the symbols from a waveform generated by generate_waveform.""" + symbols = np.zeros(self._n_waveforms, dtype=complex) + scl = self.waveform_fun(0.0) + for n in range(0, self._n_waveforms): + symbols[n] = waveform[self.center_idx[n]] / scl + return symbols diff --git a/Helpers/__init__.py b/Helpers/__init__.py index f9262cd..69f2ae4 100644 --- a/Helpers/__init__.py +++ b/Helpers/__init__.py @@ -18,6 +18,7 @@ # Sander Wahls (TU Delft) 2019 from Helpers.NFSpectrum import NFSpectrum +from Helpers.SumOfShiftedWaveforms import SumOfShiftedWaveforms def checked_get(obj, attr_name, base_class): """Auxiliary function that checks whether the object obj has set an diff --git a/Helpers/widgets.py b/Helpers/widgets.py index 14a80c3..e2dfa19 100644 --- a/Helpers/widgets.py +++ b/Helpers/widgets.py @@ -16,6 +16,7 @@ # # Contributors: # Marius Brehler (TU Dortmund) 2018-2019 +# Marius Brehler 2019 import ipywidgets as widgets from ipywidgets import Layout @@ -157,7 +158,7 @@ def update_alpha(choosen_amplification,noise_widget,path_average_widget): path_average_widget.value = True elif choosen_amplification == "Raman": from scipy.io import loadmat - matfile = loadmat('../data/gainprofile_40steps.mat') + matfile = loadmat('../Data/gainprofile_40steps.mat') selected_alpha = matfile['gain'].flatten() noise_widget.disabled = True noise_widget.value = False @@ -281,4 +282,4 @@ def create_hbox_linkHeader(): hbox_linkH = widgets.HBox([table_header_5_widget]) - return hbox_linkH \ No newline at end of file + return hbox_linkH diff --git a/Links/DDFSplitStep.py b/Links/DDFSplitStep.py new file mode 100644 index 0000000..9a1030b --- /dev/null +++ b/Links/DDFSplitStep.py @@ -0,0 +1,170 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Sander Wahls (TU Delft) 2018-2019 +# Marius Brehler (TU Dortmund) 2018-2019 +# Marius Brehler 2019 + +from Links import BaseLink +from Links._DDFssprop import _DDFssprop +import numpy as np +from scipy.constants import Planck +from Links.DDF_profile import Get_Beta_Gamma_Profile + +class DDFSplitStep(BaseLink): + """Simulates a link with one or several single mode dispersion decreasing fiber spans connected by + EDF amplifiers using a split step method based on SSPROP.""" + + def __init__(self, dt, dz, nz, + alpha=0.0, beta2=1.0, gamma=-1.0, + verbose=False, n_spans=1, + post_boost=False, noise=False, noise_figure=3, + center_frequency = 193.1e12): + """Constructor. + + Parameters + ---------- + + dt : float + Time step in s. + dz : float + Spatial step in m. + nz : int + Number of spatial steps for one fiber span. + alpha : float or numpy.array(float) + Fiber loss coefficient in 1/m. It is possible to pass a vector of + length nz in order to specify an individual loss coefficient for + each of the segments of a span. Noise-free Raman amplification can + be implemented using this feature. + beta2 : float + Fiber dispersion coefficient in s**2/m. + gamma : float + Fiber nonlinearity coefficient in (W m)**(-1). + verbose : bool + Set to True for diagnostic ouputs. + n_spans : int + Number of fiber spans in the link. + post_boost : bool + Accumulated fiber loss is compensated with a gain at the end of each + fiber span if True. Requires alpha to be a scalar. EDFA + amplification can be implemented using this feature. + noise : bool + Add amplified spontaneous emission (ASE) noise at the end of each + fiber span. EDFA amplification can be implemented using this + feature. + noise_figure : float + Noise figure in dB for determining the ASE noise if noise==True. + center_frequency : float + Center frequency in Hz for determining the ASE noise if noise==True. + """ + BaseLink.__init__(self) + self._dt = dt + self._dz = dz + self._nz = nz + self._alpha = alpha *np.log(10)*0.1 + self._gain = self._dz * self._nz * alpha + self._beta2 = beta2 + self._gamma = gamma + self._verbose = verbose + self._n_spans = n_spans + self._post_boost = post_boost + self._noise = noise + self._noise_figure = noise_figure + self._center_frequency = center_frequency + self._span_length = nz*dz + self._n_spans = n_spans + self._BETA2 = None + self._GAMMA = None + self._D_z = None + + def _ASE_noise_power(self): + '''Returns the amplified spontaneous emission (ASE) power using the + equations 7.2.11 and 7.2.15 in the 4th edition of "Fiber-Optic + Communication Systems" by G. P. Agrawal (Wiley 2010).''' + + if self._noise == False: + return 0.0 + G = 10**(self._gain/10.0) + Fn = 10**(self._noise_figure/10.0) + n_sp = (G*Fn - 1.0)/2.0/(G - 1.0) + return np.max([0, n_sp * Planck * self._center_frequency * (G - 1)]) + + def _ASE_noise(self, n_samples): + '''Generates a vector of white Gaussian noise whose variance is the ASE + noise power times the simulation bandwidth.''' + + wgn = np.sqrt(0.5)*(np.random.randn(n_samples) + 1j*np.random.randn(n_samples)) + + # The expected power of c*wgn, where c>0 is t.b.d., is + # + # P1 = E[|c*wgn[0]|**2+...+|c*wgn[n_samples-1]|^2]/n_samples = c^2 + # + # We want this to be equal to the integral of S_ASE(f) = ASE_noise_power + # over the simulation bandwidth 1/dt, i.e., + + P2 = self._ASE_noise_power() / self._dt + # Solving P1=c^2=P2 for c leads to + + c = np.sqrt(P2) + return c*wgn + + def _ASE_noise_sanity_check(self): + '''Plots an averaged periodogram of the ASE noise.''' + import matplotlib.pyplot as plt + fig = plt.figure() + N = 2**18 + K = 100 + X2 = np.zeros(N) + for k in range(0, K): + X2 += np.abs(np.fft.fft(self._ASE_noise(N)))**2 * self._dt / N + X2 /= (K - 1) + f = np.fft.fftfreq(N, 1/self._dt) + plt.semilogy(f, X2) + plt.title("The PSD should be flat at {}".format(self._ASE_noise_power())) + fig.show() + + def transmit(self, input): + # Docstring is inherited from base class. + profile = Get_Beta_Gamma_Profile(self._alpha, self._beta2, self._gamma, self._dz, self._nz) + self._BETA2 = profile['BETA2'] + self._GAMMA = profile['GAMMA'] + + if self._n_spans == 1: + uu = _DDFssprop(input, self._dt, self._dz, self._nz, self._alpha, self._BETA2, self._GAMMA) + + if self._post_boost == True: + if self._alpha.size != 1: + raise Exception('alpha array not supported together with boost') + + uu *= np.exp(self._gain*np.log(10)*0.05) + uu += self._ASE_noise(np.size(uu)) + + else: + uu = input + for span_i in range(0,self._n_spans): + uu = _DDFssprop(uu, self._dt, self._dz, self._nz, self._alpha, self._BETA2, self._GAMMA) + + if self._verbose: + print("Finished span",span_i+1) + + if self._post_boost == True: + if self._alpha.size != 1: + raise Exception('alpha array not supported together with boost') + uu *= np.exp(self._gain*np.log(10)*0.05) + uu += self._ASE_noise(np.size(uu)) + + return uu diff --git a/Links/DDF_profile.py b/Links/DDF_profile.py new file mode 100644 index 0000000..e1ed2cc --- /dev/null +++ b/Links/DDF_profile.py @@ -0,0 +1,94 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Sander Wahls (TU Delft) 2018-2019 +# Shrinivas Chimmalgi (TU Delft) 2019 +# Vinod Bajaj (TU Delft) 2019 + +import numpy as np +def Get_Beta_Gamma_Profile(alpha, beta20, gamma0, dz, nz): + ''' Returns the group velocity dispersion parameter and nonlinear parameter axial profile (vector) ie z dependent values such that + for given attenuation coefficient the fiber model is exactly solvable. + The function solves a cubic equation to find suitable values of GVM paramter and nonlinear parameter. + For details: see + Bajaj, Vinod, et al. "Exact NFDM transmission in the presence of fiber-loss." Journal of Lightwave Technology 38.11 (2020): 3051-3058. + + USAGE + + profile = Get_Beta_Gamma_Profile(alpha, beta20, gamma0, dz, nz) + + + INPUT + + alpha - power loss coefficient, ie, P=P0*exp(-alpha*z) + betap - group velocity dispersion at the starting end of the fiber + gamma - nonlinearity coefficient at the starting end of the fiber + dz - propagation stepsize + nz - number of steps to take, ie, ztotal = dz*nz + + OUTPUT + + profile - a dictionary with fiber profile parameters + with keys: + 'Rad_z' : fiber radius along the length + 'BETA2' : group velocity dispersion parameter along the length + 'GAMMA' : nonlinear parameter along the length + 'D_z' : normalized dispersion parameter profile along the length (normalized w.r.t start-end dispersion value) + 'R_z' : normalized nonlinear parameter profile along the length (normalized w.r.t start-end nonlinear parameter value) + 'avg_D_z' : average normalized dispersion parameter (used for modified NFT) + 'avg_R_z' : average normalized nonlinear parameter + ''' + alpha = alpha/2 # convert to field attenuation. + c = 3e8 # light speed m/s + lambda0 = 1.55e-6 # center wavelength + kappa = lambda0**2*1e-6/(2*np.pi*c) + BETA20 = abs(beta20) + GAMMA0 = abs(gamma0) + R0 = (BETA20/kappa + 20)/8 # core radius in micro meter + n2 = GAMMA0*(lambda0*np.pi*(R0*1e-6)**2)/(2*np.pi) + #dz = dz + #z = np.arange(0,self.n_steps_per_span)*dz + a = BETA20*R0**2 + Rad_z = np.zeros(nz) + for i in range(0,nz): + p = [8*kappa,-20*kappa,0,-a*np.exp(-2*alpha*dz*i)] + Rad_roots = np.roots(p) + ### select the one with lowest imaginary part , then use it's real part + #Rad_img = abs(Rad_roots.imag) + #j = 0 + #if Rad_imag[j]>=Rad_imag[j+1]: + # j = j+1 + #if Rad_imag[j]>=Rad_imag[j+2]: + # j = j+2 + Rad_z[i] = Rad_roots[0].real + + profile = {} + profile['Rad_z'] = Rad_z + profile['BETA2'] = -kappa*(8*np.array(Rad_z)-20) + profile['GAMMA'] = np.divide((2*np.pi*n2), (lambda0*np.pi*(Rad_z*1e-6)**2)) + profile['D_z'] = np.divide(profile['BETA2'], profile['BETA2'][0]) + profile['R_z'] = np.divide(profile['GAMMA'], profile['GAMMA'][0]) + profile['avg_D_z'] = np.mean(profile['D_z']) + profile['avg_R_z'] = np.mean(profile['R_z']) + + # self._BETA2 = -kappa*(8*np.array(Rad_z)-20) + # self._GAMMA = np.divide((2*np.pi*n2),(lambda0*np.pi*(Rad_z*1e-6)**2)) + # self._D_z = np.divide(self._BETA2,self._BETA2[0]) + # self._R_z = np.divide(self._GAMMA,self._GAMMA[0]) + # self._avg_D_z = np.mean(self._D_z) + # self._avg_R_z = np.mean(self._R_z) + return profile diff --git a/Links/SMFSplitStep.py b/Links/SMFSplitStep.py index 7b4dc4e..435e203 100644 --- a/Links/SMFSplitStep.py +++ b/Links/SMFSplitStep.py @@ -17,9 +17,12 @@ # Contributors: # Sander Wahls (TU Delft) 2018-2019 # Marius Brehler (TU Dortmund) 2018-2019 +# Marius Brehler 2019 from Links import BaseLink from Links._ssprop import _ssprop +from Links._DDFssprop import _DDFssprop +from Links.DDF_profile import Get_Beta_Gamma_Profile import numpy as np from scipy.constants import Planck @@ -32,7 +35,7 @@ class SMFSplitStep(BaseLink): polarization to the ASE noise are not included. """ - def __init__(self, dt, dz, nz, + def __init__(self, dt, dz, nz, fiber_type='SSMF', alpha=0.0, beta2=1.0, gamma=-1.0, verbose=False, n_spans=1, post_boost=False, noise=False, noise_figure=3, @@ -48,6 +51,7 @@ def __init__(self, dt, dz, nz, Spatial step in m. nz : int Number of spatial steps for one fiber span. + fiber_type : 'SSMF' standard single mode fiber (default) or 'DDF' dispersion decreasing fiber. alpha : float or numpy.array(float) Fiber loss coefficient in 1/m. It is possible to pass a vector of length nz in order to specify an individual loss coefficient for @@ -78,7 +82,7 @@ def __init__(self, dt, dz, nz, self._dt = dt self._dz = dz self._nz = nz - self._alpha = alpha *np.log(10)*0.05 + self._alpha = alpha *np.log(10)*0.1 self._gain = self._dz * self._nz * alpha self._beta2 = beta2 self._gamma = gamma @@ -90,6 +94,11 @@ def __init__(self, dt, dz, nz, self._center_frequency = center_frequency self._span_length = nz*dz self._n_spans = n_spans + self._fiber_type = fiber_type + if fiber_type == 'DDF': + self._BETA2 = None + self._GAMMA = None + self._D_z = None def _ASE_noise_power(self): '''Returns the amplified spontaneous emission (ASE) power using the @@ -99,6 +108,11 @@ def _ASE_noise_power(self): if self._noise == False: return 0.0 G = 10**(self._gain/10.0) + if isinstance(self._gain, np.ndarray): + assert(len(self._gain) == 1) # no individual loss coeffients per + # span allowed with noise enabled + # at the moment + G = G[0] Fn = 10**(self._noise_figure/10.0) n_sp = (G*Fn - 1.0)/2.0/(G - 1.0) return np.max([0, n_sp * Planck * self._center_frequency * (G - 1)]) @@ -140,15 +154,19 @@ def _ASE_noise_sanity_check(self): def transmit(self, input): # Docstring is inherited from base class. + # Depending on the fiber type choose the fiber propagation solver either the standard ssprop or DDF ssprop. + if self._fiber_type == 'DDF': + profile = Get_Beta_Gamma_Profile(self._alpha, self._beta2, self._gamma, self._dz, self._nz) + self._BETA2 = profile['BETA2'] + self._GAMMA = profile['GAMMA'] + solve_propagation = _DDFssprop + disp_nl_profile = (self._BETA2, self._GAMMA) + else: + solve_propagation = _ssprop + disp_nl_profile = ([0., 0., -self._beta2], -self._gamma) if self._n_spans == 1: - uu = _ssprop(input, - self._dt, - self._dz, - self._nz, - self._alpha, - [0., 0., -self._beta2], - -self._gamma) + uu = solve_propagation(input, self._dt, self._dz, self._nz, self._alpha, disp_nl_profile[0], disp_nl_profile[1]) if self._post_boost == True: if self._alpha.size != 1: @@ -160,13 +178,8 @@ def transmit(self, input): else: uu = input for span_i in range(0,self._n_spans): - uu = _ssprop(uu, - self._dt, - self._dz, - self._nz, - self._alpha, - [0., 0., -self._beta2], - -self._gamma) + uu = solve_propagation(uu, self._dt, self._dz, self._nz, self._alpha, disp_nl_profile[0], disp_nl_profile[1] ) + if self._verbose: print("Finished span",span_i+1) diff --git a/Links/_DDFssprop.py b/Links/_DDFssprop.py new file mode 100644 index 0000000..93d627c --- /dev/null +++ b/Links/_DDFssprop.py @@ -0,0 +1,86 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Sander Wahls (TU Delft) 2018-2019 +# Shrinivas Chimmalgi (TU Delft) 2019 +# Vinod Bajaj (TU Delft) 2019 +import numpy as np +from warnings import warn +import matplotlib.pyplot as plt + +def _DDFssprop(u0,dt,dz,nz,alpha,D,R): + """This function solves the nonlinear Schrodinger equation for + pulse propagation in an optical fiber using the split-step + Fourier method. + It has feature to include length dependent variation in the + group velocity dispersion and nonlinear parameter varies over length. + + The following effects are included in the model: group velocity + dispersion (GVD), loss, and self-phase + modulation (gamma). + + USAGE + + u1 = ssprop(u0,dt,dz,nz,alpha,D,R) + + + INPUT + + u0 - starting field amplitude (vector) + dt - time step + dz - propagation stepsize + nz - number of steps to take, ie, ztotal = dz*nz + alpha - power loss coefficient, ie, P=P0*exp(-alpha*z) + D - dispersion parameter over length of span (vector : size 1 X nz ) + R - nonlinear parameter over length of span (vector : size 1 X nz ) + + + OUTPUT + + u1 - field at the output + + NOTES The dimensions of the input and output quantities can + be anything, as long as they are self consistent. E.g., if + |u|^2 has dimensions of Watts and dz has dimensions of + meters, then gamma should be specified in W^-1*m^-1. + Similarly, if dt is given in picoseconds, and dz is given in + meters, then beta(n) should have dimensions of ps^(n-1)/m. + + AUTHORS + + This function is a Python port of the script ssprop.m by + Thomas E. Murphy (tem@umd.edu), which comes with the SSPROP + package (https://www.photonics.umd.edu/software/ssprop/). + The port has been written by Sander Wahls and has been + modified by Marius Brehler. + """ + + nt = np.size(u0) + assert nt%2 == 0 + w = 2*np.pi*np.hstack((np.arange(0, nt/2), np.arange(-nt/2, 0)))/(dt*nt) + field_attenuation = alpha/2 + u1 = u0 + for iz in range(0, nz): + beta = D[iz] # get beta_2 value at that location + gamma = R[iz] # get gamma value at that location + halfstep = np.exp(1j*beta*w**2/2*dz/2) + + uhalf = np.fft.ifft(halfstep*np.fft.fft(u1)) + uv = uhalf*np.exp(1j*gamma*(np.abs(u1)**2*dz)) + uv = np.fft.ifft(halfstep*np.fft.fft(uv)) + u1 = uv*np.exp(-field_attenuation*dz) + return u1 diff --git a/Links/__init__.py b/Links/__init__.py index 8dd1698..43b7133 100644 --- a/Links/__init__.py +++ b/Links/__init__.py @@ -18,4 +18,5 @@ # Sander Wahls (TU Delft) 2019 from Links.BaseLink import BaseLink -from Links.SMFSplitStep import SMFSplitStep \ No newline at end of file +from Links.SMFSplitStep import SMFSplitStep +from Links.DDFSplitStep import DDFSplitStep \ No newline at end of file diff --git a/Links/_ssprop.py b/Links/_ssprop.py index 45314ed..b22af32 100644 --- a/Links/_ssprop.py +++ b/Links/_ssprop.py @@ -1,3 +1,4 @@ +# Copyright 2019, Marius Brehler # Copyright 2018, Sander Wahls (TU Delft) # Copyright 2018, Marius Brehler (TU Dortmund) # Copyright 2006, Thomas E. Murphy @@ -73,11 +74,7 @@ def _ssprop(u0,dt,dz,nz,alpha,betap,gamma,maxiter=4,tol=1e-5): w = 2*np.pi*np.hstack((np.arange(0, nt/2), np.arange(-nt/2, 0)))/(dt*nt) halfstep = 0 - if np.size(alpha) == 1: - halfstep = -alpha/2 - elif np.size(alpha) == nz: - haltstep = 0 - else: + if (np.size(alpha) != 1) & (np.size(alpha) != nz): raise Exception('alpha is neither constant nor is the vector of length nz') for ii in range(0, np.size(betap)): diff --git a/Modulators/BaseModulator.py b/Modulators/BaseModulator.py index 12d616c..9ef9ce8 100644 --- a/Modulators/BaseModulator.py +++ b/Modulators/BaseModulator.py @@ -31,7 +31,7 @@ class BaseModulator(ABC): During construction, any modulator class should set the attributes - - _norm_dt + - _normalized_dt - _n_samples - _n_symbols_per_block diff --git a/Modulators/CarrierWaveforms.py b/Modulators/CarrierWaveforms.py index fb6aac6..8fa0fbb 100644 --- a/Modulators/CarrierWaveforms.py +++ b/Modulators/CarrierWaveforms.py @@ -15,7 +15,7 @@ # 02111-1307 USA # # Contributors: -# Sander Wahls (TU Delft) 2018 +# Sander Wahls (TU Delft) 2018, 2021 # Shrinivas Chimmalgi (TU Delft) 2018 import numpy as np @@ -48,6 +48,8 @@ def flat_top(xivec, T0): return vals def raised_cosine(xivec, roll_off_factor, T0): + """Raised-cosine carrier. T0 is the desired carrier spacing (i.e., the + distance between two consecutive zeros of the waveform).""" vals = np.sinc(xivec/T0) if np.isscalar(vals): xivec = np.array([xivec]) diff --git a/Modulators/ContSpecModulator.py b/Modulators/ContSpecModulator.py index 5881076..ec9bc72 100644 --- a/Modulators/ContSpecModulator.py +++ b/Modulators/ContSpecModulator.py @@ -27,6 +27,7 @@ from Modulators import BaseModulator from Helpers import NFSpectrum from Helpers import next_pow2 +from Helpers import SumOfShiftedWaveforms class ContSpecModulator(BaseModulator): """This modulator embeds symbols in the continuous spectrum using a classic @@ -43,7 +44,8 @@ def __init__(self, contspec_type, percent_precompensation = 0.0, power_control_factor = 1.0, - use_power_normalization_map = False): + use_power_normalization_map = False, + matched_filter_percentage = 0): """Constructor. Parameters @@ -86,6 +88,13 @@ def __init__(self, in Yangzhang et. al (Proc. OFC'17) to the nonlinear Fourier spectrum it has generated before passing it on to the inverse NFT. Currently only implemented for reflection coefficients (contspec_type="b/a"). + matched_filter_percentage : float + A number between 0 and 100. If >0, the modulator does not only use + b(xi) with xi being the center of the carrier of interest to retrieve + the corresponding symbol, but always considers + matched_filter_percentage percent of the neighboring frequencies in + a matched filter fashion. If matched_filter_percentage==0 (default), + only the center frequency is considered. """ # Save some given parameters for later use. @@ -99,6 +108,7 @@ def __init__(self, self._percent_precompensation = percent_precompensation self._power_control_factor = power_control_factor self._use_power_normalization_map = use_power_normalization_map + self._matched_filter_percentage = matched_filter_percentage # Choose number of time domain samples such that the user requirements # on the time step are fulfilled. @@ -136,11 +146,11 @@ def __init__(self, # Generate some time and nonlinear frequency grids for later use. self._t = self._T[0] + np.arange(0, self.n_samples)*self.normalized_dt - self._xi = self._XI[0] + np.arange(0, self._M)*self._dxi - self._carrier_centers = (np.arange(0, n_symbols_per_block)-(n_symbols_per_block-1)/2.0)*carrier_spacing - self._carrier_center_idx = np.zeros(n_symbols_per_block, dtype=int) - for n in range(0, n_symbols_per_block): - self._carrier_center_idx[n] = np.argmin(abs(self._xi - self._carrier_centers[n])) + xi = self._XI[0] + np.arange(0, self._M)*self._dxi + self._sum_shifted_carriers = SumOfShiftedWaveforms(carrier_waveform_fun, + carrier_spacing, + n_symbols_per_block, + xi) # Set the right type of continuos spectrum for the NFT routines. @@ -154,23 +164,11 @@ def __init__(self, raise Exception("Unknown contspec_type '"+str(self._contspec_type)+"'") self._opts_fwd.discspec_type = 3 # skip computation of discrete spectrum - # Determine interval [xi_min, xi_max] outside of which the carrier - # waveform is effectively zero. This is essential to speed the - # modulation process when the number of carriers is getting larger. - - vals = self._carrier_waveform_fun(self._xi) - tol = 10 * np.finfo(vals[0]).eps * np.max(np.abs(vals)) - idx = np.argwhere(np.abs(vals)>tol) - i1 = idx[0][0] - i2 = idx[-1][0] - self._xi_min = self._xi[i1] - self._xi_max = self._xi[i2] - def _new_nfspec(self): """Generates an empty NFSpectrum object with the right contspec_type and plotting hints.""" nfspec = NFSpectrum(self._contspec_type, "none") - nfspec.xi = self._xi + nfspec.xi = self._sum_shifted_carriers.grid nfspec.xi_plot_range = np.array([-1.5, 1.5])*self.n_symbols_per_block/2*self._carrier_spacing return nfspec @@ -183,13 +181,7 @@ def modulate(self, symbols): nc = np.size(symbols) assert nc == self.n_symbols_per_block nfspec = self._new_nfspec() - nfspec.cont = np.zeros(self._M, dtype=complex) - for n in range(0, nc): - shifted_xi = self._xi - self._carrier_centers[n] - idx = np.logical_and(shifted_xi>=self._xi_min, shifted_xi<=self._xi_max) - nfspec.cont[idx] = nfspec.cont[idx] + symbols[n]*self._carrier_waveform_fun(shifted_xi[idx]) - for n in range(0, nc): - i = self._carrier_center_idx[n] + nfspec.cont = self._sum_shifted_carriers.generate_waveform(symbols) # Apply the power control factor and, if requested, the power # normalization map @@ -202,15 +194,16 @@ def modulate(self, symbols): # Pre-equalize for the phase change induced by the channel - nfspec.cont *= np.exp(-2.0j*self._xi**2*self._normalized_distance*(self._percent_precompensation/100)) + xi = self._sum_shifted_carriers.grid + nfspec.cont *= np.exp(-2.0j*xi**2*self._normalized_distance*(self._percent_precompensation/100)) # Generate corresponding time domain signal on the extended grid # specified construction rdict = nsev_inverse_wrapper(self._M, nfspec.cont, - self._XI[0], - self._XI[1], + xi[0], + xi[-1], 0, [], [], @@ -264,17 +257,34 @@ def demodulate(self, q): # Remove any remaining phase changes due to the channel that have not # already been removed at the transmitter - nfspec.cont *= np.exp(-2.0j*self._xi**2*self._normalized_distance + xi = self._sum_shifted_carriers.grid + nfspec.cont *= np.exp(-2.0j*xi**2*self._normalized_distance *(1 - self._percent_precompensation/100)) # Recover the symbols symbols = np.zeros(self.n_symbols_per_block, dtype=complex) scl = np.abs(self._carrier_waveform_fun(0.0)) * self._power_control_factor - for n in range(0, self.n_symbols_per_block): - symbols[n] = nfspec.cont[self._carrier_center_idx[n]] - if self._use_power_normalization_map: - symbols[n] = np.sqrt(np.log(np.abs(symbols[n])**2 + 1.0)) * np.exp(1j*np.angle(symbols[n])) - symbols[n] /= scl + if self._matched_filter_percentage>0: + percentage_averaging = self._matched_filter_percentage # % of carrier spacing. It specifies range of 'xi' to average for symbol detection + xi_avg_range = 0.5*self._carrier_spacing*percentage_averaging/100 # xi range for averaging + idx_diff = int(xi_avg_range/self._dxi) # get the index difference + carrier_shape_filter = self._carrier_waveform_fun(xi[int(self._M/2)-idx_diff:int(self._M/2)+idx_diff])# get carrier shape for the filter + carrier_center_idx = self._sum_shifted_carriers._center_idx + for n in range(0, self.n_symbols_per_block): + idx_l = carrier_center_idx[n] - idx_diff # lower index of the xi average range; + idx_u = carrier_center_idx[n] + idx_diff # upper index of the xi average range; + matched_filtered_data = np.multiply(nfspec.cont[idx_l:idx_u],carrier_shape_filter) + symbols[n] = np.mean(matched_filtered_data)/np.mean(carrier_shape_filter) + if self._use_power_normalization_map: + symbols[n] = np.sqrt(np.log(np.abs(symbols[n])**2 + 1.0)) * np.exp(1j*np.angle(symbols[n])) + symbols[n] /= scl + + else: + symbols = self._sum_shifted_carriers.extract_symbols(nfspec.cont) + for n in range(0, np.size(symbols)): + if self._use_power_normalization_map: + symbols[n] = np.sqrt(np.log(np.abs(symbols[n])**2 + 1.0)) * np.exp(1j*np.angle(symbols[n])) + symbols[n] /= self._power_control_factor return symbols, nfspec diff --git a/Modulators/TimeDomainPulseShapingModulator.py b/Modulators/TimeDomainPulseShapingModulator.py new file mode 100644 index 0000000..32135ac --- /dev/null +++ b/Modulators/TimeDomainPulseShapingModulator.py @@ -0,0 +1,86 @@ +# This file is part of NFDMLab. +# +# NFDMLab is free software; you can redistribute it and/or +# modify it under the terms of the version 2 of the GNU General +# Public License as published by the Free Software Foundation. +# +# NFDMLab is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with NFDMLab; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA +# 02111-1307 USA +# +# Contributors: +# Sander Wahls (TU Delft) 2019, 2021 + +import numpy as np +from Modulators import BaseModulator +from Helpers import NFSpectrum, next_pow2 +from Helpers import SumOfShiftedWaveforms + +class TimeDomainPulseShapingModulator(BaseModulator): + """This is a conventional pulse shaping modulator that works in the time + domain. Symbols are embedded in superposed shifted copies of a base waveform + (e.g., a raised cosine).""" + + def __init__(self, + pulse_fun, + pulse_spacing, + requested_normalized_dt, + n_symbols_per_block, + n_guard_symbols, + make_n_samples_pow2=False): + """Constructor. + + Parameters + ---------- + pulse_fun : function + A function that maps numpy.array(float) to numpy.array(complex). The + inputs of this function are vectors of time points t. + The outputs are vectors with the values of the base waveform at + the corresponding t. + pulse_spacing : float + Spacing between two consecutive copies of the base waveform. + required_normalized_dt : float + The modulator will choose the number of time domain samples such + that the time step, in normalized units, does not exceed this value. + n_symbols_per_block : int + Number of shifted copies of the base waveform. + n_guard_symbols : int + Guard symbols are zero symbols that are added to the transmission to + avoid interference between consecutive blocks of data during + transmission. The duration of one block generated by this modulator + is pulse_duration*(n_symbols_per_block + n_guard_symbols). + make_n_samples_pow2 : bool + If True, the modular will choose the sampling interval such that the + number of samples generated by the modulator is a power of two. + """ + + self._n_samples = int(pulse_spacing/requested_normalized_dt * (n_symbols_per_block + n_guard_symbols)) + T1 = self.n_samples/2*requested_normalized_dt + T0 = -T1 + if make_n_samples_pow2: + self._n_samples = next_pow2(self._n_samples) + t = np.linspace(T0, T1, self.n_samples) + self._sum_pulses = SumOfShiftedWaveforms(pulse_fun, + pulse_spacing, + n_symbols_per_block, + t) + self._normalized_dt = (T1 - T0)/self.n_samples + self._n_symbols_per_block = n_symbols_per_block + + def modulate(self, symbols): + # Docstring is inherited from base class. + nc = np.size(symbols) + assert nc == self.n_symbols_per_block + q_tx = self._sum_pulses.generate_waveform(symbols) + return q_tx, NFSpectrum('none', 'none') + + def demodulate(self, q_rx): + # Docstring is inherited from base class. + symbols = self._sum_pulses.extract_symbols(q_rx) + return symbols, NFSpectrum('none', 'none') diff --git a/Modulators/__init__.py b/Modulators/__init__.py index 7f74035..f8a1434 100644 --- a/Modulators/__init__.py +++ b/Modulators/__init__.py @@ -15,8 +15,9 @@ # 02111-1307 USA # # Contributors: -# Sander Wahls (TU Delft) 2019 +# Sander Wahls (TU Delft) 2019, 2021 from Modulators.BaseModulator import BaseModulator from Modulators.ContSpecModulator import ContSpecModulator -from Modulators.DiscSpecModulator import DiscSpecModulator \ No newline at end of file +from Modulators.DiscSpecModulator import DiscSpecModulator +from Modulators.TimeDomainPulseShapingModulator import TimeDomainPulseShapingModulator diff --git a/QualityAssessment/ErrorVectorMagnitude.py b/QualityAssessment/ErrorVectorMagnitude.py index 0383e98..c1c51fe 100644 --- a/QualityAssessment/ErrorVectorMagnitude.py +++ b/QualityAssessment/ErrorVectorMagnitude.py @@ -37,7 +37,17 @@ def in_dB(self, correct_symbols, recovered_symbols): ------- float 10*log10(P_error / P_reference) + + References + EVM calculations reference + 1. "dB or not dB?", Application Note 1MA98 , Rohde and Schwarz Gmbh, page-29 + 2. M. Vigilante, E. McCune and P. Reynaert, "To EVM or Two EVMs?: An Answer to the Question," in IEEE Solid-State Circuits Magazine, vol. 9, no. 3, pp. 36-39, Summer 2017. + from reference 2, it is clear that relation between E_rms (dB), E_max (dB) for different modualtion formats (on page 34) hold only if we use 20*log10 definition. + EVM (in dB) = 20 log10 (EVM (in percentage)/100) + + """ + num = np.linalg.norm(correct_symbols - recovered_symbols) den = np.linalg.norm(correct_symbols) return 20*np.log10(num/den) # The *2 corrects the sqrt in linalg.norm