pam4 signal modulation

src/single-core-regen/generate_signal.py

@@ -2,7 +2,9 @@ import configparser
 from datetime import datetime
 import hashlib
 from pathlib import Path
+import time
 from matplotlib import pyplot as plt
+import scipy
 import numpy as np
 
 from rich import inspect
@@ -29,14 +31,17 @@ alpha = 0.2
 D = 17
 S = 0
 birefsteps = 1
-birefseed = 0xC0FFEE
+; birefseed = 0xC0FFEE
 
 [signal]
 modulation = "pam"
 mod_order = 4
-seed = 0xC0FFEE
+mod_depth = 0.5
+laser_power = 0
+; seed = 0xC0FFEE
 pulse_shape = "gauss_rz"
 fwhm = 0.33
+; jitter_seed = 0xC0FFEE
 
 [script]
 data_dir = "{str((Path.home()/".pypho"/"data"))}"
@@ -67,6 +72,73 @@ def get_config(config_file=None):
             #     conf[section][key] = config[section][key].strip('"')
     return conf
 
+class pam_generator:
+    def __init__(self, glova, mod_order=None, mod_depth=0.5, pulse_shape='gauss', fwhm=0.33, seed=None) -> None:
+        self.glova = glova
+        self.pulse_shape = pulse_shape
+        self.modulation_depth = mod_depth
+        self.mod_order = mod_order
+        self.fwhm = fwhm
+        self.seed = seed
+
+    def __call__(self, E, symbols, max_jitter=0):
+        symbols_x = symbols[0]/(self.mod_order or np.max(symbols[0]))
+        symbols_y = symbols[1]/(self.mod_order or np.max(symbols[1]))
+
+        
+        diffs_x = np.diff(symbols_x, prepend=symbols_x[0])
+        diffs_y = np.diff(symbols_y, prepend=symbols_y[0])
+
+        max_jitter = int(round(max_jitter*self.glova.sps))
+        digital_x = self.generate_digital_signal(diffs_x, max_jitter)
+        digital_y = self.generate_digital_signal(diffs_y, max_jitter)
+
+        digital_x = np.pad(digital_x, (0, self.glova.sps//2), 'constant', constant_values=(0,0))
+        digital_y = np.pad(digital_y, (0, self.glova.sps//2), 'constant', constant_values=(0,0))
+
+        if self.pulse_shape == 'gauss':
+            wavelet = self.gauss(oversampling=6)
+        else:
+            raise ValueError(f"Unknown pulse shape: {self.pulse_shape}")
+
+        # create analog signal of diff of symbols
+        E_x = np.convolve(digital_x, wavelet)
+        E_y = np.convolve(digital_y, wavelet)
+
+        # convert to pam and set modulation depth (scale and move up such that 1 stays at 1)
+        E_x = np.cumsum(E_x)*self.modulation_depth + (1-self.modulation_depth)
+        E_y = np.cumsum(E_y)*self.modulation_depth + (1-self.modulation_depth)
+
+        # cut off the wavelet tails
+        E_x = E_x[self.glova.sps//2+len(wavelet)//2-1:-len(wavelet)//2]
+        E_y = E_y[self.glova.sps//2+len(wavelet)//2-1:-len(wavelet)//2]
+
+        # modulate the laser
+        np.multiply(E[0]['E'][0],E_x, out=E[0]['E'][0])
+        np.multiply(E[0]['E'][1],E_y, out=E[0]['E'][1])
+
+        return E
+    
+    def generate_digital_signal(self, symbols, max_jitter=0):
+        rs = np.random.RandomState(self.seed)
+        signal = np.zeros(self.glova.nos*self.glova.sps)
+        for index in range(self.glova.nos):
+            jitter = max_jitter != 0 and rs.randint(-max_jitter, max_jitter)
+            signal_index = index*self.glova.sps + jitter
+            if signal_index < 0:
+                continue
+            if signal_index >= len(signal):
+                continue
+            signal[signal_index] = symbols[index]
+        return signal
+    
+    def gauss(self, oversampling=1):
+        sample_points = np.linspace(-oversampling*self.glova.sps, oversampling*self.glova.sps, oversampling*2*self.glova.sps, endpoint=True)
+        sigma = self.fwhm/(1*np.sqrt(2*np.log(2)))*self.glova.sps
+        pulse = 1/(sigma*np.sqrt(2*np.pi))*np.exp(-np.square(sample_points)/(2*np.square(sigma)))
+        return pulse
+
+
 def initialize_fiber_and_data(config, input_data_override=None):
     py_glova = pypho.setup(
         nos = config['glova']['nos'],
@@ -84,23 +156,21 @@ def initialize_fiber_and_data(config, input_data_override=None):
     if input_data_override is not None:
         c_data.E_in = input_data_override
     else:
+        config['signal']['seed'] = config['signal'].get('seed',(int(time.time()*1000))%2**32)
+        config['signal']['jitter_seed'] = config['signal'].get('jitter_seed', (int(time.time()*1000))%2**32)
         symbolsrc = pypho.symbols(py_glova, py_glova.nos, pattern='ones', seed=config['signal']['seed'])
-        esigsrc = pypho.signalsrc(py_glova, pulseshape=config['signal']['pulse_shape'], fwhm=config['signal']['fwhm'])
+        laser = pypho.lasmod(py_glova, power=config['signal']['laser_power'], Df=0, theta=np.pi/4)
-        sig = pypho.lasmod(py_glova, power=0, Df=0, theta=np.pi/4)
+        modulator = pam_generator(py_glova, mod_depth=config['signal']['mod_depth'], pulse_shape=config['signal']['pulse_shape'], fwhm=config['signal']['fwhm'], seed=config['signal']['jitter_seed'])
-        modulator = pypho.arbmod(py_glova)
 
         symbols_x = symbolsrc(pattern='random', p1=config['signal']['mod_order'])
         symbols_y = symbolsrc(pattern='random', p1=config['signal']['mod_order'])
+        symbols_x[:3] = 0
+        symbols_y[:3] = 0
 
-        ones = symbolsrc(pattern='ones')
+        cw = laser()
-        esig = esigsrc(bitsequence=ones)*np.sqrt(2)
+        cw[0]['E']*=np.sqrt(2)
 
-        source_signal = sig(esig)
+        source_signal = modulator(E=cw, symbols=(symbols_x, symbols_y))
-
-        constpts_x = [np.linspace(1/config['signal']['mod_order'], 1, config['signal']['mod_order'], endpoint=True, dtype=np.complex128)]
-        constpts_y = [np.linspace(1/config['signal']['mod_order'], 1, config['signal']['mod_order'], endpoint=True, dtype=np.complex128)]
-
-        source_signal = modulator(E=source_signal, symbols = (symbols_x, symbols_y), constpoints= (constpts_x, constpts_y))
 
         c_data.E_in = source_signal[0]['E']
 
@@ -112,8 +182,10 @@ def initialize_fiber_and_data(config, input_data_override=None):
         D=config['fiber']['d'],
         S=config['fiber']['s'],
     )
-    py_fiber.birefarray = pypho.birefringence_segment.create_pmd_fibre(py_fiber.l, py_fiber.l/config['fiber']['birefsteps'], 0, config['fiber']['birefseed'])
+    if config['fiber']['birefsteps'] > 0:
-    c_params = pypho.cfiber.ParamsWrapper.from_fiber(py_fiber)
+        config['fiber']['birefseed'] = config['fiber'].get('birefseed', (int(time.time()*1000))%2**32)
+        py_fiber.birefarray = pypho.birefringence_segment.create_pmd_fibre(py_fiber.l, py_fiber.l/config['fiber']['birefsteps'], 0, config['fiber']['birefseed'])
+    c_params = pypho.cfiber.ParamsWrapper.from_fiber(py_fiber, max_step=1e3)
     c_fiber = pypho.cfiber.FiberWrapper(c_data, c_params, c_glova)
 
     return c_fiber, c_data
@@ -143,14 +215,17 @@ def save_data(data, config):
         f"D = {config['fiber']['d']}",
         f"S = {config['fiber']['s']}",
         f"birefsteps = {config['fiber']['birefsteps']}",
-        f"birefseed = {config['fiber']['birefseed']}",
+        f"birefseed = {config['fiber'].get('birefseed', 'not set')}",
         "",
         "[signal]",
         f'modulation = "{config['signal']['modulation']}"',
         f"mod_order = {config['signal']['mod_order']}",
-        f"seed = {config['signal']['seed']}",
+        f"mod_depth = {config['signal']['mod_depth']}",
+        f"seed = {config['signal'].get('seed', 'not set')}",
         f'pulse_shape = "{config["signal"]["pulse_shape"]}"',
         f"fwhm = {config['signal']['fwhm']}",
+        f"max_jitter = {config['signal']['max_jitter']}",
+        f"jitter_seed = {config['signal'].get('jitter_seed', 'not set')}",
         "",
         "[script]",
         f'data_dir = "{config["script"]["data_dir"]}"',
@@ -180,37 +255,65 @@ def save_data(data, config):
     print("Saved config to", data_dir / lookup_file)
     print("Saved data to", save_dir / f"{config_hash}.npy")
 
+def length_loop(config):
+    # loop over lengths
+    # lengths = [100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 90000, 100000]
+    ranges = [1, 10, 100, 1000, 10000, 100000]
+    lengths = [list(range(length_range,10*length_range, length_range)) for length_range in ranges]
+    lengths = [length for length_range in lengths for length in length_range]
 
+    lengths = sorted(lengths)
+    input_override = None
+    for lind, length in enumerate(lengths):
+        # print(f"\nGenerating data for fiber length {length}")
+        if lind > 0:
+            # set the length to the difference between the current and previous length -> incremental
+            length = lengths[lind] - lengths[lind-1]
+        print(f"\nGenerating data for fiber length {lengths[lind]}m -> {length}m")
+        config['fiber']['length'] = length
+        # set the input data to the output data of the previous run
+        cfiber, cdata = initialize_fiber_and_data(config, input_data_override=input_override)
+        
+        if lind == 0:
+            cdata_orig = cdata
+        
+        energy_in = np.sum(np.abs(cdata.E_in[0])**2 + np.abs(cdata.E_in[1])**2)/(cfiber.glova.symbolrate*cfiber.glova.sps)
+        print(f"Energy in:  {energy_in} J")
+        
+        cfiber()
+
+        energy_out = np.sum(np.abs(cdata.E_out[0])**2 + np.abs(cdata.E_out[1])**2)/(cfiber.glova.symbolrate*cfiber.glova.sps)
+        print(f"Energy out: {energy_out} J")
+
+        input_override = cdata.E_out
+        cdata.E_in = cdata_orig.E_in
+        config['fiber']['length'] = lengths[lind]
+        save_data(cdata, config)
 
 if __name__ == "__main__":
     config = get_config()
 
-    # loop over lengths
+    length_loop(config)
-    # lengths = [1000, 2000, 4000, 8000, 16000, 32000, 64000, 128000]
+
-    lengths = [1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 90000, 100000]
-    lengths = sorted(lengths)
-    input_override = None
-    for lind, length in enumerate(lengths):
-        print(f"\nGenerating data for fiber length {length}")
-        # if lind > 0:
-        #     # set the length to the difference between the current and previous length -> incremental
-        #     length = lengths[lind] - lengths[lind-1]
-        # print(f"\nGenerating data for fiber length {lengths[lind]}m -> {length}m")
-        config['fiber']['length'] = length
-        # set the input data to the output data of the previous run
-        cfiber, cdata = initialize_fiber_and_data(config, input_data_override=input_override)
-        cfiber()
-        # input_override = cdata.E_out
-        save_data(cdata, config)
     # cfiber, cdata = initialize_fiber_and_data(config)
     # cfiber()
+
+    # energy_in = np.sum(np.abs(cdata.E_in[0])**2 + np.abs(cdata.E_in[1])**2)/(cfiber.glova.symbolrate*cfiber.glova.sps)
+    # energy_out = np.sum(np.abs(cdata.E_out[0])**2 + np.abs(cdata.E_out[1])**2)/(cfiber.glova.symbolrate*cfiber.glova.sps)
+    # print(f"Energy in:  {energy_in} J")
+    # print(f"Energy out: {energy_out} J")
     # save_data(cdata, config)
 
-    # fig, axs = plt.subplots(2, 2, sharex=True, sharey=True)
+    # fig, axs = plt.subplots(2, 1, sharex=True, sharey=True)
     # xax = np.linspace(0, cfiber.glova.nos, cfiber.glova.nos*cfiber.glova.sps)-0.5
-    # axs[0,0].plot(xax,np.abs(cdata.E_in[0]))
+    # axs[0].plot(xax,np.abs(cdata.E_in[0])**2)
-    # axs[1,0].plot(xax,np.abs(cdata.E_in[1]))
+    # axs[0].plot(xax,np.abs(cdata.E_out[0])**2)
-    # axs[0,1].plot(xax,np.abs(cdata.E_out[0]))
+    # axs[0].set_title("x")
-    # axs[1,1].plot(xax,np.abs(cdata.E_out[1]))
+    # axs[0].grid()
+
+    # axs[1].plot(xax,np.abs(cdata.E_in[1])**2)
+    # axs[1].plot(xax,np.abs(cdata.E_out[1])**2)
+    # axs[1].set_title("y")
+    # axs[1].grid()
 
     # plt.show()

src/single-core-regen/signal_generation.ini

@@ -11,8 +11,8 @@ flags = "FFTW_PATIENT"
 nthreads = 32
 
 [fiber]
-length = 1000
+length = 100000
-gamma = 1.14
+gamma = 0
 alpha = 0.2
 D = 17
 ; 1700 ps/nm/km @ 1 km length is equivalent to 17 ps/nm/km @ 100 km length but the simulation is way faster
@@ -23,9 +23,13 @@ birefsteps = 0
 [signal]
 modulation = "pam"
 mod_order = 4
-; seed = 0xC0FFEE
+mod_depth = 1
-pulse_shape = "gauss_rz"
+laser_power = 0
+seed = 0xC0FFEE
+pulse_shape = "gauss"
 fwhm = 0.33
+max_jitter=0.02
+; jitter_seed = 0xC0FFEE
 
 [script]
 data_dir = "/home/suppl/.pypho/data"