From 9093fb245217c2142b722a106a0e47f28db934fb Mon Sep 17 00:00:00 2001 From: Magnus Fischer Date: Mon, 2 Feb 2026 12:12:33 +0000 Subject: [PATCH] Revert "Merge branch 'main' of https://cau-git.rz.uni-kiel.de/nt/mitarbeiter/silas/imdd_simulation" This reverts commit 798a0ca3b30cc7e6e4aac1d42b7ca34f8343c87d --- Classes/00_signals/Signal.m | 6 +- Classes/02_optical/DP_Fiber.m | 19 +- Classes/02_optical/Optical_Demultiplex.m | 45 ++- Classes/02_optical/Optical_Multiplex.m | 45 +-- Classes/02_optical/dp_fiber_lib/CNLSE_plain.m | 189 ++++------ Classes/02_optical/dp_fiber_lib/canUseGPU.m | 20 - .../02_optical/dp_fiber_lib/getNLstepsize.m | 40 +- Classes/02_optical/dp_fiber_lib/lin_step.m | 64 ++-- Classes/04_DSP/Equalizer/EQ.m | 48 +++ Classes/04_DSP/Equalizer/ML_MLSE.m | 2 +- Classes/04_DSP/Timing_Recovery.m | 47 --- Classes/DataBaseHandler/DBHandler.m | 13 +- Functions/EQ_structures/ml_mlse.m | 15 +- .../Dissertation/mach_zehnder_modulator.m | 219 ----------- Functions/mat2tikz_improved.m | 23 -- Libs/wesanderson_colors/WesPalette.m | 127 ++----- .../minimal_example_wespalette.m | 4 +- projects/FSO_transmission/first_analysis.m | 252 ------------- .../final/Copy_of_FIGURE_WAVELENGTH.m | 26 +- .../final/FIGURE_introduction.m | 156 +++----- .../Auswertung_JLT/run_dsp_from_db.m | 5 +- .../ML_based_MLSE/analyze_filter_length.m | 16 +- .../theoretic_channel_evaluation.m | 60 +-- projects/Messung_Zürich/minimal_matched.m | 60 --- projects/WDM/WDM_model_10km_queue.m | 2 +- test/bayesopt_ffe_tuning.m | 219 ----------- test/bitwise_demapping_pam6.m | 43 --- test/duobinary_minimal_example.m | 3 + test/gpu_cpu_comparison.m | 349 ------------------ test/gpu_processing_dpfiber.m | 265 ------------- test/gpu_processing_test.m | 70 ---- 31 files changed, 385 insertions(+), 2067 deletions(-) delete mode 100644 Classes/02_optical/dp_fiber_lib/canUseGPU.m delete mode 100644 Classes/04_DSP/Timing_Recovery.m delete mode 100644 Functions/Theory/Dissertation/mach_zehnder_modulator.m delete mode 100644 Functions/mat2tikz_improved.m delete mode 100644 projects/FSO_transmission/first_analysis.m delete mode 100644 projects/Messung_Zürich/minimal_matched.m delete mode 100644 test/bayesopt_ffe_tuning.m delete mode 100644 test/bitwise_demapping_pam6.m delete mode 100644 test/gpu_cpu_comparison.m delete mode 100644 test/gpu_processing_dpfiber.m delete mode 100644 test/gpu_processing_test.m diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index f09a03e..e06c41f 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -172,9 +172,11 @@ classdef Signal hold on; if isempty(options.color) - plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1); + % plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1); + plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1); else - plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color); + % plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color); + plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Color',options.color); end % 2 c) % - xlabel if not already here: time in readable format (1 ms and not 1e-3 s) diff --git a/Classes/02_optical/DP_Fiber.m b/Classes/02_optical/DP_Fiber.m index c1f08ea..a5f3dce 100644 --- a/Classes/02_optical/DP_Fiber.m +++ b/Classes/02_optical/DP_Fiber.m @@ -24,8 +24,6 @@ classdef DP_Fiber SS_dzmax % [m] max dz (adaptive SSFM) SS_dzmin % [m] min dz (adaptive SSFM) n_waveplates % number of PMD waveplates - useGPU % GPU acceleration: true, false, or 'auto' (default) - useSingle % Use single precision on GPU (default: false) % ---- Internal state (persistent between calls) ---- state % struct mirroring legacy 'state' @@ -58,8 +56,6 @@ classdef DP_Fiber options.SS_dzmax = 2e4 % m options.SS_dzmin = 100 % m options.n_waveplates = 100 - options.useGPU = 'auto' % 'auto', true, or false - options.useSingle = false % single precision GPU end % Copy provided options into properties @@ -212,7 +208,7 @@ classdef DP_Fiber % Frequency-dependent PMD phase term (legacy form) st.brf.db0 = (R.rand(st.wave_plates,1)*2*pi - pi) * brf_multiplier; st.brf.db1 = sqrt(3*pi/8)*(st.dgd/obj.fa)/st.wave_plates .* st.omega; - st.brf.simdgd = 0; + st.brf.simdgd = 0; % cumsum used in legacy only for debug; keep compatibility variable: ~cumsum(st.brf.db0); % no-op to mirror legacy path @@ -232,18 +228,7 @@ classdef DP_Fiber x_in = signal_in(:,1).'; y_in = signal_in(:,2).'; - % Determine GPU usage - if ischar(obj.useGPU) || isstring(obj.useGPU) - if strcmpi(obj.useGPU, 'auto') - gpuFlag = []; % Let CNLSE_plain auto-detect - else - error('DP_Fiber:InvalidGPU', 'useGPU must be true, false, or ''auto'''); - end - else - gpuFlag = logical(obj.useGPU); - end - - [x_out, y_out, obj.state] = CNLSE_plain(x_in, y_in, obj.state, gpuFlag, obj.useSingle); + [x_out, y_out, obj.state] = CNLSE_plain(x_in, y_in, obj.state); obj.state.propagated_length = obj.state.propagated_length + obj.state.L; diff --git a/Classes/02_optical/Optical_Demultiplex.m b/Classes/02_optical/Optical_Demultiplex.m index 8e272d2..7d7c39e 100644 --- a/Classes/02_optical/Optical_Demultiplex.m +++ b/Classes/02_optical/Optical_Demultiplex.m @@ -42,7 +42,7 @@ classdef Optical_Demultiplex < handle function signalclasses_out = process(obj, signalclass_in) - % ---- Infer wavelength: either given or from input total signal + % ---- Infer wavelength: either given or from input total signal if isempty(obj.wavelengthplan) obj.wavelengthplan = signalclass_in.lambda; %meter else @@ -81,7 +81,7 @@ classdef Optical_Demultiplex < handle obj signal_in end - + w = obj.fs_out ./ obj.fs_in ; blocklen_in = length(signal_in); blocklen_out = w*blocklen_in; @@ -119,31 +119,30 @@ classdef Optical_Demultiplex < handle N = size(lo,1); C = size(lo,2); - % ---- VECTORIZED: Process all channels in parallel ---- - % Batched FFT operates on each column simultaneously on GPU + x_envelopes = zeros(N, C, 'like', signal_in); + y_envelopes = zeros(N, C, 'like', signal_in); - % Extract polarization signals - s1 = signal_in(:,1); % X polarization [N×1] - s2 = signal_in(:,2); % Y polarization [N×1] + s1 = signal_in(:,1); + s2 = signal_in(:,2); - % Broadcast signal to all channels and multiply with LO - % s1, s2 are [N×1], lo is [N×C] → result is [N×C] - x_mixed = att .* s1 .* lo; % [N×C] - y_mixed = att .* s2 .* lo; % [N×C] - - % Batched FFT: each column computed in parallel - x_freq = fft(x_mixed); % [N×C] - y_freq = fft(y_mixed); % [N×C] - - % Apply filter (H is [N×1], broadcasts across columns) - x_filtered = x_freq .* H; % [N×C] - y_filtered = y_freq .* H; % [N×C] - - % Batched IFFT - x_envelopes = ifft(x_filtered); % [N×C] - y_envelopes = ifft(y_filtered); % [N×C] + % Reusable work buffers (avoid reallocations) + wrk_time = zeros(N,1, 'like', signal_in); + wrk_freq = zeros(N,1, 'like', signal_in); + for c = 1:C + % ---- X branch ---- + wrk_time(:) = att .* s1 .* lo(:,c); % N×1 + wrk_freq(:) = fft(wrk_time); % N×1 + wrk_freq(:) = wrk_freq .* H; % N×1 + x_envelopes(:,c) = ifft(wrk_freq); % N×1 + % ---- Y branch ---- + wrk_time(:) = att .* s2 .* lo(:,c); + wrk_freq(:) = fft(wrk_time); + wrk_freq(:) = wrk_freq .* H; + y_envelopes(:,c) = ifft(wrk_freq); + end + end end end diff --git a/Classes/02_optical/Optical_Multiplex.m b/Classes/02_optical/Optical_Multiplex.m index 7fc8a33..1d722b0 100644 --- a/Classes/02_optical/Optical_Multiplex.m +++ b/Classes/02_optical/Optical_Multiplex.m @@ -1,10 +1,10 @@ classdef Optical_Multiplex < handle % Takes a cell array of signals - % returns a total field signal - % WDM spacing is given in wavelength plan OR via delta_F - + % returns a total field signal + % WDM spacing is given in wavelength plan OR via delta_F + % The grid is stored in the output signal -> the demux will ideally - % look this up and use this as the demux frequencies... + % look this up and use this as the demux frequencies... % signal_cell = {Opt_sig_1, Opt_sig_2}; % Opt_sig_wdm = Optical_Multiplex("fs_in",Opt_sig.fs,"fs_out",4*Opt_sig.fs,... @@ -117,7 +117,7 @@ classdef Optical_Multiplex < handle % adapt frequency shifts to match the FFT grid! Find nearest grid point [glitch(o),pos] = min(abs( freqaxis-obj.df_T(o) )); obj.df_T(o) = freqaxis(pos); - + polrots = [polrots, data_in{o}.polrot]; end @@ -139,36 +139,25 @@ classdef Optical_Multiplex < handle x_envelopes = NaN([blocklen_out N]); y_envelopes = x_envelopes; - - % ---- OPTIMIZED: Pre-compute all LO phases as [blocklen_out × N] matrix ---- - time_idx = (0:blocklen_out-1).'; % [blocklen_out × 1] - lo_phases = mod(2*pi * time_idx * obj.df_T / obj.fs_out, 2*pi); % [blocklen_out × N] - lo_all = cos(lo_phases) + 1i*sin(lo_phases); % [blocklen_out × N] - - % Collect all resampled signals first (still requires loop due to cell array) - x_signals = zeros(blocklen_out, N); - y_signals = zeros(blocklen_out, N); - + for o = 1:N - data_in_resampled = data_in{o}.resample("fs_out", obj.fs_out); - x_signals(:, o) = data_in_resampled.signal(:, 1); - y_signals(:, o) = data_in_resampled.signal(:, 2); + + pha = mod(2*pi*(0:blocklen_out-1)*obj.df_T(o)/obj.fs_out,2*pi).'; + lo = cos(pha)+1i*sin(pha); + data_in_resampled = data_in{o}.resample("fs_out",obj.fs_out); + + res_env = ifft(fft(data_in_resampled.signal(:,1)).*H); + x_envelopes(:,o) = att.*res_env.*lo; + + res_env = ifft(fft(data_in_resampled.signal(:,2)).*H); + y_envelopes(:,o) = att.*res_env.*lo; + end - % ---- VECTORIZED: Batched FFT/IFFT for all channels ---- - % Apply filter to all channels at once - x_filtered = ifft(fft(x_signals) .* H); % [blocklen_out × N] - y_filtered = ifft(fft(y_signals) .* H); % [blocklen_out × N] - - % Apply attenuation and LO shift to all channels - x_envelopes = att .* x_filtered .* lo_all; % [blocklen_out × N] - y_envelopes = att .* y_filtered .* lo_all; % [blocklen_out × N] - data_out = data_in_resampled; data_out.signal = [sum(x_envelopes,2), sum(y_envelopes,2)]; data_out.lambda = obj.lambda_T; data_out.polrot = polrots; - end end diff --git a/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m b/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m index 463adca..cfa7294 100644 --- a/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m +++ b/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m @@ -1,72 +1,41 @@ -function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state,useGPU,useSingle) +function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state) + -% GPU auto-detection if not specified -if nargin < 4 || isempty(useGPU) - useGPU = canUseGPU(); -end -if nargin < 5 || isempty(useSingle) - useSingle = false; -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% pre calculations -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % pre calculations + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -state.common_beta=struct('X',0,'Y',0); + state.common_beta=struct('X',0,'Y',0); -for n=1:2 - % get current polarization name and contrary one - curPol = state.polNames{n}; + for n=1:2 + % get current polarization name and contrary one + curPol = state.polNames{n}; - % extend linear transfer function depending on beta values for the current polarization - % Was ist der Sinn dieser komischen beta notation? zB. state.beta.X = [0.3142 0 -9.1105e-28 5.1068e-41] - for n_beta = 1:length(state.beta.(curPol)) - state.common_beta.(curPol) = state.common_beta.(curPol) + state.beta.(curPol)(n_beta) * (state.omega).^(n_beta-1) / factorial(n_beta-1); + % extend linear transfer function depending on beta values for the current polarization + % Was ist der Sinn dieser komischen beta notation? zB. state.beta.X = [0.3142 0 -9.1105e-28 5.1068e-41] + for n_beta = 1:length(state.beta.(curPol)) + state.common_beta.(curPol) = state.common_beta.(curPol) + state.beta.(curPol)(n_beta) * (state.omega).^(n_beta-1) / factorial(n_beta-1); + end + + %opt_out_struct.(curPol)=opt_in_struct.(curPol).envelope; end - %opt_out_struct.(curPol)=opt_in_struct.(curPol).envelope; -end + beta_const = state.beta.('X')(1); + beta_1 = state.beta.('X')(2); + beta_2 = state.beta.('X')(3); + beta_3 = state.beta.('X')(4); + deltaomega = state.omega; + beta_x = beta_const + beta_1 * deltaomega + 1/2 * beta_2 * deltaomega.^2 + 1/6 *beta_3 * deltaomega.^3; + +% opt_in_x = gpuArray(opt_in_x); +% opt_in_y = gpuArray(opt_in_y); -beta_const = state.beta.('X')(1); -beta_1 = state.beta.('X')(2); -beta_2 = state.beta.('X')(3); -beta_3 = state.beta.('X')(4); -deltaomega = state.omega; -beta_x = beta_const + beta_1 * deltaomega + 1/2 * beta_2 * deltaomega.^2 + 1/6 *beta_3 * deltaomega.^3; - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% GPU Transfer (if enabled) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -if useGPU - % Convert to single precision if requested (faster on most GPUs) - if useSingle - opt_in_x = gpuArray(single(opt_in_x)); - opt_in_y = gpuArray(single(opt_in_y)); - state.common_beta.X = gpuArray(single(state.common_beta.X)); - state.common_beta.Y = gpuArray(single(state.common_beta.Y)); - state.brf.db1 = gpuArray(single(state.brf.db1)); - state.brf.db0 = gpuArray(single(state.brf.db0)); - for k = 1:numel(state.brf.matR) - state.brf.matR{k} = gpuArray(single(state.brf.matR{k})); - end - else - opt_in_x = gpuArray(opt_in_x); - opt_in_y = gpuArray(opt_in_y); - state.common_beta.X = gpuArray(state.common_beta.X); - state.common_beta.Y = gpuArray(state.common_beta.Y); - state.brf.db1 = gpuArray(state.brf.db1); - state.brf.db0 = gpuArray(state.brf.db0); - for k = 1:numel(state.brf.matR) - state.brf.matR{k} = gpuArray(state.brf.matR{k}); - end - end -end - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Split Step Method -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Split Step Method + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [opt_out_x,opt_out_y] = split_step_loop(state.L,opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin,... % state.lin_z_test,state.corr_length,state.n_plates_done,state.missing_dz,state.brf,state.common_beta,.... @@ -75,33 +44,35 @@ end % [opt_out_x,opt_out_y] = split_step_loop_mex(state.L,opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin,... % state.lin_z_test,state.corr_length,state.n_plates_done,state.missing_dz,state.brf,state.common_beta,.... % state.chi,state.manakov,state.beat_len); + + % get nonlinear step size + [state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); + %[state.dz] = getNLstepsize_original(state,opt_out_struct); + + state.n_step = 0; + state.z_prop = 0; + state.test_dz = []; + state.powers = []; -% get nonlinear step size -[state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); -%[state.dz] = getNLstepsize_original(state,opt_out_struct); + tic -state.n_step = 0; -state.z_prop = 0; -state.test_dz = []; -state.powers = []; + while state.z_prop < state.L -while state.z_prop < state.L + % reduce step length (dz) if we are to overshoot the fiber length + % (L) in the next step + if state.z_prop + state.dz > state.L + state.dz = state.L - state.z_prop; + end + + % update step number (n) + state.n_step=state.n_step+1; - % reduce step length (dz) if we are to overshoot the fiber length - % (L) in the next step - if state.z_prop + state.dz > state.L - state.dz = state.L - state.z_prop; - end + % append current step length to logbook (dzs) + state.dzs(state.n_step)=state.dz; - % update step number (n) - state.n_step=state.n_step+1; - - % append current step length to logbook (dzs) - state.dzs(state.n_step)=state.dz; - - - % half linear step - [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... + + % half linear step + [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,... state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X]... @@ -111,12 +82,12 @@ while state.z_prop < state.L state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X); - % complete nonlinear step - - [opt_in_x,opt_in_y] = nl_step(opt_in_x,opt_in_y, state.dz, state.gamma, state.chi, state.manakov, state.beat_len ,state.alpha_lin.X, state.alpha_lin.Y); - - % half linear step - [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... + % complete nonlinear step + + [opt_in_x,opt_in_y] = nl_step(opt_in_x,opt_in_y, state.dz, state.gamma, state.chi, state.manakov, state.beat_len ,state.alpha_lin.X, state.alpha_lin.Y); + + % half linear step + [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,... state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X]... @@ -126,34 +97,24 @@ while state.z_prop < state.L state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X); - % get nonlinear step size - [state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); - %[state.dz] = getNLstepsize_original(state,opt_out_struct); + % get nonlinear step size + [state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); + %[state.dz] = getNLstepsize_original(state,opt_out_struct); - -end - - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% GPU Gather (if enabled) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -if useGPU - opt_out_x = gather(opt_in_x); - opt_out_y = gather(opt_in_y); - - % Gather state arrays back to CPU - state.common_beta.X = gather(state.common_beta.X); - state.common_beta.Y = gather(state.common_beta.Y); - state.brf.db1 = gather(state.brf.db1); - state.brf.db0 = gather(state.brf.db0); - for k = 1:numel(state.brf.matR) - state.brf.matR{k} = gather(state.brf.matR{k}); + end -else - opt_out_x = opt_in_x; - opt_out_y = opt_in_y; -end + + toc + + + + + opt_out_x = (opt_in_x); + opt_out_y = (opt_in_y); + +% opt_out_x = gather(opt_in_x); +% opt_out_y = gather(opt_in_y); + end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/canUseGPU.m b/Classes/02_optical/dp_fiber_lib/canUseGPU.m deleted file mode 100644 index 84666bb..0000000 --- a/Classes/02_optical/dp_fiber_lib/canUseGPU.m +++ /dev/null @@ -1,20 +0,0 @@ -function canUse = canUseGPU() -%CANUSEGPU Check if a compatible GPU is available for parallel computing -% Returns true if MATLAB Parallel Computing Toolbox is available and -% a CUDA-capable GPU is detected. - -canUse = false; - -% Check if Parallel Computing Toolbox is installed -if ~license('test', 'Distrib_Computing_Toolbox') - return; -end - -% Check for GPU device -try - gpu = gpuDevice(); - canUse = gpu.DeviceSupported; -catch - canUse = false; -end -end diff --git a/Classes/02_optical/dp_fiber_lib/getNLstepsize.m b/Classes/02_optical/dp_fiber_lib/getNLstepsize.m index fd42eba..a70edb9 100644 --- a/Classes/02_optical/dp_fiber_lib/getNLstepsize.m +++ b/Classes/02_optical/dp_fiber_lib/getNLstepsize.m @@ -1,25 +1,25 @@ function [rDZ] = getNLstepsize(ux,uy,gamma,dzmin,dzmax,dphimax,alpha_lin) -% gather() handles gpuArray inputs - ensures scalar is on CPU for comparisons -maxPow = gather(max(gamma.*max(real(ux).^2+imag(ux).^2+real(uy).^2+imag(uy).^2))); - -Leff = dphimax/maxPow; -alpha_lin = max([alpha_lin.X alpha_lin.Y]); -nl_att_len_ratio = alpha_lin*Leff; - -if nl_att_len_ratio >= 1 - rDZ = dzmax; -else - if alpha_lin == 0 - step = Leff; + + maxPow = max(gamma.*max(real(ux).^2+imag(ux).^2+real(uy).^2+imag(uy).^2)); + + Leff = dphimax/maxPow; + alpha_lin = max([alpha_lin.X alpha_lin.Y]); + nl_att_len_ratio = alpha_lin*Leff; + + if nl_att_len_ratio >= 1 + rDZ = dzmax; else - %effective length? - step = -1/alpha_lin*log(1-nl_att_len_ratio); - end - - rDZ = min([step dzmax]); - rDZ = max([rDZ dzmin]); -end - + if alpha_lin == 0 + step = Leff; + else + %effective length? + step = -1/alpha_lin*log(1-nl_att_len_ratio); + end + + rDZ = min([step dzmax]); + rDZ = max([rDZ dzmin]); + end + end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/lin_step.m b/Classes/02_optical/dp_fiber_lib/lin_step.m index cdb3562..eea8a1e 100644 --- a/Classes/02_optical/dp_fiber_lib/lin_step.m +++ b/Classes/02_optical/dp_fiber_lib/lin_step.m @@ -22,37 +22,37 @@ n_plates_left = n_plates - n_plates_done; % compute last plate size ( if it fits, it should be 0) if missing_dz > aStepSize - + last_plate = aStepSize; missing_dz = missing_dz-aStepSize; plate_sizes = last_plate; plate_numbers = n_plates; - + else - + last_plate = aStepSize - missing_dz - (n_plates_left-1)*corr_length; - + if missing_dz == 0 missing_dz = []; end - + %build vector of plate lengths with missing plate part from prev. %iterartion , then some normal plates and finally a fraction of a plate %to fit into the step length plate_sizes = [missing_dz corr_length*ones(1,n_plates_left-1) last_plate]; - + if n_plates_done == 0 plate_numbers =[(n_plates_done+1):(n_plates-1) n_plates]; else plate_numbers = [n_plates_done (n_plates_done+1):(n_plates-1) n_plates]; % not wrking yet end - - %remember for next step + + %remember for next step missing_dz = corr_length - last_plate; - + end -plate_steps = repmat(n_step,1,length(plate_sizes)); %#ok +plate_steps = repmat(n_step,1,length(plate_sizes)); % %figure;stem(plate_sizes); @@ -66,39 +66,50 @@ test_plate_numbers = [test_plate_numbers, plate_numbers]; opt_x=fft(opt_x); opt_y=fft(opt_y); -% Note: db1, db0, common_beta_x, common_beta_y are already on GPU when -% useGPU=true (transferred in CNLSE_plain) -db1 = brf.db1; -db0 = brf.db0; +% db1 = gpuArray(brf.db1); +% db0 = gpuArray(brf.db0); +% common_beta_x = gpuArray(common_beta_x); +% common_beta_y = gpuArray(common_beta_y); + +db1 = (brf.db1); +db0 = (brf.db0); +common_beta_x = (common_beta_x); +common_beta_y = (common_beta_y); % process every waveplate with given sizes in plate_sizes for n=1:length(plate_sizes) dz = plate_sizes(n); - + + % figure(87);subplot(2,1,1);plot(real(x(900:1150)));subplot(2,1,2);plot(real(y(900:1150))); + % MOV1=[MOV1 getframe(87)]; + % extract rotation matrix from pre calculated matrices matR = brf.matR{plate_numbers(n)}; - + % transform to eigenvalue of of fiber segment tOpt.X = conj(matR(1,1))*opt_x + conj(matR(2,1))*opt_y; tOpt.Y = conj(matR(1,2))*opt_x + conj(matR(2,2))*opt_y; - + % calculate statistical delta beta for pmd delta_beta = 0.5*(db1+db0(n))/corr_length; - + % build transfer function with delta beta + %common.beta = beta1+beta2*omega^2 + %accumulate delta beta for log... brf.simdgd = brf.simdgd + (db1(length(db1)/2+1)+db0(n))/corr_length; - + h.X = exp(-1j*(common_beta_x-delta_beta)*dz); h.Y = exp(-1j*(common_beta_y+delta_beta)*dz); - + % delta_beta has to be added to the transfer function + % process with transfer function tOpt.X = h.X.*tOpt.X ; tOpt.Y = h.Y.*tOpt.Y ; - + % rotate back opt_x = matR(1,1)*tOpt.X + matR(1,2)*tOpt.Y; opt_y = matR(2,1)*tOpt.X + matR(2,2)*tOpt.Y; - + end lin_z_test = lin_z_test + sum(plate_sizes,2); @@ -106,8 +117,9 @@ lin_z_test = lin_z_test + sum(plate_sizes,2); %update the number of processed plates so far n_plates_done = n_plates_done + n_plates_left; -% attenuate the signal each linear step with alpha -rOpt_x=ifft(exp(-alpha_lin_x*aStepSize/2).*opt_x); -rOpt_y=ifft(exp(-alpha_lin_y*aStepSize/2).*opt_y); +% attanuate the signal each linear state with alpha +% ( 0.2dB = 4.6052e-05 ) +rOpt_x=ifft(exp(-alpha_lin_x*aStepSize/2).*opt_x); % /2 not sure why (have to find it in formulas) +rOpt_y=ifft(exp(-alpha_lin_y*aStepSize/2).*opt_y); % but not relevant for now -end +end \ No newline at end of file diff --git a/Classes/04_DSP/Equalizer/EQ.m b/Classes/04_DSP/Equalizer/EQ.m index a03d66b..f9d58d2 100644 --- a/Classes/04_DSP/Equalizer/EQ.m +++ b/Classes/04_DSP/Equalizer/EQ.m @@ -10,6 +10,10 @@ classdef EQ < handle training_length %Number of training symbols training_loops %Number of loops through sequence for training mode ideal_dfe %Error free DFE decisions + weighted_DFE %Weighted DFE on/off + weighted_DFE_mode %Weighted DFE mode + weighted_DFE_d_min %d_min threshold parameter for weighted DFE mode 1 + weighted_DFE_I_mode %[a_s, b_s, I_max]-parameters for the weighted DFE DB_aim %Aim at duobinary output sequence @@ -68,6 +72,10 @@ classdef EQ < handle options.training_length = 1024 %Number of training symbols options.training_loops = 1 %Number of loops through sequence for training mode options.ideal_dfe = 0 %Error free DFE decisions + options.weighted_DFE = 0; + options.weighted_DFE_mode = 'R1'; + options.weighted_DFE_d_min = 0.5; + options.weighted_DFE_I_mode = [5,0.5,0.6]; options.DB_aim %Aim at duobinary output sequence @@ -438,6 +446,43 @@ classdef EQ < handle dd_out(k) = constellation_in_(dd_idx); end + % Implementation of a weighted DFE in + % order to prevent error propagation. + % For further details, study [1], chapter 3.2.2 - + % Modifications of DFE + + if obj.weighted_DFE + % define new constellations + const = unique(ref_in); + + % determine reliability factor gamma_k + if output_vec(m) > min(const) && output_vec(m) < max(const) + gamma_k = 1 - abs(output_vec(m) - dd_out(k)); + else + gamma_k = 1; + end + + % select mode + if strcmp(obj.weighted_DFE_mode,'R1') + if gamma_k >= obj.weighted_DFE_d_min + f_gamma_k = 1; + else + f_gamma_k = 0; + end + elseif strcmp(obj.weighted_DFE_mode,'R2') + f_gamma_k = gamma_k; + elseif strcmp(obj.weighted_DFE_mode,'I1') + nom = 1-exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1)); + denom = 1+exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1)); + f_gamma_k = (1/2)*((nom/denom) - 1); + elseif strcmp(obj.weighted_DFE_mode,'I2') + nom = 1-exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1)); + denom = 1+exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1)); + f_gamma_k = (obj.weighted_DFE_I_mode(3)/2)*((nom/denom) - 1); + end + % calculate weighted output + output_vec(m) = f_gamma_k.*dd_out(k)+(1-f_gamma_k).*output_vec(m); + end if obj.Nb(1) > 0 dd_DFE(2:end) = dd_DFE(1:end-1); @@ -739,3 +784,6 @@ classdef EQ < handle end end +% References +% [1] T. J. Wettlin, “Experimental Evaluation of Advanced Digital Signal Processing for Intra-Datacenter Systems using Direct-Detection,” 2023. [Online]. Available: https://nbn-resolving.org/urn:nbn:de:gbv:8:3-2023-00703-8 + diff --git a/Classes/04_DSP/Equalizer/ML_MLSE.m b/Classes/04_DSP/Equalizer/ML_MLSE.m index 168b837..9965525 100644 --- a/Classes/04_DSP/Equalizer/ML_MLSE.m +++ b/Classes/04_DSP/Equalizer/ML_MLSE.m @@ -597,7 +597,7 @@ classdef ML_MLSE < handle obj.S = obj.nSym; obj.Nf = obj.order * obj.sps; obj.nStates = obj.S^obj.L; - obj.nFeasible = obj.nStates * obj.S; %feasible state transitions + obj.nFeasible = obj.nStates * obj.S; % --- Trellis mapping obj.trellis_states = reshape(obj.constellation,1,[]); diff --git a/Classes/04_DSP/Timing_Recovery.m b/Classes/04_DSP/Timing_Recovery.m deleted file mode 100644 index 23cc331..0000000 --- a/Classes/04_DSP/Timing_Recovery.m +++ /dev/null @@ -1,47 +0,0 @@ -classdef Timing_Recovery < handle - - properties(Access=public) - timing_error_detector - sps - damping_factor - normalized_loop_bandwidth - detector_gain - end - - methods(Access=public) - function obj = Timing_Recovery(options) - arguments(Input) - - options.timing_error_detector = 'Gardner'; - options.sps = 2; - options.damping_factor = 1.0; - options.normalized_loop_bandwidth = 0.005; - options.detector_gain = 1; - - end - - fn = fieldnames(options); - for n = 1:numel(fn) - obj.(fn{n}) = options.(fn{n}); - end - - - - end - - function [data_out,timing_error] = process(obj, data_in) - - timing_synchronization = comm.SymbolSynchronizer( ... - "TimingErrorDetector", obj.timing_error_detector, ... - "SamplesPerSymbol", obj.sps, ... - "DampingFactor", obj.damping_factor, ... - "NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ... - "DetectorGain", obj.detector_gain); - - data_out = data_in; - [data_out.signal,timing_error] = timing_synchronization(data_in.signal); - - end - end -end - diff --git a/Classes/DataBaseHandler/DBHandler.m b/Classes/DataBaseHandler/DBHandler.m index b7ba026..5569746 100644 --- a/Classes/DataBaseHandler/DBHandler.m +++ b/Classes/DataBaseHandler/DBHandler.m @@ -23,12 +23,6 @@ classdef DBHandler < handle arguments options.dataBase = ""; % Default value for pathToDB if not provided options.type = "mysql"; - options.server = ""; - options.port = 3306; - options.user = ""; - options.password = ""; - - end % Assign values to class properties based on input arguments @@ -52,13 +46,12 @@ classdef DBHandler < handle obj.conn = database( ... string(obj.dataBase), ... % Database name - options.user, ... % Username - options.password, ... % Password (or getSecret) + "silas", ... % Username + "silas", ... % Password (or getSecret) "Vendor", "MySQL", ... - "Server", options.server, ... + "Server", "134.245.243.254", ... "PortNumber", 3306, ... "JDBCDriverLocation", "C:\Users\Silas\Documents\mysql-connector-j-9.3.0\mysql-connector-j-9.3.0.jar"); - end catch e diff --git a/Functions/EQ_structures/ml_mlse.m b/Functions/EQ_structures/ml_mlse.m index 04fe317..84ff5cb 100644 --- a/Functions/EQ_structures/ml_mlse.m +++ b/Functions/EQ_structures/ml_mlse.m @@ -54,6 +54,12 @@ for k = 1:numel(fn) end json_str = jsonencode(eq_small); +fn = fieldnames(eq_); +for k = 1:numel(fn) + if issparse(eq_.(fn{k})) + eq_.(fn{k}) = full(eq_.(fn{k})); + end +end ml_mlse_results.config.eq = jsonencode(eq_); ml_mlse_results.config.equalizer_structure = int32(equalizer_structure.ml_mlse); ml_mlse_results.config.comment = 'function: ML-based MLSE'; @@ -87,6 +93,7 @@ end function [bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, precode_mode, M, eth_style) % Calculate BER based on precoding mode mapper = PAMmapper(M, 0, "eth_style", eth_style); +skip_front = 150; switch precode_mode case db_mode.no_db @@ -101,11 +108,11 @@ switch precode_mode tx_bits_precoded = mapper.demap(tx_symbols_precoded); rx_bits = mapper.demap(eq_signal_hd_precoded); - [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1); % B) Just determine BER rx_bits = mapper.demap(eq_signal_hd); - [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1); case db_mode.db_precoded % Data is precoded on TX side @@ -113,12 +120,12 @@ switch precode_mode eq_signal_hd_decoded = Duobinary().encode(eq_signal_hd, "M", M); eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M); rx_bits_decoded = mapper.demap(eq_signal_hd_decoded); - [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1); % B) Omit the Coding by comparing with demapped TX symbol sequence tx_bits_demapped = mapper.demap(tx_symbols); rx_bits = mapper.demap(eq_signal_hd); - [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1); + [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1); end end diff --git a/Functions/Theory/Dissertation/mach_zehnder_modulator.m b/Functions/Theory/Dissertation/mach_zehnder_modulator.m deleted file mode 100644 index f9eca58..0000000 --- a/Functions/Theory/Dissertation/mach_zehnder_modulator.m +++ /dev/null @@ -1,219 +0,0 @@ - - -% Parameters -c0 = physconst('lightspeed'); % [m/s] -lambda0 = 1310e-9; % [m] -omega0 = 2*pi*c0/lambda0; - -L = 5e-3; % [m] effective phase section length (set as needed) -n_eff = 2.2; % [-] effective index (set as needed) - -E0 = 1; % field amplitude (arbitrary) -Vpi = 3.2; % [V] half-wave voltage (your V_pi) - -% Drive -f0 = 1e9; % [Hz] -fs = 200e9; % [Hz] -Nper = 2; % number of periods -Vpp = 0.6*Vpi; % [V] peak-to-peak of v_drive(t) - -biasV = 1.1; % [V] differential bias added to v_drive - -% Time axis + differential drive voltage v_drive(t) -T = Nper/f0; -t = (0:1/fs:T-1/fs).'; - - -if 1 - % SINE - v_drive = biasV + (Vpp/2)*sin(2*pi*f0*t); % v_drive(t) (peak = Vpp/2) - -else - - % --- Generate PAM-4 Sequence --- - symbols = linspace(-0.5, 0.5, 4); - num_symbols = 12; % Increased slightly for better visual - rng(44); - random_data = symbols(randi(4, 1, num_symbols)); - - % Create time axis (Note: T is your period from the sine code) - sps = round(T * fs); - t = (0:1/fs:(num_symbols*T)-1/fs).'; - - % Upsample to rectangular waveform - v_pam = repelem(random_data, sps).'; - - % Apply swing and bias: Resulting range is [biasV-Vpp/2, biasV+Vpp/2] - v_drive_rect = biasV + (v_pam * Vpp); - - % --- Round the edges --- - filter_span = round(sps/1.5); % Increased span for smoother "rounding" - window = gausswin(filter_span); - window = window / sum(window); - - % Apply filter (using 'same' to keep vector length, but be aware of edge transients) - v_drive = conv(v_drive_rect, window, 'same'); - -end - - -% Analytic -v_ = linspace(-1,2, 2001); -% Field transfer function (amplitude) -Field_mzm_analytic = cos((pi/2)*v_); - -% Power transfer function (intensity) -P_mzm_analytic = Field_mzm_analytic.^2; - -% Imbalance factor in YOUR notation: -rho = 1; - -% Push-pull branch voltages (consistent with v_drive = v1 - v2) -v1 = +0.5*v_drive; % arm 1 -v2 = -0.5*v_drive; % arm 2 - -% Phases phi1, phi2 -phi1 = pi * v1 / Vpi; -phi2 = pi * v2 / Vpi; - -% Fields: E_in and E_out (exactly your Eq. (mzm_e_field)) -E_in = E0 .* exp(1i*omega0*t); - -common_phase = exp(-1i * (omega0*L*n_eff/c0)); % exp(-j*omega0*L*n_eff/c0) - -E_out = E0 .* exp(1i*omega0*t) .* common_phase .* 0.5 .* ... - ( exp(-1i*phi1) + rho .* exp(-1i*phi2) ); - -% Transfer function (numerical): E_out/E_in -H_num = E_out ./ E_in; - -% Power (normalized) -Pnorm_num = abs(H_num).^2; % since |E_out/E_in|^2 - -% Ideal TF (analytic) for comparison (rho=1, push-pull) -H_ideal = common_phase .* cos( (pi/2) * (v_drive./Vpi) ); - -Pnorm_ideal = abs(H_ideal).^2; -Pnorm_math = cos( (pi/2) * (v_drive./Vpi) ).^2; - - -set(groot, 'defaultLegendInterpreter', 'tex'); -set(groot, 'defaultAxesTickLabelInterpreter', 'tex'); -set(groot, 'defaultTextInterpreter', 'tex'); - -% Normalized voltage axis (multiples of Vpi) -v_norm = v_drive./Vpi; - -colfield = [0,0,0]; %is black -colpow = linspecer(2); -colpow = colpow(1,:); -colvdrive = linspecer(2); -colvdrive = colvdrive(2,:); - -%% SIGNAL IN -figure(1); clf -plot(v_norm,t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on; -ylabel('t [ns]'); xlabel('v_{drive}(t)/V_\pi'); -title('Drive voltage (normalized)'); -xlim([min(v_) max(v_)]); -% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\linear_casee\mzm_input_signal.tex'); - - -%% IN/OUT (static transfer) — normalized x-axis + analytic curve -if 0 -figure(2); clf -plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF -plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color',colpow); hold on;% analytic power TF -% show input time signal -plot(v_norm,-1+t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on; -% show output time signal -plot(2+t*1e9, Pnorm_num, 'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colvdrive); hold on; -plot(2+t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on; -scatter(v_norm, Pnorm_num, 12, '.', 'LineWidth', 1,'MarkerEdgeColor',colvdrive); -scatter(biasV./Vpi,(cos((pi/2)*biasV./Vpi)^2),10,'Marker','o'); -line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--'); -line([max(v_drive) max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--'); -xline([min(v_norm) max(v_norm)]) - -grid on; -xlabel('v_{drive}(t)/V_\pi'); ylabel('|E_{out}/E_{in}|^2'); -% legend -xlim([min(v_) max(v_)+1]); -ylim([-1 1]); - -% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm.tex'); -end -%% - -figure(3); clf -plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF -plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color',colpow); hold on;% analytic power TF - -scatter(v_norm, Pnorm_num, 12, '.', 'LineWidth', 1,'MarkerEdgeColor',colvdrive); -scatter(biasV./Vpi,(cos((pi/2)*biasV./Vpi)^2),10,'Marker','o'); -line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--'); -line([max(v_drive) max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--'); -xline([min(v_norm) max(v_norm)]) - -grid on; -xlabel('v_{drive}(t)/V_\pi'); ylabel('|E_{out}/E_{in}|^2'); -% legend -xlim([min(v_) max(v_)]); -ylim([-1 1]); - -% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm_tramsfer_function_matlab.tex'); - - -%% - -figure(4); clf -% plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF -plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color','black'); hold on;% analytic power TF -input_dots = linspace(min(v_drive),max(v_drive),4)./Vpi; -% input_dots = unique(v_drive_rect)./Vpi; -output_dots = (cos((pi/2)*input_dots).^2); -scatter(input_dots,output_dots,'Marker','x','LineWidth',1,'MarkerEdgeColor','black'); -scatter(input_dots,zeros(size(input_dots)),'Marker','^','LineWidth',2,'MarkerEdgeColor','black'); -% scatter(ones(size(input_dots)),output_dots,'Marker','<','LineWidth',2,'MarkerEdgeColor','black'); - -for i = 1:numel(input_dots) - % Draw the dashed projection lines - line([input_dots(i), input_dots(i)], [output_dots(i), 0], 'linewidth', 0.5, 'color', 'black', 'linestyle', '--', 'handlevisibility', 'off'); - line([input_dots(i), 1], [output_dots(i), output_dots(i)], 'linewidth', 0.5, 'color', 'black', 'linestyle', '--', 'handlevisibility', 'off'); - - % Add the level annotation boxes near the output (y-axis) - % Adjust the '1.05' to move the box further right or 'output_dots(i)' for height - j = 3-(i-1)*2; - text(1, output_dots(i), sprintf('Level %d', j), ... - 'FontSize', 8, ... - 'EdgeColor', 'black', ... - 'BackgroundColor', 'white', ... - 'Margin', 2); -end - -xlim([0,1.5]); -ylim([0,1]) - -% line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), 0],'linewidth',0.5,'color','black','linestyle','--'); -% line([max(v_drive), max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), 0],'linewidth',0.5,'color','black','linestyle','--'); - -% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\linear_casee\mzm_tf.tex'); -% xticks(sort(input_dots)); -% yticks(sort(output_dots)); -grid off - - -%% -% % FIELD TF (only field here; do not mix power into this figure) -figure(5); clf -% plot(t*1e9, real(H_num), 'LineWidth', 1.0); hold on; -% plot(t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on; -plot(t*1e9, Pnorm_num, 'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colpow); hold on; -grid on; -xlabel('t [ns]'); ylabel('Re\{E_{out}/E_{in}\}'); -legend -yticks(sort(output_dots)); -% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\linear_casee\mzm_output_signal.tex'); - - - diff --git a/Functions/mat2tikz_improved.m b/Functions/mat2tikz_improved.m deleted file mode 100644 index a36b1d1..0000000 --- a/Functions/mat2tikz_improved.m +++ /dev/null @@ -1,23 +0,0 @@ -function mat2tikz_improved(filename) -arguments - % Default to the path in your example if no argument is provided - filename (1,1) string = 'C:\Users\Silas\Documents\Dissertation\00_Examples\tikz\textfig.tikz'; -end -cleanfigure; -matlab2tikz(char(filename), ... - 'width','\fwidth', ... - 'height','\fheight', ... - 'showInfo',false, ... - 'extraAxisOptions',{ ... - 'legend style={font=\footnotesize}', ... - 'xlabel style={font=\color{white!15!black},font=\small},',... - 'ylabel style={font=\color{white!15!black},font=\small},',... - 'legend columns=1', ... - 'every axis/.append style={font=\scriptsize}',... - 'legend columns=1',... - 'legend style={at={(0.02,0.98)},font=\footnotesize,draw=black!60,rounded corners=2pt,inner sep=1pt,fill=white,column sep=6pt,anchor= north west}',... - 'legend style={at={(0.02,0.98)},draw=white!0!white,font=\scriptsize,inner sep=0.1pt,fill=white,column sep=1pt,anchor= north west}',... - 'every axis/.append style={font=\scriptsize}',... - }); - -end \ No newline at end of file diff --git a/Libs/wesanderson_colors/WesPalette.m b/Libs/wesanderson_colors/WesPalette.m index 13e976f..7cd1dcb 100644 --- a/Libs/wesanderson_colors/WesPalette.m +++ b/Libs/wesanderson_colors/WesPalette.m @@ -1,13 +1,10 @@ classdef WesPalette % WESPALETTE Wes Anderson color palettes with auto-completion % Usage: - % cmap = WesPalette.Zissou1.rgb() % full palette - % cmap = WesPalette.Zissou1.rgb(3) % 3 colors (discrete default) - % cmap = WesPalette.Zissou1.rgb(12,"discrete") % any n, no interpolation - % cmap = WesPalette.Zissou1.rgb(256,"continuous") % smooth colormap (Lab interpolation) - % - % Requires: - % - colorspace.m (Pascal Getreuer) on MATLAB path for "continuous" mode + % cmap = WesPalette.Zissou1.rgb() + % cmap = WesPalette.Zissou1.rgb(3) + + % https://github.com/karthik/wesanderson?tab=readme-ov-file enumeration BottleRocket1 @@ -37,33 +34,21 @@ classdef WesPalette end methods - function cmap = rgb(obj, n, mode) + function cmap = rgb(obj, n) % Return palette as Nx3 RGB colormap [0–1] - % - % n : number of requested colors (optional) - % mode : "discrete" (default) or "continuous" - base_hex = obj.hex(); - base_rgb = WesPalette.hex2rgb(base_hex); + hex = obj.hex(); - if nargin < 2 || isempty(n) - cmap = base_rgb; - return; - end - if nargin < 3 || isempty(mode) - mode = "discrete"; - end - mode = lower(string(mode)); + rgb = hex2rgb(hex); - validateattributes(n, {'numeric'}, {'scalar','integer','positive'}, mfilename, 'n'); - if mode ~= "discrete" && mode ~= "continuous" - error('mode must be "discrete" or "continuous".'); - end - - if mode == "discrete" - cmap = WesPalette.sample_discrete(base_rgb, n); + if nargin == 2 + if n > size(rgb,1) + error('Requested %d colors, but only %d available.', ... + n, size(rgb,1)) + end + cmap = rgb(1:n,:); else - cmap = WesPalette.interpolate_continuous_lab(base_rgb, n); + cmap = rgb; end end end @@ -71,116 +56,78 @@ classdef WesPalette methods (Access = private) function hex = hex(obj) % Internal HEX storage + switch obj case WesPalette.BottleRocket1 hex = {'#A42820','#5F5647','#9B110E','#3F5151','#4E2A1E','#550307','#0C1707'}; + case WesPalette.BottleRocket2 hex = {'#FAD510','#CB2314','#273046','#354823','#1E1E1E'}; + case {WesPalette.Rushmore1, WesPalette.Rushmore} hex = {'#E1BD6D','#EABE94','#0B775E','#35274A','#F2300F'}; + case WesPalette.Royal1 hex = {'#899DA4','#C93312','#FAEFD1','#DC863B'}; + case WesPalette.Royal2 hex = {'#9A8822','#F5CDB4','#F8AFA8','#FDDDA0','#74A089'}; + case WesPalette.Zissou1 hex = {'#3B9AB2','#78B7C5','#EBCC2A','#E1AF00','#F21A00'}; + case WesPalette.Zissou1Continuous hex = {'#3A9AB2','#6FB2C1','#91BAB6','#A5C2A3','#BDC881', ... '#DCCB4E','#E3B710','#E79805','#EC7A05','#EF5703','#F11B00'}; + case WesPalette.Darjeeling1 hex = {'#FF0000','#00A08A','#F2AD00','#F98400','#5BBCD6'}; + case WesPalette.Darjeeling2 hex = {'#ECCBAE','#046C9A','#D69C4E','#ABDDDE','#000000'}; + case WesPalette.Chevalier1 hex = {'#446455','#FDD262','#D3DDDC','#C7B19C'}; + case WesPalette.FantasticFox1 hex = {'#DD8D29','#E2D200','#46ACC8','#E58601','#B40F20'}; + case WesPalette.Moonrise1 hex = {'#F3DF6C','#CEAB07','#D5D5D3','#24281A'}; + case WesPalette.Moonrise2 hex = {'#798E87','#C27D38','#CCC591','#29211F'}; + case WesPalette.Moonrise3 hex = {'#85D4E3','#F4B5BD','#9C964A','#CDC08C','#FAD77B'}; + case WesPalette.Cavalcanti1 hex = {'#D8B70A','#02401B','#A2A475','#81A88D','#972D15'}; + case WesPalette.GrandBudapest1 hex = {'#F1BB7B','#FD6467','#5B1A18','#D67236'}; + case WesPalette.GrandBudapest2 hex = {'#E6A0C4','#C6CDF7','#D8A499','#7294D4'}; + case WesPalette.IsleofDogs1 hex = {'#9986A5','#79402E','#CCBA72','#0F0D0E','#D9D0D3','#8D8680'}; + case WesPalette.IsleofDogs2 hex = {'#EAD3BF','#AA9486','#B6854D','#39312F','#1C1718'}; + case WesPalette.FrenchDispatch hex = {'#90D4CC','#BD3027','#B0AFA2','#7FC0C6','#9D9C85'}; + case WesPalette.AsteroidCity1 hex = {'#0A9F9D','#CEB175','#E54E21','#6C8645','#C18748'}; + case WesPalette.AsteroidCity2 hex = {'#C52E19','#AC9765','#54D8B1','#B67C3B','#175149','#AF4E24'}; + case WesPalette.AsteroidCity3 hex = {'#FBA72A','#D3D4D8','#CB7A5C','#5785C1'}; end end end - - methods (Static, Access = private) - function rgb = hex2rgb(hex) - % hex: cellstr like {'#RRGGBB', ...} - if isstring(hex), hex = cellstr(hex); end - n = numel(hex); - rgb = zeros(n,3); - for i = 1:n - h = char(hex{i}); - if startsWith(h,'#'), h = h(2:end); end - if numel(h) ~= 6 - error('Invalid HEX color: %s', hex{i}); - end - rgb(i,1) = hex2dec(h(1:2))/255; - rgb(i,2) = hex2dec(h(3:4))/255; - rgb(i,3) = hex2dec(h(5:6))/255; - end - end - - function cmap = sample_discrete(base_rgb, n) - % No interpolation; allow any n by sampling/repeating. - k = size(base_rgb,1); - - if n <= k - idx = round(linspace(1, k, n)); % spread across palette - idx = max(1, min(k, idx)); - cmap = base_rgb(idx,:); - else - reps = floor(n / k); - rmd = mod(n, k); - cmap = [repmat(base_rgb, reps, 1); base_rgb(1:rmd,:)]; - end - end - - function cmap = interpolate_continuous_lab(base_rgb, n) - % Smooth interpolation in Lab using colorspace(). - % Requires colorspace.m by Pascal Getreuer on MATLAB path. - - k = size(base_rgb,1); - if k == 1 - cmap = repmat(base_rgb, n, 1); - return; - end - - % Convert to Lab, interpolate each channel, convert back - lab = colorspace('Lab<-RGB', base_rgb); - - t_base = linspace(0, 1, k); - t_new = linspace(0, 1, n); - - lab_new = zeros(n,3); - for c = 1:3 - lab_new(:,c) = interp1(t_base, lab(:,c), t_new, 'linear'); - end - - rgb_new = colorspace('RGB<-Lab', lab_new); - - % Clamp to displayable gamut - cmap = min(max(rgb_new, 0), 1); - end - end end diff --git a/Libs/wesanderson_colors/minimal_example_wespalette.m b/Libs/wesanderson_colors/minimal_example_wespalette.m index f308f27..9fcd99b 100644 --- a/Libs/wesanderson_colors/minimal_example_wespalette.m +++ b/Libs/wesanderson_colors/minimal_example_wespalette.m @@ -5,8 +5,8 @@ y2 = 1e0 ./ (1 + exp(-0.4*(x-12))); % NLPN y3 = 1e-6 * 10.^(0.45*x); % RP on gamma y4 = 1e-2 * 10.^(0.18*(x-8)); % RP on beta2 -cmap = WesPalette.AsteroidCity1; -% cmap = linspecer(4); +cmap = WesPalette.AsteroidCity1.rgb(4); +cmap = linspecer(4); figure1=figure(202998);clf;hold on lw = 0.8; ms = 4; plot(x,y1,'LineWidth',lw,'Color',cmap(1,:),'Marker','o','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',[1,1,1],'MarkerSize',ms); diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m deleted file mode 100644 index 9a327ab..0000000 --- a/projects/FSO_transmission/first_analysis.m +++ /dev/null @@ -1,252 +0,0 @@ - -base = "C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC"; -mode = 0; %0 oder 1 -M = 2; - -all_files = dir(fullfile(base, "**/*.mat")); - -if M == 2 - tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat"); - filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat"); -elseif M == 4 - tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat"); - filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat"); -end - -if mode == 1 - [f, p] = uigetfile(fullfile(base, "**/*.mat")); - if f~=0 - filename = fullfile(p,f); - end -end - - -datas = load(filename); - -%% -str = filename; -M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once')); -assert(M==M_); -fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once')); -fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once')); -I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f'); -rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f'); -L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f'); -pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once')); -rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once')); -mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once')); - -%% -% Tx data - -Bits = Informationsignal(tx_data.tx_data,"fs",fsym); -Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym); - -mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there -PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme - -Symbols_ = PM.map(Bits) .* PM.scaling; -assert(isequal(Symbols.signal,Symbols_.signal)); - -Bits_ = PM.demap(Symbols); -[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal); -assert(ber == 0); - -%% For comparison, apply pulsef on Tx Symbols -Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff); -Digi_sig_tx_compare = Pform.process(Symbols); - -Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); -Digi_sync = Pform.process(Symbols); - -MF = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff); -Rx_sig_compare = MF.process(Digi_sig_tx_compare); - -%% - -% Rx Data -traceData = datas.tr.lastData(2).trace.ch3; - -%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin -scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin; -demystified = isequal(traceData.YData,scoperead_volts); -assert(demystified); - -timesig_compare = [0:1:datas.tr.lastData(1).trace.ch3.Points-1] ./ fs; -timesig = datas.tr.lastData(1).trace.ch3.XData; -assert(isequal(traceData.YData,scoperead_volts)); - -Scope_sig = Electricalsignal(traceData.YData,"fs",fs); - -%% -% 1) matched filter -% pulse is symmetric, hence we can use pulsef directly as matched filter. -% It feels off (bit I think correct) that the fsym is now the output freq.!! -% -> output 2 sps to omit timing recovery!? -apply_matched_filter = 1; -k = 1; -if apply_matched_filter - - Pform = Pulseformer("fsym",fsym,"fdac",k*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); - Rx_matched = Pform.process(Scope_sig); -else - - Rx_matched = Filter('filtdegree',4,"f_cutoff",fsym*0.5,"fs",Scope_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scope_sig); - Rx_matched = Rx_matched.resample("fs_out",k*fsym); -end -Rx_matched.spectrum(); - - -%% - -coefficients = arburg(Rx_matched.signal,25); - -figure() -[h,w] = freqz(1,coefficients,Rx_matched.length,"whole",Rx_matched.fs); -h = h/max(abs(h)); -hold on -w_ = (w - Rx_matched.fs/2); -plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(coefficients,2)), ' '],'LineWidth',2); - -%% Timing Rec -apply_timing_rec = 1; -if apply_timing_rec - [Rx_symbolsync, timing_error] = Timing_Recovery("timing_error_detector",'Gardner (non-data-aided)','sps',k,'damping_factor',0.1,'normalized_loop_bandwidth',0.1,'detector_gain',2.7).process(Rx_matched); - figure();plot(timing_error); - Rx_symbolsync.fs = fsym; - % Tsynch - [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0); - Rx_synced = Rx_synced_cell{1}; - sps = 1; -else - % Tsynch - [~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0); - Rx_synced = Rx_synced_cell{1}; - Rx_synced = Rx_synced.resample("fs_out",2*fsym); - sps = 2; -end - -if M == 2 - ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature -elseif M == 4 - ber_in_paper = 10^(-2.5); -end - - -%% - -Rx_synced = Rx_synced_cell{1}; - - -% -------------------- FFE -------------------- -% requires some more digging what is going on :-) -eq_ffe = EQ("Ne",[50, 1, 1],"Nb",[2,0,0], ... - "training_length",512,"training_loops",5,"dd_loops",5, ... - "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... - "FFEmu",0,"plotfinal",0,"ideal_dfe",1); - -vars = logspace(-4,-3,36); - -parfor i = 1:numel(vars) - - len_tr = 4096; - mu_ffe1 = 0.01;% mus(i);%0.0001; - mu_ffe2 = 0.0008; - mu_ffe3 = 0.001; - mu_dc = 0.005; - mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; - mu_dfe = vars(i); - duob_mode = db_mode.no_db; - - % requires some more digging what is going on :-) - eq_ffe_1 = EQ("Ne",[150, 1, 0],"Nb",[50,0,0], ... - "training_length",len_tr,"training_loops",5,"dd_loops",5, ... - "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... - "FFEmu",0,"plotfinal",0,"ideal_dfe",0); - - eq_ffe_2 = FFE("epochs_tr",1,"epochs_dd",vars(i),"len_tr",4096,"mu_dd",vars(i),"mu_tr",vars(i),"order",999,"sps",1,"decide",0, "adaption",adaption_method.nlms,"dd_mode",0); - % eq_ffe_2 = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-5,"dfe_mu_dd",mus(i),"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",50,"dfe_order",10,"sps",1,"decide",1); - - - ffe_results = ffe(eq_ffe_1,M,Rx_synced,Symbols,Bits, ... - "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ... - "eth_style_symbol_mapping",mapping_style); - - ffe_results.metrics.BER - bers(i) = ffe_results.metrics.BER; -end - -figure(); -plot(vars,bers); -yline(ber_in_paper) -beautifyBERplot(); -% -fprintf('Paper: %.1e \n \n',ber_in_paper); -ffe_results.metrics.print("description",'FFE'); -fprintf('FFE: %.1e \n',ffe_results.metrics.BER); - - -%% -------------------- VNLE + MLSE -------------------- -len_tr = 4096; -mu_ffe1 = 0.0001;% mus(i);%0.0001; -mu_ffe2 = 0.0008; -mu_ffe3 = 0.001; -mu_dc = 0.005; -mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; -mu_dfe = 0.0004; -duob_mode = db_mode.no_db; - -pf_ncoeffs = 4; - -vars = 1:7; -bers = zeros(size(vars)); -parfor i = 1:numel(vars) - - eqv = EQ("Ne",[200, 1, 0],"Nb",[2, 0, 0], ... - "training_length",len_tr,"training_loops",5,"dd_loops",5, ... - "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ... - "FFEmu",0,"plotfinal",0,"ideal_dfe",1); - pf_ncoeffs = vars(i); - pf = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels); - [vnle_results, mlse_results] = vnle_postfilter_mlse(eqv, pf, mlse_, M, Rx_synced, Symbols, Bits, ... - "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style); - fprintf('Paper: %.1e \n \n',ber_in_paper); - vnle_results.metrics.print("description",'VNLE'); - mlse_results.metrics.print("description",'MLSE'); - bers(i) = mlse_results.metrics.BER; -end -%% -figure();hold on -plot(vars,bers_ffe,'DisplayName','FFE [200,0,0] + PF + MLSE'); -plot(vars,bers_vnle,'DisplayName','VNLE [200,1,0] + PF + MLSE'); -plot(vars,bers_vnledfe,'DisplayName','VNLE [200,1,0] + DFE [2] + PF + MLSE'); -plot(vars,bers_vnledfe_ideal,'DisplayName','VNLE [200,1,0] + ideal DFE [2] + PF + MLSE'); -yline(ber_in_paper); -yline([2e-2, 4.85e-3, 3.8e-3, 2,2e-4],'LineWidth',2,'Color',[0.8,0.8,0.8],'LineStyle',':','HandleVisibility','off'); -ylim([1e-5,0.1]); -beautifyBERplot(); - -%% -------------------- DB target -------------------- -mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels); - -eq_ = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ... - "K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); - -dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ... - "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style); - -fprintf('Paper: %.1e \n \n',ber_in_paper); -dbt_results.metrics.print("description",'Duobinary'); - - -%% - -%ML-based MLSE (L=2) -mu_ml = 0.01; training_epochs = 100; -ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ... - "len_tr",len_tr,"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",1, ... - "traceback_depth",128,"L",3,"delta",4,"adaptive_mu",0); - -[ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode); -ml_mlse_results.metrics.print("description",'ML '); \ No newline at end of file diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m index 986fcdb..21f881f 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m @@ -6,7 +6,7 @@ db = DBHandler("dataBase", "labor_highspeed", "type", database_type); pam_levels = [4, 6, 8]; % three tiles bitrate_set = 360e9; -fiberL = 2; +fiberL = 10; fields = [ db.getTableFieldNames('power_state_info'); @@ -96,7 +96,7 @@ end %% ============================================================ % PLOT — 1×3 (PAM-4, PAM-6, PAM-8) % ============================================================ -fig = figure(9112); clf; +fig = figure(9110); clf; tiledlayout(1,3,'TileSpacing','compact','Padding','compact'); lw = 1.8; @@ -223,19 +223,19 @@ ylabel(''); end -% pos = 1e3.*[2.7770 1.2017 1.4000 0.3200]; -% set(fig, 'Position', pos); +pos = 1e3.*[2.7770 1.2017 1.4000 0.3200]; +set(fig, 'Position', pos); %% === EXPORT === -% outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz'; -% matlab2tikz(outfile, ... -% 'width','\fwidth', ... -% 'height','\fheight', ... -% 'showInfo',false, ... -% 'extraAxisOptions',{ ... -% 'legend style={font=\footnotesize}', ... -% 'legend columns=1' ... -% }); +outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz'; +matlab2tikz(outfile, ... + 'width','\fwidth', ... + 'height','\fheight', ... + 'showInfo',false, ... + 'extraAxisOptions',{ ... + 'legend style={font=\footnotesize}', ... + 'legend columns=1' ... + }); diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m index 147c7d9..2213c8a 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m @@ -1,139 +1,83 @@ -tablename = 'C:\Users\Silas\Documents\latex\JLT_400G_submission\HighSpeedExperiments_oneandonly_csv.csv'; -% Returns a Table -data = readtable(tablename,"Delimiter",';','DecimalSeparator',','); - - %% ============================================================ % PLOT % ============================================================ -%% 1. DATA EXTRACTION & SETUP -%% 1. DATA EXTRACTION & SETUP -raw_M = data.M; -raw_baud = data.BaudRate; -raw_net = data.NetRate; -raw_codes = string(data.ZoteroCode); -raw_names = string(data.Name); -raw_band = string(data.Band); - -% Filter Valid Data -target_M = [2, 4, 6, 8]; -validIdx = ismember(raw_M, target_M) & ~isnan(raw_baud) & ~isnan(raw_net); - -Mvals = raw_M(validIdx); -baud = raw_baud(validIdx); -netrate = raw_net(validIdx); -codes = raw_codes(validIdx); -names = raw_names(validIdx); -bands = raw_band(validIdx); - -pam_list = target_M; -colors = flip(cbrewer2('SET1',4)); - -%% 2. PLOT (For Visual Check only) figure; hold on; -ms = 20; -lw = 0.5; +ms = 32; % scatter size +lw = 0.8; % line width + +for k = 1:4 % PAM-2/4/6/8 -for k = 1:length(pam_list) M = pam_list(k); idxPam = (Mvals == M); - + + % Extract for this PAM x = baud(idxPam); y = netrate(idxPam); - b = bands(idxPam); n = names(idxPam); + + % Get color for this PAM format col = colors(k,:); - + + % ----- LEGEND FLAG (only add one entry per PAM) ----- firstLegend = true; - + + % ---- PLOT ALL POINTS (marker based on publication) ---- for i = 1:sum(idxPam) - % Marker Logic - ms = 20; - if strcmpi(b(i), 'O') - marker = 'o'; - elseif strcmpi(b(i), 'C') - marker = 'd'; - else - marker = 's'; - end - if strcmpi(n(i), 'THIS WORK') - marker = 'pentagram'; - ms = 100; - end - - % Plot Scatter + % marker selection by publication + pubIdx = find(pub_list == n(i), 1); + marker = markerlist{mod(pubIdx-1, nMarkers) + 1}; + if firstLegend - scatter(x(i), y(i), ms, 'Marker', marker, ... - 'MarkerEdgeColor', col, 'MarkerFaceColor', col, ... + h = scatter(x(i), y(i), ms, ... + 'Marker', marker, ... + 'MarkerEdgeColor', col, ... + 'MarkerFaceColor', col, ... 'DisplayName', sprintf('PAM-%d', M)); firstLegend = false; else - scatter(x(i), y(i), ms, 'Marker', marker, ... - 'MarkerEdgeColor', col, 'MarkerFaceColor', col, ... + h = scatter(x(i), y(i), ms, ... + 'Marker', marker, ... + 'MarkerEdgeColor', col, ... + 'MarkerFaceColor', col, ... 'HandleVisibility','off'); end + + % ====== CUSTOM DATATIP CONTENT ====== + dt = h.DataTipTemplate; + dt.DataTipRows(1).Label = 'Baud rate'; + dt.DataTipRows(2).Label = 'Net rate'; + + % Add publication name + dt.DataTipRows(end+1) = dataTipTextRow('Publication', n(i)); + + end - - % Fit lines - if length(x) >= 3 - [p, S, mu] = polyfit(x, y, 2); - xfit = linspace(min(x), max(x), 200); - yfit = polyval(p, xfit, S, mu); - plot(xfit, yfit, '-', 'LineWidth', lw, 'Color', col, 'HandleVisibility', 'off'); + + % ---- Fit (PAM-specific) ---- + valid = ~isnan(x) & ~isnan(y); + if sum(valid) >= 3 + p = polyfit(x(valid), y(valid), 2); + xfit = linspace(min(x(valid)), max(x(valid)), 200); + yfit = polyval(p, xfit); + + plot(xfit, yfit, ':', ... + 'LineWidth', lw, ... + 'Color', col, ... + 'HandleVisibility', 'off'); % do NOT add to legend end end -grid on; box on; +grid on; xlabel('Baud rate [GBd]'); ylabel('Net rate [Gb/s]'); -% title('Check Command Window for TikZ Code'); -% legend('Location','northwest'); -%% 3. GENERATE TIKZ ANNOTATION CODE -% This prints the manual \draw commands to the console +legend('Location','northwest'); +set(gca,'FontSize',11); -%% GENERATE TIKZ ANNOTATION CODE -% This prints the manual \draw commands to the console - -%% GENERATE TIKZ ANNOTATION CODE (Colored Borders + Tiny Font) -%% GENERATE TIKZ ANNOTATION CODE (No Arrow, Close Text) -fprintf('\n\n%% ===========================================================\n'); -fprintf('%% COPY THE FOLLOWING LINES INTO YOUR .TEX FILE \n'); -fprintf('%% (Paste them just before \\end{axis})\n'); -fprintf('%% ===========================================================\n\n'); - -for i = 1:length(baud) - bx = baud(i); - by = netrate(i); - key = codes(i); - M_val = Mvals(i); - - % --- PLACEMENT LOGIC --- - if M_val == 8 - % PAM-8: Place Top-Left - % 'south east' anchor means the text's bottom-right corner touches the coordinate - % shift moves it slightly up and left to clear the marker - anchorStr = 'south east'; - shiftStr = 'shift={(-3pt, 3pt)}'; - else - % Others: Place Bottom-Right - % 'north west' anchor means the text's top-left corner touches the coordinate - % shift moves it slightly down and right - anchorStr = 'north west'; - shiftStr = 'shift={(3pt, -3pt)}'; - end - - % --- PRINT COMMAND --- - % Uses \node directly at the coordinate (axis cs:...) - fprintf('\\node[anchor=%s, %s, font=\\tiny, fill=white, inner sep=1pt] at (axis cs:%.2f, %.2f) {\\cite{%s}};\n', ... - anchorStr, shiftStr, bx, by, key); -end -fprintf('\n') %% === EXPORT === -outfile = 'C:\Users\Silas\Documents\latex\JLT_400G_submission\media\matlab2tikz\highspeedresults_test.tikz'; - +outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\highspeedresults.tikz'; matlab2tikz(outfile, ... 'width','\fwidth', ... 'height','\fheight', ... diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m index 67a2ff8..ab31451 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m @@ -10,11 +10,10 @@ if dsp_options.mode == "load_run_id" if experiment == "highspeed_2024" - dsp_options.database_type = "mysql"; + dsp_options.database_type = 'mysql'; dsp_options.dataBase = 'labor_highspeed';%'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\silas_labor_newdsp_newstructure.db'; dsp_options.storage_path = 'Z:\2024\sioe_labor\'; - db = DBHandler("dataBase", [dsp_options.dataBase],... - "type", dsp_options.database_type,"server","192.168.178.192","user","silas","password","silas"); + db = DBHandler("dataBase", [dsp_options.dataBase], "type", dsp_options.database_type); elseif experiment == "mpi_ecoc_2025" diff --git a/projects/ML_based_MLSE/analyze_filter_length.m b/projects/ML_based_MLSE/analyze_filter_length.m index 1b31d6b..616b715 100644 --- a/projects/ML_based_MLSE/analyze_filter_length.m +++ b/projects/ML_based_MLSE/analyze_filter_length.m @@ -5,9 +5,9 @@ M = 4; randkey = 1; % --- Parameter sweep -order_range = 5:5:50; % FFE order -delta_range = 0:5:20; % delta -SNR_dB = 30; +order_range = 2:3:11; % FFE order +delta_range = 0:2:4; % delta +SNR_dB = 20; % --- Prepare bit sequence order_bits = 19; @@ -21,14 +21,10 @@ Symbols = PAMmapper(M,0).map(Bits); Symbols.fs = 200e9; % --- Channel (minimal ISI + AWGN) -h = abs([0.3 0.9 0.3]); h = h/norm(h); - -% h = [1 -1.67085330039878 1.17918163282514 -0.805210559745616 0.571564213123367 -0.296337147529674 0.00649773445209780 0.0854177610195952 -0.0576009020965258 0.0520994427061551 -0.0624586034913656 0.0553280962699552 -0.00705582559925755 -0.0336399056707792 0.0706903719452810 -0.0334124287931977 0.0131699455037966 0.0587431373842994 -0.0515902976066452 0.00647904355473619 0.0137506750904990 -0.0547974515885928 0.00994735499340592 -0.0135513582534086 -0.00463322575007739 0.0277311946101940]; -% h = h/norm(h); +h = [0.3 0.9 0.3]; h = h/norm(h); symbols_filt = Symbols.filter(h,1); symbols_noi = symbols_filt; symbols_noi.signal = awgn(symbols_filt.signal,SNR_dB,'measured'); -symbols_noi.spectrum(); % --- Generate all parameter pairs [O,D] = ndgrid(order_range, delta_range); @@ -42,7 +38,7 @@ ce_vec = nan(size(pairs,1),1); ce_training = nan(size(pairs,1),training_len); % --- Parallel loop over parameter pairs -for k = 1:size(pairs,1) +parfor k = 1:size(pairs,1) order_k = pairs(k,1); delta_k = pairs(k,2); @@ -93,7 +89,7 @@ end beautifyBERplot ylabel('BER'); xlabel('Filter Order [N]'); title('BER vs. Filter order'); -% ylim([1e-4, 0.1]); +ylim([1e-4, 0.1]); yline(3.8e-3,'HandleVisibility','off'); yline(2.2e-4,'HandleVisibility','off'); diff --git a/projects/ML_based_MLSE/theoretic_channel_evaluation.m b/projects/ML_based_MLSE/theoretic_channel_evaluation.m index 3eb6f09..3e1efcc 100644 --- a/projects/ML_based_MLSE/theoretic_channel_evaluation.m +++ b/projects/ML_based_MLSE/theoretic_channel_evaluation.m @@ -48,36 +48,36 @@ for i = 1:numel(SNR_dB) symbols_noi = symbols_filt; symbols_noi.signal = awgn(symbols_filt.signal, SNR_dB(i), 'measured'); % AWGN with given SNR - % Sequence Est L=5 - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',h); - mlse_.DIR = h; - [y_mlse] = mlse_.process(symbols_noi,Symbols); - mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); - [~, ~, ber_mlse_l5(i), ~] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); - fprintf('MLSE L5: %.2e \n',ber_mlse_l5(i)); - - % 2nd Approach - mu_lms = 0.0005; - pf_ncoeffs = 1; - eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",16,"sps",1,"dd_mode",1,"adaption_technique","lms"); - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - - % FFE - [y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols); - - Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe); - [~, ~, ber_ffe(i), ~] = calc_ber(Eq_bits.signal, Bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); - fprintf('FFE: %.2e \n',ber_ffe(i)); - - % Postfilter - [y_white,~] = pf_.process(y_ffe, ffe_noise); - - % Sequence Est - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',pf_.coefficients); - [y_mlse] = mlse_.process(y_white,Symbols); - mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); - [~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); - fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i)); + % % Sequence Est L=5 + % mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',h); + % mlse_.DIR = h; + % [y_mlse] = mlse_.process(symbols_noi,Symbols); + % mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); + % [~, ~, ber_mlse_l5(i), ~] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + % fprintf('MLSE L5: %.2e \n',ber_mlse_l5(i)); + % + % % 2nd Approach + % mu_lms = 0.0005; + % pf_ncoeffs = 1; + % eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",16,"sps",1,"dd_mode",1,"adaption_technique","lms"); + % pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + % + % % FFE + % [y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols); + % + % Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe); + % [~, ~, ber_ffe(i), ~] = calc_ber(Eq_bits.signal, Bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); + % fprintf('FFE: %.2e \n',ber_ffe(i)); + % + % % Postfilter + % [y_white,~] = pf_.process(y_ffe, ffe_noise); + % + % % Sequence Est + % mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',pf_.coefficients); + % [y_mlse] = mlse_.process(y_white,Symbols); + % mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); + % [~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + % fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i)); % ML-base MLSE L=2 adaptive_mu = 0; diff --git a/projects/Messung_Zürich/minimal_matched.m b/projects/Messung_Zürich/minimal_matched.m deleted file mode 100644 index a5dfa33..0000000 --- a/projects/Messung_Zürich/minimal_matched.m +++ /dev/null @@ -1,60 +0,0 @@ -L = 4; % Oversampling factor -fsym = 1e9/4; -rollOff = 0.5; % Pulse shaping roll-off factor -htx = rcosdesign(rollOff,16, L, 'sqrt'); -hrx = conj(fliplr(htx)); - -signal = zeros(1000,1); -signal(500) = 1; -Digi_sig = Informationsignal(signal,"fs",fsym); - -% Digi_sig = Digi_sig.resample("fs_out",L*fsym); -% -% Digi_sig_tx = Digi_sig.filter(htx,1); -% Digi_sig_rx = Digi_sig_tx.filter(hrx,1); -% -% figure;plot(Digi_sig_tx.signal);hold on;plot(Digi_sig_rx.signal),plot(Digi_sig.signal); -% -% -% -% -pulsef = Pulseformer("alpha",1,"matched",0,"fdac",Digi_sig.fs*L,"fsym",fsym,"pulse","rrc","pulselength",16); -Digi_sig = pulsef.process(Digi_sig); - -pulsef = Pulseformer("alpha",1,"matched",1,"fdac",Digi_sig.fs,"fsym",fsym,"pulse","rrc","pulselength",16); -Digi_sig_matched = pulsef.process(Digi_sig); - -% Filter: -htx = rcosdesign(rollOff,16, L, 'sqrt'); -% Note half of the target delay is used, because when combined -% to the matched filter, the total delay will be achieved. -hrx = conj(fliplr(htx)); - -figure -plot(htx) -title('Transmit Filter') -xlabel('Index') -ylabel('Amplitude') - -figure -plot(hrx) -title('Rx Filter (Matched Filter)') -xlabel('Index') -ylabel('Amplitude') - -p = conv(htx,hrx); - -figure -plot(p) -title('Combined Tx-Rx = Raised Cosine') -xlabel('Index') -ylabel('Amplitude') - -% And let's highlight the zero-crossings -zeroCrossings = NaN*ones(size(p)); -zeroCrossings(1:L:end) = 0; -zeroCrossings((rcDelay)*L + 1) = NaN; % Except for the central index -hold on -plot(zeroCrossings, 'o') -legend('RC Pulse', 'Zero Crossings') -hold off \ No newline at end of file diff --git a/projects/WDM/WDM_model_10km_queue.m b/projects/WDM/WDM_model_10km_queue.m index 8aae2ea..a41cb02 100644 --- a/projects/WDM/WDM_model_10km_queue.m +++ b/projects/WDM/WDM_model_10km_queue.m @@ -312,7 +312,7 @@ for realiz = 1:s.num_realiz Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... - "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",1,"useSingle",1).process(Opt_sig_wdm_fib); + "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib); % -------- Evaluate at intermediate distance (enqueue jobs) -------- if eval_ptr <= nEval && seg == eval_seg(eval_ptr) diff --git a/test/bayesopt_ffe_tuning.m b/test/bayesopt_ffe_tuning.m deleted file mode 100644 index adedb3d..0000000 --- a/test/bayesopt_ffe_tuning.m +++ /dev/null @@ -1,219 +0,0 @@ -%% Bayesian Optimization for FFE Parameter Tuning -% This script uses bayesopt to find optimal mu_dd and mu_tr values -% that minimize BER for the FFE equalizer. - -clear; clc; - -%% Setup - Same as gpu_processing_dpfiber.m -s.wavelengthplan = calcWavelengthPlan(4, 400e9, 1310); -link_length = 10; -s.pmd = 0.1; -s.gamma = 0.0023; - -s.M = 4; -fsym = 112e9; -fdac = 2*fsym; -fadc = 120000000000; -s.random_key = 1; - -% Laser / Modulator -vbias_rel = 0.5; -u_pi = 4.6; -vbias = -vbias_rel*u_pi; -laser_linewidth = 0e6; - -duob_mode = db_mode.no_db; -rcalpha = 0.05; -Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); - -s.chirpalpha = 0; -s.p_launch = 3; -s.p = "co"; - -N = numel(s.wavelengthplan); - -switch s.p - case "co" - pol_rot = 100.*ones(1,N); - d_local = 0; -end - -f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9); -margin = 25e12; -f_span = (max(f_plan)+margin)-(min(f_plan)-margin); -f_nyq = f_span/2; - -kover = 4; -upsample_required = f_nyq./(fdac*kover/2); -upsample_pow = 2^nextpow2(upsample_required); - -s.f_opt = fdac*kover*upsample_pow; -s.f_opt_nyq = s.f_opt/2; - -s.rop = -8; % Fixed ROP for optimization - -%% Generate TX signals (run once) -fprintf('Generating TX signals...\n'); -for l = 1:N - [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ... - "fsym",fsym,"M",s.M,"order",15,"useprbs",0, ... - "fs_out",fdac, ... - "applyclipping",0,"clipfactor",1.5, ... - "applypulseform",1,"pulseformer",Pform, ... - "randkey",s.random_key+l, ... - "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ... - ).process(); - - Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ... - "filterType",filtertypes.gaussian,"active",true); - - El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ... - "bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ... - "H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig); - - El_sig = El_sig.normalize("mode","oneone"); - scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); - El_sig = El_sig .* scaling; - - Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ... - "lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ... - "linewidth",laser_linewidth,"randomkey",s.random_key+l,"alpha",s.chirpalpha).process(El_sig); - - signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out); -end - -%% WDM mux + launch -Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ... - "lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell); - -Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... - "amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm); - -%% Fiber propagation -segment_length = 1; -nSegments = link_length/segment_length; -nSegments = round(nSegments); - -zdw = 1310; -randomize_D = true; -Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key); - -Opt_sig_wdm_fib = Opt_sig_wdm; -fprintf('Running fiber propagation...\n'); -for seg = 1:nSegments - fprintf('Segment %d/%d\n', seg, nSegments); - Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... - "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... - "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... - "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",true).process(Opt_sig_wdm_fib); -end - -%% Pre-process to get Rx_sig (do demux once) -fprintf('Pre-processing receiver chain...\n'); -l = 1; % Use channel 1 for optimization - -Opt_sig_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ... - "fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib); - -Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... - "amplification_db",s.rop).process(Opt_sig_demux{l}); - -PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ... - "nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx); - -rx_bwl = 100e9; -PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ... - "filterType",filtertypes.butterworth,"active",true).process(PD_sig); - -Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); -Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc, ... - "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth, ... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0, ... - "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig); - -Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym); - -[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0); -Rx_sig = Scpe_cell{1}; -Rx_sig = Rx_sig.normalize("mode","rms"); - -fprintf('Receiver pre-processing complete. Ready for optimization.\n\n'); - -%% Define the objective function for bayesopt -function ber = ffe_objective(params, Rx_sig, Symbols_l, Tx_bits_l, M, duob_mode) -mu_dd = params.mu_dd; -mu_tr = params.mu_tr; - -try - eq_ffe = FFE("epochs_tr", 5, "epochs_dd", 2, "len_tr", 2^13, ... - "mu_dd", mu_dd, "mu_tr", mu_tr, ... - "order", 50, "sps", 2, "decide", 0, ... - "adaption", adaption_method.nlms, "dd_mode", 1); - - ffe_results = ffe(eq_ffe, M, Rx_sig, Symbols_l, Tx_bits_l, ... - "precode_mode", duob_mode, ... - 'showAnalysis', 0, ... - "postFFE", [], ... - "eth_style_symbol_mapping", 0); - - ber = ffe_results.metrics.BER; - - if ber == 0 - ber = 1e-10; - end - - if ~isfinite(ber) - ber = 0.5; - end - - fprintf(' mu_dd=%.4e, mu_tr=%.4e -> BER=%.4e\n', mu_dd, mu_tr, ber); - -catch ME - fprintf(' mu_dd=%.4e, mu_tr=%.4e -> FAILED (%s)\n', mu_dd, mu_tr, ME.message); - ber = 0.5; -end -end - -%% Define optimizable variables -mu_dd_var = optimizableVariable('mu_dd', [1e-5, 0.1], 'Transform', 'log'); -mu_tr_var = optimizableVariable('mu_tr', [1e-5, 0.1], 'Transform', 'log'); - -%% Run Bayesian Optimization -fprintf('========== Starting Bayesian Optimization ==========\n'); -fprintf('Optimizing mu_dd and mu_tr to minimize BER\n'); -fprintf('Search range: mu_dd=[1e-5, 0.1], mu_tr=[1e-5, 0.1]\n\n'); - -objective_fn = @(params) ffe_objective(params, Rx_sig, Symbols{l}, Tx_bits{l}, s.M, duob_mode); - -results = bayesopt(objective_fn, [mu_dd_var, mu_tr_var], ... - 'MaxObjectiveEvaluations', 30, ... - 'AcquisitionFunctionName', 'expected-improvement-plus', ... - 'IsObjectiveDeterministic', false, ... - 'ExplorationRatio', 0.5, ... - 'Verbose', 1, ... - 'PlotFcn', []); - -%% Display Results -fprintf('\n========== FFE Optimization Complete ==========\n'); -fprintf('Best FFE parameters found:\n'); -fprintf(' mu_dd = %.6e\n', results.XAtMinObjective.mu_dd); -fprintf(' mu_tr = %.6e\n', results.XAtMinObjective.mu_tr); -fprintf(' BER = %.6e\n', results.MinObjective); - -%% Verify with optimal parameters -fprintf('\nVerifying optimal FFE parameters...\n'); -best_mu_dd = results.XAtMinObjective.mu_dd; -best_mu_tr = results.XAtMinObjective.mu_tr; - -eq_ffe_best = FFE("epochs_tr", 5, "epochs_dd", 2, "len_tr", 2^13, ... - "mu_dd", best_mu_dd, "mu_tr", best_mu_tr, ... - "order", 50, "sps", 2, "decide", 0, ... - "adaption", adaption_method.nlms, "dd_mode", 1); - -ffe_results_best = ffe(eq_ffe_best, s.M, Rx_sig, Symbols{l}, Tx_bits{l}, ... - "precode_mode", duob_mode, ... - 'showAnalysis', 1, ... - "postFFE", [], ... - "eth_style_symbol_mapping", 0); - -fprintf('\nFinal FFE BER with optimal parameters: %.6e\n', ffe_results_best.metrics.BER); diff --git a/test/bitwise_demapping_pam6.m b/test/bitwise_demapping_pam6.m deleted file mode 100644 index 6ea0b3d..0000000 --- a/test/bitwise_demapping_pam6.m +++ /dev/null @@ -1,43 +0,0 @@ - - -%%%%% SETTINGS %%%%%% -useprbs = 1; -M = 8; -randkey = 1; -fsym = 112e9; -viewresults = 0; - -%%%%% Mapping %%%%% -M = 6; -data = 0:M-1; -bitpersymbol = log2(M); - -s = RandStream('twister','Seed',1); -bitpattern = randi(s,[0 1], 2^18, 1); -bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); - -bits_tx = Informationsignal(bitpattern); - -symbols = PAMmapper(M,0).map(bits); - -pam6transitions = combvec(PAMmapper(M,0).levels,PAMmapper(M,0).levels)'; - -pam6bits = PAMmapper(6,0,"eth_style",0).demap(reshape(pam6transitions',[],1)./sqrt(10)); -symbols_rx = PAMmapper(M,0).map(pam6bits).*sqrt(10); - -pam6bits = reshape(pam6bits',5,[])'; - -figure; hold on -scatter(pam6transitions(:,1), pam6transitions(:,2), 'x', 'LineWidth', 1); -n = size(pam6transitions,1); -labels = cellstr(char(pam6bits + '0')); % -> N x 1 cell array of char rows -text(pam6transitions(:,1), pam6transitions(:,2), labels, ... - 'HorizontalAlignment','left', 'VerticalAlignment','bottom'); - - - -bits_rx = PAMmapper(M,0).demap(symbols); - -[~,error_num,ber,error_pos] = calc_ber(bits_tx.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); - -PAMmapper(8,0).showBitMapping \ No newline at end of file diff --git a/test/duobinary_minimal_example.m b/test/duobinary_minimal_example.m index c39183e..d61c044 100644 --- a/test/duobinary_minimal_example.m +++ b/test/duobinary_minimal_example.m @@ -29,6 +29,9 @@ para.skip =0; para.bruijn = 0; para.reset_prms = 0; para.method = 1; +para.pcs = 0; +para.shape_para = 0; +para.rng_num = 0; data_in = []; global loop; diff --git a/test/gpu_cpu_comparison.m b/test/gpu_cpu_comparison.m deleted file mode 100644 index 02d40ba..0000000 --- a/test/gpu_cpu_comparison.m +++ /dev/null @@ -1,349 +0,0 @@ -%% GPU vs CPU Comparison Test for DP_Fiber -% This script runs the fiber simulation with and without GPU acceleration -% and compares the numerical results. - -% clear; clc; - -%% Setup (same as gpu_processing_dpfiber.m but simplified) -s.wavelengthplan = calcWavelengthPlan(4, 400e9, 1310); -link_length = 10; -s.pmd = 0.1; -s.gamma = 0.0023; - -s.M = 4; -fsym = 112e9; -fdac = 2*fsym; -fadc = 120000000000; -s.random_key = 1; - -% Laser / Modulator -vbias_rel = 0.5; -u_pi = 4.6; -vbias = -vbias_rel*u_pi; -laser_linewidth = 0e6; - -% DB Stuff -duob_mode = db_mode.no_db; - -rcalpha = 0.05; -Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); - -s.chirpalpha = 0; -s.p_launch = 3; -s.p = "co"; - -N = numel(s.wavelengthplan); - -switch s.p - case "co" - pol_rot = 100.*ones(1,N); - d_local = 0; - case "pair" - pol_rot = repmat([100,100,0,0],1,N/4); - d_local = 0; - case "alt" - pol_rot = repmat([100,0,100,0],1,N/4); - d_local = 0; - case "seg" - pol_rot = 100.*ones(1,N); - d_local = 3; - otherwise - error('Unknown fwm_mitigation_technique: %s', string(s.p)); -end - -f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9); -margin = 5e12; -f_span = (max(f_plan)+margin)-(min(f_plan)-margin); -f_nyq = f_span/2; - -kover = 4; -upsample_required = f_nyq./(fdac*kover/2); -upsample_pow = 2^nextpow2(upsample_required); - -s.f_opt = fdac*kover*upsample_pow; -s.f_opt_nyq = s.f_opt/2; - -%% TX per channel -for l = 1:N - [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ... - "fsym",fsym,"M",s.M,"order",15,"useprbs",0, ... - "fs_out",fdac, ... - "applyclipping",0,"clipfactor",1.5, ... - "applypulseform",1,"pulseformer",Pform, ... - "randkey",s.random_key+l, ... - "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ... - ).process(); - - Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ... - "filterType",filtertypes.gaussian,"active",true); - - El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ... - "bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ... - "H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig); - - clear Digi_sig - - El_sig = El_sig.normalize("mode","oneone"); - scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); - El_sig = El_sig .* scaling; - - Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ... - "lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ... - "linewidth",laser_linewidth,"randomkey",s.random_key+l,"alpha",s.chirpalpha).process(El_sig); - - clear El_sig - - signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out); - - clear Eml_out Lp_awg -end - -disp('Signal generated for all channels.'); - -%% WDM mux + launch -Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ... - "lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell); - -Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... - "amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm); - -% Save input for both runs -Opt_sig_input = Opt_sig_wdm; - -segment_length = 1; -nSegments = link_length/segment_length; -if abs(nSegments - round(nSegments)) > 1e-12 - error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length); -end -nSegments = round(nSegments); - -zdw = 1310; -randomize_D = true; - -if nSegments > 0 - Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key); -else - Dvec = []; -end - -%% Run WITHOUT GPU -fprintf('\n========== Running WITHOUT GPU (CPU) ==========\n'); -Opt_sig_cpu = Opt_sig_input; - -tic; -for seg = 1:nSegments - fprintf('CPU Segment %d/%d \n',seg, nSegments); - - Opt_sig_cpu = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... - "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... - "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... - "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",false).process(Opt_sig_cpu); -end -time_cpu = toc; -fprintf('CPU Time: %.3f seconds\n', time_cpu); - -%% Run WITH GPU (Double Precision) -fprintf('\n========== Running WITH GPU (Double Precision) ==========\n'); -Opt_sig_gpu_double = Opt_sig_input; - -tic; -for seg = 1:nSegments - fprintf('GPU-Double Segment %d/%d \n',seg, nSegments); - - Opt_sig_gpu_double = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... - "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... - "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... - "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",false).process(Opt_sig_gpu_double); -end -time_gpu_double = toc; -fprintf('GPU Double Time: %.3f seconds\n', time_gpu_double); - -%% Run WITH GPU (Single Precision) -fprintf('\n========== Running WITH GPU (Single Precision) ==========\n'); -Opt_sig_gpu_single = Opt_sig_input; - -tic; -for seg = 1:nSegments - fprintf('GPU-Single Segment %d/%d \n',seg, nSegments); - - Opt_sig_gpu_single = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... - "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... - "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... - "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",true).process(Opt_sig_gpu_single); -end -time_gpu_single = toc; -fprintf('GPU Single Time: %.3f seconds\n', time_gpu_single); - -%% Compare Results -fprintf('\n========== Numerical Comparison ==========\n'); - -sig_cpu = Opt_sig_cpu.signal; -sig_gpu_double = Opt_sig_gpu_double.signal; -sig_gpu_single = Opt_sig_gpu_single.signal; - -% Check dimensions -fprintf('CPU signal size: [%d x %d]\n', size(sig_cpu,1), size(sig_cpu,2)); -fprintf('GPU Double signal size: [%d x %d]\n', size(sig_gpu_double,1), size(sig_gpu_double,2)); -fprintf('GPU Single signal size: [%d x %d]\n', size(sig_gpu_single,1), size(sig_gpu_single,2)); - -% CPU vs GPU Double -fprintf('\n--- CPU vs GPU Double ---\n'); -diff_cpu_double = abs(sig_cpu - sig_gpu_double); -max_diff_cpu_double = max(diff_cpu_double(:)); -mean_diff_cpu_double = mean(diff_cpu_double(:)); -rel_diff_cpu_double = max_diff_cpu_double / max(abs(sig_cpu(:))); -fprintf('Max absolute difference: %.6e\n', max_diff_cpu_double); -fprintf('Mean absolute difference: %.6e\n', mean_diff_cpu_double); -fprintf('Max relative difference: %.6e\n', rel_diff_cpu_double); - -% CPU vs GPU Single -fprintf('\n--- CPU vs GPU Single ---\n'); -diff_cpu_single = abs(sig_cpu - sig_gpu_single); -max_diff_cpu_single = max(diff_cpu_single(:)); -mean_diff_cpu_single = mean(diff_cpu_single(:)); -rel_diff_cpu_single = max_diff_cpu_single / max(abs(sig_cpu(:))); -fprintf('Max absolute difference: %.6e\n', max_diff_cpu_single); -fprintf('Mean absolute difference: %.6e\n', mean_diff_cpu_single); -fprintf('Max relative difference: %.6e\n', rel_diff_cpu_single); - -% GPU Double vs GPU Single -fprintf('\n--- GPU Double vs GPU Single ---\n'); -diff_double_single = abs(sig_gpu_double - sig_gpu_single); -max_diff_double_single = max(diff_double_single(:)); -mean_diff_double_single = mean(diff_double_single(:)); -rel_diff_double_single = max_diff_double_single / max(abs(sig_gpu_double(:))); -fprintf('Max absolute difference: %.6e\n', max_diff_double_single); -fprintf('Mean absolute difference: %.6e\n', mean_diff_double_single); -fprintf('Max relative difference: %.6e\n', rel_diff_double_single); - -% Check tolerances -fprintf('\n--- Tolerance Check ---\n'); -tol_double = 1e-10; -tol_single = 1e-5; % Single precision has ~7 significant digits - -if max_diff_cpu_double < tol_double - fprintf('✓ CPU vs GPU Double: EQUIVALENT (diff < %.0e)\n', tol_double); -else - fprintf('✗ CPU vs GPU Double: DIFFER beyond tolerance (%.0e)\n', tol_double); -end - -if max_diff_cpu_single < tol_single - fprintf('✓ CPU vs GPU Single: ACCEPTABLE (diff < %.0e)\n', tol_single); -else - fprintf('⚠ CPU vs GPU Single: Precision loss detected (diff = %.2e, tol = %.0e)\n', max_diff_cpu_single, tol_single); -end - -% Performance comparison -fprintf('\n========== Performance Summary ==========\n'); -fprintf('CPU Time: %.3f s\n', time_cpu); -fprintf('GPU Double Time: %.3f s\n', time_gpu_double); -fprintf('GPU Single Time: %.3f s\n', time_gpu_single); -fprintf('\n'); -fprintf('Speedup (GPU Double vs CPU): %.2fx\n', time_cpu/time_gpu_double); -fprintf('Speedup (GPU Single vs CPU): %.2fx\n', time_cpu/time_gpu_single); -fprintf('Speedup (GPU Single vs GPU Double): %.2fx\n', time_gpu_double/time_gpu_single); - -%% ========== BER Comparison ========== -% Process each fiber output through simplified receiver to check if -% single-precision affects actual BER performance - -fprintf('\n========== BER Comparison ==========\n'); -fprintf('Processing signals through receiver chain...\n'); - -% Receiver parameters -rop = -7; % Received optical power [dBm] -len_tr = 4096; % Training length -mu_dc = 0.005; -mu_ffe = [0.0001 0.0008 0.001]; -mu_dfe = 0.0004; - -% Helper function to process through receiver and get BER -function ber = process_receiver(Opt_sig_fib, l, Symbols, Tx_bits, ... - fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode) - -% Demux single channel -Opt_sig_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ... - "fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_fib); - -% ROP amplifier -Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... - "amplification_db",rop).process(Opt_sig_demux{l}); - -% Photodiode -PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ... - "nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx); - -% Low-pass filter -rx_bwl = 100e9; -PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ... - "filterType",filtertypes.butterworth,"active",true).process(PD_sig); - -% Scope -Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); -Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc, ... - "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth, ... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0, ... - "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig); - -% Resample to 2 sps -Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym); - -% Time sync -[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0); -Rx_sig = Scpe_cell{1}; -Rx_sig = Rx_sig.normalize("mode","rms"); - -% FFE Equalizer -ffe_order = [50, 0, 0]; -eq_ffe = EQ("Ne",ffe_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ... - "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); - -ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols{l},Tx_bits{l}, ... - "precode_mode",duob_mode, ... - 'showAnalysis',0, ... - "postFFE",[], ... - "eth_style_symbol_mapping",0); - -ber = ffe_results.metrics.BER; -end - -% Process each mode for channel 1 -l = 4; % Use first channel for comparison - -fprintf('Processing CPU result...\n'); -ber_cpu = process_receiver(Opt_sig_cpu, l, Symbols, Tx_bits, ... - fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode); - -fprintf('Processing GPU Double result...\n'); -ber_gpu_double = process_receiver(Opt_sig_gpu_double, l, Symbols, Tx_bits, ... - fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode); - -fprintf('Processing GPU Single result...\n'); -ber_gpu_single = process_receiver(Opt_sig_gpu_single, l, Symbols, Tx_bits, ... - fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode); - -% Display BER results -fprintf('\n========== BER Results (Channel %d, ROP = %d dBm) ==========\n', l, rop); -fprintf('CPU: BER = %.4e\n', ber_cpu); -fprintf('GPU Double: BER = %.4e\n', ber_gpu_double); -fprintf('GPU Single: BER = %.4e\n', ber_gpu_single); - -fprintf('\n--- BER Comparison ---\n'); -if ber_cpu == 0 && ber_gpu_double == 0 && ber_gpu_single == 0 - fprintf('✓ All BERs are zero (no errors detected)\n'); -else - ber_diff_double = abs(ber_cpu - ber_gpu_double); - ber_diff_single = abs(ber_cpu - ber_gpu_single); - fprintf('|BER_cpu - BER_gpu_double| = %.4e\n', ber_diff_double); - fprintf('|BER_cpu - BER_gpu_single| = %.4e\n', ber_diff_single); - - if ber_diff_double < 1e-6 && ber_diff_single < 1e-6 - fprintf('✓ BER differences are negligible\n'); - elseif ber_diff_single > ber_diff_double * 10 - fprintf('⚠ Single precision shows measurable BER impact\n'); - else - fprintf('✓ BER differences within acceptable range\n'); - end -end - -fprintf('\n========== Test Complete ==========\n'); diff --git a/test/gpu_processing_dpfiber.m b/test/gpu_processing_dpfiber.m deleted file mode 100644 index 317833e..0000000 --- a/test/gpu_processing_dpfiber.m +++ /dev/null @@ -1,265 +0,0 @@ - - -s.wavelengthplan = calcWavelengthPlan(16, 400e9, 1310); -N = numel(s.wavelengthplan); -link_length = 10; -s.pmd = 0.1;%0.1; -s.gamma = 0.0023; - -s.M = 4; -fsym = 112e9; -fdac = 2*fsym; -fadc = 120000000000; -s.random_key = 1; - -% Laser / s.Modulator -vbias_rel = 0.5; -u_pi = 4.6; -vbias = -vbias_rel*u_pi; -laser_linewidth = 0e6; - -% DB Stuff -duob_mode = db_mode.no_db; - -rcalpha = 0.05; -Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); - -s.chirpalpha = 0; - -s.p_launch = 3; -s.p = "co"; - -switch s.p - case "co" - pol_rot = 100.*ones(1,N); - d_local = 0; - case "pair" - pol_rot = repmat([100,100,0,0],1,N/4); - d_local = 0; - case "alt" - pol_rot = repmat([100,0,100,0],1,N/4); - d_local = 0; - case "seg" - pol_rot = 100.*ones(1,N); - d_local = 3; - otherwise - error('Unknown fwm_mitigation_technique: %s', string(s.p)); -end - -f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9); -margin = 25e12; % some THz left and right -f_span = (max(f_plan)+margin)-(min(f_plan)-margin); -f_nyq = f_span/2; - -kover = 4; -upsample_required = f_nyq./(fdac*kover/2); -upsample_pow = 2^nextpow2(upsample_required); - -s.f_opt = fdac*kover*upsample_pow; -s.f_opt_nyq = s.f_opt/2; - -s.rop = -10:1:0; - -profile on - -%% ---------- TX per channel ---------- -for l = 1:N - - [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ... - "fsym",fsym,"M",s.M,"order",17,"useprbs",0, ... - "fs_out",fdac, ... - "applyclipping",0,"clipfactor",1.5, ... - "applypulseform",1,"pulseformer",Pform, ... - "randkey",s.random_key+l, ... - "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ... - ).process(); - - Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ... - "filterType",filtertypes.gaussian,"active",true); - - El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ... - "bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ... - "H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig); - - % Digi_sig not needed after AWG - clear Digi_sig - - % Electrical Driver Amplifier - El_sig = El_sig.normalize("mode","oneone"); - scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); - El_sig = El_sig .* scaling; - - % E/O Conversion - Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ... - "lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ... - "linewidth",laser_linewidth,"randomkey",s.random_key+l,"alpha",s.chirpalpha).process(El_sig); - - % El_sig not needed after EML - clear El_sig - - signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out); - - % Eml_out not needed after pol controller - clear Eml_out Lp_awg -end - -disp('Signal generated for all channels.'); - -%% ---------- WDM mux + launch ---------- -Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ... - "lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell); - -Opt_sig_wdm.spectrum(); - -Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... - "amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm); - -Opt_sig_wdm_fib = Opt_sig_wdm; - -segment_length = 1; -nSegments = link_length/segment_length; -if abs(nSegments - round(nSegments)) > 1e-12 - error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length); -end -nSegments = round(nSegments); - -zdw = 1310; -randomize_D = true; - -% Guard for 0 km: avoid calling getDispersionVector(0,...) if it doesn't support it -if nSegments > 0 - Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key); -else - Dvec = []; -end - -for seg = 1:nSegments - - fprintf('Segment %d/%d \n',seg, nSegments); - - Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... - "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... - "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... - "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",1).process(Opt_sig_wdm_fib); - -end - -profile off -profile viewer - -Opt_sig_wdm_fib.spectrum(); - -%% ========== BER Evaluation ========== -fprintf('\n========== BER Evaluation ==========\n'); -fprintf('Processing signals through receiver chain...\n'); - -% Receiver parameters -len_tr = 4096; % Training length -mu_dc = 0.005; -mu_ffe = [0.0001 0.0008 0.001]; -mu_dfe = 0.0004; - -% Helper function to process through receiver and get BER -function ber = process_receiver(Opt_sig_fib, l, Symbols, Tx_bits, ... - fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode) - -% Demux single channel -Opt_sig_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ... - "fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_fib); - -% ROP amplifier -Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... - "amplification_db",rop).process(Opt_sig_demux{l}); - -% Photodiode -PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ... - "nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx); - -% Low-pass filter -rx_bwl = 100e9; -PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ... - "filterType",filtertypes.butterworth,"active",true).process(PD_sig); - -% Scope -Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); -Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc, ... - "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth, ... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0, ... - "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig); - -% Resample to 2 sps -Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym); - -% Time sync -[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0); -Rx_sig = Scpe_cell{1}; -Rx_sig = Rx_sig.normalize("mode","rms"); - -% FFE Equalizer -ffe_order = [50, 0, 0]; -eq_ffe = EQ("Ne",ffe_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ... - "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); -eq_ffe = FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",6.624e-05,"mu_tr",0.058136,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1); - -ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols{l},Tx_bits{l}, ... - "precode_mode",duob_mode, ... - 'showAnalysis',0, ... - "postFFE",[], ... - "eth_style_symbol_mapping",0); - -ber = ffe_results.metrics.BER; -end - -% Process each ROP value for selected channels using parfor -ber_results = zeros(length(s.rop), N); - -% Flatten loop for parfor: iterate over all (ROP, channel) combinations -num_rop = length(s.rop); -rop_vals = s.rop; -ber_flat = zeros(num_rop * N, 1); - -parfor idx = 1:(num_rop * N) - % Convert linear index to (ri, l) subscripts - ri = ceil(idx / N); - l = mod(idx - 1, N) + 1; - - fprintf('ROP %d dBm, Channel %d/%d\n', rop_vals(ri), l, N); - ber_flat(idx) = process_receiver(Opt_sig_wdm_fib, l, Symbols, Tx_bits, ... - fdac, kover, upsample_pow, fsym, fadc, rop_vals(ri), s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode); -end - -% Reshape back to [num_rop × N] matrix -ber_results = reshape(ber_flat, [N, num_rop]).'; - - -%% Display BER Results -fprintf('\n========== BER Results ==========\n'); -fprintf('ROP [dBm] | '); -for l = 1:N - fprintf('Ch%d | ', l); -end -fprintf('\n'); - -for ri = 1:length(s.rop) - fprintf('%8d | ', s.rop(ri)); - for l = 1:N - fprintf('%.2e | ', ber_results(ri, l)); - end - fprintf('\n'); -end - -% Plot BER vs ROP -figure; -semilogy(s.rop, mean(ber_results, 2), '-o', 'LineWidth', 2); -hold on; -for l = 1:N - semilogy(s.rop, ber_results(:, l), '--', 'LineWidth', 1); -end -hold off; -xlabel('ROP [dBm]'); -ylabel('BER'); -title('BER vs Received Optical Power (GPU Single Precision)'); -legend(['Mean', arrayfun(@(x) sprintf('Ch%d', x), 1:N, 'UniformOutput', false)]); -grid on; - -fprintf('\n========== Test Complete ==========\n'); diff --git a/test/gpu_processing_test.m b/test/gpu_processing_test.m deleted file mode 100644 index b08c03f..0000000 --- a/test/gpu_processing_test.m +++ /dev/null @@ -1,70 +0,0 @@ -% 1. Setup Data OUTSIDE the timer -d = gpuDevice; -N = 10000; - -fprintf('Preparing data...\n'); -A_cpu_double = rand(N, N); % Create double on CPU -A_cpu_single = single(A_cpu_double); % Create single on CPU - -% Warmup run (wakes up the GPU from idle state) -A_warm = gpuArray.rand(1000, 1000, 'single'); -B_warm = A_warm * A_warm; -wait(d); - -fprintf('------------------------------------------------\n'); - -% TEST 1: Double Precision (The "Slow" way) -% We move data to GPU first so we only measure calculation time -A_gpu_double = gpuArray(A_cpu_double); -wait(d); % Ensure transfer is done before starting timer - -fprintf('Running DOUBLE precision test... '); -tic; -B_gpu = A_gpu_double * A_gpu_double; -wait(d); % FORCE MATLAB TO WAIT FOR GPU -time_double = toc; -fprintf('Done.\n'); -fprintf('Double Precision Time: %.4f seconds\n', time_double); - -% TEST 2: Single Precision (The "Fast" way) -A_gpu_single = gpuArray(A_cpu_single); -wait(d); % Ensure transfer is done - -fprintf('Running SINGLE precision test... '); -tic; -B_gpu = A_gpu_single * A_gpu_single; -wait(d); % FORCE MATLAB TO WAIT FOR GPU -time_single = toc; -fprintf('Done.\n'); -fprintf('Single Precision Time: %.4f seconds\n', time_single); - -% Calculate Speedup -fprintf('------------------------------------------------\n'); -fprintf('Speedup Factor using single precision: %.2fx\n', time_double / time_single); - -% -% Create 10,000 small matrices (10x10) stacked in a 3D array -A_stack = gpuArray.rand(10, 10, 10000, 'single'); -B_stack = gpuArray.rand(10, 10, 10000, 'single'); - -% BAD: Looping (GPU overhead kills you) -tic; -for i=1:10000 - C(:,:,i) = A_stack(:,:,i) * B_stack(:,:,i); -end -wait(d); -loop_time = toc; - -% GOOD: Pagefun (Executes all 10,000 mults simultaneously) -tic; -C_stack = pagefun(@mtimes, A_stack, B_stack); -wait(d); -pagefun_time = toc; - -fprintf('------------------------------------------------\n'); -fprintf('Speedup Factor using pagefun: %.2fx\n', loop_time / pagefun_time); - - -fprintf('Total VRAM: %.2f GB\n', d.TotalMemory / 1e9); -fprintf('Available VRAM: %.2f GB\n', d.AvailableMemory / 1e9); -fprintf('Usage: %.1f%%\n', 100 * (1 - d.AvailableMemory / d.TotalMemory)); \ No newline at end of file