From 9d63457030eae95631f5168842334019e66dfc95 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Fri, 26 Sep 2025 14:25:52 +0200 Subject: [PATCH] ECOC Theroy stuff --- Classes/04_DSP/Sequence Detection/MLSE.m | 6 +- projects/ECOC_2025/load_signal_standalone.m | 42 +-- .../plots_from_database/plot_mpi_trial.m | 86 +++--- .../theory/analytic_mpi_evaluation.m | 19 +- .../ECOC_2025/theory/coherence_length_plot.m | 38 +-- .../Auswertung_JLT/ber_vs_rate.m | 4 +- projects/IMDD_base_system/simulation_bwl.m | 4 +- projects/WDM/WDM_auswertung.m | 251 ++++++++++++++++++ projects/WDM/WDM_model.m | 122 ++++----- 9 files changed, 397 insertions(+), 175 deletions(-) create mode 100644 projects/WDM/WDM_auswertung.m diff --git a/Classes/04_DSP/Sequence Detection/MLSE.m b/Classes/04_DSP/Sequence Detection/MLSE.m index e51d61b..c25654c 100644 --- a/Classes/04_DSP/Sequence Detection/MLSE.m +++ b/Classes/04_DSP/Sequence Detection/MLSE.m @@ -54,9 +54,9 @@ classdef MLSE < handle function [VITERBI_ESTIMATION_SYMBOLS,LLR_maxlogmap,GMI] = process_(obj,data_in,data_ref) - debug = 0; + debug = 1; - trellis_state_mode = 2; % General: States should match the target states of the prev. EQ (EQ's job was to reduce the error between signal and the target) + trellis_state_mode = 0; % General: States should match the target states of the prev. EQ (EQ's job was to reduce the error between signal and the target) % 0 = use provided states (MUST provide the correct states); % 1 = normalize to = 1 rms; % 2 = use target symbols; @@ -65,7 +65,7 @@ classdef MLSE < handle trellis_exclusion = 0; % PAM-6 only (only if data is NOT precoded!) - scale_mode = 2; % scale_mode: + scale_mode = 0; % scale_mode: % 0 = no scaling, % 1 = RMS→scale MODEL, % 2 = MMSE/time-corr→scale MODEL, diff --git a/projects/ECOC_2025/load_signal_standalone.m b/projects/ECOC_2025/load_signal_standalone.m index f58c77f..e54816a 100644 --- a/projects/ECOC_2025/load_signal_standalone.m +++ b/projects/ECOC_2025/load_signal_standalone.m @@ -1,34 +1,20 @@ +db = DBHandler("type","mysql","dataBase",'labor'); + +fp = QueryFilter(); +% fp.where('Runs', 'run_id','EQUALS', 987); +M = 4; +fp.where('Runs', 'pam_level','EQUALS', M); +fp.where('Runs', 'symbolrate','EQUALS', 112e9); +fp.where('Runs', 'fiber_length','EQUALS', 0); +fp.where('Runs', 'is_mpi','EQUALS', 1); +fp.where('Runs', 'interference_path_length','EQUALS', 70); +% fp.where('Runs', 'loop_id','GREATER_THAN', 11); +fp.where('Runs', 'sir','EQUALS',20); -savePath = 'Z:\2025\ECOC Silas\ecoc_2025\'; -databasePath = 'C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; -database_name = 'ecoc2025_loops.db'; -db = DBHandler("type","mysql"); -% db = DBHandler("pathToDB", [databasePath, database_name],"type","sqlite"); -filterParams = db.tables; -% filterParams.Configurations = struct('run_id', run_id); -filterParams.Configurations = struct( ... - 'symbolrate', 112e9, ... %[224,336,360,390,420,448] - 'fiber_length', 0, ... - 'db_mode', '"no_db"', ... - 'interference_attenuation', 4, ... - 'interference_path_length', 300, ... - 'is_mpi', 1, ... - 'pam_level', 4, ... - 'wavelength', 1310, ... - 'precomp_amp', [], ... - 'signal_attenuation', [], ... - 'v_awg', [], ... - 'v_bias', [] ... - ); - -selectedFields = {'Runs.run_id','Runs.tx_bits_path','Runs.tx_signal_path','Runs.tx_symbols_path','Runs.rx_sync_path','Runs.rx_raw_path',... - 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','Measurements.power_rop','Configurations.v_bias',... - 'Configurations.interference_attenuation', 'Measurements.power_mpi_interference'}; - -[dataTable,sql_query] = db.queryDB(filterParams, selectedFields); +[dataTable,sql_query] = db.queryDB(fp, db.getTableFieldNames('Runs')); [~, uniqueIdx] = unique(dataTable.run_id); % Get unique run_id indices @@ -55,7 +41,7 @@ for i = 1:size(dataTable,1) Scpe_sig_raw = load([savePath, char(dataTable_.rx_raw_path(1))]); Scpe_sig_raw = Scpe_sig_raw.Scpe_sig_raw; - % Scpe_sig_raw.plot("displayname",['Scope Signal (Run ID: ',num2str(dataTable_.run_id)],"fignum",dataTable_.run_id,"clear",0); + Scpe_sig_raw.plot("displayname",['Scope Signal (Run ID: ',num2str(dataTable_.run_id)],"fignum",dataTable_.run_id,"clear",0); Scpe_sig_resampled = Scpe_sig_raw.resample("fs_in",Scpe_sig_raw.fs,"fs_out",fsym); diff --git a/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m index 297f17b..68ac774 100644 --- a/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m +++ b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m @@ -1,49 +1,57 @@ -% basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; -% database = DBHandler("pathToDB",[basePath,'silas_labor.db']); -database = DBHandler("type",'mysql','dataBase','labor'); -filterParams = database.tables; -filterParams.Runs.loop_id = 209; -% filterParams.Configurations = struct( ... -% 'symbolrate', 112e9, ... %[224,336,360,390,420,448] -% 'fiber_length', 0, ... -% 'db_mode', '"no_db"', ... -% 'interference_attenuation', [], ... -% 'interference_path_length', 1000, ... -% 'is_mpi', 1, ... -% 'pam_level', 4, ... -% 'wavelength', 1310, ... -% 'precomp_amp', [], ... -% 'signal_attenuation', [], ... -% 'v_awg', [], ... -% 'v_bias', [] ... -% ); +% dsp_options.database_type = 'mysql'; +% dsp_options.dataBase = 'labor'; +% dsp_options.storage_path = 'Z:\2025\ECOC Silas\ecoc_2025\'; +% database = DBHandler("dataBase", [dsp_options.dataBase], "type", dsp_options.database_type); +% filterParams = database.tables; +% filterParams.Runs.loop_id = 209; +% % filterParams.Configurations = struct( ... +% % 'symbolrate', 112e9, ... %[224,336,360,390,420,448] +% % 'fiber_length', 0, ... +% % 'db_mode', '"no_db"', ... +% % 'interference_attenuation', [], ... +% % 'interference_path_length', 1000, ... +% % 'is_mpi', 1, ... +% % 'pam_level', 4, ... +% % 'wavelength', 1310, ... +% % 'precomp_amp', [], ... +% % 'signal_attenuation', [], ... +% % 'v_awg', [], ... +% % 'v_bias', [] ... +% % ); +% +% % if 1 +% % % filterParams.EqualizerParameters.dc_buffer_len = 1; +% % filterParams.EqualizerParameters.ffe_buffer_len = 1; +% % filterParams.EqualizerParameters.smoothing_buffer_len = 4096; +% % filterParams.EqualizerParameters.smoothing_buffer_update = 224; +% % filterParams.EqualizerParameters.DCmu = 0; +% % end +% a = database.getTableFieldNames('Runs'); +% b = database.getTableFieldNames('Results'); +% c = database.getTableFieldNames('EqualizerParameters'); +% d = [a;b;c]; +% +% [dataTable,~] = database.queryDB(filterParams, d); +% +% selectedFields = {'Configurations.run_id' 'Runs.loop_id' 'Runs.date_of_run' 'Runs.rx_raw_path' 'Runs.bitrate' 'Runs.v_bias' 'Runs.v_awg' 'Runs.precomp_amp' 'Runs.symbolrate' 'Runs.pam_level'... +% 'Runs.db_mode' 'Runs.rop_attenuation' 'Runs.is_mpi' 'Runs.interference_attenuation' 'Runs.interference_path_length' 'Runs.signal_attenuation' ... +% 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'EqualizerParameters.dc_buffer_len' 'EqualizerParameters.ffe_buffer_len' 'EqualizerParameters.smoothing_buffer_len' 'EqualizerParameters.smoothing_buffer_update' 'EqualizerParameters.DCmu' 'Measurements.power_pd_in' ... +% 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.EVM' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'}; -% if 1 -% % filterParams.EqualizerParameters.dc_buffer_len = 1; -% filterParams.EqualizerParameters.ffe_buffer_len = 1; -% filterParams.EqualizerParameters.smoothing_buffer_len = 4096; -% filterParams.EqualizerParameters.smoothing_buffer_update = 224; -% filterParams.EqualizerParameters.DCmu = 0; -% end -a = database.getTableFieldNames('Runs'); -b = database.getTableFieldNames('Results'); -c = database.getTableFieldNames('EqualizerParameters'); -d = [a;b;c]; +db = DBHandler("type","mysql","dataBase",'labor'); -[dataTable,~] = database.queryDB(filterParams, d); - -selectedFields = {'Configurations.run_id' 'Runs.loop_id' 'Runs.date_of_run' 'Runs.rx_raw_path' 'Runs.bitrate' 'Runs.v_bias' 'Runs.v_awg' 'Runs.precomp_amp' 'Runs.symbolrate' 'Runs.pam_level'... - 'Runs.db_mode' 'Runs.rop_attenuation' 'Runs.is_mpi' 'Runs.interference_attenuation' 'Runs.interference_path_length' 'Runs.signal_attenuation' ... - 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'EqualizerParameters.dc_buffer_len' 'EqualizerParameters.ffe_buffer_len' 'EqualizerParameters.smoothing_buffer_len' 'EqualizerParameters.smoothing_buffer_update' 'EqualizerParameters.DCmu' 'Measurements.power_pd_in' ... - 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.EVM' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'}; - -[dataTable_raw,sql_query] = database.queryDB(filterParams, selectedFields); +fp = QueryFilter(); +% fp.where('Runs', 'loop_id','EQUALS', 209); +% fp.where('Runs', 'sir','EQUALS', 21); +fp.where('Runs', 'pam_level','EQUALS', 4); +fn = [db.getTableFieldNames('Runs');db.getTableFieldNames('Results')]; +[dataTable,sql_query] = db.queryDB(fp,fn); %% -dataTable_clean = dataTable_raw; +dataTable_clean = dataTable; dataTable_clean.SIR = -7 - round(dataTable_clean.power_mpi_interference); dataTable_clean.NGMI = dataTable_clean.GMI ./ log2(dataTable_clean.pam_level); dataTable_clean = cleanUpTable(dataTable_clean); diff --git a/projects/ECOC_2025/theory/analytic_mpi_evaluation.m b/projects/ECOC_2025/theory/analytic_mpi_evaluation.m index 9df4f5b..460f24f 100644 --- a/projects/ECOC_2025/theory/analytic_mpi_evaluation.m +++ b/projects/ECOC_2025/theory/analytic_mpi_evaluation.m @@ -1,13 +1,13 @@ % This script is used to evaluate Fig. 1b) in the paper "Adaptive Removal of Multipath Interference in Short Reach 112 GBd PAM-4 IM/DD Systems" %% Parameters -df = 1e6;%150e3; % Laser linewidth [Hz] +df = 1e6; % Laser linewidth [Hz] SIR_dB = 20; % Interference attenuation [dB] alpha = 10^(-SIR_dB/20); % Interference attenuation [linear] n_fiber = 1.467; % Refractive index c = physconst('lightspeed'); % [m/s] -L = linspace(0,400,40); % Interference delay [m] +L = linspace(0,250,50); % Interference delay [m] tau = n_fiber./c.*L; % Interference time (= tau) [s] tau_c = 1/(pi*df); % laser coherence time [s] @@ -22,7 +22,7 @@ N = round(Tsim*fs); % number of samples for each realization max_delay_samples = round(max(tau)*fs); % largest delay that is evaluated (based on max. Interference delay) phase_noise_std = sqrt(2*pi*df/fs); % standard dev. phase noise -num_realizations = 10; % number of parallel runs +num_realizations = 50; % number of parallel runs monte_carlo_variance = zeros(num_realizations, length(L)); parfor r = 1:num_realizations @@ -48,7 +48,9 @@ avg_of_mc_variances = mean(monte_carlo_variance, 1); std_of_mc_variances = std(monte_carlo_variance, 0, 1); %% Analytic variance -analytic_variance = 2*alpha^2 * (1 - exp(-2*pi*df.*tau)).^2; +L_ = linspace(0,250,500); % Interference delay [m] +tau_ = n_fiber./c.*L_; +analytic_variance = 2*alpha^2 * (1 - exp(-2*pi*df.*tau_)).^2; %% Plot cols = [0.3467 0.5360 0.6907 @@ -62,20 +64,21 @@ hold on; plot(L, avg_of_mc_variances, 'LineWidth',2, 'DisplayName','Simulation','Color',cols(1,:),'LineStyle','-'); errorbar(L, avg_of_mc_variances,std_of_mc_variances, 'LineWidth',0.7,'LineStyle','none', 'DisplayName','Simulation','Color',cols(1,:),'HandleVisibility','off'); -plot(L, analytic_variance, 'LineWidth',2, 'DisplayName','Analytic','Color',cols(2,:),'LineStyle','-'); +plot(L_, analytic_variance, 'LineWidth',2, 'DisplayName','Analytic','Color',cols(2,:),'LineStyle','-'); xticks(coherence_length_multiples.*L_c); -xticklabels(round(coherence_length_multiples.*L_c)); +xticklabels(round(coherence_length_multiples.*L_c,1)); norm_to_coherence_len = 1; if norm_to_coherence_len xticklabels(coherence_length_multiples); - xlabel('$\tau_c$', 'FontSize',12); + xlabel('$n \cdot L_c$', 'FontSize',12); +else + xlabel('Interference Delay [m]', 'FontSize',12); end xline(L_c.*coherence_length_multiples, 'LineWidth',1.5, 'DisplayName','Coh. Length','HandleVisibility','off','Color',[0.7,0.7,0.7],'LineStyle','-'); xlim([0,L(end)]); yline(var_sat, '-.k','LineWidth',1.5, 'DisplayName','Saturation: 2$\alpha ^2$'); -xlabel('Interference Delay [m]', 'FontSize',12); grid on; ylabel('Intensity Variance', 'FontSize',12); title(sprintf('MPI Variance; %d MHz; SIR: %d dB',df.*1e-6,SIR_dB), 'FontSize',14); diff --git a/projects/ECOC_2025/theory/coherence_length_plot.m b/projects/ECOC_2025/theory/coherence_length_plot.m index e6d040d..b793687 100644 --- a/projects/ECOC_2025/theory/coherence_length_plot.m +++ b/projects/ECOC_2025/theory/coherence_length_plot.m @@ -1,32 +1,32 @@ %% Parameters -df = linspace(100e3,50e6,10000); % Laser FWHM linewidth [Hz] +df = linspace(1,50e6,10000); % Laser FWHM linewidth [Hz] n_fiber = 1.467; % Fiber group index c = 3e8; % Speed of light [m/s] % Compute coherence length (1/e of mean-fringe decay) tau_c = 1./(pi*df); -L_c = (c/n_fiber) .* tau_c; % Coherence length [m] +L_c = (c.* tau_c/n_fiber) ; % Coherence length [m] %% Plot figure('Color','w'); loglog(df/1e6, L_c, 'LineWidth',2,'LineStyle','-'); % linewidth in MHz -xticks([0.1, 1, 10, 50]); -yticks([1, 10, 100, 1000]); -yticklabels({'1','10','100','1000'}) +% xticks([0.1, 1, 10, 50]); +% yticks([1, 10, 100, 1000]); +% yticklabels({'1','10','100','1000'}) grid on; box on; -xlabel('Laser linewidth [MHz]','FontSize',12,'Interpreter','none'); -ylabel('Coherence length [m]','FontSize',12,'Interpreter','none'); -title('Coherence Length vs. Laser Linewidth','FontSize',14,'Interpreter','none'); +xlabel('Laser linewidth [MHz]','FontSize',12,'Interpreter','latex'); +ylabel('Coherence length [m]','FontSize',12,'Interpreter','latex'); +title('Coherence Length vs. Laser Linewidth','FontSize',14,'Interpreter','latex'); %% Annotate some key points -hold on; -freqs = [150e3, 1e6, 10e6, 50e6]; % [Hz] -for f = freqs - x = f/1e6; - y = (c/n_fiber) * (1/(pi*f)); - scatter(x,y,'Marker','x','LineWidth',1,'MarkerEdgeColor','black'); - - text(x*1.1,y, sprintf('%.2f MHz', f/1e6), ... - 'FontSize',10,'HorizontalAlignment','left'); - -end +% hold on; +% freqs = [150e3, 1e6, 10e6, 50e6]; % [Hz] +% for f = freqs +% x = f/1e6; +% y = (c/n_fiber) * (1/(pi*f)); +% scatter(x,y,'Marker','x','LineWidth',1,'MarkerEdgeColor','black'); +% +% text(x*1.1,y, sprintf('%.2f MHz', f/1e6), ... +% 'FontSize',10,'HorizontalAlignment','left'); +% +% end diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/ber_vs_rate.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/ber_vs_rate.m index d39a24c..1c39b25 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/ber_vs_rate.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/ber_vs_rate.m @@ -5,7 +5,7 @@ db = DBHandler("dataBase", [dataBase], "type", database_type); fp = QueryFilter(); % fp.where('Runs', 'run_id','EQUALS', 987); -M = 8; +M = 4; fp.where('Runs', 'pam_level','EQUALS', M); % fp.where('Runs', 'bitrate','LESS_THAN', 310e9); fp.where('Runs', 'fiber_length','EQUALS', 2); @@ -13,7 +13,7 @@ fp.where('Runs', 'is_mpi','EQUALS', 0); % fp.where('Runs', 'interference_path_length','EQUALS', 1000); % fp.where('Runs', 'loop_id','GREATER_THAN', 11); % fp.where('Runs', 'sir','EQUALS',18); -fp.where('Runs', 'wavelength','EQUALS', 1293); +fp.where('Runs', 'wavelength','EQUALS', 1310); % fp.where('Runs', 'db_mode','EQUALS', 0); fp.where('Runs', 'rop_attenuation','EQUALS', 0); diff --git a/projects/IMDD_base_system/simulation_bwl.m b/projects/IMDD_base_system/simulation_bwl.m index 40ccbb2..79fc4c7 100644 --- a/projects/IMDD_base_system/simulation_bwl.m +++ b/projects/IMDD_base_system/simulation_bwl.m @@ -1,6 +1,6 @@ %%% Run parameters % TX -M = 6; +M = 4; m = floor(log2(M)*10)/10; fsym = 224e9; @@ -48,7 +48,7 @@ cols = linspecer(6); rop = [-6]; bwl = [0.5:0.1:1.5]; fsym = [120:8:256].*1e9; -fsym =150e9; +fsym =210e9; ber_vnle = []; ber_mlse = []; diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m new file mode 100644 index 0000000..32c6cee --- /dev/null +++ b/projects/WDM/WDM_auswertung.m @@ -0,0 +1,251 @@ +figure(); hold on; +cols = cbrewer2('set2',N); % one color per wavelength (Ch) + +try + rop = res.settings.rop; % 12 points + wavelengthplan = res.settings.wavelengthplan; +catch + wavelengthplan = [1295,1305,1315,1325]; + % wavelengthplan = calcWavelengthPlan(16,400e9,1310); + rop = -8.25:0.75:0; +end + +N = length(wavelengthplan); + +fec = 2.2e-4; +fec = 3.8e-3; +Sffe = cell(1,N); +Svnle = cell(1,N); +Smlse = cell(1,N); +Sdbt = cell(1,N); + +% Choose your quantile band. For your old style, use 0.04/0.99: +qLow = 0.0; % lower quantile (e.g., 0.04 for old script) +qHigh = 1; % upper quantile (e.g., 0.99 for old script) +cols = linspecer(N); % one color per wavelength (Ch) +cols = cbrewer2('set1',N); + +for l = 1:N + % Slice 12x50 cell arrays + ffe_cells = reshape(squeeze(res.ffe(l,:,:)),length(rop),[]); + vnle_cells = reshape(squeeze(res.vnle(l,:,:)),length(rop),[]); + mlse_cells = reshape(squeeze(res.mlse(l,:,:)),length(rop),[]); + dbt_cells = reshape(squeeze(res.dbt(l,:,:)),length(rop),[]); + + [Sffe{l}, noX_ffe] = fecCrossings(rop, ffe_cells, fec); + + [Svnle{l}, noX_ffe] = fecCrossings(rop, vnle_cells, fec); + + [Smlse{l}, noX_ffe] = fecCrossings(rop, mlse_cells, fec); + + [Sdbt{l}, noX_ffe] = fecCrossings(rop, dbt_cells, fec); + + % Extract BER matrices using only complete realizations (12/12 ROP filled) + ffe_mat = extractCompleteBER(ffe_cells); % 12 x K_ffe + vnle_mat = extractCompleteBER(vnle_cells); % 12 x K_vnle + mlse_mat = extractCompleteBER(mlse_cells); % 12 x K_mlse + mlse_alpha_mat = extractCompleteAlphas(mlse_cells); % 12 x K_mlse + dbt_mat = extractCompleteBER(dbt_cells); % 12 x K_dbt + + showLegend = 1; % one legend entry per technique + + % Plot shaded band + mean line with boundedline + % plotBandMeanBL(rop, ffe_mat, cols(l,:), sprintf('FFE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--s', showLegend); + % scatter(Sffe,fec.*ones(size(Sffe)),20,'v','MarkerFaceColor','black'); + + plotBandMeanBL(rop, vnle_mat, cols(l,:), sprintf('VNLE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--x', showLegend); + + % plotBandMeanBL(rop, mlse_mat, cols(l,:), sprintf('VNLE+PF+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '-o', showLegend); + + % plotBandMeanBL(rop, dbt_mat, cols(l,:), sprintf('DBt.+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--v', showLegend); + + set(gca,'XScale','linear','YScale','log','TickLabelInterpreter','latex','FontSize',11); + yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5); + +end + +ylabel('BER'); +xlabel('ROP'); +title('BER vs. ROP'); +xlim([min(rop) max(rop)]); +ylim([1e-5 0.3]); +grid on; +legend show; + + +S_cell = Sdbt; +S_cell =Smlse; +S_cell = {Svnle,Smlse,Sdbt}; +S_cell = {Svnle}; +figure(5); hold on; +for i = 1:length(S_cell) + % Pad to rectangular matrix: rows = realizations, cols = wavelengths + Kmax = max(cellfun(@numel, S_cell{i})); + S_mat = NaN(Kmax, N); + for l = 1:N + k = numel(S_cell{i}{l}); + if k > 0 + S_mat(1:k, l) = S_cell{i}{l}; + end + end + + % --- Violin plot over wavelengths (columns) --- + + cols=linspecer(3); + catLabels = arrayfun(@(nm) sprintf('%d nm', nm), wavelengthplan, 'UniformOutput', false); + vs = violinplot(S_mat, catLabels, ... + 'ViolinColor', cols(i,:), ... + 'ViolinAlpha', 0.10, ... + 'MarkerSize', 20, ... + 'ShowMedian', true, ... + 'EdgeColor', cols(i,:), ... + 'ShowWhiskers', false, ... + 'ShowData', true, ... + 'ShowBox', false, ... + 'Bandwidth', 0.05); + + ylim([floor(min(S_mat,[],'all')), ceil(max(S_mat,[],'all'))]) + ylim([-8 0]); + ylabel('ROP at FEC crossing'); + title(sprintf('RROP to cross BER %.2e', fec)); + grid on; box on; + +end + + + + + + + + + + +%% ================= helper ================= +function plotBandMeanBL(x, Y, color, techLabel, qLow, qHigh, lineSpec, showLegend) + % Y: (nPoints x nRealizations) + % Remove realizations that are entirely zero (like removeZeros behavior) + badCols = all(Y == 0, 1); + Y(:, badCols) = []; + + Y(Y==0) = 1e-8; + % Stats across realizations + mu = mean(Y, 2, 'omitnan'); % mean line + lo = quantile(Y, qLow, 2); % lower bound + hi = quantile(Y, qHigh, 2); % upper bound + + % Convert to asymmetric distances required by boundedline: + % b(:,1) = distance to lower side; b(:,2) = distance to upper side + b = [mu - lo, hi - mu]; + + % Call boundedline with alpha shading + [hl, hp] = boundedline(x(:), mu(:), b, lineSpec, 'alpha', ... + 'transparency', 0.18); + % Color styling + set(hl, 'Color', color, 'LineWidth', 1.4, 'MarkerSize', 4); + set(hp, 'FaceColor', color, 'HandleVisibility','off'); % patch hidden in legend + + % Single legend entry per technique (use first wavelength only) + if showLegend + set(hl, 'DisplayName', techLabel); + else + set(hl, 'HandleVisibility','off'); + end + + % Optional: outline the bounds if outlinebounds is available + if exist('outlinebounds','file') == 2 + ho = outlinebounds(hl, hp); + set(ho, 'linestyle', ':', 'color', color, 'linewidth', 1, ... + 'HandleVisibility','off'); + end +end + +function [S, noCrossingMask, Y_keep] = fecCrossings(rop, cells12xR, fec) +% cells12xR: 12xR cell array (one wavelength + scheme slice) +% each cell must be a struct with .metrics.BER +% rop: 12x1 numeric vector of ROP points +% fec: scalar FEC threshold (e.g., 3.8e-3) +% +% Outputs: +% S 1xK vector of crossing ROP per kept realization (NaN if none) +% noCrossingMask 1xK logical mask: true if no crossing for that realization +% Y_keep 12xK numeric BER matrix used for the crossing detection + + % 1) keep only complete realization columns + Y = extractCompleteBER(cells12xR); % -> 12 x K + if isempty(Y) + S = []; + noCrossingMask = []; + Y_keep = Y; + return; + end + + % 2) optionally drop realizations with mean BER > 0.1 + ok = mean(Y,1,'omitnan') <= 0.1; + Y = Y(:, ok); + if isempty(Y) + S = []; + noCrossingMask = []; + Y_keep = Y; + return; + end + + % 3) find crossings per realization + nR = size(Y,2); + S = nan(1,nR); + noCrossingMask = true(1,nR); + + rop = rop(:); % ensure column + for j = 1:nR + y = Y(:,j); + + % sign change from >fec to <=fec (first time it drops below FEC) + above = (y > fec); + idx = find(above(1:end-1) & ~above(2:end), 1, 'first'); + + if ~isempty(idx) + % linear interpolation between (x1,y1) and (x2,y2) + x1 = rop(idx); y1 = y(idx); + x2 = rop(idx+1); y2 = y(idx+1); + + if isfinite(y1) && isfinite(y2) && y2 ~= y1 + t = (fec - y1) / (y2 - y1); + S(j) = x1 + t*(x2 - x1); + noCrossingMask(j) = false; + end + end + end + + Y_keep = Y; +end + + + +function Y = extractCompleteBER(cellSlice) +% cellSlice: 12xR cell array; each cell should be a struct with .metrics.BER +% Keep only those realization columns where ALL 12 ROP entries are valid. + if isempty(cellSlice), Y = []; return; end + nR = size(cellSlice,2); + keep = false(1,nR); + for r = 1:nR + col = cellSlice(:,r); + keep(r) = all(cellfun(@(c) ~isempty(c) , col)); + end + if ~any(keep), Y = []; return; end + Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true); +end + +function Y = extractCompleteAlphas(cellSlice) +% cellSlice: 12xR cell array; each cell should be a struct with .metrics.BER +% Keep only those realization columns where ALL 12 ROP entries are valid. + if isempty(cellSlice), Y = []; return; end + nR = size(cellSlice,2); + keep = false(1,nR); + for r = 1:nR + col = cellSlice(:,r); + keep(r) = all(cellfun(@(c) ~isempty(c) , col)); + end + if ~any(keep), Y = []; return; end + Y = cellfun(@(c) c.metrics.Alpha, cellSlice(:,keep), 'UniformOutput', true); +end + diff --git a/projects/WDM/WDM_model.m b/projects/WDM/WDM_model.m index 9563ef6..faa74e6 100644 --- a/projects/WDM/WDM_model.m +++ b/projects/WDM/WDM_model.m @@ -2,23 +2,22 @@ % TX % --- FIRST LINE: evaluate settings located beside this script --- run(fullfile(fileparts(mfilename('fullpath')),'WDM_settings.m')); +s = struct; +s.num_realiz = 1; +% s.wavelengthplan = calcWavelengthPlan(16,400e9,1310); +s.wavelengthplan = [1295,1305,1315,1325]; +s.link_length = 2; +s.pmd = 0.0; +s.gamma = 0.00; -num_realiz = 50; -% wavelengthplan = calcWavelengthPlan(16,400e9,1310); -wavelengthplan = [1295,1305,1315,1325]; -link_length = 2; -pmd = 0.1; -gamma = 0.0023; - - -M = 4; -m = floor(log2(M)*10)/10; +s.M = 4; +m = floor(log2(s.M)*10)/10; fsym = 224e9; fdac = 2*fsym; fadc = 2*fsym; -random_key = 2; +s.random_key = 100; -% Laser / Modulator +% Laser / s.Modulator vbias_rel = 0.5; u_pi = 3.2; vbias = -vbias_rel*u_pi; @@ -48,8 +47,8 @@ apply_pulsef = 0; rcalpha = 0.05; Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); -N = numel(wavelengthplan); -f_plan = physconst('lightspeed')./(wavelengthplan.*1e-9); +N = numel(s.wavelengthplan); +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; @@ -58,38 +57,38 @@ upsample_required = f_nyq./(fdac*kover/2); upsample_pow = 2^nextpow2(upsample_required); upsample_ceil = ceil(upsample_required); -f_opt = fdac*kover*upsample_pow; -f_opt_nyq = f_opt/2; +s.f_opt = fdac*kover*upsample_pow; +s.f_opt_nyq = s.f_opt/2; signal_cell = {}; Symbols = {}; Tx_bits = {}; -rop = -8.25:0.75:0; +s.rop = -6:0.75:-0.75; -output_ffe = cell(length(wavelengthplan),length(rop),num_realiz); -output_vnle = cell(length(wavelengthplan),length(rop),num_realiz); -output_mlse = cell(length(wavelengthplan),length(rop),num_realiz); -output_dbt = cell(length(wavelengthplan),length(rop),num_realiz); +output_ffe = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz); +output_vnle = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz); +output_mlse = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz); +output_dbt = cell(length(s.wavelengthplan),length(s.rop),s.num_realiz); -for realiz = 1:num_realiz +for realiz = 1:s.num_realiz parfor l = 1:N [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource(... - "fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fsym",fsym,"M",s.M,"order",18,"useprbs",0,... "fs_out",fdac,... "applyclipping",0,"clipfactor",1.5,... "applypulseform",apply_pulsef,"pulseformer",Pform,... - "randkey",random_key+l+realiz,... + "randkey",s.random_key+l+realiz,... "db_precode",db_precode,"db_encode",db_encode,... "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process(); % Digi_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0); Lp_awg = Filter('filtdegree',3,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig); - % El_sig = M8199B("kover",kover).process(Digi_sig); + % El_sig = s.M8199B("kover",kover).process(Digi_sig); % El_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0); %%%%% Electrical Driver Amplifier %%%%%% @@ -97,10 +96,10 @@ for realiz = 1:num_realiz % El_sig = El_sig.setPower(1,"dBm"); % figure;histogram(El_sig.signal); - %%%%% MODULATE E/O CONVERSION %%%%% - Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",wavelengthplan(l),"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+l+realiz).process(El_sig); + %%%%% s.MODULATE 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+realiz).process(El_sig); - signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",30).process(Eml_out); + signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",100).process(Eml_out); end Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover,... @@ -119,12 +118,12 @@ for realiz = 1:num_realiz zdw = 1310; D_local = 0; %if ~=0, simulation uses "segmented fiber with d+,d-) randomize_D = true; - Dvec = getDispersionVector(nSegments, D_local, zdw, randomize_D, random_key+realiz); - for s = 1:nSegments + Dvec = getDispersionVector(nSegments, D_local, zdw, randomize_D, s.random_key+realiz); + for seg = 1:nSegments - Opt_sig_wdm_fib = DP_Fiber("L",link_length/nSegments,"D",Dvec(s),"Dpmd",pmd,"Ds",0.07,... + Opt_sig_wdm_fib = DP_Fiber("L",s.link_length/nSegments,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07,... "beat_len",10,"corr_len",100,"dz",1,"manakov",0,... - "gamma",gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01,... + "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).process(Opt_sig_wdm_fib); end @@ -133,23 +132,23 @@ for realiz = 1:num_realiz % Opt_sig_wdm_fib.move_it_spectrum("fignum",100212,"displayname",'bla'); - % Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig) + % Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",s.link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"s.gamma",0,"Dslope",0.07).process(Opt_sig) - parfor ri = 1:length(rop) + for ri = 1:length(s.rop) %%%%%% ROP %%%%%% - Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop(ri)+10*log10(N)).process(Opt_sig_wdm_fib); + Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)+10*log10(N)).process(Opt_sig_wdm_fib); Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",Opt_sig_wdm_rx.fs/upsample_pow,"fs_in",Opt_sig_wdm_rx.fs,"lambda_center",1310).process(Opt_sig_wdm_rx); PD_cell = {}; - for l = 1:N + parfor l = 1:N %%%%%% PD Square Law %%%%%% assert(fdac*kover==Opt_sig_wdm_demux{l}.fs,'Sampling Frequencies do not match! Check previous steps'); - PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",random_key+l+realiz).process(Opt_sig_wdm_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+realiz).process(Opt_sig_wdm_demux{l}); - PD_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1); + % PD_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1); %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% rx_bwl = 100e9; @@ -176,7 +175,7 @@ for realiz = 1:num_realiz % FFE 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,M,Rx_sig,Symbols{l},Tx_bits{l},... + ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols{l},Tx_bits{l},... "precode_mode",duob_mode,... 'showAnalysis',0,... "postFFE",[],... @@ -194,12 +193,12 @@ for realiz = 1:num_realiz useviterbi = 0; if useviterbi - mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + mlse_ = MLSE_viterbi("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels); else - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + mlse_ = MLSE("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels); end - [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Rx_sig, Symbols{l},Tx_bits{l}, ... + [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_, pf_, mlse_, s.M, Rx_sig, Symbols{l},Tx_bits{l}, ... "precode_mode", duob_mode,... 'showAnalysis', 0, ... "postFFE", [],... @@ -209,18 +208,17 @@ for realiz = 1:num_realiz output_mlse{l,ri,realiz} = mlse_results; - % DB tgt. useviterbi = 0; if useviterbi - mlse_db_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + mlse_db_ = MLSE_viterbi("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels); else - mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",s.M,"trellis_states",PAMmapper(s.M,0).levels); end ffe_order = [50, 5, 5]; eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"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",1); - dbt_results = duobinary_target(eq_, mlse_db_, M, Rx_sig, Symbols{l},Tx_bits{l}, ... + dbt_results = duobinary_target(eq_, mlse_db_, s.M, Rx_sig, Symbols{l},Tx_bits{l}, ... "precode_mode", duob_mode, ... 'showAnalysis', 0,... "postFFE", []); @@ -232,6 +230,7 @@ for realiz = 1:num_realiz end res = struct(); + res.settings = s; res.ffe = output_ffe; res.vnle = output_vnle; res.mlse = output_mlse; @@ -240,38 +239,13 @@ for realiz = 1:num_realiz % Save results save(fullfile(output_root, fname), 'res', '-v7.3'); fprintf('Saved results to: %s\n', fullfile(output_root, fname)); + disp(datetime('now','TimeZone','local','Format','yyyyMs.Mdd_HHmmss')); + end - - -figure();hold on; -cols = linspecer(N); -for l = 1:N - % plot(rop,mean(squeeze(ber_vnle(l,:,:)),2,'omitnan'),'Marker','*','DisplayName',sprintf('Ch: %d',wavelengthplan(l))) - % plot(rop,cellfun(@(c) c.metrics.BER, output_ffe(l,:), 'UniformOutput', true),'Marker','*','DisplayName',sprintf('Ch: %d',wavelengthplan(l)),'Color',cols(l,:),'HandleVisibility','on','LineStyle',':'); - plot(rop,cellfun(@(c) c.metrics.BER, res.vnle(l,:), 'UniformOutput', true),'Marker','x','DisplayName',sprintf('Ch: %d',wavelengthplan(l)),'Color',cols(l,:),'HandleVisibility','on','LineStyle','--') - plot(rop,cellfun(@(c) c.metrics.BER, res.mlse(l,:), 'UniformOutput', true),'Marker','o','DisplayName',sprintf('Ch: %d',wavelengthplan(l)),'Color',cols(l,:),'HandleVisibility','on','LineStyle','-') -end -yline([3.8e-3,2.2e-4],'HandleVisibility','off'); -ylabel('BER'); -xlabel('ROP') -title('BER vs. ROP'); -set(gca, 'XScale', 'linear', ... - 'YScale', 'log', ... - 'TickLabelInterpreter', 'latex', ... - 'FontSize', 11); -xlim([min(rop) max(rop)]) -ylim([1e-5 0.3]) - -% --- save as PNG --- -outname = fullfile(output_root, 'BER_vs_ROP.png'); % saves to current folder -print(gcf, outname, '-dpng', '-r300'); % 300 dpi -fprintf('Saved figure to %s\n', outname); - - function dispersion_vector = getDispersionVector(N, D, ref_zdw, randomize_ZDW, randomkey) -% MATLAB version of the Python generator shown above. +% s.MATLAB version of the Python generator shown above. % Returns an N×1 vector (ps/(nm·km)). % % D is the nominal dispersion magnitude. For D>0 the link is segmented with