diff --git a/src/single-core-regen/generate_signal.py b/src/single-core-regen/generate_signal.py index 2028f53..4295971 100644 --- a/src/single-core-regen/generate_signal.py +++ b/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']) - sig = pypho.lasmod(py_glova, power=0, Df=0, theta=np.pi/4) - modulator = pypho.arbmod(py_glova) + laser = pypho.lasmod(py_glova, power=config['signal']['laser_power'], 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']) 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') - esig = esigsrc(bitsequence=ones)*np.sqrt(2) + cw = laser() + cw[0]['E']*=np.sqrt(2) - source_signal = sig(esig) - - 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)) + source_signal = modulator(E=cw, symbols=(symbols_x, symbols_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']) - c_params = pypho.cfiber.ParamsWrapper.from_fiber(py_fiber) + if config['fiber']['birefsteps'] > 0: + 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 - # 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) + length_loop(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[1,0].plot(xax,np.abs(cdata.E_in[1])) - # axs[0,1].plot(xax,np.abs(cdata.E_out[0])) - # axs[1,1].plot(xax,np.abs(cdata.E_out[1])) + # axs[0].plot(xax,np.abs(cdata.E_in[0])**2) + # axs[0].plot(xax,np.abs(cdata.E_out[0])**2) + # axs[0].set_title("x") + # 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() diff --git a/src/single-core-regen/signal_generation.ini b/src/single-core-regen/signal_generation.ini index 9d2ae16..16348c9 100644 --- a/src/single-core-regen/signal_generation.ini +++ b/src/single-core-regen/signal_generation.ini @@ -11,8 +11,8 @@ flags = "FFTW_PATIENT" nthreads = 32 [fiber] -length = 1000 -gamma = 1.14 +length = 100000 +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 -pulse_shape = "gauss_rz" +mod_depth = 1 +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" \ No newline at end of file