new plots for Diss. Mostly AI gen. Few changes in actual codebase

This commit is contained in:
Silas Oettinghaus
2026-07-30 08:35:45 +02:00
parent 125d8508ca
commit 7a9deaeb0c
62 changed files with 6171 additions and 630 deletions

View File

@@ -403,7 +403,7 @@ classdef Signal
% Divide by 2*pi for the f/fs axis where Nyquist is 0.5.
end
p_lin = movmean(p_lin,10);
p_lin = movmean(p_lin,20);
if options.normalizeTo0dB
p_lin = p_lin ./ max(p_lin);
@@ -1042,34 +1042,23 @@ classdef Signal
if isempty(finite_eye)
finite_eye = sig(isfinite(sig));
end
amp_min = min(finite_eye);
amp_max = max(finite_eye);
amp_center = (amp_max + amp_min) / 2;
amp_span = amp_max - amp_min;
if amp_span == 0
amp_span = max(abs(amp_center),1);
end
amp_margin = 0.08 * amp_span;
maxA = amp_center + amp_span/2 + amp_margin;
minA = amp_center - amp_span/2 - amp_margin;
if ~isa(obj,'Opticalsignal') && minA < 0 && maxA > 0
targetStep = max(abs([minA maxA])) / 2;
if targetStep > 0
stepMagnitude = 10^floor(log10(targetStep));
normalizedStep = targetStep / stepMagnitude;
if normalizedStep <= 1
tickStep = stepMagnitude;
elseif normalizedStep <= 2
tickStep = 2 * stepMagnitude;
elseif normalizedStep <= 5
tickStep = 5 * stepMagnitude;
else
tickStep = 10 * stepMagnitude;
end
axisLimit = 2 * tickStep;
maxA = axisLimit;
minA = -axisLimit;
if isa(obj,'Opticalsignal')
amp_min = min(finite_eye);
amp_max = max(finite_eye);
amp_center = (amp_max + amp_min) / 2;
amp_span = amp_max - amp_min;
if amp_span == 0
amp_span = max(abs(amp_center),1);
end
amp_margin = 0.08 * amp_span;
maxA = amp_center + amp_span/2 + amp_margin;
minA = amp_center - amp_span/2 - amp_margin;
else
% Normalized electrical/digital eye display range.
% Samples outside this range are clipped in the
% histogram image so all eye plots use the same scale.
minA = -3;
maxA = 3;
end
% maxA = 0.12;
@@ -1109,13 +1098,13 @@ classdef Signal
max_ = abs(max(obj.signal(100:end-100)).^2);
elseif isa(obj,'Electricalsignal')
title(['Electrical Eye ',options.displayname])
ylabel("Voltage in V");
ylabel("Normalized amplitude");
yTickValues = linspace(maxA,minA,5);
min_ = min(obj.signal(100:end-100));
max_ = abs(max(obj.signal(100:end-100)));
else
title(['Digital Eye ',options.displayname])
ylabel("Digital Signal Amplitude");
ylabel("Normalized amplitude");
yTickValues = linspace(maxA,minA,5);
min_ = min(obj.signal(100:end-100));
max_ = abs(max(obj.signal(100:end-100)));

View File

@@ -212,6 +212,81 @@ classdef ChannelFreqResp < handle
end
function Target = apply(obj,Target,options)
% Apply the measured magnitude-only channel response.
arguments
obj
Target
options.fileName = ''
options.loadPath = ''
options.maxampdb double = 0
end
if isempty(obj.H)
obj.load("fileName",options.fileName,"loadPath",options.loadPath);
end
fstarget = Target.fs;
N_target = numel(Target.signal);
% Ignore measured DC attenuation when setting the response
% reference level.
H_measured = abs(obj.H(:)).';
H_measured_valid = isfinite(H_measured);
H_measured_dB = 20*log10(H_measured);
idx_old = (obj.faxis > 0) & ...
(obj.faxis < fstarget/2) & H_measured_valid;
H_reference_dB = max(H_measured_dB(idx_old));
H_measured_dB = H_measured_dB - H_reference_dB + options.maxampdb;
H_measured = 10.^(H_measured_dB/20);
H_measured(~H_measured_valid) = 0;
% Positive-frequency FFT bins, excluding DC and Nyquist.
N_positive = floor((N_target - 1)/2);
fnew = (1:N_positive)*fstarget/N_target;
f_old = obj.faxis(idx_old);
H_old = H_measured(idx_old);
H_interp = interp1(f_old,H_old,fnew,'linear');
% Match the edge handling used by the precompensation path.
nH = find(isfinite(H_interp),1,'first');
if isempty(nH)
error('ChannelFreqResp:InterpolationFailed', ...
'The measured response does not overlap the target frequency grid.');
end
H_interp(1:nH-1) = H_interp(nH);
nH = find(isfinite(H_interp),1,'last');
H_interp(nH+1:end) = H_interp(nH);
% Real, zero-phase, conjugate-symmetric FFT response.
if mod(N_target,2) == 0
H_apply = [H_interp(1), H_interp, 0, fliplr(H_interp)];
else
H_apply = [H_interp(1), H_interp, fliplr(H_interp)];
end
H_apply = H_apply(:);
% Re-anchor the final interpolated response using only the
% positive-frequency bins, so DC attenuation is ignored.
H_apply_dB = 20*log10(abs(H_apply));
positive_bins = 2:(N_positive + 1);
H_apply_dB = H_apply_dB - max(H_apply_dB(positive_bins)) + options.maxampdb;
H_apply = 10.^(H_apply_dB/20);
obj.H_apply = H_apply;
% Apply the channel in the frequency domain and preserve shape.
signal_shape = size(Target.signal);
signal_in = real(Target.signal(:));
signal_out = real(ifft(fft(signal_in).*H_apply));
Target.signal = reshape(signal_out,signal_shape);
end
function plot(obj)
figure();

View File

@@ -45,6 +45,8 @@ classdef ML_MLSE_DUOBINARY < handle
valid
valid_to_idx
valid_from_idx
incoming_edge_idx
incoming_from_idx
w
% Fast lookup
@@ -124,6 +126,21 @@ classdef ML_MLSE_DUOBINARY < handle
end
[obj.valid_to_idx,obj.valid_from_idx] = find(obj.valid);
% Precompute compact incoming-transition lookup tables. Every
% trellis state has obj.S incoming transitions, so compare-select
% can operate on valid edges only instead of an nStates-by-nStates
% matrix for every symbol.
obj.incoming_edge_idx = zeros(obj.S, obj.nStates);
obj.incoming_from_idx = zeros(obj.S, obj.nStates, 'uint32');
incoming_count = zeros(obj.nStates, 1);
for edge_idx = 1:numel(obj.valid_to_idx)
to_idx = obj.valid_to_idx(edge_idx);
slot = incoming_count(to_idx) + 1;
obj.incoming_edge_idx(slot, to_idx) = edge_idx;
obj.incoming_from_idx(slot, to_idx) = obj.valid_from_idx(edge_idx);
incoming_count(to_idx) = slot;
end
% --- Initialize weights
if isempty(obj.w) || any(size(obj.w) ~= [obj.Nf+1,obj.nFeasible])
% obj.w = randn(obj.Nf+1,obj.nFeasible);
@@ -148,6 +165,8 @@ classdef ML_MLSE_DUOBINARY < handle
% TRAINING
% ==============================================================
fprintf('\n--- Training mode ---\n');
obj.ber = nan(1, obj.epochs_tr);
obj.ce = nan(1, obj.epochs_tr);
obj.equalize(X.signal, D.signal, obj.mu_tr, obj.epochs_tr, obj.len_tr, true);
obj.e_tr = obj.e;
@@ -155,6 +174,7 @@ classdef ML_MLSE_DUOBINARY < handle
% DECISION-DIRECTED / TESTING
% ==============================================================
fprintf('--- Decision-directed / detection mode ---\n');
obj.ber_dd = nan(1, obj.epochs_dd);
[y, y_vit] = obj.equalize(X.signal, D.signal, obj.mu_dd, obj.epochs_dd, X.length, false);
X_viterbi = X;
@@ -173,9 +193,18 @@ classdef ML_MLSE_DUOBINARY < handle
for epoch = 1:epochs
pm = zeros(obj.nStates,1);
pred = zeros(nSymbols,obj.nStates,'uint32');
pm_sto = nan(obj.nStates,nSymbols,'like',pm);
needTraceback = ~training || epoch == epochs;
if needTraceback
pred = zeros(nSymbols,obj.nStates,'uint32');
end
if debug && showPlots
pm_sto = nan(obj.nStates,nSymbols,'like',pm);
end
CE_accum = 0;
if training
CE_symbol = zeros(nSymbols,1);
CE_smooth = zeros(nSymbols,1);
end
start_sample = 1;
end_sample = N;
@@ -292,13 +321,24 @@ classdef ML_MLSE_DUOBINARY < handle
% ===================================================================
% DECODING MODE (Viterbi only)
% ===================================================================
% Compare-Select (always executed)
vmat=inf(obj.nStates,obj.nStates);
vmat(obj.valid)=v_tilde;
[pm_next,pred(symbol,:)]=min(vmat,[],2);
% Compare-select over valid incoming transitions only.
incoming_metrics = v_tilde(obj.incoming_edge_idx);
[pm_next, predecessor_slot] = min(incoming_metrics,[],1);
pm_next = pm_next.';
pm_next=pm_next-min(pm_next);
pm=pm_next;
pm_sto(:,symbol)=pm;
if needTraceback
linear_idx = predecessor_slot + (0:obj.nStates-1) .* obj.S;
pred(symbol,:) = obj.incoming_from_idx(linear_idx);
end
if debug && showPlots
pm_sto(:,symbol)=pm;
end
end
if training && ~needTraceback
obj.ce(epoch)=CE_accum/symbol;
continue
end
% --- Traceback

View File

@@ -14,6 +14,21 @@ classdef clr
% Paired colormap
Paired = struct( ...
'lblue', [0.6510, 0.8078, 0.8902], ...
'dblue', [0.1216, 0.4706, 0.7059], ...
'mblue', [0.3863, 0.6392, 0.7981], ...
'lgreen', [0.6980, 0.8745, 0.5412], ...
'dgreen', [0.2000, 0.6275, 0.1725], ...
'mgreen', [0.4490, 0.7510, 0.3569], ...
'lred', [0.9843, 0.6039, 0.6000], ...
'dred', [0.8902, 0.1020, 0.1098], ...
'mred', [0.9373, 0.3530, 0.3549], ...
'lorange', [0.9922, 0.7490, 0.4353], ...
'dorange', [1.0000, 0.4980, 0.0000], ...
'morange', [0.9961, 0.6235, 0.2176], ...
'llila', [0.7922, 0.6980, 0.8392], ...
'dlila', [0.4157, 0.2392, 0.6039], ...
'mlila', [0.6040, 0.4686, 0.7216], ...
'lightblue', [0.6510, 0.8078, 0.8902], ...
'blue', [0.1216, 0.4706, 0.7059], ...
'lightgreen', [0.6980, 0.8745, 0.5412], ...

View File

@@ -104,13 +104,13 @@ function [db_results] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx
[bits_db,errors_db,ber_db,a] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
burst_db = count_error_bursts(a, 40);
cols = linspecer(8);
figure();hold on;
stem(1:40,burst_db,'LineWidth',1,'Color',cols(4,:),'Marker','_','DisplayName','w/o diff. precoder');
stem(1:40,burst_db_precoded,'LineWidth',1,'Color',cols(3,:),'Marker','.','LineStyle','-','DisplayName','w diff. precoder');
xlabel('Bit Error Burst Length')
ylabel('Occurence')
set(gca, 'yscale', 'log');
% cols = linspecer(8);
% figure();hold on;
% stem(1:40,burst_db,'LineWidth',1,'Color',cols(4,:),'Marker','_','DisplayName','w/o diff. precoder');
% stem(1:40,burst_db_precoded,'LineWidth',1,'Color',cols(3,:),'Marker','.','LineStyle','-','DisplayName','w diff. precoder');
% xlabel('Bit Error Burst Length')
% ylabel('Occurence')
% set(gca, 'yscale', 'log');
end
% M = numel(unique(tx_symbols.signal));
@@ -118,6 +118,8 @@ function [db_results] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx
[bits_db,errors_db,ber_db,errorIndice_db] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
[snr_db_target, snr_db_target_lvl] = calc_snr(db_ref_sequence.signal, eq_noise.signal);
alpha = arburg(eq_noise.signal,1);%pf_.coefficients(2);
alpha = alpha(2);
@@ -152,6 +154,8 @@ function [db_results] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx
db_results.metrics.AIR = air_mlse;
db_results.metrics.MLSE_dir = mlse_.DIR;
db_results.metrics.Alpha = alpha;
db_results.metrics.SNR = snr_db_target;
db_results.metrics.SNR_level = snr_db_target_lvl;
% Create DB results structure
db_results.config = Equalizerstruct();
@@ -177,7 +181,10 @@ function [db_results] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx
showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise after Equalization');
fprintf('DB tgt BER: %.2e \n',ber_db);
if options.decoding_mode == db_decoder.sequencedetection
fprintf('DB tgt BER: %.2e \n',ber_db); %not relevant for memoryless, wont work
end
fprintf('DB tgt BER precoded: %.2e \n',ber_db_diff_precoded);
figure(341); clf;
tx_symbols_uncoded = Duobinary().decode(db_ref_sequence);

View File

@@ -50,6 +50,11 @@ eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd);
[tx_symbols_pr,~] = pf_.process(tx_symbols, eq_noise);
% Calc SNR of Partial Response Signal
pr_noise = tx_symbols_pr-mlse_sig_sd;
[snr_partial_response, ~] = calc_snr(tx_symbols_pr.signal, pr_noise.signal);
if 0 %tx_symbols.fs > 190e9
if pf_.ncoeff == 1
if pf_.coefficients(2) < 0
@@ -83,7 +88,7 @@ mlse_sig_hd = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).qua
%% Calculate performance metrics
% VNLE metrics
[snr_vnle, snr_vnle_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal);
[snr_full_response, snr_vnle_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal);
% [gmi_vnle] = calc_air(eq_signal_sd, tx_symbols, "skip_front", 10000, "skip_end", 10000);
%calculate bitwise GMI
@@ -131,7 +136,7 @@ ffe_results.metrics.numBits = numbits.vnle;
ffe_results.metrics.numBitErr = errors.vnle;
ffe_results.metrics.BER_precoded = bers.vnle_precoded;
ffe_results.metrics.numBitErr_precoded = errors.vnle_precoded;
ffe_results.metrics.SNR = snr_vnle;
ffe_results.metrics.SNR = snr_full_response;
ffe_results.metrics.SNR_level = snr_vnle_lvl;
ffe_results.metrics.STD = std_vnle_total;
ffe_results.metrics.STD_level = std_vnle_lvl;
@@ -165,7 +170,7 @@ mlse_results.metrics.numBits = numbits.mlse;
mlse_results.metrics.numBitErr = errors.mlse;
mlse_results.metrics.BER_precoded = bers.mlse_precoded;
mlse_results.metrics.numBitErr_precoded = errors.mlse_precoded;
% mlse_results.metrics.SNR = NaN;
mlse_results.metrics.SNR = snr_partial_response;
mlse_results.metrics.GMI = gmi_mlse;
mlse_results.metrics.AIR = air_mlse;
% mlse_results.metrics.EVM = NaN;

View File

@@ -225,6 +225,7 @@ else
ml_mlse_db_results.config.comment = 'function: ML-based MLSE; duobinary encoded';
ml_mlse_db_results.recipe_config = collectRecipeConfig("ml_mlse_db", ml_mlse_db_equalizer, p, options);
output.mlmlse_db_package = ml_mlse_db_results;
ml_mlse_db_results.metrics.print;
end
if p.run_mlse_db
@@ -275,10 +276,10 @@ p = struct();
p.run_ffe = false;
p.run_vnle = false;
p.run_dfe = false;
p.run_vnle_mlse = false;
p.run_dbtgt = false;
p.run_ml_mlse = true; % non-encoded and precoded branches
p.run_ml_mlse_db = true; % db_encoded: ML-based MLSE
p.run_vnle_mlse = true;
p.run_dbtgt = true;
p.run_ml_mlse = false; % non-encoded and precoded branches
p.run_ml_mlse_db = false; % db_encoded: ML-based MLSE
p.run_mlse_db = true; % db_encoded: conventional MLSE
@@ -321,7 +322,7 @@ p.decoding_mode = [];
p.ml_mlse_mu_tr = 0.03;
p.ml_mlse_mu_dd = 0.03;
p.ml_mlse_epochs_tr = 100;
p.ml_mlse_epochs_tr = 150;
p.ml_mlse_epochs_dd = 1;
p.ml_mlse_len_tr = [];
p.ml_mlse_order = 11;
@@ -376,20 +377,22 @@ end
function dbtgt_results = runDuobinaryTarget(eq_dbtgt, mlse_db, ...
Scpe_sig, Symbols, Tx_bits, options, p)
if isempty(p.decoding_mode)
dbtgt_results = duobinary_target(eq_dbtgt, mlse_db, options.M, ...
Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
"showAnalysis", options.debug_plots, ...
"postFFE", []);
else
dbtgt_results = duobinary_target(eq_dbtgt, mlse_db, options.M, ...
Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
"showAnalysis", options.debug_plots, ...
"postFFE", [], ...
"decoding_mode", p.decoding_mode);
end
if isempty(p.decoding_mode)
dbtgt_results = duobinary_target(eq_dbtgt, mlse_db, options.M, ...
Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
"showAnalysis", options.debug_plots, ...
"postFFE", []);
else
dbtgt_results = duobinary_target(eq_dbtgt, mlse_db, options.M, ...
Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", options.duob_mode, ...
"showAnalysis", options.debug_plots, ...
"postFFE", [], ...
"decoding_mode", p.decoding_mode);
end
end
function mlse = buildMlse(M, duobMode, p, mode)

View File

@@ -30,7 +30,7 @@ end
% Ensure the figure is ready before calling spectrum
eq_noise = eq_noise - mean(eq_noise.signal);
eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 0,"color",options.color,"fft_length",4096);
eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 0,"color",options.color,"fft_length",4096,"normalizeToDC",0);
if ~isnan(options.postfilter_taps)
% Hold on to the figure for further plotting

View File

@@ -5,7 +5,8 @@ arguments
options.fignum (1,1) double = NaN % Default to NaN if not provided
options.displayname (1,:) char = '' % Default to an empty string if not provided
options.ref_symbol_uncoded = []
options.nbins (1,1) double = 1000
options.nbins (1,1) double = 200
options.legendLabels = []
end
if isa(eq_signal,'Signal')
@@ -21,6 +22,7 @@ end
eq_signal = eq_signal(:);
ref_symbols = ref_symbols(:);
ref_symbol_uncoded = options.ref_symbol_uncoded(:);
legend_labels = string(options.legendLabels);
assert(numel(eq_signal) == numel(ref_symbols), ...
'showLevelHistogram:LengthMismatch', ...
@@ -66,11 +68,16 @@ end
for lvl = 1:numel(constellation)
intermediate = received_sd(lvl,:);
cnt(lvl) = round(numel(intermediate(~isnan(intermediate)))./numel(eq_signal),3).*100;
if isempty(legend_labels)
display_name = ['Lvl ',num2str(lvl),' ; ',num2str(cnt(lvl)),' '];
else
display_name = char(legend_labels(lvl));
end
hold on
warning off
histogram(received_sd(lvl,:),options.nbins, ...
"EdgeAlpha",0, ...
"DisplayName",['Lvl ',num2str(lvl),' ; ',num2str(cnt(lvl)),' '], ...
"DisplayName",display_name, ...
"FaceColor",lvlcol(lvl,:), ...
"Normalization","pdf");
warning on
@@ -90,7 +97,7 @@ end
rm_idx = mid + (-1:2); % two before mid and two after (2x2 removal centered)
rm_idx = max(1,min(size(lvlcol,1),rm_idx));
lvlcol(rm_idx,:) = [];
lvlcol = lvlcol(1:numel(sir_group_labels),:);
% lvlcol = lvlcol(1:numel(sir_group_labels),:);
for db_lvl = 1:numel(db_constellation)
db_mask = ref_symbols == db_constellation(db_lvl);
@@ -99,12 +106,20 @@ end
cnt = round(nnz(ref_symbol_uncoded == mapped_class)./numel(eq_signal),3).*100;
if db_lvl == findFirstMappedDbLevel(ref_symbols,ref_symbol_uncoded,db_constellation,mapped_class)
display_name = ['p(y|x_',num2str(class_idx),') ; ',num2str(cnt),' '];
display_name = ['p(y|x_',num2str(class_idx),')'];
if isempty(legend_labels)
display_name = ['p(y|x_',num2str(class_idx),') ; ',num2str(cnt),' '];
display_name = ['p(y|x_',num2str(class_idx),')'];
else
display_name = char(legend_labels(class_idx));
end
handle_visibility = "on";
else
display_name = ['p(y|x_',num2str(class_idx),') ; ',num2str(cnt),' '];
display_name = ['p(y|x_',num2str(class_idx),')'];
if isempty(legend_labels)
display_name = ['p(y|x_',num2str(class_idx),') ; ',num2str(cnt),' '];
display_name = ['p(y|x_',num2str(class_idx),')'];
else
display_name = char(legend_labels(class_idx));
end
handle_visibility = "off";
end
@@ -123,7 +138,7 @@ end
end
end
xlim([-3 3]);
legend
legend("Interpreter", "latex");
grid on
% view([90 -90]);

View File

@@ -3,7 +3,11 @@ function [Bits, Symbols, Scpe_cell, found_sync] = loadAndSyncRunSignals(dataTabl
%
% Inputs:
% dataTable - one-row table with run metadata and signal file paths
% options - struct with storage_path, start_occurence and max_occurences
% options - scalar struct with the following fields:
% storage_path - root directory containing the signal files.
% start_occurence - first synchronized occurrence to return (default: 1).
% max_occurences - maximum number of synchronized occurrences to return
% (default: 1).
%
% Outputs:
% Bits - transmitted bit reference
@@ -11,6 +15,14 @@ function [Bits, Symbols, Scpe_cell, found_sync] = loadAndSyncRunSignals(dataTabl
% Scpe_cell - synchronized received signal occurrences
% found_sync - true when a valid synchronization was found
arguments
dataTable
options
% options.storage_path
% options.start_occurence = 1
% options.max_occurences = 1
end
found_sync = 0;
tempLocalStorage = 1;
Scpe_cell = {};
@@ -130,7 +142,7 @@ end
function Scpe_cell = selectSyncedOccurrences(Scpe_cell, options)
available_occurences = length(Scpe_cell);
start_occurence = floor(getOption(options, 'start_occurence', 1));
max_occurences = floor(getOption(options, 'max_occurences', available_occurences));
max_occurences = floor(getOption(options, 'max_occurences', 1));
if available_occurences < 1
warning('loadAndSyncRunSignals:NoSyncedOccurrences', ...

View File

@@ -29,15 +29,21 @@ function [snr_all, snr_per_level] = calc_snr(tx_signal, eq_noise)
% Get the unique amplitude levels in the transmitted signal
levels = unique(tx_signal);
snr_per_level = zeros(numel(levels),1);
% Preallocate an array to store the SNR for each unique level
% Loop over each unique level to compute the SNR for that level
for i = 1:length(levels)
% Find indices where tx_signal equals the current level
idx = (tx_signal == levels(i));
try
% Preallocate an array to store the SNR for each unique level
% Loop over each unique level to compute the SNR for that level
for i = 1:length(levels)
% Find indices where tx_signal equals the current level
idx = (tx_signal == levels(i));
% Compute the SNR for these indices
snr_per_level(i) = snr(tx_signal(idx), eq_noise(idx));
% Compute the SNR for these indices
snr_per_level(i) = snr(tx_signal(idx), eq_noise(idx));
end
catch ME
% Handle any errors that occur during SNR calculation
% warning('Error calculating SNR for level %d: %s', levels(i), ME.message);
end
end

View File

@@ -0,0 +1,24 @@
function channelplan_nm = calcWavelengthPlan(N, df_hz, center_nm)
vec = [N:-1:1]-(N/2+0.5);
channelplan_hz = nm2hz(center_nm) + (vec * df_hz) ;
channelplan_nm = hz2nm(channelplan_hz);
% dfft = 1.367053998632946e+07; %hz
% a = diff(channelplan_hz)./2./dfft;
end
function hz = nm2hz(nm)
wavelen_in_m = nm.* 1e-9;
hz = (299792458 ./ wavelen_in_m); %frequency in Terahertz
end
function nm = hz2nm(hz)
m = (299792458 ./ hz); %wavelength in meter
nm = m .* 10^9;
end

View File

@@ -0,0 +1,4 @@
function dw = df2dw(df,f)
%delta frequency to delta wavelength
dw = 299792458 .* df ./ f.^2;
end

View File

@@ -0,0 +1,4 @@
function df = dw2df(dw,w)
%delta wavelength to delta frequency
df = 299792458 .* dw ./ w.^2;
end

View File

@@ -0,0 +1,7 @@
function positions = getSweepWavelengths(N, df_hz, StartNm)
positions = zeros(1, N);
for i = 1:N
position = (299792458 / (299792458 / (StartNm * 1e-9) + (i - 1) * df_hz)) * 1e9;
positions(i) = position;
end
end

View File

@@ -0,0 +1,4 @@
function nm = hz2nm(hz)
m = (299792458 ./ hz); %wavelength in meter
nm = m .* 10^9;
end

View File

@@ -0,0 +1,4 @@
function hz = nm2hz(nm)
wavelen_in_m = nm.* 1e-9;
hz = (299792458 ./ wavelen_in_m); %frequency in Hertz
end

View File

@@ -0,0 +1,4 @@
function thz = nm2thz(nm)
wavelen_in_m = nm.* 1e-9;
thz = (299792458 ./ wavelen_in_m) .*1e-12; %frequency in Terahertz
end

View File

@@ -0,0 +1,4 @@
function nm = thz2nm(thz)
freq_in_hz = thz.* 1e12;
nm = (299792458 ./ freq_in_hz) .*1e9; %wavelength in nanometer
end

View File

@@ -16,6 +16,10 @@ matlab2tikz(char(filename), ...
'height', '\fheight', ...
'showInfo', false, ...
'extraAxisOptions', { ...
'clip=true',...
'clip marker paths=true',...
'grid style={dashed, draw=gray80, line width=0.25pt}',...
'minor grid style={densely dotted, draw=gray80, line width=0.25pt}',...
'xlabel style={font=\color{white!15!black}\small}', ...
'ylabel style={font=\color{white!15!black}\small}', ...
'scaled ticks=false', ...

View File

@@ -0,0 +1,30 @@
function outputFile = saveDirectoryInfo(directoryPath)
% Save metadata for all files and subdirectories under directoryPath.
directoryPath = string(directoryPath);
if ~isfolder(directoryPath)
error("Directory does not exist: %s", directoryPath);
end
% Recursively list all files and subdirectories
info = dir(fullfile(directoryPath, "**", "*"));
% Remove "." and ".." entries
info = info(~ismember({info.name}, {'.', '..'}));
% Create an output filename in MATLAB's current folder
[~, folderName] = fileparts(char(directoryPath));
if isempty(folderName)
folderName = 'directory';
end
timestamp = datestr(now, 'yyyymmdd_HHMMSS');
outputFile = fullfile(pwd, ...
sprintf('directory_info_%s_%s.mat', folderName, timestamp));
% Save the path and directory information
save(outputFile, 'directoryPath', 'info');
fprintf('Saved information to:\n%s\n', outputFile);
end

View File

@@ -46,6 +46,27 @@ classdef ML_MLSE_DUOBINARY_test < IMDDTestCase
txSymbols.signal(11:end-10), "AbsTol", 1e-12);
end
function intermediateTrainingEpochsKeepWeightsWithoutTracebackBer(testCase)
[~, ~, txPrecoded, txEncoded] = makePam4Fixture(512);
eq = ML_MLSE_DUOBINARY( ...
"sps", 1, ...
"order", 1, ...
"len_tr", length(txEncoded), ...
"epochs_tr", 2, ...
"epochs_dd", 1, ...
"mu_tr", 0.01, ...
"mu_dd", 0.01, ...
"adaptive_mu", false, ...
"L", 1);
eq.process(txEncoded, txPrecoded);
testCase.verifyTrue(isnan(eq.ber(1)));
testCase.verifyTrue(isfinite(eq.ber(2)));
testCase.verifyTrue(all(isfinite(eq.ce)));
end
function resultWrapperReportsBothPrecodedAndOriginalBer(testCase)
[txBits, ~, txPrecoded, txEncoded] = makePam4Fixture(32000);

Binary file not shown.

View File

@@ -5,25 +5,32 @@ fdac = 256e9;%fsym;
fadc = 256e9;
fsym = ([32:16:240].*1e9);
%fsym = 160e9;
% 1) PRBS Generation
O = 18; %order of prbs
N = 2^(O-1); %length of prbs
[~,seed] = prbs(O,1); %initialize first seed of prbs
bitpattern=[];
% %fsym = 160e9;
% % 1) PRBS Generation
% O = 18; %order of prbs
% N = 2^(O-1); %length of prbs
% [~,seed] = prbs(O,1); %initialize first seed of prbs
% bitpattern=[];
%
% for i = 1:log2(M)
% [bitpattern(:,i),seed] = prbs(O,N,seed);
% end
%
%
% if M == 6
% bitpattern = reshape(bitpattern,[],1);
% bitpattern = bitpattern(1:end-mod(length(bitpattern),5));
% end
for i = 1:log2(M)
[bitpattern(:,i),seed] = prbs(O,N,seed);
end
bits = Signalgenerator( ...
"form", signalform.prms, ...
"M", M, ...
"order", 17).process();
% symbols = PAMmapper(4, 0).map(bits);
if M == 6
bitpattern = reshape(bitpattern,[],1);
bitpattern = bitpattern(1:end-mod(length(bitpattern),5));
end
% 2 ) Build Inf. signal class
bits = Informationsignal(bitpattern);
% bits = Informationsignal(bitpattern);
% 5) AWG (lowpass, quantization, sample and hold)
kover = 8;
@@ -131,7 +138,7 @@ ylabel("Vpp in V")
legend
ylim([0.3 1.4])
%%
figure(7)
hold on
scatter(fsym.*1e-9,powerlist1,'DisplayName','M8199B','MarkerFaceColor',cols(1,:),'MarkerEdgeColor',cols(1,:),'LineWidth',2);

View File

@@ -1,5 +1,5 @@
dsp_options.storage_path = 'Z:\2024\sioe_labor\';
dsp_options.max_occurences = 1;
dsp_options.storage_path = 'W:\labdata\sioe_labor\';
dsp_options.max_occurences = 10;
database = DBHandler("dataBase", 'labor_highspeed', "type", 'mysql' );
@@ -7,34 +7,32 @@ cols = cbrewer2('BuPu',25);
cols = [cols(end-10:2:end,:)];
cols = cbrewer2('Set1',6);
fignum = 200;
fig=figure(fignum);clf;
dbmode = 0;
% 1 - PAM 4 with preemphasis
fp = QueryFilter();
M = 8;
rate = [360e9];
M = 4;
% rate = [300e9];
fp.where('Runs', 'pam_level','EQUALS', M);
fp.where('Runs', 'bitrate','EQUALS', rate);%360,390
% fp.where('Runs', 'symbolrate','EQUALS',165e9);
fp.where('Runs', 'fiber_length','EQUALS', 2);
fp.where('Runs', 'wavelength','EQUALS', 1310);
fp.where('Runs', 'db_mode','EQUALS', dbmode);
fp.where('Runs', 'rop_attenuation','EQUAL', 0);
[dataTable,~] = database.queryDB(fp, database.getTableFieldNames('Runs'));
dataTable = queryRunid(dataTable.run_id, database);
dataTable = dataTable(1,:);
% dataTable = queryRunid(dataTable.run_id, database);
fsym = dataTable.symbolrate;
M = double(dataTable.pam_level);
duob_mode = db_mode(strrep(dataTable.db_mode,'"',''));
% Load and Sync signal data from DB
[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable, dsp_options);
[Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable,dsp_options);
Scpe_sig_syncd = Scpe_cell{1};
Scpe_sig_syncd.eye(fsym,M,"fignum",rate.*1e-9*M+1,"displayname",' Eye of Signal');
% Scpe_sig_syncd.eye(fsym,M,"fignum",rate.*1e-9*M+1,"displayname",' Eye of Signal');
%%%%%% SNR CHEAT - Avges the measured signal occurences found after correlation in "tsynch" %%%%%%
average_signals = 1;
if average_signals
@@ -46,16 +44,19 @@ if average_signals
scope_mean = scope_mean ./ n;
Scpe_sig_avg.signal = scope_mean;
Scpe_sig_avg.spectrum("displayname","Scope PSD","fignum",20,"normalizeTo0dB",1);
Scpe_sig_avg.plot("displayname","Scope raw signal","fignum",27,"clear",1);
% Scpe_sig_avg.spectrum("displayname","Scope PSD","fignum",20,"normalizeTo0dB",1);
% Scpe_sig_avg.plot("displayname","Scope raw signal","fignum",27,"clear",1);
Scpe_sig_avg = Scpe_sig_avg.*1.25;
Scpe_sig_avg.eye(fsym,M,"fignum",rate.*1e-9*M,"displayname",' Eye of AVG Signal');
% Scpe_sig_avg.eye(fsym,M,"fignum",rate.*1e-9*M,"displayname",' Eye of AVG Signal');
Scpe_sig = Scpe_sig_avg.normalize("mode","rms");
end
% Preprocess signal
Scpe_sig = preprocessSignal(Scpe_sig_avg, Symbols, fsym);
% Scpe_sig = preprocessSignal(Scpe_sig_avg, Symbols, fsym);
Scpe_sig.eye(fsym,M,"fignum",M*10);
%% Eye of Preprocess signal
% Scpe_sig.eye(fsym,M,"fignum",M*10);
@@ -71,8 +72,67 @@ Scpe_sig.eye(fsym,M,"fignum",M*10);
% 'legend style={font=\footnotesize}', ...
% 'legend columns=1' ...
% } );
%% Simply FFE
%%
len_tr = 4096*2;
% combine and minimize repeated params
mu_ffe = [0.0001,0.0008,0.001];
mu_dd = 0.05;
mu_dc = 0.005;
% single ffe_order used
ffe_order = [50];
% map into p for FFE init (compact)
p.epochs_tr = 5;
p.epochs_dd = 5;
p.len_tr = 4096*2;
p.ffe_mu_tr = mu_ffe(1);
p.ffe_mu_dd = mu_dd;
p.ffe_order = ffe_order;
p.eq_sps = 2;
p.optimize_mus = true;
p.dd_mode = true;
p.ffe_adaption = 1;
p.mu_dc = mu_dc;
% Initialize FFE with mapped parameters
eq_ffe = FFE("epochs_tr", p.epochs_tr, ...
"epochs_dd", p.epochs_dd, ...
"len_tr", p.len_tr, ...
"mu_dd", p.ffe_mu_dd, ...
"mu_tr", p.ffe_mu_tr, ...
"order", p.ffe_order(1), ...
"sps", p.eq_sps, ...
"decide", false, ...
"optmize_mus", p.optimize_mus, ...
"dd_mode", p.dd_mode, ...
"adaption_technique", p.ffe_adaption, ...
"dc_tracking_mu", p.mu_dc);
mu_ffe = [0.0001, 0.0008, 0.001];
mu_dfe = 0.0004;
eq_ffe = EQ("Ne", [50,0,0], ...
"Nb", [0,0,0], ...
"training_length", p.len_tr, ...
"training_loops", p.epochs_tr, ...
"dd_loops", p.epochs_dd, ...
"K", 2, ...
"DCmu", p.mu_dc, ...
"DDmu", [mu_ffe mu_dfe], ...
"DFEmu", 0.005, ...
"FFEmu", 0, ...
"plotfinal", 0, ...
"ideal_dfe", false);
[ffe_results, equalized_signal] = ffe(eq_ffe, M, Scpe_sig, Symbols, Tx_bits, ...
"precode_mode", dbmode, ...
"showAnalysis", 0, ...
"postFFE", [], ...
"eth_style_symbol_mapping", 0);
equalized_signal.eye(fsym,M,"fignum",M*2);
%% Duobinary
if duob_mode == db_mode.no_db && M == 6 %only for PAM-6 and no duobinary precoding, otherwise leads to false sequence estimation
trellexlusion = 1;
@@ -102,178 +162,4 @@ dbt_results = duobinary_target(eq_, mlse_db_, M, Scpe_sig, Symbols, Tx_bits, ...
'showAnalysis', 1,...
"postFFE", []);
%% === FINAL FIGURE SIZE ===
% Existing figure numbers
figEye = 249;
figConst = 341;
% Find axes in the source figures
srcAxEye = findobj(figEye, 'Type', 'axes');
srcAxConst = findobj(figConst, 'Type', 'axes');
% Create new combined figure
figCombined = figure;
t = tiledlayout(figCombined, 1, 2);
t.TileSpacing = 'compact';
t.Padding = 'compact';
% ------------------------------------------------------------
% LEFT TILE: EYE DIAGRAM
% ------------------------------------------------------------
ax1 = nexttile(t, 1);
hold(ax1, 'on')
% Copy children (images, lines, patches, hist objects, etc.)
copyobj(srcAxEye.Children, ax1);
% Copy labels and title
ax1.XLabel.String = srcAxEye.XLabel.String;
ax1.YLabel.String = srcAxEye.YLabel.String;
ax1.Title.String = srcAxEye.Title.String;
% Copy axis limits
ax1.XLim = srcAxEye.XLim;
ax1.YLim = srcAxEye.YLim;
ax1.YDir = srcAxEye.YDir;
% Copy ticks + labels EXACTLY (including remapped/scaled ones)
ax1.XTick = srcAxEye.XTick;
ax1.XTickLabel = srcAxEye.XTickLabel;
ax1.YTick = srcAxEye.YTick;
ax1.YTickLabel = srcAxEye.YTickLabel;
% Copy colormap + clim (important for density eye)
colormap(ax1, colormap(srcAxEye.Parent));
ax1.CLim = srcAxEye.CLim;
% Copy any style props that matter
ax1.TickDir = srcAxEye.TickDir;
ax1.TickLength = srcAxEye.TickLength;
ax1.FontSize = srcAxEye.FontSize;
ax1.Box = srcAxEye.Box;
grid(ax1,'on');
% ------------------------------------------------------------
% RIGHT TILE: CONSTELLATION HISTOGRAM
% ------------------------------------------------------------
ax2 = nexttile(t, 2);
hold(ax2, 'on')
copyobj(srcAxConst.Children, ax2);
% Copy labels and title
ax2.XLabel.String = srcAxConst.XLabel.String;
ax2.YLabel.String = srcAxConst.YLabel.String;
ax2.Title.String = srcAxConst.Title.String;
% The histogram uses the same y-axis as the eye
% Extract mapping from eye
rawTicks = ax1.YTick;
rawLabelsCell = ax1.YTickLabel;
trueVoltages = str2double(rawLabelsCell);
% Apply true voltages to the histogram axis
ax2.XTick = flip(trueVoltages);
ax2.XTickLabel = flip(rawLabelsCell);
% Set histogram y-limits to match the actual voltages
ax2.XLim = [min(trueVoltages) max(trueVoltages)];
% Ensure eye diagram prints the same (we *do not* touch ax1.YLim)
ax1.XTickLabel = rawLabelsCell;
% Copy colormap (your histogram uses same palette)
colormap(ax2, colormap(srcAxConst.Parent));
% Style properties
ax2.TickDir = srcAxConst.TickDir;
ax2.TickLength = srcAxConst.TickLength;
ax2.FontSize = srcAxConst.FontSize;
ax2.Box = srcAxConst.Box;
grid(ax2,'on');
% ============================================================
% remove right y-axis completely
% ============================================================
ax2.XAxis.Visible = 'off'; % hides ticks + labels + axis line
% BUT we still keep the YTick positions internally for alignment:
% ax2.YTick = <values already set earlier> ;
% ============================================================
% minimize distance between the two plots
% ============================================================
t.TileSpacing = 'none'; % no space between tiles
t.Padding = 'none'; % no outer padding
% Also reduce internal padding for each axis
ax1.Position(3) = ax1.Position(3) + 0.02; % widen eye a bit
ax2.Position(1) = ax2.Position(1) - 0.02; % pull histogram closer
% Keep left axis grid visible
ax2.YGrid = 'off';
%
% =====================================================================
% FINAL POLISHING: unified visual style
% =======================================================================
% --- unified font size ---
FS = 12;
set([ax1 ax2], 'FontSize', FS);
% --- unified axis line width (outline stroke thickness) ---
LW = 1.0;
set([ax1 ax2], 'LineWidth', LW);
% --- unified tick length ---
TL = [.015 .015];
set([ax1 ax2], 'TickLength', TL);
% --- unified grid style ---
set([ax1 ax2], 'XGrid', 'on', 'YGrid', 'on');
set([ax1 ax2], 'GridLineStyle', '--');
set([ax1 ax2], 'GridAlpha', 0.2);
% --- remove right y-axis ticks and labels ---
ax2.YAxis.Visible = 'off';
% --- copy colormap + CLim from the eye to histogram (synchronize look) ---
colormap(ax1, colormap(srcAxEye.Parent));
colormap(ax2, colormap(srcAxEye.Parent));
ax2.CLim = ax1.CLim;
% --- minimal spacing between tiles ---
t.TileSpacing = 'none';
t.Padding = 'none';
% --- pull the panels together (touching boundary effect) ---
pos1 = ax1.Position;
pos2 = ax2.Position;
% Shift histogram left until the outlines touch
pos2(1) = pos1(1) + pos1(3) - 0.002; % 0.002 = fine overlap control
ax2.Position = pos2;
% Expand histogram slightly, remove white band
pos2 = ax2.Position;
pos2(3) = pos2(3) + 0.01;
ax2.Position = pos2;
% Ensure the left plot stays correct after the move
ax1.Position = pos1;
% --- enforce same visible outline ---
% For ax2, create a fake left spine (since YAxis is hidden)
ax2.Box = 'on'; % keep outline but no ticks on the right
ax1.Box = 'on';
ax2.View = [90 -90];

View File

@@ -1,4 +1,4 @@
dsp_options.storage_path = 'Z:\2024\sioe_labor\';
dsp_options.storage_path = 'W:\labdata\sioe_labor\';
dsp_options.max_occurences = 1;
% database = DBHandler("dataBase", 'labor_highspeed', "type", 'mysql' );
dsp_options.database_type = "mysql";
@@ -18,51 +18,64 @@ cols = cbrewer2('BuPu',25);
cols = [cols(end-10:2:end,:)];
cols = cbrewer2('Set1',6);
fignum = 200;
fig=figure(fignum);clf;
% fignum = 200;
% fig=figure(fignum);clf;
for dbmode = 0%length(rates)
if 1
rcalpha = 0.05;
fsym = rates/2;
pulsef = 1;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha);
rcalpha = 0.05;
fsym = rates/2;
pulsef = 1;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha);
Pamsource = PAMsource(...
"fsym",fsym,"M",4,"order",18,"useprbs",0,...
"fs_out",fdac,...
"applyclipping",0,"clipfactor",1.2,...
"applypulseform",pulsef,"pulseformer",Pform,...
"randkey",20,...
"db_precode",dbmode,"db_encode",0,...
"mrds_code",0,"mrds_blocklength",512);
Pamsource = PAMsource(...
"fsym",fsym,"M",4,"order",18,"useprbs",0,...
"fs_out",256e9,...
"applyclipping",0,"clipfactor",1.2,...
"applypulseform",pulsef,"pulseformer",Pform,...
"randkey",20,...
"duobinary_mode",dbmode,...
"mrds_code",0,"mrds_blocklength",512);
[Digi_sig,Symbols,Bits] = Pamsource.process();
[Digi_sig,Symbols,Bits] = Pamsource.process();
Digi_sig = Digi_sig.normalize("mode","rms");
Digi_sig = Digi_sig.normalize("mode","rms");
%%% 1) PLOT FULL RESPONSE SIGNAL
Digi_sig.spectrum("displayname","Full Response","fignum",fignum+dbmode,"normalizeToNyquist",0,"normalizeTo0dB",0,"color",[0.2,0.2,0.2],"linestyle",'-','addDCoffset',0,'normalizeToDC',1);
%%% 1) PLOT FULL RESPONSE SIGNAL
Digi_sig.spectrum("displayname","Full Response","fignum",fignum+dbmode,"normalizeToNyquist",0,"normalizeTo0dB",0,"color",clr.Paired.dgreen,"linestyle",'-','addDCoffset',0,'normalizeToDC',1);
%%% 2) PLOT PREEMPH. TX SIGNAL
if dbmode == 0
maxamp = -37;
precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs);
%%% 2) PLOT PREEMPH. TX SIGNAL
if dbmode == 0
maxamp = -37;
precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs);
precomp_path = "W:\labdata\sioe_labor\precomp";
precomp_fn = "lab_high_speed";
Digi_sig_pre = precomp_est.precomp(Digi_sig,'maxampdb',maxamp,'loadPath',precomp_path,'fileName',precomp_fn);
precomp_path = "W:\labdata\sioe_labor\precomp";
precomp_fn = "lab_high_speed";
Digi_sig_pre = precomp_est.precomp(Digi_sig,'maxampdb',maxamp,'loadPath',precomp_path,'fileName',precomp_fn);
Digi_sig_pre = Digi_sig_pre.resample("fs_out",fdac);
Digi_sig_pre = Digi_sig_pre.resample("fs_out",256e9);
Digi_sig_pre= Digi_sig_pre.normalize("mode","rms");
Digi_sig_pre= Digi_sig_pre.normalize("mode","rms");
Digi_sig_pre.spectrum("displayname","Strong Precomp","fignum",fignum+dbmode,"normalizeToNyquist",0,"normalizeTo0dB",0,"color",[0,0,0],"linestyle",'-.','addDCoffset',0,'normalizeToDC',1);
end
Digi_sig_pre.spectrum("displayname","Tx w/ pre-emph.","fignum",fignum+dbmode,"normalizeToNyquist",0,"normalizeTo0dB",0,"color",clr.Paired.dred,"linestyle",'-.','addDCoffset',0,'normalizeToDC',1);
elseif dbmode == 2
maxamp = -38;
precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs);
precomp_path = "W:\labdata\sioe_labor\precomp";
precomp_fn = "lab_high_speed";
Digi_sig_pre = precomp_est.precomp(Digi_sig,'maxampdb',maxamp,'loadPath',precomp_path,'fileName',precomp_fn);
Digi_sig_pre = Digi_sig_pre.resample("fs_out",256e9);
Digi_sig_pre= Digi_sig_pre.normalize("mode","rms");
Digi_sig_pre.spectrum("displayname","Strong Precomp","fignum",fignum+dbmode,"normalizeToNyquist",0,"normalizeTo0dB",0,"color",[0,0,0],"linestyle",'-.','addDCoffset',0,'normalizeToDC',1);
end
end
% 1 - PAM 4 with preemphasis
@@ -70,7 +83,7 @@ for dbmode = 0%length(rates)
M = 4;
fp.where('Runs', 'pam_level','EQUALS', M);
fp.where('Runs', 'bitrate','EQUALS', rates);%360,390
fp.where('Runs', 'fiber_length','EQUALS', 2);
fp.where('Runs', 'fiber_length','EQUALS', 10);
fp.where('Runs', 'wavelength','EQUALS', 1310); %1327.4
fp.where('Runs', 'db_mode','EQUALS', dbmode);
fp.where('Runs', 'rop_attenuation','EQUAL', 0);
@@ -89,14 +102,18 @@ for dbmode = 0%length(rates)
Scpe_sig = Scpe_cell{1};
%%% 3) PLOT DB Tgt. SIGNAL
if 1
if dbmode ~= 2
DB_Symbols = Duobinary().encode(Symbols);
DB_Symbols.spectrum("fignum",fignum+dbmode,"normalizeTo0dB",1,"displayname",'DB-Response','addDCoffset',0,'color',clr.Set1.blue,'normalizeToNyquist',0,'linestyle','--');
else
DB_Symbols = Symbols;
end
DB_Symbols.spectrum("fignum",fignum+dbmode,"normalizeTo0dB",1,"displayname",'DB response','addDCoffset',0,'color',clr.Paired.dblue,'normalizeToNyquist',0,'linestyle','-');
%%% 4) Plot RX Signal
Scpe_sig = Scpe_sig - mean(Scpe_sig.signal);
Scpe_sig.spectrum("fignum",fignum+dbmode,"normalizeTo0dB",1,"displayname",'Rx','addDCoffset',1,'color',[0,0,0],'normalizeToNyquist',0,'linestyle',':');
Scpe_sig.eye(fsym,M,"fignum",47,"displayname",' Eye of AVG Signal');
% Scpe_sig.eye(fsym,M,"fignum",47,"displayname",' Eye of AVG Signal');
% xline(Symbols.fs/2.*1e-9,'Color',cols(r,:),'HandleVisibility','off');
average_signals = 1;
@@ -113,8 +130,8 @@ for dbmode = 0%length(rates)
Symbols.spectrum("fignum",20,"normalizeTo0dB",1,"displayname",'Full Response','addDCoffset',0,'color',clr.Set1.red,'normalizeToNyquist',0,'linestyle','--');
DB_Symbols.spectrum("fignum",20,"normalizeTo0dB",1,"displayname",'DB-Response','addDCoffset',0,'color',clr.Set1.blue,'normalizeToNyquist',0,'linestyle','--');
Scpe_sig_avg.spectrum("displayname","Scope PSD","fignum",20,"normalizeTo0dB",1,"addDCoffset",5);
Scpe_sig_avg.plot("displayname","Scope raw signal","fignum",27,"clear",1);
Scpe_sig_avg.eye(fsym,M,"fignum",48,"displayname",' Eye of AVG Signal');
% Scpe_sig_avg.plot("displayname","Scope raw signal","fignum",27,"clear",1);
% Scpe_sig_avg.eye(fsym,M,"fignum",48,"displayname",' Eye of AVG Signal');
end
@@ -128,9 +145,9 @@ for dbmode = 0%length(rates)
xticks(-100:20:100);
yticks(-20:10:10);
beautifyBERplot("logscale",0,"setmarkers",0)
pos = [100.3333 991.6667 358.0000 192.6667];
set(fig, 'Position', pos);
beautifyBERplot("logscale",0,"setmarkers",0,"setcolors",0,"changemarkers",0)
% pos = [100.3333 991.6667 358.0000 192.6667];
% set(fig, 'Position', pos);
%%%%%%%%%%%%
drawnow;
@@ -164,7 +181,7 @@ for dbmode = 0%length(rates)
%%% DB TARGET
db_ref_sequence = Duobinary().encode(Symbols);
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",0);
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",0);
[eq_signal, db_noise] = eq_.process(Scpe_sig,db_ref_sequence);
% db_noise.signal = db_noise.signal - mean(db_noise.signal);
@@ -187,8 +204,8 @@ for dbmode = 0%length(rates)
yticks(-50:10:10);
beautifyBERplot("logscale",0,"setmarkers",0)
pos = [100.3333 991.6667 358.0000 192.6667];
set(fig, 'Position', pos);
% pos = [100.3333 991.6667 358.0000 192.6667];
% set(fig, 'Position', pos);
end

View File

@@ -1,160 +0,0 @@
%% ============================================================
% PARAMETERS
% ============================================================
database_type = 'mysql';
db = DBHandler("dataBase", "labor_highspeed", "type", database_type);
pam_level = 4; % FIXED for this figure
baudrates = [300e9 330e9 360e9 390e9];
fiberL = 10;
fields = [
db.getTableFieldNames('power_state_info');
db.getTableFieldNames('dashboard_ungrouped_alltime')
];
%% ============================================================
% DEFINE DSP SCHEMES
% ============================================================
curves = struct;
curves(1).name = 'VNLE';
curves(1).eq = equalizer_structure.vnle;
curves(1).color = clr.Paired.red;
curves(2).name = 'PF + MLSE';
curves(2).eq = equalizer_structure.vnle_pf_mlse;
curves(2).color = clr.Paired.green;
curves(3).name = 'DB-target + MLSE';
curves(3).eq = equalizer_structure.vnle_db_mlse;
curves(3).color = clr.Paired.blue;
curves(4).name = 'ML-based MLSE';
curves(4).eq = equalizer_structure.ml_mlse;
curves(4).color = clr.Paired.purple;
%% ============================================================
% ANALYSIS — results(b, k): b = baudrate index, k = DSP scheme index
% ============================================================
results = struct;
for b = 1:length(baudrates)
Rb = baudrates(b);
% --- query matching runs ---
fp = QueryFilter();
fp.where('Runs','pam_level','EQUALS', pam_level);
fp.where('Runs','fiber_length','EQUALS', fiberL);
fp.where('Runs','bitrate','EQUALS', Rb);
fp.where('Runs','is_mpi','EQUALS', 0);
[dataTable, ~] = db.queryDB(fp, fields);
for k = 1:numel(curves)
%% ---- DECIDE PRE-EMPH & PRECoded RULES for PAM-4 ----
pre_emph = decide_preemph(pam_level, curves(k).eq);
use_precoded = decide_precoded(pam_level, curves(k).eq);
%% ---- SETUP ANALYSIS CONFIG ----
cfg = struct;
cfg.x_axis = 'wavelength';
cfg.y_axis = 'BER';
cfg.agg = 'min';
cfg.outlier = 'none';
% cfg.group_by = {'wavelength'};
cfg.show_raw = false;
cfg.filters = struct( ...
'pam_level', pam_level, ...
'is_mpi', 0, ...
'bitrate', Rb, ...
'fiber_length', fiberL, ...
'equalizer_structure', curves(k).eq, ...
'pre_emph', pre_emph);
%% ---- RUN ANALYSIS ----
A = analyze_measurements_gpt(dataTable, cfg);
results(b,k).wavelength = A.group{1}.x;
if use_precoded
results(b,k).ber = A.group{1}.y_precoded;
else
results(b,k).ber = A.group{1}.y;
end
end
end
%% ============================================================
% PLOT — 1×4 (one tile per baudrate)
% ============================================================
fig = figure(); clf;
tiledlayout(1,4,'TileSpacing','compact','Padding','compact');
lw = 1.8;
ms = 6;
for b = 1:length(baudrates)
nexttile; hold on;
for k = 1:numel(curves)
plot(results(b,k).wavelength, results(b,k).ber, ...
'-o', ...
'Color', curves(k).color, ...
'MarkerFaceColor', curves(k).color, ...
'MarkerSize', ms, ...
'LineWidth', lw, ...
'DisplayName', curves(k).name);
end
set(gca,'YScale','log');
grid on;
xlabel('Wavelength [nm]');
ylabel('BER');
ylim([1e-4 0.1]);
title(sprintf('PAM-%d @ %.0f GBd',pam_level, baudrates(b)/1e9));
legend('Location','best');
beautifyBERplot();
end
% Optional figure size
pos = 1e3.*[0.1 0.55 1.4 0.32];
set(fig, 'Position', pos);
%% ============================================================
% DECISION LOGIC (INLINE FUNCTIONS)
% ============================================================
function pe = decide_preemph(M, eq)
% PRE-EMPH RULES:
switch M
case 4
if eq == equalizer_structure.vnle
pe = 1; % PAM4: VNLE → pre-emph on
else
pe = 0; % PAM4: all others → off
end
case {6,8}
pe = 1; % PAM6/8: all → pre-emph on
otherwise
pe = 0;
end
end
function flag = decide_precoded(M, eq)
% PRE-CODE RULES:
if eq == equalizer_structure.vnle_db_mlse
flag = 1; % Always for DB-target
elseif eq == equalizer_structure.ml_mlse && M == 4
flag = 1; % PAM4: ML-based → precoded
else
flag = 0;
end
end

View File

@@ -107,14 +107,14 @@ if 1
fprintf('Plotting: %s\n', precomp_filename);
freqresp.plot();
outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\spectrum_2.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\spectrum_2.tikz';
% matlab2tikz(outfile, ...
% 'width','\fwidth', ...
% 'height','\fheight', ...
% 'showInfo',false, ...
% 'extraAxisOptions',{ ...
% 'legend style={font=\footnotesize}', ...
% 'legend columns=1' ...
% });
end

View File

@@ -8,6 +8,7 @@ fp.where('power_state_info', 'pam_level','EQUALS', 4);
fp.where('power_state_info', 'db_mode','EQUALS', 0);
% fp.where('power_state_info', 'fiber_length','EQUALS', 1);
fp.where('power_state_info', 'is_mpi','EQUALS', 0);
fp.where('rop_attenuation', 'is_mpi','EQUALS', 0);
fields = db.getTableFieldNames('power_state_info');
[dataTable,~] = db.queryDB(fp, fields);
@@ -60,7 +61,7 @@ for fl = 1:numel(fiber_len)
%
dname = sprintf('%s; %d km',y_variable, fiber_len(fl));
h2 = plot(fl_filtered_.(x_variable), fl_filtered_.(['mean_',y_variable]), 'LineWidth', 1, 'MarkerSize', 4,'Marker','o','LineStyle','-','Color',cols(fl,:),'MarkerFaceColor',cols(fl,:),'DisplayName',dname);
h2 = scatter(fl_filtered_.(x_variable), fl_filtered_.(['mean_',y_variable]), 36, 'Marker','o', 'MarkerEdgeColor',cols(fl,:), 'MarkerFaceColor',cols(fl,:), 'DisplayName',dname);
h2.DataTipTemplate.DataTipRows(end+1) = ...
dataTipTextRow('run\_id', run_ids);
@@ -94,3 +95,28 @@ yline([4.85e-3, 2e-2],'--','LineWidth',1,'HandleVisibility','off');
posH = get(f, 'Position'); % [left, bottom, width, height]
newPos = [posH(1), posH(2), 750, 300];
set(f, 'Position', newPos);
%% Average laser power over all measurements for each wavelength
valid_rows = ~ismissing(dataTable.wavelength) & ...
~ismissing(dataTable.power_laser);
laser_measurements = dataTable(valid_rows, :);
laser_power_by_wavelength = groupsummary( ...
laser_measurements, ...
'wavelength', ...
'mean', ...
'power_mzm');
laser_power_by_wavelength = sortrows( ...
laser_power_by_wavelength, 'wavelength', 'ascend');
f_all = figure(4);
clf(f_all);
scatter(laser_power_by_wavelength.wavelength, ...
laser_power_by_wavelength.mean_power_mzm, ...
36, 'o', 'filled');
grid on;
xlabel('Wavelength in nm', 'FontSize', 12);
ylabel('Average laser power in dB', 'FontSize', 12);
title('Average laser power over all measurements', ...
'FontSize', 14, 'FontWeight', 'bold');
set(gca, 'FontSize', 11);

View File

@@ -1,5 +1,5 @@
filename = "Z:\2024\sioe\High Speed Messungen Oktober\bias_5km\PAMX_5km_20241025_204334_wh.mat";
filename = "W:\labdata\sioe_labor\bias_5km\PAMX_5km_20241025_204334_wh.mat";
a = load(filename);
wh = a.obj;
@@ -30,7 +30,7 @@ for l = 1:numel(lambda_vals)
for m = 1:numel(M_vals)
ber_ffe = wh.getStoValue('ber_ffe',v_bias_vals,awg_vpp_vals(1),precomp_amp_max_vals(1),rop_atten_vals(1),M_vals(m),lambda_vals(l));
ber = wh.getStoValue('ber_ffe',v_bias_vals,awg_vpp_vals(1),precomp_amp_max_vals(1),rop_atten_vals(1),M_vals(m),lambda_vals(l));
ber = wh.getStoValue('ber_collect',v_bias_vals,awg_vpp_vals(1),precomp_amp_max_vals(1),rop_atten_vals(1),M_vals(m),lambda_vals(l));
exfo = wh.getStoValue('exfo',v_bias_vals,awg_vpp_vals(1),precomp_amp_max_vals(1),rop_atten_vals(1),M_vals(m),lambda_vals(l));
lb = wh.getStoValue('exfo',v_bias_vals,awg_vpp_vals(1),precomp_amp_max_vals(1),rop_atten_vals(1),M_vals(m),lambda_vals(l));

View File

@@ -0,0 +1,158 @@
% Minimal all-time database plot: PAM-4 baudrate sweep versus VNLE noise.
% The dashed curves are the Burg postfilter responses estimated by
% VNLE + postfilter + MLSE.
%% SettingsC:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Advanced_DSP_for_400G_IMDD_experiments\Auswertung_JLT\final\FIGURE_EQ_NOISE_VS_BAUDRATE.m
M = 4;
bitrate = [390e9]; % [] -> all matching baudrates in the DB
dsp_options.storage_path = "W:\labdata\sioe_labor\";
dsp_options.max_occurences = 1; % standard routine: first synchronized signal
fiber_length = 10;
wavelength = 1310;
is_mpi = 0;
db_mode_filter = int32(db_mode.no_db);
rop_attenuation = 0;
vnle_order = [50 5 5];
dfe_order = [0 0 0];
training_length = 4096*2;
postfilter_order = 1;
%% Query the all-time database
db = DBHandler( ...
"dataBase", "labor_highspeed", ...
"type", "mysql", ...
"server", "192.168.178.192", ...
"user", "silas", ...
"password", "silas");
fp = QueryFilter();
fp.where('Runs', 'pam_level', SqlOperator.EQUALS, M);
fp.where('Runs', 'fiber_length', SqlOperator.EQUALS, fiber_length);
fp.where('Runs', 'wavelength', SqlOperator.EQUALS, wavelength);
fp.where('Runs', 'is_mpi', SqlOperator.EQUALS, is_mpi);
fp.where('Runs', 'db_mode', 'LESS_THAN', 2);
fp.where('Runs', 'rop_attenuation', SqlOperator.EQUALS, rop_attenuation);
% fp.where('Runs', 'symbolrate','EQUALS', 180e9);
% if ~isempty(baudrate_GBd)
% fp.where('Runs', 'symbolrate', SqlOperator.IN, baudrate_GBd*1e9);
% end
fields = [db.getTableFieldNames('dashboard_ungrouped_alltime')];
fields = unique(fields, 'stable');
[runs, ~] = db.queryDB(fp, fields);
if ~isempty(bitrate)
runs = runs(ismember(runs.bitrate, bitrate), :);
end
if isempty(runs)
error('FIGURE_EQ_NOISE_VS_BAUDRATE:NoData', ...
'No PAM-%d runs matched the selected baudrate range.', M);
end
run_ids = unique(runs.run_id, 'stable');
%% One plot: VNLE residual noise and Burg tap response
fig = figure(220);
clf(fig);
hold on;
max_fs = 0;
for runIndex = 1:numel(run_ids)
runData = queryRunid(run_ids(runIndex), db);
fsym = double(runData.symbolrate(1));
bitrate = double(runData.bitrate(1));
dsp_options.max_occurences = 20;
dsp_options.start_occurence = 2;
[Tx_bits, Symbols, Scpe_cell, found_sync] = loadAndSyncRunSignals(runData, dsp_options);
if ~found_sync || isempty(Scpe_cell)
warning('Skipping run %g: no synchronized signal found.', runData.run_id(1));
continue
end
useavg = 0;
Scpe_sig = averageScopeSignals(Scpe_cell, useavg, 1);
Scpe_sig = preprocessSignal(Scpe_sig, Symbols, fsym);
Scpe_sig = Scpe_sig.normalize("mode", "rms");
% VNLE: retain the difference between equalizer output and reference.
eq_vnle = EQ( ...
"Ne", vnle_order, ...
"Nb", dfe_order, ...
"training_length", training_length, ...
"training_loops", 5, ...
"dd_loops", 5, ...
"K", 2, ...
"DCmu", 0.001, ...
"DDmu", [0.0004 0.0004 0.0004 0.0004], ...
"DFEmu", 0.05, ...
"FFEmu", 0, ...
"plotfinal", 0, ...
"ideal_dfe", 1);
% eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ...
% "mu_dd",1e-1,"mu_tr",0.4,"order",50, ...
% "sps",2,"decide",0,"optmize_mus",1,"dd_mode",1, ...
% "adaption_technique","nlms","dc_tracking_mu",1.021e-05);
[equalized_signal, eq_noise] = eq_vnle.process(Scpe_sig, Symbols);
postfilter_order = 3;
pf = Postfilter("ncoeff", postfilter_order, "useBurg", 1);
[mlse_sig_sd,whitened_noise] = pf.process(equalized_signal, eq_noise);
% mlse = MLSE( ...
% "DIR", [0 0], ...
% "duobinary_output", 0, ...
% "M", M, ...
% "trellis_states", PAMmapper(M, 0).levels);
%
% [results, ~] = vnle_postfilter_mlse( ...
% eq_vnle, pf, mlse, M, Scpe_sig, Symbols, Tx_bits, ...
% "precode_mode", db_mode.no_db, ...
% "showAnalysis", 0, ...
% "postFFE", [], ...
% "eth_style_symbol_mapping", 0);
eq_noise = eq_noise - mean(eq_noise.signal);
fig = figure(220+runIndex);
showEQNoisePSD(eq_noise, ...
"fignum", fig.Number, ...
"displayname", sprintf('%.0f Gbps: VNLE Noise', bitrate*1e-9), ...
"postfilter_taps", pf.coefficients, ...
"colormode", "qualitative");
whitened_noise = whitened_noise - mean(whitened_noise.signal);
whitened_noise.spectrum("displayname", 'Whitened Noise', "fignum", fig.Number, "normalizeTo0dB", 0,"fft_length",4096,"normalizeToDC",0);
max_fs = max(max_fs, eq_noise.fs);
end
xlabel('Frequency in GHz');
ylabel('normalized to 0 dB');
title('Noise of soft decision signal (not MLSE)');
grid on;
grid minor;
legend('show', 'Interpreter', 'none', 'Location', 'best');
xlim([-max_fs/2 max_fs/2]*1e-9);
ylim([-20 0]);
function signalOut = averageScopeSignals(scopeCell, averageSignals, gain)
signalOut = scopeCell{1};
if ~averageSignals
return
end
commonLength = min(cellfun(@(s) numel(s.signal), scopeCell));
scopeMean = zeros(commonLength, 1);
for idx = 1:numel(scopeCell)
scopeMean = scopeMean + scopeCell{idx}.signal(1:commonLength);
end
scopeMean = scopeMean ./ numel(scopeCell);
signalOut.signal = gain .* scopeMean;
end

View File

@@ -8,7 +8,7 @@ clear; clc;
%% 1) Query data
selectedPamLevel = 8;
selectedPamLevels = [4, 6, 8];
selectedFiberLengthKm = 10;
selectedWavelengthNm = 1310;
selectedRopAttenuation = 0; % set [] to use all ROP attenuation values
@@ -33,7 +33,6 @@ db.refresh();
fp = QueryFilter();
fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm);
fp.where('Runs', 'pam_level', 'EQUALS', selectedPamLevel);
fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm);
if ~isempty(selectedRopAttenuation)
fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation);
@@ -66,6 +65,8 @@ for fieldIdx = 1:numel(numericFields)
end
end
data = data(ismember(data.pam_level, selectedPamLevels), :);
if ~ismember("precomp_amp", string(data.Properties.VariableNames))
warning("plot_best_algos:NoPrecompAmp", ...
"Runs.precomp_amp was not returned. Falling back to pre_emphasis = (db_mode == 0).");
@@ -85,6 +86,9 @@ if ismember("equalizer_structure", string(duobinaryRows.Properties.VariableNames
equalizerMask(duobinaryRows.equalizer_structure, ...
equalizer_structure.db_encoded), :);
end
duobinaryRows = duobinaryRows(datetime(duobinaryRows.date_of_processing)>datetime("2026-01-01 00:00:00"),:);
duobinaryPlotData = buildDuobinarySignalingRows(duobinaryRows);
duobinaryPlotData = duobinaryPlotData(isfinite(duobinaryPlotData.BER_plot) & ...
duobinaryPlotData.BER_plot > 0 & ...
@@ -117,70 +121,75 @@ if isempty(availableStyles)
return
end
fig = figure(); clf;
ax = axes(fig); hold(ax, "on");
fig = figure(432); clf;
t = tiledlayout(fig, 1, numel(selectedPamLevels), ...
"TileSpacing", "compact", ...
"Padding", "compact");
for styleIdx = 1:height(availableStyles)
style = availableStyles(styleIdx, :);
rowMask = bestPlotData.algorithm_key == style.algorithm_key;
if ~any(rowMask)
continue
for pamIdx = 1:numel(selectedPamLevels)
selectedPamLevel = selectedPamLevels(pamIdx);
ax = nexttile(t); hold(ax, "on");
pamMask = bestPlotData.pam_level == selectedPamLevel;
for styleIdx = 1:height(availableStyles)
style = availableStyles(styleIdx, :);
rowMask = pamMask & bestPlotData.algorithm_key == style.algorithm_key;
if ~any(rowMask)
continue
end
algoData = sortrows(bestPlotData(rowMask, :), "grossrate_Gbps");
if showRawEntries
scatter(ax, algoData.grossrate_Gbps, algoData.BER_plot, ...
9, ...
"Marker", ".", ...
"MarkerEdgeColor", style.color, ...
"MarkerFaceColor", style.color, ...
"MarkerEdgeAlpha", 0.25, ...
"MarkerFaceAlpha", 0.25, ...
"HandleVisibility", "off");
end
if showBestLine
plot(ax, algoData.grossrate_Gbps, algoData.BER_plot, ...
"LineStyle", style.lineStyle, ...
"Marker", style.marker, ...
"MarkerSize", 5, ...
"LineWidth", 1.5, ...
"Color", style.color, ...
"MarkerFaceColor", style.markerFaceColor, ...
"MarkerEdgeColor", style.color, ...
"DisplayName", style.name);
end
end
algoData = sortrows(bestPlotData(rowMask, :), "grossrate_Gbps");
yline(ax, [2.2e-4, 4.85e-3, 2e-2], ...
"LineWidth", 1, ...
"LineStyle", "--", ...
"Color", [0.25 0.25 0.25], ...
"HandleVisibility", "off");
if showRawEntries
scatter(ax, algoData.grossrate_Gbps, algoData.BER_plot, ...
9, ...
"Marker", ".", ...
"MarkerEdgeColor", style.color, ...
"MarkerFaceColor", style.color, ...
"MarkerEdgeAlpha", 0.25, ...
"MarkerFaceAlpha", 0.25, ...
"HandleVisibility", "off");
end
% title(ax, sprintf("PAM-%d, %.0f km, %.0f nm", ...
% selectedPamLevel, selectedFiberLengthKm, selectedWavelengthNm));
xlabel(ax, "Gross rate [Gb/s]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
ylim(ax, [8e-5, 0.1]);
grid(ax, "on");
box(ax, "on");
if showBestLine
plot(ax, algoData.grossrate_Gbps, algoData.BER_plot, ...
"LineStyle", style.lineStyle, ...
"Marker", style.marker, ...
"MarkerSize", 5, ...
"LineWidth", 1.5, ...
"Color", style.color, ...
"MarkerFaceColor", style.markerFaceColor, ...
"MarkerEdgeColor", style.color, ...
"DisplayName", style.name);
xticks(ax, 300:30:480);
xlim(ax, [300, 480]);
legend(ax, "Location", "northeast", "Interpreter", "none");
if exist("beautifyBERplot", "file")
beautifyBERplot("logscale", true, "setcolors", false, ...
"setmarkers", false, "changemarkers", false);
end
end
yline(ax, [2.2e-4, 4.85e-3, 2e-2], ...
"LineWidth", 1, ...
"LineStyle", "--", ...
"Color", [0.25 0.25 0.25], ...
"HandleVisibility", "off");
title(ax, sprintf("PAM-%d, %.0f km, %.0f nm", ...
selectedPamLevel, selectedFiberLengthKm, selectedWavelengthNm));
xlabel(ax, "Gross rate [Gb/s]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
ylim(ax, [1e-5, maxBerForPlot]);
grid(ax, "on");
box(ax, "on");
xTicks = unique(bestPlotData.grossrate_Gbps(isfinite(bestPlotData.grossrate_Gbps)));
if ~isempty(xTicks)
xticks(ax, xTicks);
xlim(ax, [min(xTicks), max(xTicks)]);
end
legend(ax, "Location", "best", "Interpreter", "none");
if exist("beautifyBERplot", "file")
beautifyBERplot("logscale", true, "setcolors", false, ...
"setmarkers", false, "changemarkers", false);
end
set(fig, "Position", 1e3 .* [0.1000 0.5500 0.7200 0.4200]);
set(fig, "Position", 1e3 .* [0.1070 0.5497 1.0585 0.2282]);
%% Local helpers
@@ -220,7 +229,7 @@ styles = table( ...
"VNLE + PF + MLSE"; ...
"VNLE DBt. + MLSE"; ...
"ML pre-EQ + Viterbi"; ...
"Duobinary signaling"], ...
"DBS + VNLE + MLSE"], ...
["o"; "square"; "diamond"; "^"; "v"], ...
["-"; "-"; "-"; "-"; "-"], ...
["w"; "w"; "w"; "w"; "w"], ...
@@ -322,7 +331,7 @@ value = double(enumEntry);
end
function bestData = bestBerByAlgorithmAndGrossRate(data)
groupVars = ["algorithm_key", "grossrate_Gbps"];
groupVars = ["pam_level", "algorithm_key", "grossrate_Gbps"];
groupId = findgroups(data(:, groupVars));
keepIdx = NaN(max(groupId), 1);

View File

@@ -8,7 +8,7 @@ clear; clc;
%% 1) Gather data
selectedPamLevel = 4;
selectedFiberLengthKm = 10;
selectedFiberLengthKm = 2;
selectedWavelength =1310; % set [] to use all wavelengths
selectedRopAttenuation = []; % set [] to use all ROP attenuation values
selectedIsMpi = []; % set [] to use all entries
@@ -98,7 +98,7 @@ if isempty(availableEqStyles)
return
end
fig = figure(401); clf;
fig = figure(); clf;
tiledlayout(1, height(availableEqStyles), ...
"TileSpacing", "compact", ...
"Padding", "compact");

View File

@@ -14,7 +14,9 @@ selectedWavelengthNm = 1310;
selectedRopAttenuation = 0;
selectedIsMpi = 0;
selectedDbMode = db_mode.db_encoded;
selectedEqualizerStructure = equalizer_structure.db_encoded;
selectedEqualizerStructures = [ ...
equalizer_structure.db_encoded, ...
equalizer_structure.ml_mlse];
maxBerForPlot = 0.5;
showRawEntries = false;
@@ -60,8 +62,13 @@ end
data = data(ismember(data.pam_level, selectedPamLevels), :);
if ismember("equalizer_structure", string(data.Properties.VariableNames)) && ...
~isempty(selectedEqualizerStructure)
data = data(equalizerMask(data.equalizer_structure, selectedEqualizerStructure), :);
~isempty(selectedEqualizerStructures)
eqMask = false(height(data), 1);
for eqIdx = 1:numel(selectedEqualizerStructures)
eqMask = eqMask | equalizerMask(data.equalizer_structure, ...
selectedEqualizerStructures(eqIdx));
end
data = data(eqMask, :);
end
plotData = buildDetectionMetricRows(data);
@@ -150,6 +157,9 @@ for pamIdx = 1:numel(availablePamLevels)
xlim(ax, [min(xTicks), max(xTicks)]);
end
xticks(ax, 300:30:480);
xlim(ax, [300, 480]);
legend(ax, "Location", "best", "Interpreter", "none");
if exist("beautifyBERplot", "file")
beautifyBERplot("logscale", true, "setcolors", false, ...
@@ -182,31 +192,52 @@ values = double(values);
end
function plotData = buildDetectionMetricRows(data)
baseRows = data(isfinite(data.BER), :);
dbEncodedRows = data(equalizerMask(data.equalizer_structure, ...
equalizer_structure.db_encoded), :);
baseRows = dbEncodedRows(isfinite(dbEncodedRows.BER), :);
baseRows.detection_type = repmat("VNLE + MLSE", height(baseRows), 1);
baseRows.BER_plot = baseRows.BER;
if ismember("BER_precoded", string(data.Properties.VariableNames))
memorylessRows = data(isfinite(data.BER_precoded), :);
memorylessRows = dbEncodedRows(isfinite(dbEncodedRows.BER_precoded), :);
memorylessRows.detection_type = repmat("VNLE + memoryless", ...
height(memorylessRows), 1);
memorylessRows.BER_plot = memorylessRows.BER_precoded;
plotData = [baseRows; memorylessRows];
else
warning("plot_duobinary_detection:NoPrecodedBer", ...
"BER_precoded was not returned. Plotting only VNLE + MLSE rows.");
plotData = baseRows;
"BER_precoded was not returned. Plotting only BER rows.");
memorylessRows = dbEncodedRows([], :);
end
mlMlseRows = data(equalizerMask(data.equalizer_structure, ...
equalizer_structure.ml_mlse), :);
if ismember("BER_precoded", string(data.Properties.VariableNames))
mlMlsePrecodedRows = mlMlseRows(isfinite(mlMlseRows.BER_precoded), :);
mlMlsePrecodedRows.detection_type = repmat( ...
"ML pre-EQ + Viterbi", ...
height(mlMlsePrecodedRows), 1);
mlMlsePrecodedRows.BER_plot = mlMlsePrecodedRows.BER_precoded;
else
mlMlsePrecodedRows = mlMlseRows([], :);
end
plotData = [baseRows; memorylessRows; mlMlsePrecodedRows];
end
function styles = defaultDetectionStyles()
styles = table( ...
["VNLE + MLSE"; "VNLE + memoryless"], ...
["VNLE + MLSE"; "VNLE + memoryless"], ...
["o"; "square"], ...
["-"; "--"], ...
["w"; "none"], ...
[clr.Paired.blue; clr.Paired.orange], ...
["VNLE + MLSE"; ...
"VNLE + memoryless"; ...
"ML pre-EQ + Viterbi"], ...
["VNLE + MLSE"; ...
"VNLE + memoryless"; ...
"ML pre-EQ + Viterbi"], ...
["o"; "square"; "^"], ...
["-"; "--"; "-"], ...
["w"; "none"; "w"], ...
[clr.Paired.blue; ...
clr.Paired.orange; ...
clr.Paired.purple], ...
'VariableNames', ["detection_type", "name", "marker", ...
"lineStyle", "markerFaceColor", "color"]);
end

View File

@@ -0,0 +1,384 @@
%% 400G revisit: eye comparison for PAM-4/6/8
% The native Signal.eye() method is used for the eye display. Optional TikZ
% export is configured below and uses mat2tikz_improved().
%
% db_mode_setting:
% 0 - no duobinary processing; use the copied FFE configuration
% 1 - DB-targeted VNLE; the VNLE target is Duobinary().encode(Symbols)
% 2 - DB-encoded data; equalize the stored DB-encoded Symbols with a VNLE
clc;
%% Configuration
M_values = [4];
selectedBitrate = 360e9;
selectedFiberLengthKm = 10;
selectedWavelengthNm = 1310;
selectedRopAttenuation = 0;
selectedIsMpi = 0;
db_mode_setting = 2; % Set to 0, 1, or 2.
average_signals = 1;
show_histograms = true;
export_tikz = 0;
max_occurences = 10;
histogram_figure_base = 600;
eye_figure_base = 500;
psd_figure_base = 700;
% PSD colors follow the established DSP algorithm styles. The target is
% black and the received waveform is grey.
dsp_styles = defaultAlgorithmStyles();
psd_color_target = [0, 0, 0];
psd_color_received = [0.3500, 0.3500, 0.3500];
psd_color_equalized = dsp_styles.color(dsp_styles.algorithm_key == "vnle", :);
tikz_root = 'C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\04_Experimental_Evaluation\tikz\400g';
tikz_eye_folder = fullfile(tikz_root, 'eyes');
tikz_histogram_folder = fullfile(tikz_root, 'histograms');
tikz_psd_folder = fullfile(tikz_root, 'psd');
% The averaging/gain/normalization follows FIGURE_EYES.m.
average_gain = 1.25;
% Common equalizer parameters copied from FIGURE_EYES.m.
len_tr = 4096*2;
training_loops = 5;
dd_loops = 5;
eq_K = 2;
mu_dc = 0.005;
mu_ffe = [0.0001, 0.0008, 0.001];
mu_dfe = 0.0004;
if ~ismember(db_mode_setting, [0, 1, 2])
error('PLOT_EYES_400G_REVISIT:InvalidDbMode', ...
'db_mode_setting must be 0, 1, or 2.');
end
if export_tikz
if ~exist(tikz_eye_folder, 'dir')
mkdir(tikz_eye_folder);
end
if ~exist(tikz_histogram_folder, 'dir')
mkdir(tikz_histogram_folder);
end
if ~exist(tikz_psd_folder, 'dir')
mkdir(tikz_psd_folder);
end
end
%% Database setup
dsp_options = struct();
dsp_options.storage_path = 'W:\labdata\sioe_labor\';
dsp_options.start_occurence = 1;
dsp_options.max_occurences = max_occurences;
database = DBHandler( ...
"dataBase", "labor_highspeed", ...
"type", "mysql", ...
"server", "192.168.178.192", ...
"user", "silas", ...
"password", "silas");
%% DSP phase: query, synchronize, equalize, and collect all signals
selectedRuns = [];
eyeSignals = cell(size(M_values));
rawSignals = cell(size(M_values));
eyeReferences = cell(size(M_values));
symbolRates = NaN(size(M_values));
for mIdx = 1:numel(M_values)
M = M_values(mIdx);
fp = QueryFilter();
fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm);
fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm);
fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation);
fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi);
fp.where('Runs', 'pam_level', 'EQUALS', M);
fp.where('Runs', 'db_mode', 'EQUALS', db_mode_setting);
if ~isempty(selectedBitrate)
fp.where('Runs', 'bitrate', 'EQUALS', selectedBitrate);
end
[dataTable, ~] = database.queryDB(fp, database.getTableFieldNames('Runs'));
if isempty(dataTable)
error('PLOT_EYES_400G_REVISIT:MissingRun', ...
'No run found for PAM-%d with the selected filters.', M);
end
dataTable = sortrows(dataTable, {'bitrate', 'run_id'});
dataTable = dataTable(1, :);
if isempty(selectedRuns)
selectedRuns = dataTable;
else
selectedRuns = [selectedRuns; dataTable]; %#ok<AGROW>
end
fsym = double(dataTable.symbolrate(1));
symbolRates(mIdx) = fsym;
[~, Symbols, Scpe_cell, found_sync] = ...
loadAndSyncRunSignals(dataTable, dsp_options);
if ~found_sync || isempty(Scpe_cell)
error('PLOT_EYES_400G_REVISIT:SynchronizationFailed', ...
'Could not synchronize the selected PAM-%d run.', M);
end
Scpe_sig = averageScopeSignals(Scpe_cell, average_signals, average_gain);
Scpe_sig = preprocessSignal(Scpe_sig, Symbols, fsym);
Scpe_sig = Scpe_sig.normalize("mode", "rms");
switch db_mode_setting
case 0
% Simply FFE, copied from FIGURE_EYES.m.
eq_ = EQ("Ne", [50, 0, 0], ...
"Nb", [0, 0, 0], ...
"training_length", len_tr, ...
"training_loops", training_loops, ...
"dd_loops", dd_loops, ...
"K", eq_K, ...
"DCmu", mu_dc, ...
"DDmu", [mu_ffe, mu_dfe], ...
"DFEmu", 0.005, ...
"FFEmu", 0, ...
"plotfinal", 0, ...
"ideal_dfe", false);
referenceForEye = Symbols;
case 1
% DB-targeted VNLE: train against the DB waveform generated
% from the stored (precoded but not DB-encoded) symbols.
eq_ = makeVnle(len_tr, training_loops, dd_loops, eq_K, ...
mu_dc, mu_ffe, mu_dfe);
referenceForEye = Duobinary().encode(Symbols, "M", M);
case 2
% DB-encoded data: the stored Symbols already contain the DB
% encoded reference waveform used by the VNLE.
eq_ = makeVnle(len_tr, training_loops, dd_loops, eq_K, ...
mu_dc, mu_ffe, mu_dfe);
referenceForEye = Symbols;
end
[equalized_signal, ~] = eq_.process(Scpe_sig, referenceForEye);
eyeSignals{mIdx} = equalized_signal;
rawSignals{mIdx} = Scpe_sig;
eyeReferences{mIdx} = referenceForEye;
fprintf('PAM-%d: run %d, %.3f GBd, db_mode %d\n', ...
M, dataTable.run_id(1), fsym*1e-9, db_mode_setting);
end
%% Plotting and export phase
show_histograms = 1;
showEye = 1;
showPSD = 1;
for mIdx = 1:numel(M_values)
M = M_values(mIdx);
fsym = symbolRates(mIdx);
rawSignal = rawSignals{mIdx};
equalized_signal = eyeSignals{mIdx};
referenceForEye = eyeReferences{mIdx};
if showPSD
psdFigure = psd_figure_base + mIdx + db_mode_setting;
figure(psdFigure); clf;
switch db_mode_setting
case 0
referencePsdName = sprintf('PAM-%d reference', M);
equalizedPsdName = 'FFE output';
psd_color_equalized = dsp_styles.color( ...
dsp_styles.algorithm_key == "vnle", :);
case 1
referencePsdName = 'DB-target reference';
equalizedPsdName = 'DB-targeted VNLE output';
psd_color_equalized = dsp_styles.color( ...
dsp_styles.algorithm_key == "vnle_db_mlse", :);
case 2
referencePsdName = 'DB-encoded reference';
equalizedPsdName = 'DB-encoded VNLE output';
psd_color_equalized = dsp_styles.color( ...
dsp_styles.algorithm_key == "db_encoded", :);
end
referenceForEye.spectrum( ...
"fignum", psdFigure, ...
"show_onesided", 1, ...
"fft_length", 4096, ...
"displayname", referencePsdName, ...
"color", psd_color_target);
rawSignal.spectrum( ...
"fignum", psdFigure, ...
"show_onesided", 1, ...
"fft_length", 4096, ...
"displayname", 'Received signal', ...
"color", psd_color_received);
equalized_signal.spectrum( ...
"fignum", psdFigure, ...
"show_onesided", 1, ...
"fft_length", 4096, ...
"displayname", equalizedPsdName, ...
"color", psd_color_equalized);
eq_noise = equalized_signal - referenceForEye;
showEQNoisePSD(eq_noise, ...
"fignum", 800, ...
"displayname", sprintf('%.0f Gbps: VNLE Noise', selectedBitrate), ...
"colormode", "diverging");
% PSDs are line-only, irrespective of the marker styles used by
% algorithm comparison plots.
set(findall(gca, 'Type', 'Line'), 'Marker', 'none');
figure(psdFigure);
title(sprintf('PSD comparison: %.0f GBd PAM-%d', fsym*1e-9, M), ...
'Interpreter', 'none');
legend('Location', 'best', 'Interpreter', 'none');
grid on;
box on;
if export_tikz
baudrateLabel = sprintf('%.0fGBd', fsym*1e-9);
psdFilename = fullfile(tikz_psd_folder, ...
sprintf('psd_pam_%d_%s_db_mode_%d.tikz', ...
M, baudrateLabel, db_mode_setting));
mat2tikz_improved(psdFilename);
end
end
if showEye
% Use the native Signal eye implementation for the displayed eye.
eyeFigure = eye_figure_base + mIdx+ db_mode_setting;
equalized_signal.eye(fsym, M, ...
"fignum", eyeFigure, ...
"displayname", sprintf('PAM-%d, db_mode = %d', M, db_mode_setting));
figure(eyeFigure);
set(gcf, "Position", [100 + 430*(mIdx - 1), 100, 400, 360]);
stripEyeAxes(gca);
if export_tikz
baudrateLabel = sprintf('%.0fGBd', fsym*1e-9);
eyeFilename = fullfile(tikz_eye_folder, ...
sprintf('eye_pam_%d_%s_db_mode_%d.tikz', ...
M, baudrateLabel, db_mode_setting));
mat2tikz_improved(eyeFilename);
end
end
if show_histograms || export_tikz
histogramFigure = histogram_figure_base + M + db_mode_setting;
plotHistogram(equalized_signal, referenceForEye, ...
M, histogramFigure, db_mode_setting, histogramLegendLabels(M));
if export_tikz
figure(histogramFigure);
histogramFilename = fullfile(tikz_histogram_folder, ...
sprintf('histogram_pam_%d_%s_db_mode_%d.tikz', ...
M, baudrateLabel, db_mode_setting));
mat2tikz_improved(histogramFilename,"cleanfigure",1);
end
end
end
%% Local functions
function signalOut = averageScopeSignals(scopeCell, averageSignals, gain)
signalOut = scopeCell{1};
if ~averageSignals
return
end
commonLength = min(cellfun(@(s) numel(s.signal), scopeCell));
scopeMean = zeros(commonLength, 1);
for idx = 1:numel(scopeCell)
scopeMean = scopeMean + scopeCell{idx}.signal(1:commonLength);
end
scopeMean = scopeMean ./ numel(scopeCell);
signalOut.signal = gain .* scopeMean;
end
function eq_ = makeVnle(lenTr, trainingLoops, ddLoops, eqK, muDc, muFfe, muDfe)
eq_ = EQ("Ne", [50, 5, 5], ...
"Nb", [0, 0, 0], ...
"training_length", lenTr, ...
"training_loops", trainingLoops, ...
"dd_loops", ddLoops, ...
"K", eqK, ...
"DCmu", muDc, ...
"DDmu", [muFfe, muDfe], ...
"DFEmu", 0.005, ...
"FFEmu", 0, ...
"plotfinal", 0, ...
"ideal_dfe", true);
end
function plotHistogram(equalizedSignal, referenceForEye, M, figureNumber, ...
dbModeSetting, legendLabels)
figure(figureNumber); clf;
if dbModeSetting ~= 0
uncodedReference = Duobinary().decode(referenceForEye, "M", M);
showLevelHistogram(equalizedSignal, referenceForEye, ...
"fignum", figureNumber, ...
"ref_symbol_uncoded", uncodedReference, ...
"legendLabels", legendLabels);
else
showLevelHistogram(equalizedSignal, referenceForEye, ...
"fignum", figureNumber, ...
"legendLabels", legendLabels);
end
title(sprintf('PAM-%d histogram, db\_mode = %d', M, dbModeSetting));
end
function labels = histogramLegendLabels(M)
labels = arrayfun(@(idx) sprintf('$p(z|d=d_{%d})$', idx), ...
1:M, 'UniformOutput', false);
end
function stripEyeAxes(ax)
title(ax, '');
xlabel(ax, '');
ylabel(ax, '');
set(ax, ...
'XTick', [], ...
'YTick', [], ...
'XTickLabel', [], ...
'YTickLabel', [], ...
'XMinorTick', 'off', ...
'YMinorTick', 'off');
grid(ax, 'off');
end
function styles = defaultAlgorithmStyles()
styles = table( ...
["vnle"; ...
"vnle_pf_mlse"; ...
"vnle_db_mlse"; ...
"ml_mlse"; ...
"db_encoded"], ...
[equalizer_structure.vnle; ...
equalizer_structure.vnle_pf_mlse; ...
equalizer_structure.vnle_db_mlse; ...
equalizer_structure.ml_mlse; ...
equalizer_structure.db_encoded], ...
["VNLE"; ...
"VNLE + PF + MLSE"; ...
"VNLE DBt. + MLSE"; ...
"ML pre-EQ + Viterbi"; ...
"DBS + VNLE + MLSE"], ...
["o"; "square"; "diamond"; "^"; "v"], ...
["-"; "-"; "-"; "-"; "-"], ...
["w"; "w"; "w"; "w"; "w"], ...
[clr.Paired.red; ...
clr.Paired.green; ...
clr.Paired.blue; ...
clr.Paired.purple; ...
clr.Paired.orange], ...
'VariableNames', ["algorithm_key", "eq", "name", "marker", ...
"lineStyle", "markerFaceColor", "color"]);
end

View File

@@ -0,0 +1,168 @@
%% Memoryless DB-target BER versus baudrate for PAM4/6/8
clear; clc;
warehouseFile = "C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Diss\400G_revisit\results_duobinary_eq_memoryless_2km_pam468.mat";
pamLevels = [4, 6, 8];
%% Load warehouse and matching run metadata
S = load(warehouseFile, "wh");
wh = S.wh;
db = DBHandler("dataBase", "labor_highspeed", ...
"type", "mysql", ...
"server", "192.168.178.192", ...
"user", "silas", ...
"password", "silas");
fp = QueryFilter();
fp.where('Runs', 'fiber_length', 'EQUALS', 2);
fp.where('Runs', 'wavelength', 'EQUALS', 1310);
fp.where('Runs', 'rop_attenuation', 'EQUALS', 0);
fp.where('Runs', 'is_mpi', 'EQUALS', 0);
fp.where('Runs', 'db_mode', 'EQUALS', 0);
[runTable, ~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
warehouseRunIds = warehouseIds(wh);
runTable = runTable(ismember(double(runTable.run_id), warehouseRunIds), :);
%% Extract BER precoded from dbtgt_package
warehouseRows = zeros(0, 3); % PAM level, baudrate [GBd], BER precoded
for k = 1:numel(wh.sto.dbtgt_package)
[phys, realizationResults] = wh.getPhysAndValueByLinIndex( ...
"dbtgt_package", k);
if ~isfield(phys, "run_id")
continue
end
runRow = find(double(runTable.run_id) == double(phys.run_id), 1);
if isempty(runRow)
continue
end
for realization = 1:numel(realizationResults)
package = realizationResults{realization};
berPrecoded = readMetric(package, "BER_precoded");
if isfinite(berPrecoded) && berPrecoded > 0
warehouseRows(end+1, :) = [ ...
double(runTable.pam_level(runRow)), ...
double(runTable.grossrate(runRow)) * 1e-9, ...
berPrecoded]; %#ok<AGROW>
end
end
end
warehouseCurve = minByBaudrate(warehouseRows, "BER_precoded");
%% Load database BER and BER precoded for VNLE + DB target/MLSE
selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime');
[databaseRows, ~] = db.queryDB(fp, selectedFields);
databaseRows = databaseRows(ismember(double(databaseRows.pam_level), pamLevels), :);
databaseRows = databaseRows(equalizerMask( ...
databaseRows.equalizer_structure, equalizer_structure.vnle_db_mlse), :);
berRows = databaseRows(isfinite(double(databaseRows.BER)) & ...
double(databaseRows.BER) > 0, :);
berPrecodedRows = databaseRows(isfinite(double(databaseRows.BER_precoded)) & ...
double(databaseRows.BER_precoded) > 0, :);
databaseBerCurve = minByBaudrate( ...
[double(berRows.pam_level), double(berRows.grossrate) * 1e-9, ...
double(berRows.BER)], "BER");
databaseBerPrecodedCurve = minByBaudrate( ...
[double(berPrecodedRows.pam_level), ...
double(berPrecodedRows.grossrate) * 1e-9, ...
double(berPrecodedRows.BER_precoded)], "BER_precoded");
%% Plot: three curves per PAM format
figure(470); clf;
tiledlayout(1, 3, "TileSpacing", "compact", "Padding", "compact");
for pamLevel = pamLevels
ax = nexttile; hold(ax, "on");
plotCurve(ax, warehouseCurve, pamLevel, "BER_precoded", ...
"o-", "precode + memoryless");
plotCurve(ax, databaseBerCurve, pamLevel, "BER", ...
"s-", "DBt. + MLSE");
plotCurve(ax, databaseBerPrecodedCurve, pamLevel, "BER_precoded", ...
"s--", "precode + DBt. + MLSE");
ylim([1e-4, 1e-1]);
title(ax, sprintf("PAM-%d", pamLevel));
xlabel(ax, "Baudrate [GBd]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
end
%% Local helpers
function runIds = warehouseIds(wh)
runIds = zeros(0, 1);
for k = 1:numel(wh.sto.dbtgt_package)
[phys, ~] = wh.getPhysAndValueByLinIndex("dbtgt_package", k);
if isfield(phys, "run_id")
runIds(end+1, 1) = double(phys.run_id); %#ok<AGROW>
end
end
runIds = unique(runIds);
end
function value = readMetric(package, metricName)
value = NaN;
if isempty(package) || ~isstruct(package) || ~isfield(package, "metrics")
return
end
metrics = package.metrics;
if isstruct(metrics) && isfield(metrics, metricName)
value = double(metrics.(metricName));
elseif isobject(metrics) && isprop(metrics, metricName)
value = double(metrics.(metricName));
end
end
function curve = minByBaudrate(rows, metricName)
if isempty(rows)
curve = table(zeros(0, 1), zeros(0, 1), ...
'VariableNames', ["pam_level", "baudrate_GBd"]);
curve.(metricName) = zeros(0, 1);
return
end
curve = table(rows(:, 1), rows(:, 2), rows(:, 3), ...
'VariableNames', ["pam_level", "baudrate_GBd", metricName]);
curve = groupsummary(curve, ["pam_level", "baudrate_GBd"], ...
"min", metricName);
curve.Properties.VariableNames(end) = metricName;
end
function plotCurve(ax, curve, pamLevel, metricName, style, label)
if isempty(curve) || ~ismember(metricName, string(curve.Properties.VariableNames))
return
end
rows = curve(curve.pam_level == pamLevel, :);
if isempty(rows)
return
end
rows = sortrows(rows, "baudrate_GBd");
plot(ax, rows.baudrate_GBd, rows.(metricName), style, ...
"LineWidth", 1.3, "MarkerSize", 5, "DisplayName", label);
end
function mask = equalizerMask(values, target)
if isa(values, "equalizer_structure")
mask = values == target;
elseif isnumeric(values)
mask = double(values) == double(target);
else
valuesString = string(values);
mask = valuesString == string(target) | ...
str2double(valuesString) == double(target);
end
mask = mask(:);
end

View File

@@ -0,0 +1,307 @@
%% Minimal BERp plot for MLSE postfilter orders and duobinary modes
clear; clc;
%% Configuration
warehouseFile = "C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Diss\400G_revisit\mlse_n_tap_pam4.mat";
pfValues = [1, 2, 3];
%% Load warehouse and query run metadata
loadedData = load(warehouseFile, "wh");
wh = loadedData.wh;
wh.showInfo;
db = DBHandler("dataBase", "labor_highspeed", ...
"type", "mysql", ...
"server", "192.168.178.192", ...
"user", "silas", ...
"password", "silas");
warehouseRunIds = getWarehouseRunIds(wh, "mlse_package");
runTable = queryRunsById(db, warehouseRunIds);
if isempty(runTable)
error("plot_mlse_n_tap:NoRunOverlap", ...
"The DB returned no metadata for the warehouse run_ids %s.", ...
mat2str(warehouseRunIds));
end
%% Extract minimum BERp for every run_id and pf_ncoeffs
plotData = warehouseToTable(wh, runTable, pfValues);
hasBerp = isfinite(plotData.BERp) & plotData.BERp > 0;
hasBer = isfinite(plotData.BER) & plotData.BER > 0;
if ~any(hasBerp | hasBer)
error("plot_mlse_n_tap:NoBerp", ...
"No positive BER or BER_precoded values were found in mlse_package.");
end
berpData = plotData(hasBerp, :);
berpData = groupsummary(berpData, ...
["db_mode", "pf_ncoeffs", "bitrate_Gbps"], "min", "BERp");
berData = plotData(hasBer, :);
berData = groupsummary(berData, ...
["db_mode", "pf_ncoeffs", "bitrate_Gbps"], "min", "BER");
vnleData = vnleBerTable(wh, runTable);
vnleBerData = vnleData(isfinite(vnleData.BER) & vnleData.BER > 0, :);
if ~isempty(vnleBerData)
vnleBerData = groupsummary(vnleBerData, ...
["db_mode", "bitrate_Gbps"], "min", "BER");
end
vnleBerpData = vnleData(isfinite(vnleData.BERp) & vnleData.BERp > 0, :);
if ~isempty(vnleBerpData)
vnleBerpData = groupsummary(vnleBerpData, ...
["db_mode", "bitrate_Gbps"], "min", "BERp");
end
%% Plot db_mode = 0 and db_mode = 1 in separate panels
fig = figure(461); clf;
tiledlayout(1, 2, "TileSpacing", "compact", "Padding", "compact");
markers = ["o", "square", "diamond"];
for dbMode = [1, 0]
ax = nexttile; hold(ax, "on");
colors = modeColors(dbMode);
for pfIdx = 1:numel(pfValues)
pf = pfValues(pfIdx);
berpMask = berpData.db_mode == dbMode & ...
berpData.pf_ncoeffs == pf;
if any(berpMask)
curveData = sortrows(berpData(berpMask, :), "bitrate_Gbps");
plot(ax, curveData.bitrate_Gbps, curveData.min_BERp, ...
"LineStyle", "--", ...
"LineWidth", 1.3, ...
"Marker", markers(pfIdx), ...
"MarkerSize", 5, ...
"Color", colors(pfIdx, :), ...
"DisplayName", sprintf("L = %d, precoded",pf));
end
berMask = berData.db_mode == dbMode & ...
berData.pf_ncoeffs == pf;
if any(berMask)
curveData = sortrows(berData(berMask, :), "bitrate_Gbps");
plot(ax, curveData.bitrate_Gbps, curveData.min_BER, ...
"LineStyle", "-", ...
"LineWidth", 1.3, ...
"Marker", markers(pfIdx), ...
"MarkerSize", 5, ...
"Color", colors(pfIdx, :), ...
"DisplayName", sprintf("L = %d",pf));
end
end
vnleBerMask = vnleBerData.db_mode == dbMode;
if any(vnleBerMask)
curveData = sortrows(vnleBerData(vnleBerMask, :), "bitrate_Gbps");
plot(ax, curveData.bitrate_Gbps, curveData.min_BER, ...
"LineStyle", "-", ...
"LineWidth", 1.5, ...
"Color", clr.Set1.gray, ...
"DisplayName", sprintf("VNLE"));
end
vnleBerpMask = vnleBerpData.db_mode == dbMode;
if any(vnleBerpMask)
curveData = sortrows(vnleBerpData(vnleBerpMask, :), "bitrate_Gbps");
plot(ax, curveData.bitrate_Gbps, curveData.min_BERp, ...
"LineStyle", "--", ...
"LineWidth", 1.5, ...
"Color", clr.Set1.gray, ...
"DisplayName", sprintf("VNLE precoded"));
end
yline(ax, [2.2e-4, 4.85e-3, 2e-2], ...
"LineStyle", "--", ...
"LineWidth", 0.8, ...
"Color", [0.25, 0.25, 0.25], ...
"HandleVisibility", "off");
title(ax, sprintf("db\\_mode = %d", dbMode));
xlabel(ax, "Gross rate [Gb/s]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
ylim(ax, [3e-4, 0.1]);
xlim([360, 410])
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
beautifyBERplot("setcolors",false,"setmarkers",false);
end
set(fig, "Position", [100, 500, 1050, 380]);
%% Local helpers
function colors = modeColors(dbMode)
switch dbMode
case 0
colors = [clr.Paired.dred; clr.Paired.dblue; clr.Paired.dgreen];
case 1
colors = [clr.Paired.dred; clr.Paired.dblue; clr.Paired.dgreen];
otherwise
error("plot_mlse_n_tap:UnknownDbMode", ...
"Unsupported db_mode: %d.", dbMode);
end
end
function runIds = getWarehouseRunIds(wh, storageName)
runIds = zeros(0, 1);
if ~isfield(wh.sto, storageName)
return
end
storageValues = wh.sto.(storageName);
for linIdx = 1:numel(storageValues)
[phys, ~] = wh.getPhysAndValueByLinIndex(storageName, linIdx);
if isfield(phys, "run_id")
runIds(end+1, 1) = double(phys.run_id); %#ok<AGROW>
end
end
runIds = unique(runIds);
end
function runTable = queryRunsById(db, runIds)
runTable = table();
fields = db.getTableFieldNames('Runs');
for runIdx = 1:numel(runIds)
fp = QueryFilter();
fp.where('Runs', 'run_id', 'EQUALS', runIds(runIdx));
[oneRun, ~] = db.queryDB(fp, fields);
if isempty(oneRun)
warning("plot_mlse_n_tap:MissingRunMetadata", ...
"No DB metadata found for warehouse run_id %d.", runIds(runIdx));
continue
end
if isempty(runTable)
runTable = oneRun;
else
runTable = [runTable; oneRun]; %#ok<AGROW>
end
end
if ~isempty(runTable) && ismember("db_mode", string(runTable.Properties.VariableNames))
runTable = runTable(double(runTable.db_mode) < 2, :);
end
end
function plotData = warehouseToTable(wh, runTable, pfValues)
runIdCol = zeros(0, 1);
dbModeCol = zeros(0, 1);
pfCol = zeros(0, 1);
bitrateCol = zeros(0, 1);
berCol = zeros(0, 1);
berpCol = zeros(0, 1);
storageValues = wh.sto.mlse_package;
for linIdx = 1:numel(storageValues)
[phys, storedValue] = wh.getPhysAndValueByLinIndex("mlse_package", linIdx);
if ~isfield(phys, "run_id") || ~isfield(phys, "pf_ncoeffs")
continue
end
pf = double(phys.pf_ncoeffs);
if ~ismember(pf, pfValues)
continue
end
runId = double(phys.run_id);
rowIdx = find(double(runTable.run_id) == runId, 1, "first");
if isempty(rowIdx)
continue
end
ber = minMetricValue(storedValue, "BER");
berp = minMetricValue(storedValue, "BER_precoded");
if (~isfinite(ber) || ber <= 0) && ...
(~isfinite(berp) || berp <= 0)
continue
end
runIdCol(end+1, 1) = runId; %#ok<AGROW>
dbModeCol(end+1, 1) = double(runTable.db_mode(rowIdx)); %#ok<AGROW>
pfCol(end+1, 1) = pf; %#ok<AGROW>
bitrateCol(end+1, 1) = double(runTable.bitrate(rowIdx)) .* 1e-9; %#ok<AGROW>
berCol(end+1, 1) = ber; %#ok<AGROW>
berpCol(end+1, 1) = berp; %#ok<AGROW>
end
plotData = table(runIdCol, dbModeCol, pfCol, bitrateCol, berCol, berpCol, ...
'VariableNames', ["run_id", "db_mode", "pf_ncoeffs", ...
"bitrate_Gbps", "BER", "BERp"]);
end
function plotData = vnleBerTable(wh, runTable)
plotData = table(zeros(0, 1), zeros(0, 1), zeros(0, 1), zeros(0, 1), zeros(0, 1), ...
'VariableNames', ["run_id", "db_mode", "bitrate_Gbps", "BER", "BERp"]);
if ~isfield(wh.sto, "vnle_package")
return
end
runIdCol = zeros(0, 1);
dbModeCol = zeros(0, 1);
bitrateCol = zeros(0, 1);
berCol = zeros(0, 1);
berpCol = zeros(0, 1);
storageValues = wh.sto.vnle_package;
for linIdx = 1:numel(storageValues)
[phys, storedValue] = wh.getPhysAndValueByLinIndex("vnle_package", linIdx);
if ~isfield(phys, "run_id")
continue
end
runId = double(phys.run_id);
rowIdx = find(double(runTable.run_id) == runId, 1, "first");
if isempty(rowIdx)
continue
end
ber = minMetricValue(storedValue, "BER");
berp = minMetricValue(storedValue, "BER_precoded");
if (~isfinite(ber) || ber <= 0) && ...
(~isfinite(berp) || berp <= 0)
continue
end
runIdCol(end+1, 1) = runId; %#ok<AGROW>
dbModeCol(end+1, 1) = double(runTable.db_mode(rowIdx)); %#ok<AGROW>
bitrateCol(end+1, 1) = double(runTable.bitrate(rowIdx)) .* 1e-9; %#ok<AGROW>
berCol(end+1, 1) = ber; %#ok<AGROW>
berpCol(end+1, 1) = berp; %#ok<AGROW>
end
plotData = table(runIdCol, dbModeCol, bitrateCol, berCol, berpCol, ...
'VariableNames', ["run_id", "db_mode", "bitrate_Gbps", "BER", "BERp"]);
end
function minValue = minMetricValue(value, metricName)
minValue = NaN;
if isempty(value)
return
end
if ~iscell(value)
value = {value};
end
metricValues = NaN(1, numel(value));
for idx = 1:numel(value)
package = value{idx};
if ~isstruct(package) || ~isfield(package, "metrics")
continue
end
metrics = package.metrics;
fieldName = char(metricName);
if isstruct(metrics) && isfield(metrics, fieldName)
metricValues(idx) = metrics.(fieldName);
elseif isobject(metrics) && isprop(metrics, fieldName)
metricValues(idx) = metrics.(fieldName);
end
end
metricValues = metricValues(isfinite(metricValues) & metricValues > 0);
if ~isempty(metricValues)
minValue = min(metricValues);
end
end

View File

@@ -0,0 +1,407 @@
%% Best NGMI, GMI, AIR and FEC rates over symbol rate
% The normal algorithms are reduced to the best result per algorithm,
% PAM format and symbol rate. Duobinary signaling is added as a separate
% algorithm and is restricted to the post-2026 processing results.
%
% The figure has one column per PAM format and one row per metric:
% BER, NGMI, GMI, AIR, NGMI-based SD+HD NDR, and the best BER-based NDR.
% NDR values are calculated from the measured BER/NGMI using
% TransmissionPerformance.
clear; clc;
%% 1) Query data
selectedPamLevels = [4, 6, 8];
selectedFiberLengthKm = 10;
selectedWavelengthNm = 1310;
selectedRopAttenuation = 0;
selectedIsMpi = 0;
normalDbModes = [double(db_mode.no_db), double(db_mode.db_precoded)];
duobinaryDbMode = double(db_mode.db_encoded);
showRawEntries = false;
maxBerForPlot = 0.5;
algoStyles = defaultAlgorithmStyles();
db = DBHandler( ...
"dataBase", "labor_highspeed", ...
"type", "mysql", ...
"server", "192.168.178.192", ...
"user", "silas", ...
"password", "silas");
db.refresh();
fp = QueryFilter();
fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm);
fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm);
fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation);
fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi);
selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime');
selectedFields = appendMissingFields(selectedFields, {'Runs.precomp_amp'});
[rawData, query] = db.queryDB(fp, selectedFields);
disp(query);
fprintf("Fetched %d rows.\n", height(rawData));
%% 2) Clean data and keep the requested algorithms
data = rawData;
numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ...
"symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ...
"rop_attenuation", "precomp_amp", "is_mpi", "numBits", "numBitErr", ...
"BER", "numBitErr_precoded", "BER_precoded", "GMI", "AIR", "NGMI"];
for fieldIdx = 1:numel(numericFields)
fieldName = numericFields(fieldIdx);
if ismember(fieldName, string(data.Properties.VariableNames))
data.(char(fieldName)) = numericColumn(data.(char(fieldName)));
end
end
data = data(ismember(data.pam_level, selectedPamLevels), :);
normalRows = data(ismember(data.db_mode, normalDbModes), :);
normalPlotData = buildNormalMetricRows(normalRows);
duobinaryRows = data(data.db_mode == duobinaryDbMode, :);
if ismember("equalizer_structure", string(duobinaryRows.Properties.VariableNames))
duobinaryRows = duobinaryRows( ...
equalizerMask(duobinaryRows.equalizer_structure, ...
equalizer_structure.db_encoded), :);
end
% Duobinary results before this date are not comparable to the current set.
if ismember("date_of_processing", string(duobinaryRows.Properties.VariableNames))
duobinaryRows.date_of_processing = datetime(string(duobinaryRows.date_of_processing));
duobinaryRows = duobinaryRows(datetime(duobinaryRows.date_of_processing)>datetime("2026-01-01 00:00:00"),:);
else
warning("plot_best_metrics:NoProcessingDate", ...
"date_of_processing was not returned; no duobinary date filtering was applied.");
end
duobinaryPlotData = buildDuobinaryMetricRows(duobinaryRows);
plotData = [normalPlotData; duobinaryPlotData];
if isempty(plotData)
warning("plot_best_metrics:NoRows", ...
"No rows remain after the database and PAM-format filters.");
return
end
plotData = addDerivedMetrics(plotData);
plotData = plotData(plotData.symbolrate_GBd > 0 & ...
isfinite(plotData.symbolrate_GBd), :);
fprintf("Remaining metric rows: %d\n", height(plotData));
disp(groupcounts(plotData, ["algorithm_key", "pam_level", "precode"]));
%% 3) Calculate BER-/NGMI-dependent net rates
tp = TransmissionPerformance;
grossRate = double(plotData.grossrate);
ngmi = double(plotData.NGMI);
ber = double(plotData.BER_plot);
ngmi(~isfinite(ngmi) | ngmi < 0 | ngmi > 1.05) = NaN;
ber(~isfinite(ber) | ber <= 0 | ber > maxBerForPlot) = NaN;
ndr = tp.calculateNetRate(grossRate, "NGMI", ngmi, "BER", ber);
plotData.NDR_SDHD = columnVector(ndr.SDHD.NetRate) .* 1e-9;
berBasedRates = [ ...
columnVector(ndr.STAIR.NetRate); ...
columnVector(ndr.HD.NetRate); ...
columnVector(ndr.KP4.NetRate); ...
columnVector(ndr.KP4_hamming.NetRate); ...
columnVector(ndr.O_FEC.NetRate)];
plotData.NDR_BER_BEST = max(reshape(berBasedRates, height(plotData), []), [], 2, "omitnan") .* 1e-9;
%% 4) Keep the best row for every plotted metric and group
metricDefinitions = struct( ...
"field", {"BER_plot", "NGMI", "GMI", "AIR_Gbps", "NDR_SDHD", "NDR_BER_BEST"}, ...
"label", {"BER", "NGMI", "GMI [bit/sym]", "AIR [Gb/s]", ...
"SD+HD NDR [Gb/s]", "Best BER-FEC NDR [Gb/s]"}, ...
"scale", {"log", "linear", "linear", "linear", "linear", "linear"});
bestMetricData = cell(numel(metricDefinitions), 1);
for metricIdx = 1:numel(metricDefinitions)
bestMetricData{metricIdx} = bestMetricRows(plotData, ...
metricDefinitions(metricIdx).field);
end
availableStyles = algoStyles(hasAlgorithmRows(plotData, algoStyles), :);
if isempty(availableStyles)
warning("plot_best_metrics:NoSelectedAlgorithms", ...
"None of the configured algorithm styles match the queried rows.");
return
end
%% 5) Plot one figure: metric rows x PAM columns
fig = figure(433); clf;
t = tiledlayout(fig, numel(metricDefinitions), numel(selectedPamLevels), ...
"TileSpacing", "compact", ...
"Padding", "compact");
for metricIdx = 1:numel(metricDefinitions)
metric = metricDefinitions(metricIdx);
metricData = bestMetricData{metricIdx};
for pamIdx = 1:numel(selectedPamLevels)
selectedPamLevel = selectedPamLevels(pamIdx);
ax = nexttile(t); hold(ax, "on");
pamMask = metricData.pam_level == selectedPamLevel;
for styleIdx = 1:height(availableStyles)
style = availableStyles(styleIdx, :);
rowMask = pamMask & metricData.algorithm_key == style.algorithm_key;
if ~any(rowMask)
continue
end
algoData = sortrows(metricData(rowMask, :), "symbolrate_GBd");
y = algoData.(metric.field);
valid = isfinite(y);
if ~any(valid)
continue
end
if showRawEntries
scatter(ax, algoData.symbolrate_GBd(valid), y(valid), ...
9, ...
"Marker", ".", ...
"MarkerEdgeColor", style.color, ...
"MarkerFaceColor", style.color, ...
"MarkerEdgeAlpha", 0.25, ...
"MarkerFaceAlpha", 0.25, ...
"HandleVisibility", "off");
end
plot(ax, algoData.symbolrate_GBd(valid), y(valid), ...
"LineStyle", style.lineStyle, ...
"Marker", style.marker, ...
"MarkerSize", 4, ...
"LineWidth", 1.35, ...
"Color", style.color, ...
"MarkerFaceColor", style.markerFaceColor, ...
"MarkerEdgeColor", style.color, ...
"DisplayName", style.name);
end
formatMetricAxis(ax, metric, selectedPamLevel);
if metricIdx == 1 && pamIdx == 1
legend(ax, "Location", "southwest", "Interpreter", "none");
end
end
end
title(t, sprintf("Best information metrics, %.0f km, %.0f nm", ...
selectedFiberLengthKm, selectedWavelengthNm));
set(fig, "Position", 1e3 .* [0.08 0.04 1.18 0.86]);
%% Local helpers
function fields = appendMissingFields(fields, extraFields)
fields = cellstr(fields);
extraFields = cellstr(extraFields);
for idx = 1:numel(extraFields)
if ~any(strcmp(fields, extraFields{idx}))
fields{end+1, 1} = extraFields{idx}; %#ok<AGROW>
end
end
end
function values = numericColumn(values)
if iscell(values)
values = string(values);
end
if isstring(values) || ischar(values)
values = str2double(values);
end
values = double(values);
end
function values = columnVector(values)
values = double(values(:));
end
function styles = defaultAlgorithmStyles()
styles = table( ...
["vnle"; "vnle_pf_mlse"; "vnle_db_mlse"; "ml_mlse"; "db_encoded"], ...
["VNLE"; "VNLE + PF + MLSE"; "VNLE DBt. + MLSE"; ...
"ML pre-EQ + Viterbi"; "DBS + VNLE + MLSE"], ...
["o"; "square"; "diamond"; "^"; "v"], ...
["-"; "-"; "-"; "-"; "-"], ...
["w"; "w"; "w"; "w"; "w"], ...
[clr.Paired.red; clr.Paired.green; clr.Paired.blue; ...
clr.Paired.purple; clr.Paired.orange], ...
'VariableNames', ["algorithm_key", "name", "marker", ...
"lineStyle", "markerFaceColor", "color"]);
end
function plotData = buildNormalMetricRows(data)
baseRows = data(isfinite(data.BER), :);
baseRows.precode = false(height(baseRows), 1);
baseRows.BER_plot = baseRows.BER;
baseRows.algorithm_key = algorithmKeyFromEqualizer(baseRows.equalizer_structure);
baseRows = baseRows(baseRows.algorithm_key ~= "", :);
if ismember("BER_precoded", string(data.Properties.VariableNames))
precodedRows = data(isfinite(data.BER_precoded), :);
precodedRows.precode = true(height(precodedRows), 1);
precodedRows.BER_plot = precodedRows.BER_precoded;
precodedRows.algorithm_key = algorithmKeyFromEqualizer( ...
precodedRows.equalizer_structure);
precodedRows = precodedRows(precodedRows.algorithm_key ~= "", :);
plotData = [baseRows; precodedRows];
else
warning("plot_best_metrics:NoPrecodedBer", ...
"BER_precoded was not returned; plotting only BER rows.");
plotData = baseRows;
end
end
function plotData = buildDuobinaryMetricRows(data)
plotData = data(isfinite(data.BER), :);
plotData.precode = false(height(plotData), 1);
plotData.BER_plot = plotData.BER;
plotData.algorithm_key = repmat("db_encoded", height(plotData), 1);
end
function plotData = addDerivedMetrics(plotData)
plotData.symbolrate_GBd = double(plotData.symbolrate) .* 1e-9;
plotData.grossrate_Gbps = double(plotData.grossrate) .* 1e-9;
% AIR is stored in bit/s. Reconstruct it from GMI when the stored value is
% missing or outside the physically meaningful range.
plotData.AIR_Gbps = numericOrNaN(plotData, "AIR") .* 1e-9;
gmi = numericOrNaN(plotData, "GMI");
symbolrate = double(plotData.symbolrate);
grossrate = double(plotData.grossrate);
fallbackAir = gmi .* symbolrate .* 1e-9;
useFallback = ~isfinite(plotData.AIR_Gbps) | plotData.AIR_Gbps < 0 | ...
(isfinite(grossrate) & plotData.AIR_Gbps > grossrate .* 1.05e-9);
plotData.AIR_Gbps(useFallback) = fallbackAir(useFallback);
plotData.GMI = gmi;
plotData.NGMI = numericOrNaN(plotData, "NGMI");
end
function values = numericOrNaN(data, fieldName)
if ismember(fieldName, string(data.Properties.VariableNames))
values = numericColumn(data.(char(fieldName)));
else
values = NaN(height(data), 1);
end
values = values(:);
end
function algorithmKey = algorithmKeyFromEqualizer(equalizerColumn)
if isnumeric(equalizerColumn) || islogical(equalizerColumn)
eqNumeric = double(equalizerColumn);
algorithmKey = strings(size(eqNumeric));
knownKeys = ["vnle", "vnle_pf_mlse", "vnle_db_mlse", "ml_mlse"];
knownValues = [double(equalizer_structure.vnle), ...
double(equalizer_structure.vnle_pf_mlse), ...
double(equalizer_structure.vnle_db_mlse), ...
double(equalizer_structure.ml_mlse)];
for idx = 1:numel(knownKeys)
algorithmKey(eqNumeric == knownValues(idx)) = knownKeys(idx);
end
return
end
algorithmKey = lower(string(equalizerColumn));
algorithmKey(~ismember(algorithmKey, ...
["vnle", "vnle_pf_mlse", "vnle_db_mlse", "ml_mlse"])) = "";
end
function mask = equalizerMask(equalizerColumn, eqValue)
if isnumeric(equalizerColumn) || islogical(equalizerColumn)
mask = double(equalizerColumn) == double(eqValue);
else
mask = lower(string(equalizerColumn)) == lower(string(eqValue));
end
end
function bestData = bestMetricRows(data, metricField)
valid = isfinite(data.(metricField));
if strcmp(metricField, "BER_plot")
valid = valid & data.(metricField) > 0;
elseif strcmp(metricField, "NGMI")
valid = valid & data.(metricField) >= 0 & data.(metricField) <= 1.05;
else
valid = valid & data.(metricField) >= 0;
end
candidateData = data(valid, :);
if isempty(candidateData)
bestData = candidateData;
return
end
groupVars = ["pam_level", "algorithm_key", "symbolrate_GBd"];
groupId = findgroups(candidateData(:, groupVars));
keepIdx = zeros(max(groupId), 1);
for curGroup = 1:max(groupId)
rowIdx = find(groupId == curGroup);
values = candidateData.(metricField)(rowIdx);
if strcmp(metricField, "BER_plot")
[~, localIdx] = min(values);
else
[~, localIdx] = max(values);
end
keepIdx(curGroup) = rowIdx(localIdx(1));
end
bestData = sortrows(candidateData(keepIdx, :), groupVars);
end
function keep = hasAlgorithmRows(data, algoStyles)
keep = false(height(algoStyles), 1);
for idx = 1:height(algoStyles)
keep(idx) = any(data.algorithm_key == algoStyles.algorithm_key(idx));
end
end
function formatMetricAxis(ax, metric, pamLevel)
set(ax, "FontSize", 8, "TickLabelInterpreter", "none");
xlabel(ax, "Symbol rate [GBd]");
ylabel(ax, metric.label);
xlim(ax, [95 245]);
xticks(ax, 100:20:240);
grid(ax, "on");
grid(ax, "minor");
box(ax, "on");
switch metric.field
case "BER_plot"
set(ax, "YScale", "log");
ylim(ax, [1e-5 0.2]);
yline(ax, [2.2e-4 4.85e-3 2e-2], ...
"LineWidth", 0.8, "LineStyle", ":", ...
"Color", [0.25 0.25 0.25], "HandleVisibility", "off");
case "NGMI"
ylim(ax, [0.45 1.02]);
case "GMI"
set(ax, "YScale", "linear");
ylim(ax, [0 max(3.2, log2(pamLevel) + 0.15)]);
otherwise
set(ax, "YScale", "linear");
ylim(ax, [0 500]);
yline(ax, 400, "LineWidth", 0.8, "LineStyle", "--", ...
"Color", [0.25 0.25 0.25], "HandleVisibility", "off");
end
if pamLevel == 4
title(ax, "PAM-4");
elseif pamLevel == 6
title(ax, "PAM-6");
elseif pamLevel == 8
title(ax, "PAM-8");
else
title(ax, sprintf("PAM-%d", pamLevel));
end
end

View File

@@ -0,0 +1,471 @@
%% Analysis plots for the reprocessed dsp_400g_recipe warehouse
% Source warehouse is produced by RUN_REPROCESS_BAUDRATE_FILES_DSP400G.m.
clear; clc;
%% 1) Configuration and warehouse loading
warehouseFile = "D:\baudrate_sweep_b2b\PAMX_b2b_baudrate20241024_210648_dsp400g_reprocessed_wh.mat";
selectedRopAttenForBaudrate = 0;
selectedEqForRopCurves = "mlse"; % "ffe", "mlse", or "db"
acquisitionAggregation = "min"; % "min", "median", or "mean"
fecThreshold = 1e-2;
polyfitOrderMax = 4;
showRateAsDatarate = false;
showRawRopMarkers = true;
showPolynomialFits = true;
loadedData = load(warehouseFile);
if isfield(loadedData, "wh")
wh = loadedData.wh;
elseif isfield(loadedData, "obj")
wh = loadedData.obj;
else
error("plot_reprocessed_wh:NoWarehouse", ...
"Warehouse file must contain a variable named wh or obj.");
end
wh.showInfo;
eqStyles = defaultEqStyles();
availableEqStyles = eqStyles(hasWarehouseStorage(wh, eqStyles.storage), :);
if isempty(availableEqStyles)
error("plot_reprocessed_wh:NoEqStorage", ...
"None of the configured BER storages are present in wh.sto.");
end
rawData = warehouseToTable(wh, availableEqStyles);
rawData = rawData(isfinite(rawData.ber) & rawData.ber > 0, :);
allData = aggregateAcquisitions(rawData, acquisitionAggregation);
[allData.rop_axis, ropAxisLabel] = deriveRopAxis(allData);
fprintf("Loaded %d finite acquisition BER rows from warehouse.\n", height(rawData));
fprintf("Aggregated to %d grouped BER rows using acquisitionAggregation = %s.\n", ...
height(allData), acquisitionAggregation);
disp(groupcounts(allData, ["M", "eq"]));
pamVals = sort(unique(allData.M).');
fsymVals = sort(unique(allData.fsym).');
ropAttenVals = sort(unique(allData.rop_atten).');
ropAxisVals = sort(unique(allData.rop_axis(isfinite(allData.rop_axis))).');
%% 2) BER versus baud rate at fixed ROP attenuation
fig = figure(450); clf;
tiledlayout(1, numel(pamVals), ...
"TileSpacing", "compact", ...
"Padding", "compact");
for pamIdx = 1:numel(pamVals)
pamLevel = pamVals(pamIdx);
ax = nexttile; hold(ax, "on");
for eqIdx = 1:height(availableEqStyles)
eqStyle = availableEqStyles(eqIdx, :);
rowMask = allData.M == pamLevel & ...
allData.eq == eqStyle.eq & ...
allData.rop_atten == selectedRopAttenForBaudrate;
if ~any(rowMask)
continue
end
curData = sortrows(allData(rowMask, :), "fsym_GBd");
plot(ax, curData.fsym_GBd, curData.ber, ...
"LineStyle", eqStyle.lineStyle, ...
"Marker", eqStyle.marker, ...
"MarkerSize", 5, ...
"LineWidth", 1.4, ...
"Color", eqStyle.color, ...
"MarkerFaceColor", "w", ...
"MarkerEdgeColor", eqStyle.color, ...
"DisplayName", eqStyle.name);
end
plotFecLines(ax, fecThreshold);
title(ax, sprintf("PAM-%d, ROP atten. %.1f dB", ...
pamLevel, selectedRopAttenForBaudrate));
xlabel(ax, "Symbol rate [GBd]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
ylim(ax, [1e-5, 0.5]);
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
applyBerStyle();
end
set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]);
%% 3) ROP attenuation curves for one EQ scheme
selectedRopEqStyle = availableEqStyles(availableEqStyles.eq == selectedEqForRopCurves, :);
if isempty(selectedRopEqStyle)
error("plot_reprocessed_wh:MissingSelectedEq", ...
"selectedEqForRopCurves = %s is not available in this warehouse.", ...
selectedEqForRopCurves);
end
fig = figure(451); clf;
tiledlayout(1, numel(pamVals), ...
"TileSpacing", "compact", ...
"Padding", "compact");
rateColors = rateColorMap(numel(fsymVals));
for pamIdx = 1:numel(pamVals)
pamLevel = pamVals(pamIdx);
ax = nexttile; hold(ax, "on");
for fsymIdx = 1:numel(fsymVals)
fsym = fsymVals(fsymIdx);
rowMask = allData.M == pamLevel & ...
allData.eq == selectedRopEqStyle.eq & ...
allData.fsym == fsym;
if ~any(rowMask)
continue
end
curData = sortrows(allData(rowMask, :), "rop_axis");
curveColor = rateColors(fsymIdx, :);
displayName = sprintf("%.0f GBd", fsym .* 1e-9);
if showRawRopMarkers
plot(ax, curData.rop_axis, curData.ber, ...
"LineStyle", "none", ...
"Marker", "o", ...
"MarkerSize", 3, ...
"LineWidth", 0.8, ...
"Color", curveColor, ...
"MarkerFaceColor", "w", ...
"MarkerEdgeColor", curveColor, ...
"DisplayName", displayName);
end
if showPolynomialFits
[xFit, yFit] = fitLogBerCurve(curData.rop_axis, ...
curData.ber, polyfitOrderMax);
if ~isempty(xFit)
plot(ax, xFit, yFit, ...
"LineStyle", "-", ...
"LineWidth", 1.1, ...
"Color", curveColor, ...
"HandleVisibility", "off");
end
end
end
title(ax, sprintf("PAM-%d, %s", pamLevel, selectedRopEqStyle.name));
xlabel(ax, ropAxisLabel);
ylabel(ax, "BER");
set(ax, "YScale", "log");
ylim(ax, [1e-5, 0.5]);
xlim(ax, [min(ropAxisVals), max(ropAxisVals)]);
plotFecLines(ax, fecThreshold);
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
applyBerStyle();
end
set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]);
%% 4) Required ROP attenuation at FEC threshold versus baud rate
fecData = computeFecCrossings(allData, availableEqStyles, fecThreshold, ...
polyfitOrderMax);
if isempty(fecData)
warning("plot_reprocessed_wh:NoFecCrossings", ...
"No FEC crossings were found for threshold %.3g.", fecThreshold);
else
fig = figure(452); clf;
tiledlayout(1, numel(pamVals), ...
"TileSpacing", "compact", ...
"Padding", "compact");
for pamIdx = 1:numel(pamVals)
pamLevel = pamVals(pamIdx);
ax = nexttile; hold(ax, "on");
for eqIdx = 1:height(availableEqStyles)
eqStyle = availableEqStyles(eqIdx, :);
rowMask = fecData.M == pamLevel & fecData.eq == eqStyle.eq;
if ~any(rowMask)
continue
end
curData = sortrows(fecData(rowMask, :), "fsym_GBd");
if showRateAsDatarate
xData = curData.datarate_Gbps;
xLabelText = "Datarate [Gb/s]";
else
xData = curData.fsym_GBd;
xLabelText = "Symbol rate [GBd]";
end
plot(ax, xData, curData.required_rop_axis, ...
"LineStyle", eqStyle.lineStyle, ...
"Marker", eqStyle.marker, ...
"MarkerSize", 5, ...
"LineWidth", 1.4, ...
"Color", eqStyle.color, ...
"MarkerFaceColor", "w", ...
"MarkerEdgeColor", eqStyle.color, ...
"DisplayName", eqStyle.name);
end
title(ax, sprintf("PAM-%d, BER = %.2g", pamLevel, fecThreshold));
xlabel(ax, xLabelText);
ylabel(ax, "Required " + ropAxisLabel);
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
end
set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]);
end
%% Local helpers
function styles = defaultEqStyles()
styles = table( ...
["ffe"; "mlse"; "db"], ...
["ber_ffe"; "ber_mlse"; "ber_db"], ...
["FFE"; "VNLE + PF + MLSE"; "VNLE DBt. + MLSE"], ...
["o"; "square"; "diamond"], ...
["-"; "-"; "-"], ...
[clr.Paired.red; clr.Paired.green; clr.Paired.blue], ...
'VariableNames', ["eq", "storage", "name", "marker", ...
"lineStyle", "color"]);
end
function keep = hasWarehouseStorage(wh, storageNames)
stoFields = string(fieldnames(wh.sto));
keep = ismember(storageNames, stoFields);
end
function data = warehouseToTable(wh, eqStyles)
lastIdx = wh.getLastLinIndice();
rows = cell(lastIdx * height(eqStyles), 11);
rowIdx = 0;
for eqIdx = 1:height(eqStyles)
eqStyle = eqStyles(eqIdx, :);
for linIdx = 1:lastIdx
[phys, value] = wh.getPhysAndValueByLinIndex(eqStyle.storage, linIdx);
if isempty(value) || ~isnumeric(value) || ~isscalar(value)
continue
end
fsym = double(phys.fsym);
pamLevel = double(phys.M);
ropAtten = double(phys.rop_atten);
acquisitionIdx = double(phys.acquisition_idx);
ropValue = getOptionalStoredScalarByLinIndex(wh, "rop", linIdx);
pdInValue = getOptionalStoredScalarByLinIndex(wh, "pd_in", linIdx);
rowIdx = rowIdx + 1;
rows(rowIdx, :) = { ...
fsym, ...
fsym .* 1e-9, ...
ropAtten, ...
pamLevel, ...
acquisitionIdx, ...
eqStyle.eq, ...
eqStyle.storage, ...
double(value), ...
ropValue, ...
pdInValue, ...
fsym .* floor(log2(pamLevel) * 10) / 10 .* 1e-9};
end
end
rows = rows(1:rowIdx, :);
data = cell2table(rows, 'VariableNames', ...
["fsym", "fsym_GBd", "rop_atten", "M", "acquisition_idx", ...
"eq", "storage", "ber", "rop_dBm", "pd_in_dBm", "datarate_Gbps"]);
data.eq = string(data.eq);
data.storage = string(data.storage);
end
function data = aggregateAcquisitions(rawData, aggregationMode)
groupVars = ["fsym", "fsym_GBd", "rop_atten", "M", "eq", ...
"storage", "rop_dBm", "pd_in_dBm", "datarate_Gbps"];
switch aggregationMode
case "min"
data = groupsummary(rawData, groupVars, "min", "ber");
data.ber = data.min_ber;
data = removevars(data, "min_ber");
case "median"
data = groupsummary(rawData, groupVars, "median", "ber");
data.ber = data.median_ber;
data = removevars(data, "median_ber");
case "mean"
data = groupsummary(rawData, groupVars, "mean", "ber");
data.ber = data.mean_ber;
data = removevars(data, "mean_ber");
otherwise
error("plot_reprocessed_wh:UnknownAggregation", ...
"Unknown acquisitionAggregation: %s", aggregationMode);
end
end
function value = getOptionalStoredScalarByLinIndex(wh, storageName, linIdx)
value = NaN;
if ~isfield(wh.sto, storageName)
return
end
storedValue = wh.sto.(storageName){linIdx};
if isnumeric(storedValue) && isscalar(storedValue)
value = double(storedValue);
end
end
function [ropAxis, ropAxisLabel] = deriveRopAxis(data)
if ismember("rop_dBm", string(data.Properties.VariableNames)) && ...
any(isfinite(data.rop_dBm))
ropAxis = data.rop_dBm;
ropAxisLabel = "ROP [dBm]";
elseif ismember("pd_in_dBm", string(data.Properties.VariableNames)) && ...
any(isfinite(data.pd_in_dBm))
ropAxis = data.pd_in_dBm;
ropAxisLabel = "PD input power [dBm]";
else
ropAxis = data.rop_atten;
ropAxisLabel = "ROP attenuation [dB]";
end
end
function cmap = rateColorMap(numColors)
anchors = [ ...
clr.Paired.lightblue; ...
clr.Paired.blue; ...
clr.Paired.green; ...
clr.Paired.orange; ...
clr.Paired.red; ...
clr.Paired.purple];
if numColors <= size(anchors, 1)
cmap = anchors(1:numColors, :);
return
end
xAnchor = linspace(0, 1, size(anchors, 1));
xQuery = linspace(0, 1, numColors);
cmap = interp1(xAnchor, anchors, xQuery, "linear");
end
function [xFit, yFit] = fitLogBerCurve(x, y, maxOrder)
valid = isfinite(x) & isfinite(y) & y > 0;
x = x(valid);
y = y(valid);
if numel(unique(x)) < 2
xFit = [];
yFit = [];
return
end
[x, orderIdx] = sort(x(:));
y = y(orderIdx);
fitOrder = min(maxOrder, numel(unique(x)) - 1);
coeff = polyfit(x, log10(y), fitOrder);
xFit = linspace(min(x), max(x), 300).';
yFit = 10 .^ polyval(coeff, xFit);
end
function fecData = computeFecCrossings(data, eqStyles, fecThreshold, maxOrder)
pamVals = unique(data.M).';
fsymVals = unique(data.fsym).';
rows = cell(height(eqStyles) * numel(pamVals) * numel(fsymVals), 6);
rowIdx = 0;
for eqIdx = 1:height(eqStyles)
eqStyle = eqStyles(eqIdx, :);
for pamIdx = 1:numel(pamVals)
pamLevel = pamVals(pamIdx);
for fsymIdx = 1:numel(fsymVals)
fsym = fsymVals(fsymIdx);
rowMask = data.eq == eqStyle.eq & ...
data.M == pamLevel & ...
data.fsym == fsym;
if ~any(rowMask)
continue
end
curData = sortrows(data(rowMask, :), "rop_axis");
requiredRopAxis = fecCrossingFromCurve( ...
curData.rop_axis, curData.ber, fecThreshold, maxOrder);
if ~isfinite(requiredRopAxis)
continue
end
rowIdx = rowIdx + 1;
rows(rowIdx, :) = { ...
fsym, ...
fsym .* 1e-9, ...
pamLevel, ...
eqStyle.eq, ...
requiredRopAxis, ...
fsym .* floor(log2(pamLevel) * 10) / 10 .* 1e-9};
end
end
end
if rowIdx == 0
fecData = table();
else
rows = rows(1:rowIdx, :);
fecData = cell2table(rows, 'VariableNames', ...
["fsym", "fsym_GBd", "M", "eq", "required_rop_axis", ...
"datarate_Gbps"]);
fecData.eq = string(fecData.eq);
end
end
function requiredRopAtten = fecCrossingFromCurve( ...
ropAtten, ber, fecThreshold, maxOrder)
[xFit, yFit] = fitLogBerCurve(ropAtten, ber, maxOrder);
requiredRopAtten = NaN;
if isempty(xFit)
return
end
crossingMask = isfinite(yFit) & yFit > 0;
xFit = xFit(crossingMask);
yFit = yFit(crossingMask);
if numel(xFit) < 2
return
end
delta = log10(yFit) - log10(fecThreshold);
crossingIdx = find(delta(1:end-1) .* delta(2:end) <= 0, 1, "first");
if isempty(crossingIdx)
return
end
x1 = xFit(crossingIdx);
x2 = xFit(crossingIdx + 1);
y1 = delta(crossingIdx);
y2 = delta(crossingIdx + 1);
requiredRopAtten = x1 - y1 .* (x2 - x1) ./ (y2 - y1);
end
function plotFecLines(ax, fecThreshold)
xl = xlim(ax);
h = plot(ax, xl, [fecThreshold, fecThreshold], ...
"LineStyle", "--", ...
"LineWidth", 1, ...
"Color", [0.25 0.25 0.25], ...
"HandleVisibility", "off");
h.Annotation.LegendInformation.IconDisplayStyle = "off";
end
function applyBerStyle()
if exist("beautifyBERplot", "file")
beautifyBERplot("logscale", true, "setcolors", false, ...
"setmarkers", false, "changemarkers", false);
end
end

View File

@@ -0,0 +1,377 @@
%% 400G BER over ROP: best normal algorithms plus duobinary signaling
% Normal algorithms are reduced to the best BER per ROP point across
% db_mode 0/1, pre-emphasis on/off, and BER/BER_precoded result variants.
% Duobinary signaling uses db_mode = 2 and only the sequence-detection BER
% stored in the BER field.
clear; clc;
%% 1) Query data
selectedPamLevels = [4, 6, 8];
selectedFiberLengthKm = 1;
selectedWavelengthNm = 1310;
selectedBitrateGbps = 360;
selectedIsMpi = 0; % set [] to use all entries
maxPowerPdIn = []; % set a numeric limit to enable
normalDbModes = [double(db_mode.no_db)];
duobinaryDbMode = double(db_mode.db_encoded);
maxBerForPlot = 0.5;
showRawEntries = false;
showBestLine = true;
algoStyles = defaultAlgorithmStyles();
db = DBHandler( ...
"dataBase", "labor_highspeed", ...
"type", "mysql", ...
"server", "192.168.178.192", ...
"user", "silas", ...
"password", "silas");
db.refresh();
fp = QueryFilter();
fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm);
fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm);
fp.where('Runs', 'bitrate', 'EQUALS', selectedBitrateGbps .* 1e9);
if ~isempty(selectedIsMpi)
fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi);
end
if ~isempty(maxPowerPdIn)
fp.where('Runs', 'power_pd_in', 'LESS_THAN', maxPowerPdIn);
end
selectedFields = [ ...
db.getTableFieldNames('power_state_info'); ...
db.getTableFieldNames('dashboard_ungrouped_alltime')];
selectedFields = appendMissingFields(selectedFields, ...
{'Runs.precomp_amp'; 'Runs.is_mpi'; 'Runs.power_pd_in'});
selectedFields = selectedFields(:);
[rawData, query] = db.queryDB(fp, selectedFields);
disp(query);
fprintf("Fetched %d 400G ROP result rows.\n", height(rawData));
%% 2) Clean data and build the five plotted curves
data = rawData;
numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ...
"symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ...
"rop_attenuation", "precomp_amp", "is_mpi", "power_rop", ...
"power_mzm", "power_pd_in", "voa_atten", "numBits", "numBitErr", ...
"BER", "numBitErr_precoded", "BER_precoded", "STD", "STDrx", ...
"GMI", "AIR", "NGMI", "EVM", "Alpha"];
for fieldIdx = 1:numel(numericFields)
fieldName = numericFields(fieldIdx);
if ismember(fieldName, string(data.Properties.VariableNames))
data.(char(fieldName)) = numericColumn(data.(char(fieldName)));
end
end
data = data(ismember(data.pam_level, selectedPamLevels), :);
if ~ismember("precomp_amp", string(data.Properties.VariableNames))
warning("plot_rop_best_algos:NoPrecompAmp", ...
"Runs.precomp_amp was not returned. Falling back to pre_emphasis = (db_mode == 0).");
data.pre_emphasis = data.db_mode == double(db_mode.no_db);
else
data.pre_emphasis = derivePreEmphasis(data.precomp_amp, data.db_mode);
end
normalRows = data(ismember(data.db_mode, normalDbModes), :);
normalPlotData = buildNormalMetricRows(normalRows);
normalPlotData = normalPlotData(isfinite(normalPlotData.BER_plot) & ...
normalPlotData.BER_plot > 0 & normalPlotData.BER_plot < maxBerForPlot, :);
duobinaryRows = data(data.db_mode == duobinaryDbMode, :);
if ismember("equalizer_structure", string(duobinaryRows.Properties.VariableNames))
duobinaryRows = duobinaryRows( ...
equalizerMask(duobinaryRows.equalizer_structure, ...
equalizer_structure.db_encoded), :);
end
duobinaryPlotData = buildDuobinarySignalingRows(duobinaryRows);
duobinaryPlotData = duobinaryPlotData(isfinite(duobinaryPlotData.BER_plot) & ...
duobinaryPlotData.BER_plot > 0 & ...
duobinaryPlotData.BER_plot < maxBerForPlot, :);
plotData = [normalPlotData; duobinaryPlotData];
if isempty(plotData)
warning("plot_rop_best_algos:NoRows", ...
"No rows remain after length/PAM/rate/BER filtering.");
return
end
[plotData.rop_axis, ropAxisLabel] = deriveRopAxis(plotData);
plotData.rop_axis = round(plotData.rop_axis, 4);
plotData.bitrate_Gbps = plotData.bitrate .* 1e-9;
plotData.grossrate_Gbps = plotData.grossrate .* 1e-9;
plotData = plotData(isfinite(plotData.rop_axis), :);
fprintf("Remaining candidate BER rows: %d\n", height(plotData));
disp(groupcounts(plotData, ["algorithm_key", "db_mode", "pre_emphasis", "precode"]));
bestPlotData = bestBerByAlgorithmAndRop(plotData);
fprintf("Keeping %d best-BER rows across PAM/algorithm/ROP groups.\n", ...
height(bestPlotData));
disp(groupcounts(bestPlotData, "algorithm_key"));
%% 3) Plot one 1x3 figure with five lines per PAM
availableStyles = algoStyles(hasAlgorithmRows(bestPlotData, algoStyles), :);
if isempty(availableStyles)
warning("plot_rop_best_algos:NoSelectedAlgorithms", ...
"None of the configured algorithm styles match the queried rows.");
return
end
fig = figure(432); clf;
t = tiledlayout(fig, 1, numel(selectedPamLevels), ...
"TileSpacing", "compact", ...
"Padding", "compact");
for pamIdx = 1:numel(selectedPamLevels)
selectedPamLevel = selectedPamLevels(pamIdx);
ax = nexttile(t); hold(ax, "on");
pamMask = bestPlotData.pam_level == selectedPamLevel;
for styleIdx = 1:height(availableStyles)
style = availableStyles(styleIdx, :);
rowMask = pamMask & bestPlotData.algorithm_key == style.algorithm_key;
if ~any(rowMask)
continue
end
algoData = sortrows(bestPlotData(rowMask, :), "rop_axis");
if showRawEntries
scatter(ax, algoData.rop_axis, algoData.BER_plot, ...
9, ...
"Marker", ".", ...
"MarkerEdgeColor", style.color, ...
"MarkerFaceColor", style.color, ...
"MarkerEdgeAlpha", 0.25, ...
"MarkerFaceAlpha", 0.25, ...
"HandleVisibility", "off");
end
if showBestLine
plot(ax, algoData.rop_axis, algoData.BER_plot, ...
"LineStyle", style.lineStyle, ...
"Marker", style.marker, ...
"MarkerSize", 5, ...
"LineWidth", 1.5, ...
"Color", style.color, ...
"MarkerFaceColor", style.markerFaceColor, ...
"MarkerEdgeColor", style.color, ...
"DisplayName", style.name);
end
end
yline(ax, [2.2e-4, 4.85e-3, 2e-2], ...
"LineWidth", 1, ...
"LineStyle", "--", ...
"Color", [0.25 0.25 0.25], ...
"HandleVisibility", "off");
xlabel(ax, ropAxisLabel);
ylabel(ax, "BER");
set(ax, "YScale", "log");
ylim(ax, [8e-5, 0.1]);
grid(ax, "on");
box(ax, "on");
xTicks = unique(bestPlotData.rop_axis(pamMask & ...
isfinite(bestPlotData.rop_axis)));
if ~isempty(xTicks)
if isscalar(xTicks)
xlim(ax, xTicks + [-0.5, 0.5]);
else
xlim(ax, [min(xTicks), max(xTicks)]);
end
end
legend(ax, "Location", "northeast", "Interpreter", "none");
if exist("beautifyBERplot", "file")
beautifyBERplot("logscale", true, "setcolors", false, ...
"setmarkers", false, "changemarkers", false);
end
end
set(fig, "Position", 1e3 .* [0.1070 0.5497 1.0585 0.2282]);
%% Local helpers
function fields = appendMissingFields(fields, extraFields)
fields = cellstr(fields);
extraFields = cellstr(extraFields);
for idx = 1:numel(extraFields)
if ~any(strcmp(fields, extraFields{idx}))
fields{end+1, 1} = extraFields{idx}; %#ok<AGROW>
end
end
end
function values = numericColumn(values)
if iscell(values)
values = string(values);
end
if isstring(values) || ischar(values)
values = str2double(values);
end
values = double(values);
end
function styles = defaultAlgorithmStyles()
styles = table( ...
["vnle"; ...
"vnle_pf_mlse"; ...
"vnle_db_mlse"; ...
"ml_mlse"; ...
"db_encoded"], ...
[equalizer_structure.vnle; ...
equalizer_structure.vnle_pf_mlse; ...
equalizer_structure.vnle_db_mlse; ...
equalizer_structure.ml_mlse; ...
equalizer_structure.db_encoded], ...
["VNLE"; ...
"VNLE + PF + MLSE"; ...
"VNLE DBt. + MLSE"; ...
"ML pre-EQ + Viterbi"; ...
"DBS + VNLE + MLSE"], ...
["o"; "square"; "diamond"; "^"; "v"], ...
["-"; "-"; "-"; "-"; "-"], ...
["w"; "w"; "w"; "w"; "w"], ...
[clr.Paired.red; ...
clr.Paired.green; ...
clr.Paired.blue; ...
clr.Paired.purple; ...
clr.Paired.orange], ...
'VariableNames', ["algorithm_key", "eq", "name", "marker", ...
"lineStyle", "markerFaceColor", "color"]);
end
function preEmphasis = derivePreEmphasis(precompAmp, dbMode)
preEmphasis = false(size(dbMode));
validPrecomp = isfinite(precompAmp);
preEmphasis(validPrecomp) = precompAmp(validPrecomp) > -45;
missingPrecomp = ~validPrecomp;
preEmphasis(missingPrecomp) = dbMode(missingPrecomp) == double(db_mode.no_db);
end
function plotData = buildNormalMetricRows(data)
baseRows = data(isfinite(data.BER), :);
baseRows.precode = false(height(baseRows), 1);
baseRows.BER_plot = baseRows.BER;
baseRows.algorithm_key = algorithmKeyFromEqualizer(baseRows.equalizer_structure);
baseRows = baseRows(baseRows.algorithm_key ~= "", :);
if ismember("BER_precoded", string(data.Properties.VariableNames))
precodedRows = data(isfinite(data.BER_precoded), :);
precodedRows.precode = true(height(precodedRows), 1);
precodedRows.BER_plot = precodedRows.BER_precoded;
precodedRows.algorithm_key = algorithmKeyFromEqualizer( ...
precodedRows.equalizer_structure);
precodedRows = precodedRows(precodedRows.algorithm_key ~= "", :);
plotData = [baseRows; precodedRows];
else
warning("plot_rop_best_algos:NoPrecodedBer", ...
"BER_precoded was not returned. Plotting only BER rows for normal algorithms.");
plotData = baseRows;
end
end
function plotData = buildDuobinarySignalingRows(data)
plotData = data(isfinite(data.BER), :);
plotData.precode = false(height(plotData), 1);
plotData.BER_plot = plotData.BER;
plotData.algorithm_key = repmat("db_encoded", height(plotData), 1);
end
function [ropAxis, label] = deriveRopAxis(data)
if ismember("power_mzm", string(data.Properties.VariableNames)) && ...
any(isfinite(data.power_mzm))
ropAxis = data.power_mzm;
label = "ROP [dBm]";
elseif ismember("power_pd_in", string(data.Properties.VariableNames)) && ...
any(isfinite(data.power_pd_in))
ropAxis = data.power_pd_in;
label = "PD input power [dBm]";
else
ropAxis = data.rop_attenuation;
label = "ROP attenuation [dB]";
end
end
function algorithmKey = algorithmKeyFromEqualizer(equalizerColumn)
eqNumeric = equalizerNumeric(equalizerColumn);
algorithmKey = strings(size(eqNumeric));
algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle)) = "vnle";
algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle_pf_mlse)) = ...
"vnle_pf_mlse";
algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle_db_mlse)) = ...
"vnle_db_mlse";
algorithmKey(eqNumeric == enumValue(equalizer_structure.ml_mlse)) = "ml_mlse";
end
function mask = equalizerMask(equalizerColumn, eqValue)
eqNumeric = equalizerNumeric(equalizerColumn);
mask = eqNumeric == enumValue(eqValue);
end
function eqNumeric = equalizerNumeric(equalizerColumn)
if isa(equalizerColumn, "equalizer_structure")
eqNumeric = double(equalizerColumn);
elseif isnumeric(equalizerColumn)
eqNumeric = double(equalizerColumn);
else
equalizerString = string(equalizerColumn);
eqNumeric = str2double(equalizerString);
enumNames = ["vnle", "ffe", "dfe", "vnle_pf_mlse", ...
"vnle_db_mlse", "db_encoded", "ml_mlse"];
enumValues = [ ...
enumValue(equalizer_structure.vnle), ...
enumValue(equalizer_structure.ffe), ...
enumValue(equalizer_structure.dfe), ...
enumValue(equalizer_structure.vnle_pf_mlse), ...
enumValue(equalizer_structure.vnle_db_mlse), ...
enumValue(equalizer_structure.db_encoded), ...
enumValue(equalizer_structure.ml_mlse)];
for idx = 1:numel(enumNames)
missingNumeric = isnan(eqNumeric);
eqNumeric(missingNumeric & equalizerString == enumNames(idx)) = ...
enumValues(idx);
end
end
end
function value = enumValue(enumEntry)
value = double(enumEntry);
end
function bestData = bestBerByAlgorithmAndRop(data)
groupVars = ["pam_level", "algorithm_key", "rop_axis"];
groupId = findgroups(data(:, groupVars));
keepIdx = NaN(max(groupId), 1);
for curGroup = 1:max(groupId)
rowIdx = find(groupId == curGroup);
[~, localBestIdx] = min(data.BER_plot(rowIdx));
keepIdx(curGroup) = rowIdx(localBestIdx);
end
bestData = sortrows(data(keepIdx, :), groupVars);
end
function keep = hasAlgorithmRows(data, algoStyles)
keep = false(height(algoStyles), 1);
for idx = 1:height(algoStyles)
keep(idx) = any(data.algorithm_key == algoStyles.algorithm_key(idx));
end
end

View File

@@ -0,0 +1,317 @@
%% PAM4/6/8 SNR and uncoded BER versus symbolrate
clear; clc;
warehouseFile = "C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Diss\400G_revisit\results_snr_duobinary_and_partialresponse_pam468_10km.mat";
pamLevels = [4, 6, 8];
storageNames = ["mlse_package", "vnle_package", ...
"dbtgt_package", "mlse_db_package"];
techniqueNames = ["MLSE", "VNLE", "DB target", ...
"duobinary signaling"];
%% Load warehouse and matching run metadata
S = load(warehouseFile, "wh");
wh = S.wh;
db = DBHandler("dataBase", "labor_highspeed", ...
"type", "mysql", ...
"server", "192.168.178.192", ...
"user", "silas", ...
"password", "silas");
fp = QueryFilter();
fp.where('Runs', 'fiber_length', 'EQUALS', 10);
fp.where('Runs', 'wavelength', 'EQUALS', 1310);
fp.where('Runs', 'rop_attenuation', 'EQUALS', 0);
fp.where('Runs', 'is_mpi', 'EQUALS', 0);
[runTable, ~] = db.queryDB(fp, db.getTableFieldNames('Runs'));
runTable = runTable(ismember(double(runTable.pam_level), pamLevels), :);
%% Extract package metrics
techniqueCol = strings(0, 1);
pamCol = zeros(0, 1);
dbModeCol = zeros(0, 1);
symbolrateCol = zeros(0, 1);
snrCol = zeros(0, 1);
berCol = zeros(0, 1);
for techniqueIdx = 1:numel(storageNames)
storageName = storageNames(techniqueIdx);
if ~isfield(wh.sto, storageName)
continue
end
for k = 1:numel(wh.sto.(storageName))
[phys, realizationResults] = wh.getPhysAndValueByLinIndex( ...
storageName, k);
if ~isfield(phys, "run_id") || isempty(realizationResults)
continue
end
runRow = find(double(runTable.run_id) == double(phys.run_id), 1);
if isempty(runRow)
continue
end
for realization = 1:numel(realizationResults)
package = realizationResults{realization};
snr = readMetric(package, "SNR");
ber = readMetric(package, "BER");
if ~isfinite(snr) && ~isfinite(ber)
continue
end
techniqueCol(end+1, 1) = techniqueNames(techniqueIdx); %#ok<AGROW>
pamCol(end+1, 1) = double(runTable.pam_level(runRow)); %#ok<AGROW>
dbModeCol(end+1, 1) = double(runTable.db_mode(runRow)); %#ok<AGROW>
symbolrateCol(end+1, 1) = double(runTable.symbolrate(runRow)) * 1e-9; %#ok<AGROW>
snrCol(end+1, 1) = snr; %#ok<AGROW>
berCol(end+1, 1) = ber; %#ok<AGROW>
end
end
end
data = table(techniqueCol, pamCol, dbModeCol, symbolrateCol, snrCol, berCol, ...
'VariableNames', ["technique", "pam_level", "db_mode", ...
"symbolrate_GBd", "SNR", "BER"]);
% SNR: maximum realization value at each PAM/mode/symbolrate point.
snrData = data(isfinite(data.SNR), :);
snrData = groupsummary(snrData, ...
["technique", "pam_level", "db_mode", "symbolrate_GBd"], ...
"max", "SNR");
snrData.Properties.VariableNames(end) = "SNR";
% BER: minimum positive, uncoded BER at each PAM/mode/symbolrate point.
berData = data(isfinite(data.BER) & data.BER > 0, :);
berData = groupsummary(berData, ...
["technique", "pam_level", "db_mode", "symbolrate_GBd"], ...
"min", "BER");
berData.Properties.VariableNames(end) = "BER";
%% SNR figure: one tile per technique, PAM4/6/8 together
figure(472); clf;
tiledlayout(1, numel(pamLevels), "TileSpacing", "compact", "Padding", "compact");
for pamLevel = pamLevels
ax = nexttile; hold(ax, "on");
plotMetric(ax, snrData, pamLevel, techniqueNames, "SNR");
title(ax, sprintf("PAM-%d", pamLevel));
xlabel(ax, "Symbolrate [GBd]");
ylabel(ax, "SNR [dB]");
ylim(ax, [15, 25]);
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
if exist("beautifyBERplot", "file")
beautifyBERplot("changemarkers", 0, "setcolors", false);
end
end
%% BER figure: one tile per PAM format
figure(473); clf;
tiledlayout(1, numel(pamLevels), "TileSpacing", "compact", "Padding", "compact");
for pamLevel = pamLevels
ax = nexttile; hold(ax, "on");
plotMetric(ax, berData, pamLevel, techniqueNames, "BER");
title(ax, sprintf("PAM-%d", pamLevel));
xlabel(ax, "Symbolrate [GBd]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
end
%% Shared SNR figure: one tile per technique, PAM4/6/8 together
figure(474); clf;
tiledlayout(1, 4, "TileSpacing", "compact", "Padding", "compact");
ax = nexttile; hold(ax, "on");
plotTechniqueSNR(ax, snrData, "VNLE", pamLevels, ...
[0, 1], [clr.Paired.lred; clr.Paired.dred]);
title(ax, "Full Response");
ax = nexttile; hold(ax, "on");
plotTechniqueSNR(ax, snrData, "MLSE", pamLevels, ...
[0, 1], [clr.Paired.lblue; clr.Paired.dblue]);
title(ax, "Partial Response");
ax = nexttile; hold(ax, "on");
plotTechniqueSNR(ax, snrData, "DB target", pamLevels, ...
[0, 1], [clr.Paired.lgreen; clr.Paired.dgreen]);
title(ax, "Duobinary Target");
ax = nexttile; hold(ax, "on");
plotTechniqueSNR(ax, snrData, "duobinary signaling", pamLevels, ...
2, clr.Paired.dlila);
title(ax, "Duobinary Signaling");
%% Shared BER figure: one tile per technique, PAM4/6/8 together
figure(475); clf;
tiledlayout(1, 4, "TileSpacing", "compact", "Padding", "compact");
ax = nexttile; hold(ax, "on");
plotTechniqueBER(ax, berData, "VNLE", pamLevels, ...
[0, 1], [clr.Paired.lred; clr.Paired.dred]);
title(ax, "VNLE");
ax = nexttile; hold(ax, "on");
plotTechniqueBER(ax, berData, "MLSE", pamLevels, ...
[0, 1], [clr.Paired.lblue; clr.Paired.dblue]);
title(ax, "MLSE");
ax = nexttile; hold(ax, "on");
plotTechniqueBER(ax, berData, "DB target", pamLevels, ...
[0, 1], [clr.Paired.lgreen; clr.Paired.dgreen]);
title(ax, "DBt.");
ax = nexttile; hold(ax, "on");
plotTechniqueBER(ax, berData, "duobinary signaling", pamLevels, ...
2, clr.Paired.dlila);
title(ax, "DB signaling");
%% Local helpers
function value = readMetric(package, metricName)
value = NaN;
if isempty(package) || ~isstruct(package) || ~isfield(package, "metrics")
return
end
metrics = package.metrics;
if isstruct(metrics) && isfield(metrics, metricName)
value = double(metrics.(metricName));
elseif isobject(metrics) && isprop(metrics, metricName)
value = double(metrics.(metricName));
end
end
function plotMetric(ax, data, pamLevel, techniqueNames, metricName)
markers = ["o", "s", "diamond", "^"];
lineStyles = ["-", "--", ":"];
for techniqueIdx = 1:numel(techniqueNames)
for dbMode = 0:2
rows = data(data.pam_level == pamLevel & ...
data.technique == techniqueNames(techniqueIdx) & ...
data.db_mode == dbMode, :);
if isempty(rows)
continue
end
rows = sortrows(rows, "symbolrate_GBd");
plot(ax, rows.symbolrate_GBd, rows.(metricName), ...
"LineStyle", lineStyles(dbMode + 1), ...
"Marker", markers(techniqueIdx), ...
"MarkerSize", 5, ...
"LineWidth", 1.3, ...
"Color", pamModeColor(pamLevel, dbMode), ...
"DisplayName", sprintf("%s, %s", ...
techniqueNames(techniqueIdx), modeLabel(dbMode)));
end
end
end
function plotTechniqueSNR(ax, data, techniqueName, pamLevels, dbModes, colors)
markers = ["square", "hexagram", "*"];
lineStyles = ["-", "--", ":"];
for modeIdx = 1:numel(dbModes)
dbMode = dbModes(modeIdx);
for pamIdx = 1:numel(pamLevels)
pamLevel = pamLevels(pamIdx);
rows = data(data.technique == techniqueName & ...
data.pam_level == pamLevel & data.db_mode == dbMode, :);
if isempty(rows)
continue
end
rows = sortrows(rows, "symbolrate_GBd");
plot(ax, rows.symbolrate_GBd, rows.SNR, ...
"LineStyle", lineStyles(dbMode + 1), ...
"Marker", markers(pamIdx), ...
"MarkerSize", 6, ...
"LineWidth", 1.3, ...
"Color", pamModeColor(pamLevel, dbMode), ...
"DisplayName", sprintf("%s, PAM-%d", ...
modeLabel(dbMode), pamLevel));
end
end
xlabel(ax, "Symbolrate [GBd]");
ylabel(ax, "SNR [dB]");
ylim(ax, [13, 24]);
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
end
function plotTechniqueBER(ax, data, techniqueName, pamLevels, dbModes, colors)
markers = ["square", "hexagram", "*"];
lineStyles = ["-", "--", ":"];
for modeIdx = 1:numel(dbModes)
dbMode = dbModes(modeIdx);
for pamIdx = 1:numel(pamLevels)
pamLevel = pamLevels(pamIdx);
rows = data(data.technique == techniqueName & ...
data.pam_level == pamLevel & data.db_mode == dbMode, :);
if isempty(rows)
continue
end
rows = sortrows(rows, "symbolrate_GBd");
plot(ax, rows.symbolrate_GBd, rows.BER, ...
"LineStyle", lineStyles(dbMode + 1), ...
"Marker", markers(pamIdx), ...
"MarkerSize", 6, ...
"LineWidth", 1.3, ...
"Color", pamModeColor(pamLevel, dbMode), ...
"DisplayName", sprintf("%s, PAM-%d", ...
modeLabel(dbMode), pamLevel));
end
end
xlabel(ax, "Symbolrate [GBd]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
ylim(ax, [1e-4, 1e-1]);
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
end
function color = pamModeColor(pamLevel, dbMode)
switch pamLevel
case 4
lightColor = clr.Paired.lgreen;
darkColor = clr.Paired.dgreen;
case 6
lightColor = clr.Paired.lblue;
darkColor = clr.Paired.dblue;
case 8
lightColor = clr.Paired.lred;
darkColor = clr.Paired.dred;
otherwise
error("plot_pam_comparison:UnknownPam", ...
"Unsupported PAM level %d.", pamLevel);
end
if dbMode == 0
color = lightColor;
else
color = darkColor;
end
end
function label = modeLabel(dbMode)
switch dbMode
case 0
label = "with preemphasis";
case 1
label = "no preemphasis";
case 2
label = "db encoded";
otherwise
label = sprintf("db\_mode = %d", dbMode);
end
end

View File

@@ -0,0 +1,440 @@
%% Warehouse baud-rate and ROP analysis for FFE, MLSE, and duobinary target
% Template source:
% projects/Diss/400G_revisit/auswertung_baudrate.mlx
clear; clc;
%% 1) Configuration and warehouse loading
warehouseFile = "W:\labdata\sioe_labor\baudrate_sweep_b2b\PAMX_b2b_baudrate20241024_210648_wh_final.mat";
selectedRopAttenForBaudrate = 0;
selectedEqForRopCurves = "mlse"; % "ffe", "mlse", or "db"
fecThreshold = 2.2e-2;
polyfitOrderMax = 4;
showRateAsDatarate = false;
showRawRopMarkers = true;
showPolynomialFits = true;
loadedData = load(warehouseFile);
if isfield(loadedData, "obj")
wh = loadedData.obj;
elseif isfield(loadedData, "wh")
wh = loadedData.wh;
else
error("plot_warehouse:NoWarehouse", ...
"Warehouse file must contain a variable named obj or wh.");
end
wh.showInfo;
fsymVals = double(wh.parameter.fsym.values(:).');
ropAttenVals = double(wh.parameter.rop_atten.values(:).');
pamVals = sort(double(wh.parameter.M.values(:).'));
eqStyles = defaultEqStyles();
availableEqStyles = eqStyles(hasWarehouseStorage(wh, eqStyles.storage), :);
if isempty(availableEqStyles)
error("plot_warehouse:NoEqStorage", ...
"None of the configured BER storages are present in wh.sto.");
end
allData = warehouseToTable(wh, fsymVals, ropAttenVals, pamVals, availableEqStyles);
allData = allData(isfinite(allData.ber) & allData.ber > 0, :);
fprintf("Loaded %d finite BER rows from warehouse.\n", height(allData));
disp(groupcounts(allData, ["M", "eq"]));
for eqIdx = 1:height(availableEqStyles)
if ~any(allData.eq == availableEqStyles.eq(eqIdx))
warning("plot_warehouse:NoPositiveBer", ...
"Storage %s exists, but contains no positive BER values to plot.", ...
availableEqStyles.storage(eqIdx));
end
end
%% 2) BER versus baud rate at fixed ROP attenuation
fig = figure(440); clf;
tiledlayout(1, numel(pamVals), ...
"TileSpacing", "compact", ...
"Padding", "compact");
for pamIdx = 1:numel(pamVals)
pamLevel = pamVals(pamIdx);
ax = nexttile; hold(ax, "on");
for eqIdx = 1:height(availableEqStyles)
eqStyle = availableEqStyles(eqIdx, :);
rowMask = allData.M == pamLevel & ...
allData.eq == eqStyle.eq & ...
allData.rop_atten == selectedRopAttenForBaudrate;
if ~any(rowMask)
continue
end
curData = sortrows(allData(rowMask, :), "fsym_GBd");
plot(ax, curData.fsym_GBd, curData.ber, ...
"LineStyle", eqStyle.lineStyle, ...
"Marker", eqStyle.marker, ...
"MarkerSize", 5, ...
"LineWidth", 1.4, ...
"Color", eqStyle.color, ...
"MarkerFaceColor", "w", ...
"MarkerEdgeColor", eqStyle.color, ...
"DisplayName", eqStyle.name);
end
yline(ax, fecThreshold, ...
"LineStyle", "--", ...
"LineWidth", 1, ...
"Color", [0.25 0.25 0.25], ...
"HandleVisibility", "off");
title(ax, sprintf("PAM-%d, ROP atten. %.1f dB", ...
pamLevel, selectedRopAttenForBaudrate));
xlabel(ax, "Symbol rate [GBd]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
ylim(ax, [1e-5, 0.5]);
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
applyBerStyle();
end
set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]);
%% 3) ROP attenuation curves for one EQ scheme
selectedRopEqStyle = availableEqStyles(availableEqStyles.eq == selectedEqForRopCurves, :);
if isempty(selectedRopEqStyle)
error("plot_warehouse:MissingSelectedEq", ...
"selectedEqForRopCurves = %s is not available in this warehouse.", ...
selectedEqForRopCurves);
end
fig = figure(441); clf;
tiledlayout(1, numel(pamVals), ...
"TileSpacing", "compact", ...
"Padding", "compact");
rateColors = rateColorMap(numel(fsymVals));
for pamIdx = 1:numel(pamVals)
pamLevel = pamVals(pamIdx);
ax = nexttile; hold(ax, "on");
for fsymIdx = 1:numel(fsymVals)
fsym = fsymVals(fsymIdx);
rowMask = allData.M == pamLevel & ...
allData.eq == selectedRopEqStyle.eq & ...
allData.fsym == fsym;
if ~any(rowMask)
continue
end
curData = sortrows(allData(rowMask, :), "rop_atten");
curveColor = rateColors(fsymIdx, :);
displayName = sprintf("%.0f GBd", fsym .* 1e-9);
if showRawRopMarkers
plot(ax, curData.rop_atten, curData.ber, ...
"LineStyle", "none", ...
"Marker", "o", ...
"MarkerSize", 3, ...
"LineWidth", 0.8, ...
"Color", curveColor, ...
"MarkerFaceColor", "w", ...
"MarkerEdgeColor", curveColor, ...
"DisplayName", displayName);
end
if showPolynomialFits
[xFit, yFit] = fitLogBerCurve(curData.rop_atten, ...
curData.ber, polyfitOrderMax);
if ~isempty(xFit)
plot(ax, xFit, yFit, ...
"LineStyle", "-", ...
"LineWidth", 1.1, ...
"Color", curveColor, ...
"HandleVisibility", "off");
end
end
end
yline(ax, fecThreshold, ...
"LineStyle", "--", ...
"LineWidth", 1, ...
"Color", [0.25 0.25 0.25], ...
"HandleVisibility", "off");
title(ax, sprintf("PAM-%d, %s", pamLevel, selectedRopEqStyle.name));
xlabel(ax, "ROP attenuation [dB]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
ylim(ax, [1e-5, 0.5]);
xlim(ax, [min(ropAttenVals), max(ropAttenVals)]);
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
applyBerStyle();
end
set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]);
%% 4) Required power at FEC threshold versus baud rate
fecData = computeFecCrossings(allData, availableEqStyles, fecThreshold, ...
polyfitOrderMax);
if isempty(fecData)
warning("plot_warehouse:NoFecCrossings", ...
"No FEC crossings were found for threshold %.3g.", fecThreshold);
else
fig = figure(442); clf;
tiledlayout(1, numel(pamVals), ...
"TileSpacing", "compact", ...
"Padding", "compact");
for pamIdx = 1:numel(pamVals)
pamLevel = pamVals(pamIdx);
ax = nexttile; hold(ax, "on");
for eqIdx = 1:height(availableEqStyles)
eqStyle = availableEqStyles(eqIdx, :);
rowMask = fecData.M == pamLevel & fecData.eq == eqStyle.eq;
if ~any(rowMask)
continue
end
curData = sortrows(fecData(rowMask, :), "fsym_GBd");
if showRateAsDatarate
xData = curData.datarate_Gbps;
xLabelText = "Datarate [Gb/s]";
else
xData = curData.fsym_GBd;
xLabelText = "Symbol rate [GBd]";
end
plot(ax, xData, curData.required_rop_dBm, ...
"LineStyle", eqStyle.lineStyle, ...
"Marker", eqStyle.marker, ...
"MarkerSize", 5, ...
"LineWidth", 1.4, ...
"Color", eqStyle.color, ...
"MarkerFaceColor", "w", ...
"MarkerEdgeColor", eqStyle.color, ...
"DisplayName", eqStyle.name);
end
title(ax, sprintf("PAM-%d, BER = %.2g", pamLevel, fecThreshold));
xlabel(ax, xLabelText);
ylabel(ax, "Required ROP [dBm]");
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
end
set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]);
end
%% Local helpers
function styles = defaultEqStyles()
styles = table( ...
["ffe"; "mlse"; "db"], ...
["ber_ffe"; "ber_mlse"; "ber_db"], ...
["FFE"; "MLSE"; "Duobinary target"], ...
["o"; "square"; "diamond"], ...
["-"; "-"; "-"], ...
[clr.Paired.red; clr.Paired.green; clr.Paired.blue], ...
'VariableNames', ["eq", "storage", "name", "marker", ...
"lineStyle", "color"]);
end
function keep = hasWarehouseStorage(wh, storageNames)
stoFields = string(fieldnames(wh.sto));
keep = ismember(storageNames, stoFields);
end
function data = warehouseToTable(wh, fsymVals, ropAttenVals, pamVals, eqStyles)
rows = {};
for eqIdx = 1:height(eqStyles)
eqStyle = eqStyles(eqIdx, :);
for pamIdx = 1:numel(pamVals)
pamLevel = pamVals(pamIdx);
for fsymIdx = 1:numel(fsymVals)
fsym = fsymVals(fsymIdx);
berValues = wh.getStoValue(eqStyle.storage, fsym, ...
ropAttenVals, pamLevel);
ropValues = getOptionalStoValues(wh, "rop", fsym, ...
ropAttenVals, pamLevel);
pdInValues = getOptionalStoValues(wh, "pd_in", fsym, ...
ropAttenVals, pamLevel);
berValues = berValues(:);
ropValues = ropValues(:);
pdInValues = pdInValues(:);
for ropIdx = 1:numel(ropAttenVals)
rows(end+1, :) = { ... %#ok<AGROW>
fsym, ...
fsym .* 1e-9, ...
ropAttenVals(ropIdx), ...
pamLevel, ...
eqStyle.eq, ...
eqStyle.storage, ...
berValues(ropIdx), ...
ropValues(ropIdx), ...
pdInValues(ropIdx), ...
fsym .* floor(log2(pamLevel) * 10) / 10 .* 1e-9};
end
end
end
end
data = cell2table(rows, 'VariableNames', ...
["fsym", "fsym_GBd", "rop_atten", "M", "eq", "storage", ...
"ber", "rop_dBm", "pd_in_dBm", "datarate_Gbps"]);
data.eq = string(data.eq);
data.storage = string(data.storage);
end
function values = getOptionalStoValues(wh, storageName, fsym, ropAttenVals, pamLevel)
if ismember(storageName, string(fieldnames(wh.sto)))
values = wh.getStoValue(storageName, fsym, ropAttenVals, pamLevel);
else
values = nan(size(ropAttenVals));
end
end
function cmap = rateColorMap(numColors)
anchors = [ ...
clr.Paired.lightblue; ...
clr.Paired.blue; ...
clr.Paired.green; ...
clr.Paired.orange; ...
clr.Paired.red; ...
clr.Paired.purple];
if numColors <= size(anchors, 1)
cmap = anchors(1:numColors, :);
return
end
xAnchor = linspace(0, 1, size(anchors, 1));
xQuery = linspace(0, 1, numColors);
cmap = interp1(xAnchor, anchors, xQuery, "linear");
end
function [xFit, yFit] = fitLogBerCurve(x, y, maxOrder)
valid = isfinite(x) & isfinite(y) & y > 0;
x = x(valid);
y = y(valid);
if numel(unique(x)) < 2
xFit = [];
yFit = [];
return
end
[x, orderIdx] = sort(x(:));
y = y(orderIdx);
fitOrder = min(maxOrder, numel(unique(x)) - 1);
coeff = polyfit(x, log10(y), fitOrder);
xFit = linspace(min(x), max(x), 300).';
yFit = 10 .^ polyval(coeff, xFit);
end
function fecData = computeFecCrossings(data, eqStyles, fecThreshold, maxOrder)
rows = {};
pamVals = unique(data.M).';
fsymVals = unique(data.fsym).';
for eqIdx = 1:height(eqStyles)
eqStyle = eqStyles(eqIdx, :);
for pamIdx = 1:numel(pamVals)
pamLevel = pamVals(pamIdx);
for fsymIdx = 1:numel(fsymVals)
fsym = fsymVals(fsymIdx);
rowMask = data.eq == eqStyle.eq & ...
data.M == pamLevel & ...
data.fsym == fsym;
if ~any(rowMask)
continue
end
curData = sortrows(data(rowMask, :), "rop_atten");
[requiredRopAtten, requiredRop] = fecCrossingFromCurve( ...
curData.rop_atten, curData.ber, curData.rop_dBm, ...
fecThreshold, maxOrder);
if ~isfinite(requiredRop)
continue
end
rows(end+1, :) = { ... %#ok<AGROW>
fsym, ...
fsym .* 1e-9, ...
pamLevel, ...
eqStyle.eq, ...
requiredRopAtten, ...
requiredRop, ...
fsym .* floor(log2(pamLevel) * 10) / 10 .* 1e-9};
end
end
end
if isempty(rows)
fecData = table();
else
fecData = cell2table(rows, 'VariableNames', ...
["fsym", "fsym_GBd", "M", "eq", "required_rop_atten", ...
"required_rop_dBm", "datarate_Gbps"]);
fecData.eq = string(fecData.eq);
end
end
function [requiredRopAtten, requiredRop] = fecCrossingFromCurve( ...
ropAtten, ber, rop, fecThreshold, maxOrder)
[xFit, yFit] = fitLogBerCurve(ropAtten, ber, maxOrder);
requiredRopAtten = NaN;
requiredRop = NaN;
if isempty(xFit)
return
end
crossingMask = isfinite(yFit) & yFit > 0;
xFit = xFit(crossingMask);
yFit = yFit(crossingMask);
if numel(xFit) < 2
return
end
delta = log10(yFit) - log10(fecThreshold);
crossingIdx = find(delta(1:end-1) .* delta(2:end) <= 0, 1, "first");
if isempty(crossingIdx)
return
end
x1 = xFit(crossingIdx);
x2 = xFit(crossingIdx + 1);
y1 = delta(crossingIdx);
y2 = delta(crossingIdx + 1);
requiredRopAtten = x1 - y1 .* (x2 - x1) ./ (y2 - y1);
validRop = isfinite(ropAtten) & isfinite(rop);
if nnz(validRop) >= 2
requiredRop = interp1(ropAtten(validRop), rop(validRop), ...
requiredRopAtten, "linear", "extrap");
else
requiredRop = requiredRopAtten;
end
end
function applyBerStyle()
if exist("beautifyBERplot", "file")
beautifyBERplot("logscale", true, "setcolors", false, ...
"setmarkers", false, "changemarkers", false);
end
end

View File

@@ -0,0 +1,691 @@
%% 400G BER over wavelength: best normal algorithms plus duobinary signaling
% Normal algorithms are reduced to the best BER per wavelength across
% db_mode 0/1, pre-emphasis on/off, and BER/BER_precoded result variants.
% Duobinary signaling uses db_mode = 2 and only the sequence-detection BER
% stored in the BER field.
clear; clc;
%% 1) Query data
selectedPamLevels = [4, 6, 8];
selectedFiberLengthKm = 5;
selectedBitratesGbps = 300:30:480; % set [] to use all available bitrates
selectedRopAttenuation = 0; % set [] to use all ROP attenuation values
selectedIsMpi = 0; % set [] to use all entries
normalDbModes = [double(db_mode.no_db), double(db_mode.db_precoded)];
duobinaryDbMode = double(db_mode.db_encoded);
maxBerForPlot = 0.5;
showRawEntries = false;
showBestLine = true;
showPolynomialFits = true;
polyfitOrder = 4;
fecThreshold = 2e-2;
algoStyles = defaultAlgorithmStyles();
areaResults = cell(numel(selectedPamLevels), 1);
for pamIdx = 1:numel(selectedPamLevels)
selectedPamLevel = selectedPamLevels(pamIdx);
db = DBHandler( ...
"dataBase", "labor_highspeed", ...
"type", "mysql", ...
"server", "192.168.178.192", ...
"user", "silas", ...
"password", "silas");
db.refresh();
fp = QueryFilter();
fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm);
fp.where('Runs', 'pam_level', 'EQUALS', selectedPamLevel);
if ~isempty(selectedRopAttenuation)
fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation);
end
if ~isempty(selectedIsMpi)
fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi);
end
selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime');
selectedFields = appendMissingFields(selectedFields, ...
{'Runs.precomp_amp'; 'Runs.is_mpi'});
selectedFields = selectedFields(:);
[rawData, query] = db.queryDB(fp, selectedFields);
disp(query);
fprintf("Fetched %d wavelength-sweep result rows.\n", height(rawData));
%% 2) Clean data and build the five plotted curves
data = rawData;
numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ...
"symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ...
"rop_attenuation", "precomp_amp", "is_mpi", "numBits", "numBitErr", ...
"BER", "numBitErr_precoded", "BER_precoded", "STD", "STDrx", ...
"GMI", "AIR", "NGMI", "EVM", "Alpha"];
for fieldIdx = 1:numel(numericFields)
fieldName = numericFields(fieldIdx);
if ismember(fieldName, string(data.Properties.VariableNames))
data.(char(fieldName)) = numericColumn(data.(char(fieldName)));
end
end
data.bitrate_Gbps = data.bitrate .* 1e-9;
if ~isempty(selectedBitratesGbps)
data = data(ismember(round(data.bitrate_Gbps, 6), selectedBitratesGbps), :);
end
if ~ismember("precomp_amp", string(data.Properties.VariableNames))
warning("plot_wavelength_best_algos:NoPrecompAmp", ...
"Runs.precomp_amp was not returned. Falling back to pre_emphasis = (db_mode == 0).");
data.pre_emphasis = data.db_mode == double(db_mode.no_db);
else
data.pre_emphasis = derivePreEmphasis(data.precomp_amp, data.db_mode);
end
normalRows = data(ismember(data.db_mode, normalDbModes), :);
normalPlotData = buildNormalMetricRows(normalRows);
normalPlotData = normalPlotData(isfinite(normalPlotData.BER_plot) & ...
normalPlotData.BER_plot > 0 & normalPlotData.BER_plot < maxBerForPlot, :);
duobinaryRows = data(data.db_mode == duobinaryDbMode, :);
if ismember("equalizer_structure", string(duobinaryRows.Properties.VariableNames))
duobinaryRows = duobinaryRows( ...
equalizerMask(duobinaryRows.equalizer_structure, ...
equalizer_structure.db_encoded), :);
end
duobinaryPlotData = buildDuobinarySignalingRows(duobinaryRows);
duobinaryPlotData = duobinaryPlotData(isfinite(duobinaryPlotData.BER_plot) & ...
duobinaryPlotData.BER_plot > 0 & ...
duobinaryPlotData.BER_plot < maxBerForPlot, :);
plotData = [normalPlotData; duobinaryPlotData];
if isempty(plotData)
warning("plot_wavelength_best_algos:NoRows", ...
"No rows remain after length/PAM/rate/BER filtering.");
continue
end
plotData.bitrate_Gbps = plotData.bitrate .* 1e-9;
plotData.grossrate_Gbps = plotData.grossrate .* 1e-9;
fprintf("Remaining candidate BER rows: %d\n", height(plotData));
disp(groupcounts(plotData, ["algorithm_key", "db_mode", "pre_emphasis", "precode"]));
bestPlotData = bestBerByAlgorithmAndWavelength(plotData);
fprintf("Keeping %d best-BER rows across algorithm/wavelength groups.\n", ...
height(bestPlotData));
disp(groupcounts(bestPlotData, "algorithm_key"));
%% 3) Plot one tile per algorithm
availableStyles = algoStyles(hasAlgorithmRows(bestPlotData, algoStyles), :);
if isempty(availableStyles)
warning("plot_wavelength_best_algos:NoSelectedAlgorithms", ...
"None of the configured algorithm styles match the queried rows.");
continue
end
fig = figure(432 + pamIdx); clf;
t = tiledlayout(fig, 1, height(availableStyles), ...
"TileSpacing","compact", ...
"Padding", "compact");
availableBitratesGbps = unique(bestPlotData.bitrate_Gbps(isfinite( ...
bestPlotData.bitrate_Gbps))).';
bitrateMarkers = bitrateMarkerSet(numel(availableBitratesGbps));
for styleIdx = 1:height(availableStyles)
style = availableStyles(styleIdx, :);
ax = nexttile(t); hold(ax, "on");
bitrateColors = sequentialColors(style.color, numel(availableBitratesGbps), ...
style.algorithm_key);
for bitrateIdx = 1:numel(availableBitratesGbps)
bitrateGbps = availableBitratesGbps(bitrateIdx);
rowMask = bestPlotData.algorithm_key == style.algorithm_key & ...
bestPlotData.bitrate_Gbps == bitrateGbps;
if ~any(rowMask)
continue
end
algoData = sortrows(bestPlotData(rowMask, :), "wavelength");
curveColor = bitrateColors(bitrateIdx, :);
if showRawEntries
scatter(ax, algoData.wavelength, algoData.BER_plot, ...
9, ...
"Marker", ".", ...
"MarkerEdgeColor", curveColor, ...
"MarkerFaceColor", curveColor, ...
"MarkerEdgeAlpha", 0.25, ...
"MarkerFaceAlpha", 0.25, ...
"HandleVisibility", "off");
end
if showBestLine
plot(ax, algoData.wavelength, algoData.BER_plot, ...
"LineStyle", style.lineStyle, ...
"Marker", bitrateMarkers(bitrateIdx), ...
"MarkerSize", 2.5, ...
"LineWidth", 1.0, ...
"Color", curveColor, ...
"MarkerFaceColor", style.markerFaceColor, ...
"MarkerEdgeColor", curveColor, ...
"HandleVisibility", "off");
end
if showPolynomialFits
[xFit, yFit] = fitLogBerCurve(algoData.wavelength, ...
algoData.BER_plot, polyfitOrder);
if ~isempty(xFit)
plot(ax, xFit, yFit, ...
"LineStyle", ":", ...
"LineWidth", 1.1, ...
"Color", curveColor, ...
"HandleVisibility", "off");
end
end
end
% title(ax, style.name, "Interpreter", "none");
xlabel(ax, "Wavelength [nm]");
ylabel(ax, "BER");
set(ax, "YScale", "log");
if selectedPamLevel == 4
ylim(ax, [1e-6, 0.3]);
elseif selectedPamLevel == 6
ylim(ax, [1e-4, 0.3]);
elseif selectedPamLevel == 8
ylim(ax, [1e-4, 0.3]);
end
grid(ax, "on");
box(ax, "on");
if styleIdx > 1
ylabel(ax, "");
yticklabels(ax, {});
end
xTicks = unique(bestPlotData.wavelength(isfinite(bestPlotData.wavelength)));
xTicks = xTicks(1:2:end);
xTicks_own = [1300,1310,1320];
set(ax, "XTick", xTicks_own);
if ~isempty(xTicks)
xlim(ax, [min(xTicks)-1, max(xTicks)+1]);
end
plotFecLines(ax, fecThreshold);
if exist("beautifyBERplot", "file")
beautifyBERplot("logscale", true, "setcolors", false, ...
"setmarkers", false, "changemarkers", false);
end
text(ax, 0.02, 0.02, style.name, ...
"Units", "normalized", ...
"HorizontalAlignment", "left", ...
"VerticalAlignment", "bottom", ...
"BackgroundColor", "white", ...
"EdgeColor", [0.60 0.60 0.60], ...
"LineWidth", 0.5, ...
"Margin", 2, ...
"FontSize", 8, ...
"Interpreter", "none", ...
"Clipping", "on");
end
set(fig, "Position", 1e3 .* [0.1070 0.5497 1.4113 0.3253]);
wavelengthTikzPath = sprintf( ...
"C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/04_Experimental_Evaluation/tikz/400g/10km_wavelength_compare_pam%d.tikz", ...
selectedPamLevel);
% mat2tikz_improved(wavelengthTikzPath, "cleanfigure", 0);
areaData = computePermissibleWavelengthAreas(bestPlotData, availableStyles, ...
fecThreshold, polyfitOrder);
if isempty(areaData)
warning("plot_wavelength_best_algos:NoPermissibleAreas", ...
"No measured wavelength samples are at or below the BER threshold %.3g.", ...
fecThreshold);
else
areaResults{pamIdx} = areaData;
end
end
%% 4) Permissible wavelength area versus gross rate for all PAM levels
areaFig = figure(436); clf;
areaLayout = tiledlayout(areaFig, 1, 3, ...
"TileSpacing", "compact", ...
"Padding", "compact");
for pamIdx = 1:numel(selectedPamLevels)
areaAx = nexttile(areaLayout); hold(areaAx, "on");
areaData = areaResults{pamIdx};
if isempty(areaData)
grid(areaAx, "on");
box(areaAx, "on");
else
for styleIdx = 1:height(algoStyles)
style = algoStyles(styleIdx, :);
rowMask = areaData.algorithm_key == style.algorithm_key;
if ~any(rowMask)
continue
end
algoAreaData = sortrows(areaData(rowMask, :), "grossrate_Gbps");
plot(areaAx, algoAreaData.grossrate_Gbps, ...
algoAreaData.permissible_wavelength_area_nm, ...
"LineStyle", style.lineStyle, ...
"Marker", style.marker, ...
"MarkerSize", 5, ...
"LineWidth", 1.4, ...
"Color", style.color, ...
"MarkerFaceColor", style.markerFaceColor, ...
"MarkerEdgeColor", style.color, ...
"DisplayName", style.name, ...
"HandleVisibility", "on");
end
end
ylim([5 35]);
xlabel(areaAx, "Gross rate [Gb/s]");
if pamIdx == 1
ylabel(areaAx, "Permissible wavelength area [nm]");
else
ylabel(areaAx, "");
% Ensure yticklabels refers to the function, not a variable
if exist("yticklabels", "var")
clear yticklabels
end
yticks(areaAx,[5:5:35]);
yticklabels(areaAx, "");
end
grid(areaAx, "on");
box(areaAx, "on");
text(areaAx, 0.2, 0.06, sprintf("PAM%d", selectedPamLevels(pamIdx)), ...
"Units", "normalized", ...
"HorizontalAlignment", "left", ...
"VerticalAlignment", "bottom", ...
"BackgroundColor", "white", ...
"EdgeColor", [0.60 0.60 0.60], ...
"LineWidth", 0.5, ...
"Margin", 2, ...
"FontSize", 8, ...
"Interpreter", "none", ...
"Clipping", "on");
end
legend
set(areaFig, "Position", 1e3 .* [0.3500 0.3500 0.7000 0.8000]);
% mat2tikz_improved( ...
% "C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/04_Experimental_Evaluation/tikz/400g/10km_permissible_wavelength.tikz", ...
% "cleanfigure", 0);
%% Local helpers
function fields = appendMissingFields(fields, extraFields)
fields = cellstr(fields);
extraFields = cellstr(extraFields);
for idx = 1:numel(extraFields)
if ~any(strcmp(fields, extraFields{idx}))
fields{end+1, 1} = extraFields{idx}; %#ok<AGROW>
end
end
end
function values = numericColumn(values)
if iscell(values)
values = string(values);
end
if isstring(values) || ischar(values)
values = str2double(values);
end
values = double(values);
end
function styles = defaultAlgorithmStyles()
styles = table( ...
["vnle"; ...
"vnle_pf_mlse"; ...
"vnle_db_mlse"; ...
"ml_mlse"; ...
"db_encoded"], ...
[equalizer_structure.vnle; ...
equalizer_structure.vnle_pf_mlse; ...
equalizer_structure.vnle_db_mlse; ...
equalizer_structure.ml_mlse; ...
equalizer_structure.db_encoded], ...
["VNLE"; ...
"VNLE + PF + MLSE"; ...
"VNLE DBt. + MLSE"; ...
"ML pre-EQ + Viterbi"; ...
"DBS + VNLE + MLSE"], ...
["o"; "square"; "diamond"; "^"; "v"], ...
["-"; "-"; "-"; "-"; "-"], ...
["w"; "w"; "w"; "w"; "w"], ...
[clr.Paired.red; ...
clr.Paired.green; ...
clr.Paired.blue; ...
clr.Paired.purple; ...
clr.Paired.orange], ...
'VariableNames', ["algorithm_key", "eq", "name", "marker", ...
"lineStyle", "markerFaceColor", "color"]);
end
function preEmphasis = derivePreEmphasis(precompAmp, dbMode)
preEmphasis = false(size(dbMode));
validPrecomp = isfinite(precompAmp);
preEmphasis(validPrecomp) = precompAmp(validPrecomp) > -45;
missingPrecomp = ~validPrecomp;
preEmphasis(missingPrecomp) = dbMode(missingPrecomp) == double(db_mode.no_db);
end
function plotData = buildNormalMetricRows(data)
baseRows = data(isfinite(data.BER), :);
baseRows.precode = false(height(baseRows), 1);
baseRows.BER_plot = baseRows.BER;
baseRows.algorithm_key = algorithmKeyFromEqualizer(baseRows.equalizer_structure);
baseRows = baseRows(baseRows.algorithm_key ~= "", :);
if ismember("BER_precoded", string(data.Properties.VariableNames))
precodedRows = data(isfinite(data.BER_precoded), :);
precodedRows.precode = true(height(precodedRows), 1);
precodedRows.BER_plot = precodedRows.BER_precoded;
precodedRows.algorithm_key = algorithmKeyFromEqualizer( ...
precodedRows.equalizer_structure);
precodedRows = precodedRows(precodedRows.algorithm_key ~= "", :);
plotData = [baseRows; precodedRows];
else
warning("plot_wavelength_best_algos:NoPrecodedBer", ...
"BER_precoded was not returned. Plotting only BER rows for normal algorithms.");
plotData = baseRows;
end
end
function plotData = buildDuobinarySignalingRows(data)
plotData = data(isfinite(data.BER), :);
plotData.precode = false(height(plotData), 1);
plotData.BER_plot = plotData.BER;
plotData.algorithm_key = repmat("db_encoded", height(plotData), 1);
end
function algorithmKey = algorithmKeyFromEqualizer(equalizerColumn)
eqNumeric = equalizerNumeric(equalizerColumn);
algorithmKey = strings(size(eqNumeric));
algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle)) = "vnle";
algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle_pf_mlse)) = ...
"vnle_pf_mlse";
algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle_db_mlse)) = ...
"vnle_db_mlse";
algorithmKey(eqNumeric == enumValue(equalizer_structure.ml_mlse)) = "ml_mlse";
end
function mask = equalizerMask(equalizerColumn, eqValue)
eqNumeric = equalizerNumeric(equalizerColumn);
mask = eqNumeric == enumValue(eqValue);
end
function eqNumeric = equalizerNumeric(equalizerColumn)
if isa(equalizerColumn, "equalizer_structure")
eqNumeric = double(equalizerColumn);
elseif isnumeric(equalizerColumn)
eqNumeric = double(equalizerColumn);
else
equalizerString = string(equalizerColumn);
eqNumeric = str2double(equalizerString);
enumNames = ["vnle", "ffe", "dfe", "vnle_pf_mlse", ...
"vnle_db_mlse", "db_encoded", "ml_mlse"];
enumValues = [ ...
enumValue(equalizer_structure.vnle), ...
enumValue(equalizer_structure.ffe), ...
enumValue(equalizer_structure.dfe), ...
enumValue(equalizer_structure.vnle_pf_mlse), ...
enumValue(equalizer_structure.vnle_db_mlse), ...
enumValue(equalizer_structure.db_encoded), ...
enumValue(equalizer_structure.ml_mlse)];
for idx = 1:numel(enumNames)
missingNumeric = isnan(eqNumeric);
eqNumeric(missingNumeric & equalizerString == enumNames(idx)) = ...
enumValues(idx);
end
end
end
function value = enumValue(enumEntry)
value = double(enumEntry);
end
function bestData = bestBerByAlgorithmAndWavelength(data)
groupVars = ["algorithm_key", "bitrate_Gbps", "wavelength"];
groupId = findgroups(data(:, groupVars));
keepIdx = NaN(max(groupId), 1);
for curGroup = 1:max(groupId)
rowIdx = find(groupId == curGroup);
[~, localBestIdx] = min(data.BER_plot(rowIdx));
keepIdx(curGroup) = rowIdx(localBestIdx);
end
bestData = sortrows(data(keepIdx, :), groupVars);
end
function [xFit, yFit, fitModel] = fitLogBerCurve(x, y, maxOrder)
valid = isfinite(x) & isfinite(y) & y > 0;
x = x(valid);
y = y(valid);
[x, orderIdx] = sort(x(:));
y = y(orderIdx);
[x, uniqueIdx] = unique(x);
y = y(uniqueIdx);
if numel(x) < 2
xFit = [];
yFit = [];
fitModel = [];
return
end
% Keep the requested fourth-order fit when enough wavelength points exist.
% With fewer than five unique points, use the highest identifiable order.
fitOrder = min(maxOrder, numel(x) - 1);
[coefficients, ~, mu] = polyfit(x, log10(y), fitOrder);
fitModel.coefficients = coefficients;
fitModel.mu = mu;
xFit = linspace(min(x), max(x), 300).';
yFit = 10 .^ polyval(coefficients, (xFit - mu(1)) ./ mu(2));
end
function areaData = computePermissibleWavelengthAreas(data, algoStyles, ...
fecThreshold, maxOrder)
grossrateValues = unique(data.grossrate_Gbps(isfinite( ...
data.grossrate_Gbps))).';
wavelengthGrid = unique(data.wavelength(isfinite(data.wavelength))).';
rows = cell(height(algoStyles) * numel(grossrateValues), 5);
rowIdx = 0;
for styleIdx = 1:height(algoStyles)
style = algoStyles(styleIdx, :);
algorithmMask = data.algorithm_key == style.algorithm_key;
for grossrateIdx = 1:numel(grossrateValues)
grossrateGb = grossrateValues(grossrateIdx);
rowMask = algorithmMask & data.grossrate_Gbps == grossrateGb;
if ~any(rowMask)
continue
end
curveData = sortrows(data(rowMask, :), "wavelength");
[wavelengthMin, wavelengthMax] = permissibleWavelengthRange( ...
curveData.wavelength, curveData.BER_plot, fecThreshold, ...
maxOrder, wavelengthGrid);
if any(~isfinite([wavelengthMin, wavelengthMax]))
continue
end
rowIdx = rowIdx + 1;
rows(rowIdx, :) = { ...
style.algorithm_key, ...
grossrateGb, ...
wavelengthMin, ...
wavelengthMax, ...
wavelengthMax - wavelengthMin};
end
end
if rowIdx == 0
areaData = table();
else
rows = rows(1:rowIdx, :);
areaData = cell2table(rows, "VariableNames", ...
["algorithm_key", "grossrate_Gbps", ...
"wavelength_min_allowed_nm", "wavelength_max_allowed_nm", ...
"permissible_wavelength_area_nm"]);
areaData.algorithm_key = string(areaData.algorithm_key);
end
end
function [wavelengthMin, wavelengthMax] = permissibleWavelengthRange(x, y, ...
fecThreshold, maxOrder, wavelengthGrid)
wavelengthMin = NaN;
wavelengthMax = NaN;
valid = isfinite(x) & isfinite(y) & y > 0;
x = x(valid);
y = y(valid);
if isempty(x)
return
end
[x, orderIdx] = sort(x(:));
y = y(orderIdx);
[x, uniqueIdx] = unique(x);
y = y(uniqueIdx);
wavelengthGrid = sort(wavelengthGrid(:));
if isempty(wavelengthGrid)
return
end
globalRange = [wavelengthGrid(1), wavelengthGrid(end)];
edgeTolerance = max(1e-9, 1e-9 .* max(abs(globalRange)));
missingLeftEdge = x(1) > globalRange(1) + edgeTolerance;
missingRightEdge = x(end) < globalRange(2) - edgeTolerance;
% If all observed points are below the FEC, use the complete sweep span
% instead of requiring polynomial roots outside the observed range.
if all(y <= fecThreshold)
wavelengthMin = globalRange(1);
wavelengthMax = globalRange(2);
return
end
[~, ~, fitModel] = fitLogBerCurve(x, y, maxOrder);
if isempty(fitModel)
return
end
xRange = [x(1), x(end)];
crossingCoefficients = fitModel.coefficients;
crossingCoefficients(end) = crossingCoefficients(end) - log10(fecThreshold);
rootsAtThreshold = roots(crossingCoefficients);
realRoots = real(rootsAtThreshold(abs(imag(rootsAtThreshold)) < ...
1e-7 .* max(1, abs(real(rootsAtThreshold)))));
realRoots = realRoots .* fitModel.mu(2) + fitModel.mu(1);
realRoots = sort(realRoots(realRoots >= xRange(1) & realRoots <= xRange(2)));
if isempty(realRoots)
return
end
% Include the observed boundaries in the interval search. This lets each
% side be handled independently: a side that remains below FEC uses the
% corresponding observed edge, while a side that crosses FEC uses its root.
intervalEdges = [xRange(1); realRoots(:); xRange(2)];
intervalWidths = diff(intervalEdges);
intervalIsAllowed = false(size(intervalWidths));
thresholdLogBer = log10(fecThreshold);
for intervalIdx = 1:numel(intervalEdges)-1
midpoint = mean(intervalEdges(intervalIdx:intervalIdx + 1));
normalizedMidpoint = (midpoint - fitModel.mu(1)) ./ fitModel.mu(2);
intervalIsAllowed(intervalIdx) = polyval(fitModel.coefficients, ...
normalizedMidpoint) <= thresholdLogBer;
end
allowedIntervals = find(intervalIsAllowed);
if isempty(allowedIntervals)
return
end
[~, widestAllowedIdx] = max(intervalWidths(allowedIntervals));
selectedInterval = allowedIntervals(widestAllowedIdx);
wavelengthMin = intervalEdges(selectedInterval);
wavelengthMax = intervalEdges(selectedInterval + 1);
% A missing terminal wavelength sample is treated as below FEC when the
% permissible interval reaches that side of the observed curve.
if selectedInterval == 1 && missingLeftEdge
wavelengthMin = globalRange(1);
end
if selectedInterval == numel(intervalWidths) && missingRightEdge
wavelengthMax = globalRange(2);
end
end
function keep = hasAlgorithmRows(data, algoStyles)
keep = false(height(algoStyles), 1);
for idx = 1:height(algoStyles)
keep(idx) = any(data.algorithm_key == algoStyles.algorithm_key(idx));
end
end
function plotFecLines(ax, fecLevels)
xl = xlim(ax);
for idx = 1:numel(fecLevels)
h = plot(ax, xl, [fecLevels(idx), fecLevels(idx)], ...
"LineWidth", 1, ...
"LineStyle", "--", ...
"Color", [0.25 0.25 0.25], ...
"HandleVisibility", "off");
h.Annotation.LegendInformation.IconDisplayStyle = "off";
end
end
function markers = bitrateMarkerSet(n)
markerOptions = ["o", "square", "diamond", "^", "v", ">", "<"];
if n <= numel(markerOptions)
markers = markerOptions(1:n);
else
markers = markerOptions(mod(0:n-1, numel(markerOptions)) + 1);
end
end
function colors = sequentialColors(baseColor, n, algorithmKey)
if n <= 1
colors = baseColor;
return
end
lightBlend = linspace(0.72, 0.00, n).';
darkBlend = linspace(0.00, 0.35, n).';
if algorithmKey == "db_encoded"
lightBlend = linspace(0.62, 0.00, n).';
darkBlend = linspace(0.00, 0.12, n).';
end
colors = zeros(n, 3);
for idx = 1:n
color = (1 - lightBlend(idx)) .* baseColor + lightBlend(idx) .* [1 1 1];
color = (1 - darkBlend(idx)) .* color;
colors(idx, :) = color;
end
colors = flip(colors);
end

View File

@@ -0,0 +1,85 @@
%% Replot opt_filt_analysis.fig grouped by PAM level
clear; close all; clc;
figPath = fullfile(fileparts(mfilename("fullpath")), "opt_filt_analysis.fig");
figIn = openfig(figPath, "invisible");
srcAx = findall(figIn, "Type", "axes");
srcAx = srcAx(1);
srcLines = flipud(findall(srcAx, "Type", "line"));
% Get every plotted datapoint and its display name.
data = arrayfun(@(h) struct( ...
"XData", h.XData, ...
"YData", h.YData, ...
"DisplayName", string(h.DisplayName)), srcLines);
displayNames = string({data.DisplayName})';
disp(displayNames);
selectedEqType = "ffe"; % choose "ffe" or "mlse"
eqtype = ["ffe", "mlse"];
eqGroup = repmat("", size(displayNames));
for k = 1:numel(eqtype)
eqGroup(contains(lower(displayNames), eqtype(k))) = eqtype(k);
end
if ~ismember(selectedEqType, eqtype)
error("selectedEqType must be either 'ffe' or 'mlse'.");
end
keep = eqGroup == selectedEqType;
srcLines = srcLines(keep);
data = data(keep);
displayNames = displayNames(keep);
pamLevels = ["pam4", "pam6", "pam8"];
pamGroup = repmat("", size(displayNames));
for k = 1:numel(pamLevels)
pamGroup(contains(lower(displayNames), pamLevels(k))) = pamLevels(k);
end
% Group by the first part of the display name.
groupLabels = ["A", "B", "C", "D", "E"];
groupGroup = repmat("", size(displayNames));
nameLower = lower(displayNames);
groupGroup(startsWith(nameLower, "120 gb") | ...
startsWith(nameLower, "144 gb") | startsWith(nameLower, "170 gb")) = "A";
groupGroup(startsWith(nameLower, "optfil")) = "B";
groupGroup(startsWith(nameLower, "thormax")) = "C";
groupGroup(startsWith(nameLower, "thormax optfil")) = "E";
groupGroup(startsWith(nameLower, "thormax optfil ohne fl")) = "D";
groupColors = [clr.Paired.red; clr.Paired.blue; clr.Paired.green; ...
clr.Paired.orange; clr.Paired.purple];
groupStyles = {"-", "--", ":", ":", ":"};
groupMarkers = {"o", "s", "^", "d", "v"};
figOut = figure("Name", "opt_filt_analysis by " + upper(selectedEqType) + " and PAM");
layout = tiledlayout(figOut, 1, 3, "TileSpacing", "compact", "Padding", "compact");
for k = 1:numel(pamLevels)
ax = nexttile(layout);
hold(ax, "on");
for lineIdx = find(pamGroup == pamLevels(k))'
newLine = copyobj(srcLines(lineIdx), ax);
groupIdx = find(groupLabels == groupGroup(lineIdx), 1);
newLine.Color = groupColors(groupIdx, :);
newLine.Marker = groupMarkers{groupIdx};
newLine.LineWidth = 1;
newLine.LineStyle = groupStyles{groupIdx};
newLine.MarkerSize = 3;
newLine.DisplayName = groupGroup(lineIdx);
end
ax.XScale = srcAx.XScale;
ax.YScale = srcAx.YScale;
title(ax, upper(pamLevels(k)));
xlabel(ax, srcAx.XLabel.String);
ylabel(ax, srcAx.YLabel.String);
grid(ax, "on");
legend(ax, "show", "Interpreter", "none", "Location", "best");
beautifyBERplot("setcolors", false, "setmarkers", false);
end
close(figIn);

View File

@@ -0,0 +1,430 @@
%% Reprocess baudrate-sweep MAT files with dsp_400g_recipe
% File naming convention:
% PAMX_b2b_baudrate20241024_210648PAM_4_fsym_100_bits.mat
% PAMX_b2b_baudrate20241024_210648PAM_4_fsym_100_symbols.mat
% PAMX_b2b_baudrate20241024_210648PAM_4_fsym_100_rop_0_5_rx_signal.mat
clear; clc;
%% Configuration
dataDir = "D:\baudrate_sweep_b2b";
filePrefix = "PAMX_b2b_baudrate20241024_210648";
outputFile = fullfile(dataDir, filePrefix + "_dsp400g_reprocessed_wh.mat");
originalWarehouseFile = fullfile(dataDir, filePrefix + "_wh_final.mat");
pamLevels = [4, 6, 8];
fsymGBd = 100:6:196;
ropAtten = 0:0.5:7;
% ropAtten = 0:0.5:7;
useParallel = true;
batchSize = 24; % warehouse writes happen after each batch
maxJobs = inf; % use a small number for smoke tests
saveEveryBatches = 1;
debugPlots = false;
storeDspOutput = false; % full output objects can make the warehouse very large
copyOriginalPowerMetadata = true;
recipeParams = struct( ...
"run_ffe", true, ...
"run_vnle_mlse", true, ...
"run_dbtgt", true, ...
"run_ml_mlse", true, ...
"run_ml_mlse_db", false, ...
"run_mlse_db", false, ...
"plot_input_signal", false, ...
"plot_output_signals", false);
%% Discover all file/acquisition jobs
jobs = buildJobList(dataDir, filePrefix, pamLevels, fsymGBd, ropAtten);
if isempty(jobs)
error("reprocess:NoJobs", "No processable RX acquisition jobs were found.");
end
if isfinite(maxJobs)
jobs = jobs(1:min(numel(jobs), maxJobs));
end
maxAcquisitionIdx = max([jobs.acquisition_idx]);
fprintf("Discovered %d acquisition jobs, max acquisition index S{%d}.\n", ...
numel(jobs), maxAcquisitionIdx);
powerLookup = loadOriginalPowerLookup(originalWarehouseFile, ...
copyOriginalPowerMetadata);
%% Build output warehouse
params = struct();
params.fsym = fsymGBd .* 1e9;
params.rop_atten = ropAtten;
params.M = pamLevels;
params.acquisition_idx = 1:maxAcquisitionIdx;
wh = DataStorage(params);
wh.addStorage("ber_ffe");
wh.addStorage("ber_mlse");
wh.addStorage("ber_db");
wh.addStorage("ber_ffe_precoded");
wh.addStorage("ber_mlse_precoded");
wh.addStorage("ber_db_precoded");
wh.addStorage("rop");
wh.addStorage("pd_in");
wh.addStorage("rx_file");
wh.addStorage("status");
wh.addStorage("error_message");
if storeDspOutput
wh.addStorage("dsp_output");
end
parallelEnabled = useParallel && ensureParallelPool();
if parallelEnabled
pool = gcp("nocreate");
fprintf("Processing with parfor on %d workers.\n", pool.NumWorkers);
else
fprintf("Processing serially.\n");
end
%% Process jobs in parallel batches, write warehouse serially
numBatches = ceil(numel(jobs) / batchSize);
for batchIdx = 1:numBatches
firstJob = (batchIdx - 1) * batchSize + 1;
lastJob = min(batchIdx * batchSize, numel(jobs));
batchJobs = jobs(firstJob:lastJob);
batchResults = cell(numel(batchJobs), 1);
fprintf("Batch %d/%d: jobs %d-%d of %d\n", ...
batchIdx, numBatches, firstJob, lastJob, numel(jobs));
if parallelEnabled
parfor localIdx = 1:numel(batchJobs)
batchResults{localIdx} = processOneJob( ...
batchJobs(localIdx), recipeParams, debugPlots, storeDspOutput);
end
else
for localIdx = 1:numel(batchJobs)
batchResults{localIdx} = processOneJob( ...
batchJobs(localIdx), recipeParams, debugPlots, storeDspOutput);
end
end
for localIdx = 1:numel(batchResults)
wh = writeResultToWarehouse(wh, batchResults{localIdx}, ...
storeDspOutput, powerLookup);
end
if mod(batchIdx, saveEveryBatches) == 0 || batchIdx == numBatches
save(outputFile, "wh", "recipeParams", "dataDir", "filePrefix", ...
"jobs", "storeDspOutput", "copyOriginalPowerMetadata", ...
"originalWarehouseFile", "-v7.3");
fprintf("Saved checkpoint to %s\n", outputFile);
end
end
fprintf("Saved reprocessed warehouse to %s\n", outputFile);
wh.showInfo;
%% Local helpers
function jobs = buildJobList(dataDir, filePrefix, pamLevels, fsymGBd, ropAtten)
jobs = struct( ...
"M", {}, ...
"fsym_GBd", {}, ...
"fsym", {}, ...
"rop_atten", {}, ...
"acquisition_idx", {}, ...
"bits_file", {}, ...
"symbols_file", {}, ...
"rx_file", {});
for M = pamLevels
for fsymGb = fsymGBd
fsym = fsymGb .* 1e9;
bitsFile = fullfile(dataDir, sprintf("%sPAM_%d_fsym_%d_bits.mat", ...
filePrefix, M, fsymGb));
symbolsFile = fullfile(dataDir, sprintf("%sPAM_%d_fsym_%d_symbols.mat", ...
filePrefix, M, fsymGb));
if ~isfile(bitsFile) || ~isfile(symbolsFile)
warning("reprocess:MissingReference", ...
"Skipping PAM-%d %.0f GBd: missing bits or symbols file.", ...
M, fsymGb);
continue
end
for rop = ropAtten
rxFile = fullfile(dataDir, sprintf("%sPAM_%d_fsym_%d_rop_%s_rx_signal.mat", ...
filePrefix, M, fsymGb, ropToken(rop)));
if ~isfile(rxFile)
warning("reprocess:MissingRx", "Missing RX file: %s", rxFile);
continue
end
numAcquisitions = acquisitionCount(rxFile);
for acqIdx = 1:numAcquisitions
jobs(end+1) = struct( ... %#ok<AGROW>
"M", M, ...
"fsym_GBd", fsymGb, ...
"fsym", fsym, ...
"rop_atten", rop, ...
"acquisition_idx", acqIdx, ...
"bits_file", bitsFile, ...
"symbols_file", symbolsFile, ...
"rx_file", rxFile);
end
end
end
end
end
function numAcq = acquisitionCount(rxFile)
info = whos("-file", rxFile, "S");
if isempty(info)
warning("reprocess:MissingS", "RX file has no variable S: %s", rxFile);
numAcq = 0;
else
numAcq = prod(info.size);
end
end
function result = processOneJob(job, recipeParams, debugPlots, storeDspOutput)
result = emptyJobResult(job);
try
bitsData = load(job.bits_file, "Bits");
symbolsData = load(job.symbols_file, "Symbols");
ScpeSigRaw = loadRxSignal(job.rx_file, job.acquisition_idx);
dataTable = makeRecipeDataTable(job);
fprintf("PAM-%d, %.0f GBd, ROP atten %.1f dB, S{%d}\n", ...
job.M, job.fsym_GBd, job.rop_atten, job.acquisition_idx);
dspOut = dsp_400g_recipe(ScpeSigRaw, symbolsData.Symbols, bitsData.Bits, ...
"fsym", job.fsym, ...
"M", job.M, ...
"duob_mode", db_mode.no_db, ...
"dataTable", dataTable, ...
"userParameters", recipeParams, ...
"debug_plots", debugPlots);
result.ber_ffe = extractBer(dspOut, "ffe_package", "BER");
result.ber_mlse = extractBer(dspOut, "mlse_package", "BER");
result.ber_db = extractBer(dspOut, "dbtgt_package", "BER");
result.ber_ffe_precoded = extractBer(dspOut, "ffe_package", "BER_precoded");
result.ber_mlse_precoded = extractBer(dspOut, "mlse_package", "BER_precoded");
result.ber_db_precoded = extractBer(dspOut, "dbtgt_package", "BER_precoded");
result.status = "ok";
if storeDspOutput
result.dsp_output = dspOut;
end
catch err
result.status = "failed";
result.error_message = string(err.message);
warning("reprocess:DspFailed", ...
"DSP failed for PAM-%d %.0f GBd ROP %.1f dB S{%d}: %s", ...
job.M, job.fsym_GBd, job.rop_atten, job.acquisition_idx, err.message);
end
end
function result = emptyJobResult(job)
result = struct( ...
"M", job.M, ...
"fsym", job.fsym, ...
"fsym_GBd", job.fsym_GBd, ...
"rop_atten", job.rop_atten, ...
"acquisition_idx", job.acquisition_idx, ...
"rx_file", job.rx_file, ...
"ber_ffe", NaN, ...
"ber_mlse", NaN, ...
"ber_db", NaN, ...
"ber_ffe_precoded", NaN, ...
"ber_mlse_precoded", NaN, ...
"ber_db_precoded", NaN, ...
"status", "not_run", ...
"error_message", "", ...
"dsp_output", []);
end
function wh = writeResultToWarehouse(wh, result, storeDspOutput, powerLookup)
idx = {result.fsym, result.rop_atten, result.M, result.acquisition_idx};
[ropValue, pdInValue] = lookupOriginalPower(powerLookup, ...
result.fsym, result.rop_atten, result.M);
wh.addValueToStorage(result.ber_ffe, "ber_ffe", idx{:});
wh.addValueToStorage(result.ber_mlse, "ber_mlse", idx{:});
wh.addValueToStorage(result.ber_db, "ber_db", idx{:});
wh.addValueToStorage(result.ber_ffe_precoded, "ber_ffe_precoded", idx{:});
wh.addValueToStorage(result.ber_mlse_precoded, "ber_mlse_precoded", idx{:});
wh.addValueToStorage(result.ber_db_precoded, "ber_db_precoded", idx{:});
wh.addValueToStorage(ropValue, "rop", idx{:});
wh.addValueToStorage(pdInValue, "pd_in", idx{:});
wh.addValueToStorage(char(result.rx_file), "rx_file", idx{:});
wh.addValueToStorage(char(result.status), "status", idx{:});
wh.addValueToStorage(char(result.error_message), "error_message", idx{:});
if storeDspOutput
wh.addValueToStorage(result.dsp_output, "dsp_output", idx{:});
end
end
function powerLookup = loadOriginalPowerLookup(originalWarehouseFile, enabled)
powerLookup = struct("enabled", false, "data", table());
if ~enabled
return
end
if ~isfile(originalWarehouseFile)
warning("reprocess:MissingOriginalWarehouse", ...
"Original warehouse not found, storing NaN for rop/pd_in: %s", ...
originalWarehouseFile);
return
end
loadedData = load(originalWarehouseFile);
if isfield(loadedData, "wh")
originalWh = loadedData.wh;
elseif isfield(loadedData, "obj")
originalWh = loadedData.obj;
else
warning("reprocess:NoOriginalWarehouseObject", ...
"Original warehouse file contains neither wh nor obj: %s", ...
originalWarehouseFile);
return
end
requiredStorages = ["rop", "pd_in"];
availableStorages = string(fieldnames(originalWh.sto));
if ~all(ismember(requiredStorages, availableStorages))
warning("reprocess:MissingPowerStorage", ...
"Original warehouse has no complete rop/pd_in storage. Storing NaN.");
return
end
fsymVals = double(originalWh.parameter.fsym.values(:).');
ropAttenVals = double(originalWh.parameter.rop_atten.values(:).');
pamVals = double(originalWh.parameter.M.values(:).');
rows = cell(numel(fsymVals) * numel(ropAttenVals) * numel(pamVals), 5);
rowIdx = 0;
for pamLevel = pamVals
for fsym = fsymVals
for ropAtten = ropAttenVals
rowIdx = rowIdx + 1;
rows(rowIdx, :) = { ...
fsym, ...
ropAtten, ...
pamLevel, ...
getOriginalScalar(originalWh, "rop", fsym, ropAtten, pamLevel), ...
getOriginalScalar(originalWh, "pd_in", fsym, ropAtten, pamLevel)};
end
end
end
powerLookup.enabled = true;
powerLookup.data = cell2table(rows, 'VariableNames', ...
["fsym", "rop_atten", "M", "rop", "pd_in"]);
fprintf("Loaded %d original rop/pd_in metadata rows from %s\n", ...
height(powerLookup.data), originalWarehouseFile);
end
function value = getOriginalScalar(originalWh, storageName, fsym, ropAtten, pamLevel)
value = NaN;
try
rawValue = originalWh.getStoValue(storageName, fsym, ropAtten, pamLevel);
catch
return
end
if isnumeric(rawValue) && isscalar(rawValue)
value = double(rawValue);
end
end
function [ropValue, pdInValue] = lookupOriginalPower(powerLookup, fsym, ropAtten, pamLevel)
ropValue = NaN;
pdInValue = NaN;
if ~powerLookup.enabled || isempty(powerLookup.data)
return
end
match = powerLookup.data.fsym == fsym & ...
powerLookup.data.rop_atten == ropAtten & ...
powerLookup.data.M == pamLevel;
if ~any(match)
return
end
firstMatch = find(match, 1, "first");
ropValue = powerLookup.data.rop(firstMatch);
pdInValue = powerLookup.data.pd_in(firstMatch);
end
function sig = loadRxSignal(rxFile, acquisitionIdx)
rxData = load(rxFile, "S");
if ~iscell(rxData.S) || isempty(rxData.S)
error("RX file %s does not contain a nonempty cell array S.", rxFile);
end
if acquisitionIdx > numel(rxData.S)
error("Requested acquisition S{%d}, but %s only contains %d acquisitions.", ...
acquisitionIdx, rxFile, numel(rxData.S));
end
sig = rxData.S{acquisitionIdx};
end
function token = ropToken(rop)
token = strrep(sprintf("%.1f", rop), ".", "_");
token = regexprep(token, "_0$", "");
end
function dataTable = makeRecipeDataTable(job)
dataTable = table();
dataTable.run_id = makeRunId(job.M, job.fsym, job.rop_atten, job.acquisition_idx);
dataTable.pam_level = job.M;
dataTable.symbolrate = job.fsym;
dataTable.bitrate = job.fsym .* log2(job.M);
dataTable.grossrate = dataTable.bitrate;
dataTable.fiber_length = 0;
dataTable.wavelength = 1310;
dataTable.rop_attenuation = job.rop_atten;
dataTable.acquisition_idx = job.acquisition_idx;
dataTable.rx_file = string(job.rx_file);
end
function runId = makeRunId(M, fsym, rop, acquisitionIdx)
runId = M .* 1e10 + round(fsym .* 1e-6) .* 1e2 + ...
round(rop .* 10) .* 10 + acquisitionIdx;
end
function ber = extractBer(dspOut, packageName, metricName)
ber = NaN;
if ~isfield(dspOut, packageName)
return
end
pkg = dspOut.(packageName);
if ~isfield(pkg, "metrics")
return
end
if isobject(pkg.metrics) && isprop(pkg.metrics, metricName)
ber = pkg.metrics.(metricName);
elseif isstruct(pkg.metrics) && isfield(pkg.metrics, metricName)
ber = pkg.metrics.(metricName);
end
end
function ok = ensureParallelPool()
ok = false;
try
pool = gcp("nocreate");
if ~isempty(pool) && contains(string(class(pool)), "ThreadPool")
delete(pool);
pool = [];
end
if isempty(pool)
pool = parpool("local");
end
ok = ~isempty(pool);
catch err
warning("reprocess:NoParallelPool", ...
"Parallel pool unavailable, falling back to serial processing: %s", ...
err.message);
end
end

View File

@@ -0,0 +1,97 @@
% === 400G DSP settings ===
dsp_options = struct();
dsp_options.mode = "run_id";
dsp_options.recipe = @dsp_400g_recipe;
dsp_options.append_to_db = false;
% dsp_options.append_mpi_reduction_db = false;
dsp_options.start_occurence = 1;
dsp_options.max_occurences = 3;
dsp_options.debug_plots = false;
dsp_options.database_type = "mysql";
dsp_options.dataBase = "labor_highspeed";
if ismac
dsp_options.storage_path = "/Volumes/media/labdata/sioe_labor";
else
dsp_options.storage_path = "W:\labdata\sioe_labor";
end
dsp_options.server = "192.168.178.192";
dsp_options.port = 3306;
dsp_options.user = "silas";
dsp_options.password = "silas";
db = DBHandler("dataBase", [dsp_options.dataBase], ...
"type", dsp_options.database_type, ...
"server", dsp_options.server, ...
"user", dsp_options.user, ...
"password", dsp_options.password);
%% Load normal Signal w/o preemphasis
fp = QueryFilter();
fp.where('Runs','fiber_length','EQUALS', 10);
fp.where('Runs','wavelength','EQUALS', 1310);
fp.where('Runs','bitrate','EQUALS', 420e9);
fp.where('Runs','pam_level','EQUALS', 4);
fp.where('Runs','rop_attenuation','EQUALS', 0);
fp.where('Runs','is_mpi','EQUALS', 0);
fp.where('Runs', 'db_mode','EQUALS', 0);
fields = db.getTableFieldNames('Runs');
[dataTable, query] = db.queryDB(fp, fields);
[~, Symbols_preemph, Scpe_cell_preemph, ~] = loadAndSyncRunSignals(dataTable(1,:), dsp_options);
ScopeSignal = Scpe_cell_preemph{1};
ScopeSignal_preemph = preprocessSignal(ScopeSignal, Symbols_preemph, Symbols_preemph.fs);
%% Load normal Signal w/o preemphasis
fp = QueryFilter();
fp.where('Runs','fiber_length','EQUALS', 10);
fp.where('Runs','wavelength','EQUALS', 1310);
fp.where('Runs','bitrate','EQUALS', 420e9);
fp.where('Runs','pam_level','EQUALS', 4);
fp.where('Runs','rop_attenuation','EQUALS', 0);
fp.where('Runs','is_mpi','EQUALS', 0);
fp.where('Runs', 'db_mode','EQUALS', 1);
fields = db.getTableFieldNames('Runs');
[dataTable, query] = db.queryDB(fp, fields);
[~, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable(1,:), dsp_options);
ScopeSignal = Scpe_cell{1};
ScopeSignal_no_preemph = preprocessSignal(ScopeSignal, Symbols, Symbols.fs);
%% Duobinary
fp = QueryFilter();
fp.where('Runs','fiber_length','EQUALS', 10);
fp.where('Runs','wavelength','EQUALS', 1310);
fp.where('Runs','bitrate','EQUALS', 420e9);
fp.where('Runs','pam_level','EQUALS', 4);
fp.where('Runs','rop_attenuation','EQUALS', 0);
fp.where('Runs','is_mpi','EQUALS', 0);
fp.where('Runs', 'db_mode','EQUALS', 2);
fields = db.getTableFieldNames('Runs');
[dataTable, query] = db.queryDB(fp, fields);
[~, Symbols_db, Scpe_cell_db, found_sync] = loadAndSyncRunSignals(dataTable(1,:), dsp_options);
ScopeSignal = Scpe_cell_db{1};
ScopeSignal_DB = preprocessSignal(ScopeSignal, Symbols_db, Symbols_db.fs);
%%
ScopeSignal_no_preemph.spectrum("displayname",'Full Response w/o preemphasis','fignum',2,'normalizeTo0dB',0,'color',clr.Paired.dblue);
ScopeSignal_preemph.spectrum("displayname",'Full Response w/ preemphasis','fignum',2,'normalizeTo0dB',0,'color',clr.Paired.dgreen);
ScopeSignal_DB.spectrum("displayname",'DB Response w/ preemphasis','fignum',2,'normalizeTo0dB',0,'color',clr.Paired.dorange);

View File

@@ -0,0 +1,115 @@
rates = [420e9];
rcalpha = 0.05;
fsym = rates/2;
apply_pulsef = 1;
M = 4;
%% Normal Tx Signal
duob_mode = db_mode.no_db;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha);
[Digi_sig,~,~] = PAMsource(...
"fsym",fsym,"M",M,"order",21,"useprbs",0,...
"fs_out",256e9,...
"applyclipping",0,"clipfactor",1.5,...
"applypulseform",apply_pulsef,"pulseformer",Pform,...
"randkey",1,...
'duobinary_mode',duob_mode,...
"mrds_code",0,"mrds_blocklength",512).process();
Digi_sig= Digi_sig.normalize("mode","oneone");
Digi_sig = Digi_sig-mean(Digi_sig.signal);
maxamp = -28; %optimized for DB!
precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs);
precomp_path = "W:\labdata\sioe_labor\precomp";
precomp_fn = "lab_high_speed";
Digi_sig_pre = precomp_est.precomp(Digi_sig,'maxampdb',maxamp,'loadPath',precomp_path,'fileName',precomp_fn);
Digi_sig_pre.signal = reshape(Digi_sig_pre.signal,[],1);
Digi_sig_pre = Digi_sig_pre.resample("fs_out",256e9);
Digi_sig_pre= Digi_sig_pre.normalize("mode","oneone");
Digi_sig_pre = Digi_sig_pre-mean(Digi_sig_pre.signal);
% precomp_est.plot
Digi_sig_rx = precomp_est.apply(...
Digi_sig,'maxampdb',0,'loadPath',precomp_path,'fileName',precomp_fn);
Digi_sig_pre_rx = precomp_est.apply(...
Digi_sig_pre,'maxampdb',0,'loadPath',precomp_path,'fileName',precomp_fn);
Digi_sig.spectrum(...
"displayname",'Tx w/o pre-emphasis',...
"fignum",2,"normalizeTo0dB",0,"color",[0,0,0]);
Digi_sig_pre.spectrum(...
"displayname",'Tx w/ pre-emphasis',...
"fignum",2,"normalizeTo0dB",0,"color",clr.Paired.lblue);
Digi_sig_rx.spectrum(...
"displayname",'Rx w/o pre-emphasis',...
"fignum",2,"normalizeTo0dB",0,"color",clr.Paired.dred);
Digi_sig_pre_rx.spectrum(...
"displayname",'Rx w/ pre-emphasis',...
"fignum",2,"normalizeTo0dB",0,"color",clr.Paired.dblue);
%% Duobinary Encoded Tx Signal
duob_mode = db_mode.db_encoded;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha);
[Digi_sig_DB,Symbols,Tx_bits] = PAMsource(...
"fsym",fsym,"M",M,"order",21,"useprbs",0,...
"fs_out",256e9,...
"applyclipping",0,"clipfactor",1.5,...
"applypulseform",apply_pulsef,"pulseformer",Pform,...
"randkey",1,...
'duobinary_mode',duob_mode,...
"mrds_code",0,"mrds_blocklength",512).process();
Digi_sig_DB= Digi_sig_DB.normalize("mode","oneone");
maxamp = -38; %optimized for DB!
precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig_DB.fs);
precomp_path = "W:\labdata\sioe_labor\precomp";
precomp_fn = "lab_high_speed";
Digi_sig_DB_pre = precomp_est.precomp(Digi_sig_DB,'maxampdb',maxamp,'loadPath',precomp_path,'fileName',precomp_fn);
Digi_sig_DB_pre.signal = reshape(Digi_sig_DB_pre.signal,[],1);
Digi_sig_DB_pre = Digi_sig_DB_pre.resample("fs_out",256e9);
Digi_sig_DB_pre = Digi_sig_DB_pre.normalize("mode","oneone");
%%
AWG_ = M8199B("kover",4);
% AWG_ = AWG("fdac",256e9,"f_cutoff",fsym,"lpf_active",0,"kover",4,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1);
El_sig = AWG_.process(Digi_sig);
El_sig_preemph = AWG_.process(Digi_sig_pre);
El_sig_channel = AWG_.process(Digi_sig_channel);
El_sig.spectrum("displayname",'Full Response w/o preemphasis','fignum',1,'normalizeTo0dB',0,'color',clr.Paired.lblue);
El_sig_preemph.spectrum("displayname",'Full Response w/ preemphasis','fignum',1,'normalizeTo0dB',0,'color',clr.Paired.dblue);
El_sig_channel.spectrum("displayname",'Full Response w/ measured channel','fignum',1,'normalizeTo0dB',0,'color',clr.Paired.blue);
El_sig_DB_preemph = AWG_.process(Digi_sig_DB_pre);
El_sig_DB = AWG_.process(Digi_sig_DB);
El_sig_DB.spectrum("displayname",'DB Response w/o preemphasis','fignum',1,'normalizeTo0dB',0,'color',clr.Paired.lorange);
El_sig_DB_preemph.spectrum("displayname",'DB Response w/ preemphasis','fignum',1,'normalizeTo0dB',0,'color',clr.Paired.dorange);

Binary file not shown.

View File

@@ -5,9 +5,11 @@ dsp_options.recipe = @dsp_400g_recipe;
dsp_options.append_to_db = false;
% dsp_options.append_mpi_reduction_db = false;
dsp_options.start_occurence = 1;
dsp_options.max_occurences = 1;
dsp_options.max_occurences = 3;
dsp_options.debug_plots = false;
dsp_options.database_type = "mysql";
dsp_options.dataBase = "labor_highspeed";
@@ -34,15 +36,15 @@ maxRunIds = 1; % keep small until the recipe settings are sett
fp = QueryFilter();
fp.where('Runs','fiber_length','EQUALS', 10);
fp.where('Runs','wavelength','EQUALS', 1310);
fp.where('Runs','bitrate','EQUALS', 330e9);
fp.where('Runs','pam_level','EQUALS', 4);
% fp.where('Runs','bitrate','LESS_THAN', 480e9);
% fp.where('Runs','pam_level','EQUALS', 4);
fp.where('Runs','rop_attenuation','EQUALS', 0);
fp.where('Runs','is_mpi','EQUALS', 0);
fp.where('Runs', 'db_mode','EQUALS', 2);
% fp.where('Runs', 'db_mode','EQUALS', 2);
fields = db.getTableFieldNames('Runs');
[dataTable, query] = db.queryDB(fp, fields);
disp(query);
% disp(query);
dataTable = sortrows(dataTable, {'bitrate', 'run_id'});
@@ -60,21 +62,12 @@ fprintf("Selected %d run_id(s): %s\n", numel(run_ids), mat2str(run_ids));
%% Parameter sweep
dsp_options.userParameters = struct();
dsp_options.userParameters.run_ml_mlse_db = true;
dsp_options.userParameters.run_mlse_db = true;
% Enable/disable equalizer branches.
% dsp_options.userParameters.run_ffe = false;
% dsp_options.userParameters.run_vnle = false;
% dsp_options.userParameters.run_dfe = false;
% dsp_options.userParameters.run_vnle_mlse = true;
% dsp_options.userParameters.run_dbtgt = true;
% Examples for parameter loops. DataStorage expands every vector-valued field.
% dsp_options.userParameters.len_tr = 4096*2;
% dsp_options.userParameters.pf_ncoeffs = 1;
% dsp_options.userParameters.pf_ncoeffs = [1,2,3];
% dsp_options.userParameters.decoding_mode = [db_decoder.memoryless,db_decoder.sequencedetection];
dsp_options.userParameters.decoding_mode = [db_decoder.sequencedetection];
% dsp_options.userParameters.pf_ncoeffs = [1, 2, 3];
% dsp_options.userParameters.mu_dc = [0, 1e-5, 1e-4];
% dsp_options.userParameters.run_ml_mlse_db = [false, true];
@@ -93,7 +86,7 @@ fprintf("-> [ %d run_id(s) x %d userParam combination(s) = %d job(s) ] x %d real
%% Run
[results, wh] = submitJobs(run_ids, dsp_options, processingMode.serial, ...
[results, wh] = submitJobs(run_ids, dsp_options, processingMode.parallel, ...
"wh", wh, ...
"waitbar", true);
@@ -102,6 +95,7 @@ fprintf("-> [ %d run_id(s) x %d userParam combination(s) = %d job(s) ] x %d real
printBerSummary(wh);
plotBerVsBitrateQuick(wh, dataTable);
plotMlseBerPrecodedVsBitrateByPf(wh, dataTable);
function printBerSummary(wh)
storageNames = fieldnames(wh.sto);
@@ -182,7 +176,7 @@ for storageIdx = 1:numel(storageNames)
end
modeData = sortrows(plotData(rowMask, :), "bitrate_Gbps");
summaryData = groupsummary(modeData, "bitrate_Gbps", "median", "BER");
summaryData = groupsummary(modeData, "bitrate_Gbps", "min", "BER");
summaryData = sortrows(summaryData, "bitrate_Gbps");
color = quickPlotColor(modeIdx);
marker = markers(1 + mod(storageIdx - 1, numel(markers)));
@@ -196,7 +190,7 @@ for storageIdx = 1:numel(storageNames)
"MarkerEdgeAlpha", 0.25, ...
"HandleVisibility", "off");
plot(ax, summaryData.bitrate_Gbps, summaryData.median_BER, ...
plot(ax, summaryData.bitrate_Gbps, summaryData.min_BER, ...
"LineStyle", lineStyles(metricIdx), ...
"Marker", marker, ...
"MarkerSize", 5, ...
@@ -279,6 +273,108 @@ if ~isempty(plotData)
end
end
function plotMlseBerPrecodedVsBitrateByPf(wh, dataTable)
storageName = "mlse_package";
if ~isfield(wh.sto, storageName)
fprintf("No %s storage available for BERp plot.\n", storageName);
return
end
runIds = dataTable.run_id(:);
bitrates = dataTable.bitrate(:);
pfCol = zeros(0, 1);
bitrateCol = zeros(0, 1);
berpCol = zeros(0, 1);
storageValues = wh.sto.(storageName);
for linIdx = 1:numel(storageValues)
[phys, storedValue] = wh.getPhysAndValueByLinIndex(storageName, linIdx);
if ~isfield(phys, "pf_ncoeffs")
continue
end
berp = extractMinMetricValue(storedValue, "BER_precoded");
if ~isfinite(berp) || berp <= 0
continue
end
runId = resolveRunId(phys, runIds);
bitrate = resolveBitrate(runId, runIds, bitrates);
if ~isfinite(bitrate)
continue
end
pfCol(end+1, 1) = double(phys.pf_ncoeffs); %#ok<AGROW>
bitrateCol(end+1, 1) = double(bitrate) * 1e-9; %#ok<AGROW>
berpCol(end+1, 1) = double(berp); %#ok<AGROW>
end
if isempty(berpCol)
fprintf("No BERp values available for %s.\n", storageName);
return
end
plotData = table(pfCol, bitrateCol, berpCol, ...
'VariableNames', ["pf_ncoeffs", "bitrate_Gbps", "BERp"]);
summaryData = groupsummary(plotData, ["pf_ncoeffs", "bitrate_Gbps"], ...
"min", "BERp");
fig = figure(403); clf;
ax = axes(fig); hold(ax, "on");
pfValues = [1, 2, 3];
markers = ["o", "square", "diamond"];
colors = lines(numel(pfValues));
for pfIdx = 1:numel(pfValues)
pfValue = pfValues(pfIdx);
rowMask = summaryData.pf_ncoeffs == pfValue;
if ~any(rowMask)
continue
end
curveData = sortrows(summaryData(rowMask, :), "bitrate_Gbps");
plot(ax, curveData.bitrate_Gbps, curveData.min_BERp, ...
"LineWidth", 1.3, ...
"Marker", markers(pfIdx), ...
"MarkerSize", 5, ...
"Color", colors(pfIdx, :), ...
"DisplayName", sprintf("pf\\_ncoeffs = %d", pfValue));
end
xlabel(ax, "Bitrate [Gb/s]");
ylabel(ax, "min BERp");
title(ax, "MLSE BERp vs bitrate");
set(ax, "YScale", "log");
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
set(fig, "Position", [1200, 450, 560, 360]);
end
function minValue = extractMinMetricValue(value, metricName)
minValue = NaN;
if isempty(value)
return
end
if ~iscell(value)
value = {value};
end
metricValues = NaN(1, numel(value));
for packageIdx = 1:numel(value)
package = value{packageIdx};
if ~isstruct(package) || ~isfield(package, "metrics")
continue
end
metricValues(packageIdx) = readMetricValue(package.metrics, metricName);
end
metricValues = metricValues(isfinite(metricValues) & metricValues > 0);
if ~isempty(metricValues)
minValue = min(metricValues);
end
end
function metricRows = extractBerMetricRows(values)
metricNames = strings(0, 1);
berValues = zeros(0, 1);

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

@@ -0,0 +1,274 @@
%AMPLIFIER_GAIN_ANALYSIS Compare the Aeon SOA and Thorlabs PDFA gain.
%
% The first figure shows gain versus wavelength for the pump levels that
% are available in the wavelength-sweep warehouses. The second figure
% shows gain versus pump level at -10 dBm input power and 1310 nm. The
% third and fourth figures show output power and OSNR versus wavelength.
% The fifth figure estimates the optical noise figure from the ASE level.
%
% The 2 dB setup correction follows amplifier_input_output_curve.m:
% gain = Pout - (Pin - 2 dB) = Pout - Pin + 2 dB.
scriptDir = fileparts(mfilename("fullpath"));
repoDir = fileparts(fileparts(scriptDir));
warehouseClassDir = fullfile(repoDir, "Classes", "Warehouse_class", "classes");
if ~contains(string(path), warehouseClassDir)
addpath(warehouseClassDir);
end
inputPowerDb = -10;
inputPowerCorrectionDb = 2;
correctedInputPowerDb = inputPowerDb - inputPowerCorrectionDb;
wavelengthForPumpSweepNm = 1310;
osaResolutionNm = 0.1;
% Wavelength-dependent measurements: pump values are 50, 75, and 100 %.
lambdaSweepAeon = load(fullfile(scriptDir, "aeon_soa_measurement_lambda_plaser_pump.mat"), "wh");
lambdaSweepThorlabs = load(fullfile(scriptDir, "thorlabs_pdfa_measurement_lambda_plaser_pump.mat"), "wh");
whAeonLambda = lambdaSweepAeon.wh;
whThorlabsLambda = lambdaSweepThorlabs.wh;
% Pump-level measurements: pump values are available from 0 to 100 % in 5 % steps.
pumpSweepAeon = load(fullfile(scriptDir, "aeon_soa_measurement_pump_level_sweep.mat"), "wh");
pumpSweepThorlabs = load(fullfile(scriptDir, "thorlabs_pdfa_measurement_pump_level_sweep.mat"), "wh");
whAeonPump = pumpSweepAeon.wh;
whThorlabsPump = pumpSweepThorlabs.wh;
if ~ismember(inputPowerDb, whAeonLambda.parameter.laserpower.values)
error("Input power %g dBm is not present in the wavelength-sweep warehouse.", inputPowerDb);
end
if ~ismember(wavelengthForPumpSweepNm, whAeonPump.parameter.lambda.values)
error("Wavelength %g nm is not present in the pump-sweep warehouse.", wavelengthForPumpSweepNm);
end
% The warehouse uses physical values as query arguments. Vector queries
% return one row per requested physical value.
wavelengthsNm = whAeonLambda.parameter.lambda.values;
pumpLevelsLambda = whAeonLambda.parameter.pump.values;
pumpLevelsSweep = whAeonPump.parameter.pump.values;
%% Gain versus wavelength
figure("Name", "Amplifier gain versus wavelength", "Color", "w");
tiledlayout(1, 2, "TileSpacing", "compact", "Padding", "compact");
lambdaWarehouses = {whAeonLambda, whThorlabsLambda};
amplifierNames = ["Aeon SOA", "Thorlabs PDFA"];
for amplifierIndex = 1:numel(lambdaWarehouses)
nexttile;
hold on;
colors = lines(numel(pumpLevelsLambda));
wh = lambdaWarehouses{amplifierIndex};
for pumpIndex = 1:numel(pumpLevelsLambda)
pumpLevel = pumpLevelsLambda(pumpIndex);
signalOutputDb = wh.getStoValue("psig_osa", inputPowerDb, wavelengthsNm, pumpLevel);
totalOutputDb = wh.getStoValue("psig_total", inputPowerDb, wavelengthsNm, pumpLevel);
signalGainDb = signalOutputDb - correctedInputPowerDb;
totalGainDb = totalOutputDb - correctedInputPowerDb;
plot(wavelengthsNm, signalGainDb, "-", "Color", colors(pumpIndex, :), ...
"LineWidth", 1.5, ...
"DisplayName", sprintf("Pump %g%%: signal", pumpLevel));
plot(wavelengthsNm, totalGainDb, "--", "Color", colors(pumpIndex, :), ...
"LineWidth", 1.2, ...
"DisplayName", sprintf("Pump %g%%: total", pumpLevel));
end
title(amplifierNames(amplifierIndex), "Interpreter", "none");
xlabel("Wavelength [nm]", "Interpreter", "none");
ylabel("Gain [dB]", "Interpreter", "none");
grid on;
box on;
legend("Location", "best", "Interpreter", "none");
end
sgtitle(sprintf("Gain versus wavelength at P_{in} = %g dBm (2 dB setup correction)", inputPowerDb), ...
"Interpreter", "tex");
%% Gain versus pump level at -10 dBm input power
figure("Name", "Amplifier gain versus pump level", "Color", "w");
hold on;
colors = lines(numel(lambdaWarehouses));
pumpWarehouses = {whAeonPump, whThorlabsPump};
for amplifierIndex = 1:numel(pumpWarehouses)
wh = pumpWarehouses{amplifierIndex};
signalOutputDb = wh.getStoValue("psig_osa", inputPowerDb, wavelengthForPumpSweepNm, pumpLevelsSweep);
totalOutputDb = wh.getStoValue("psig_total", inputPowerDb, wavelengthForPumpSweepNm, pumpLevelsSweep);
signalGainDb = signalOutputDb - correctedInputPowerDb;
totalGainDb = totalOutputDb - correctedInputPowerDb;
plot(pumpLevelsSweep, signalGainDb, "-", "Color", colors(amplifierIndex, :), ...
"LineWidth", 1.6, ...
"DisplayName", amplifierNames(amplifierIndex) + " - signal");
plot(pumpLevelsSweep, totalGainDb, "--", "Color", colors(amplifierIndex, :), ...
"LineWidth", 1.3, ...
"DisplayName", amplifierNames(amplifierIndex) + " - total");
end
title(sprintf("Gain versus pump level at P_{in} = %g dBm, %g nm", ...
inputPowerDb, wavelengthForPumpSweepNm), "Interpreter", "tex");
xlabel("Pump level [percent]", "Interpreter", "none");
ylabel("Gain [dB]", "Interpreter", "none");
grid on;
box on;
legend("Location", "best", "Interpreter", "none");
%% Output power versus wavelength
figure("Name", "Amplifier output power versus wavelength", "Color", "w");
tiledlayout(1, 2, "TileSpacing", "compact", "Padding", "compact");
for amplifierIndex = 1:numel(lambdaWarehouses)
nexttile;
hold on;
colors = lines(numel(pumpLevelsLambda));
wh = lambdaWarehouses{amplifierIndex};
legendHandles = gobjects(numel(pumpLevelsLambda), 1);
for pumpIndex = 1:numel(pumpLevelsLambda)
pumpLevel = pumpLevelsLambda(pumpIndex);
signalOutputDb = wh.getStoValue("psig_osa", inputPowerDb, wavelengthsNm, pumpLevel);
totalOutputDb = wh.getStoValue("psig_total", inputPowerDb, wavelengthsNm, pumpLevel);
plot(wavelengthsNm, signalOutputDb, "-", "Color", colors(pumpIndex, :), ...
"LineWidth", 1.5, ...
"DisplayName", sprintf("Pump %g%%: signal", pumpLevel));
plot(wavelengthsNm, totalOutputDb, "--", "Color", colors(pumpIndex, :), ...
"LineWidth", 1.2, ...
"DisplayName", sprintf("Pump %g%%: total", pumpLevel));
end
title(amplifierNames(amplifierIndex), "Interpreter", "none");
xlabel("Wavelength [nm]", "Interpreter", "none");
ylabel("P_{out} [dBm]", "Interpreter", "tex");
grid on;
box on;
legend("Location", "best", "Interpreter", "none");
end
sgtitle(sprintf("Output power versus wavelength at P_{in} = %g dBm", inputPowerDb), ...
"Interpreter", "tex");
%% OSNR versus wavelength
figure("Name", "Amplifier OSNR versus wavelength", "Color", "w");
tiledlayout(1, 2, "TileSpacing", "compact", "Padding", "compact");
for amplifierIndex = 1:numel(lambdaWarehouses)
nexttile;
hold on;
colors = lines(numel(pumpLevelsLambda));
wh = lambdaWarehouses{amplifierIndex};
for pumpIndex = 1:numel(pumpLevelsLambda)
pumpLevel = pumpLevelsLambda(pumpIndex);
% osnr_osa is calculated from spectrum_osa during measurement.
osnrResults = wh.getStoValue("osnr_osa", inputPowerDb, wavelengthsNm, pumpLevel);
osnrDb = cellfun(@(result) result.corrected_dB, osnrResults);
plot(wavelengthsNm, osnrDb, "-", "Color", colors(pumpIndex, :), ...
"LineWidth", 1.5, ...
"DisplayName", sprintf("Pump %g%%", pumpLevel));
end
title(amplifierNames(amplifierIndex), "Interpreter", "none");
xlabel("Wavelength [nm]", "Interpreter", "none");
ylabel("OSNR [dB]", "Interpreter", "none");
grid on;
box on;
legend("Location", "best", "Interpreter", "none");
end
sgtitle(sprintf("OSNR versus wavelength at P_{in} = %g dBm", inputPowerDb), ...
"Interpreter", "tex");
%% Noise figure versus wavelength
% The stored pase_osa values are OSA power readings per resolution
% bandwidth. The measurement script configured the OSA to 0.1 nm.
hPlanck = 6.62607015e-34;
speedOfLight = 299792458;
lambdaMeters = wavelengthsNm(:) .* 1e-9;
opticalFrequencyHz = speedOfLight ./ lambdaMeters;
opticalNoiseBandwidthHz = speedOfLight ./ lambdaMeters.^2 .* (osaResolutionNm * 1e-9);
quantumNoisePowerW = hPlanck .* opticalFrequencyHz .* opticalNoiseBandwidthHz;
figure("Name", "Amplifier noise figure versus wavelength", "Color", "w");
tiledlayout(1, 2, "TileSpacing", "compact", "Padding", "compact");
for amplifierIndex = 1:numel(lambdaWarehouses)
nexttile;
hold on;
colors = lines(numel(pumpLevelsLambda));
wh = lambdaWarehouses{amplifierIndex};
for pumpIndex = 1:numel(pumpLevelsLambda)
pumpLevel = pumpLevelsLambda(pumpIndex);
signalOutputDb = wh.getStoValue("psig_osa", inputPowerDb, wavelengthsNm, pumpLevel);
aseOutputDb = wh.getStoValue("pase_osa", inputPowerDb, wavelengthsNm, pumpLevel);
gainDb = signalOutputDb - correctedInputPowerDb;
linearGain = 10.^(gainDb ./ 10);
asePowerW = 10.^((aseOutputDb - 30) ./ 10);
% Exact estimate including the amplified input shot-noise term.
noiseFactor = 1 ./ linearGain + ...
asePowerW ./ (linearGain .* quantumNoisePowerW);
noiseFigureDb = 10 .* log10(noiseFactor);
legendHandles(pumpIndex) = plot(wavelengthsNm(:), noiseFigureDb, "-", ...
"Color", colors(pumpIndex, :), ...
"LineWidth", 1.5, ...
"DisplayName", sprintf("Pump %g%%", pumpLevel));
end
title(amplifierNames(amplifierIndex), "Interpreter", "none");
xlabel("Wavelength [nm]", "Interpreter", "none");
ylabel("Noise figure [dB]", "Interpreter", "none");
grid on;
box on;
legend(legendHandles, "Location", "best", "Interpreter", "none");
end
sgtitle(sprintf("Estimated noise figure versus wavelength at P_{in} = %g dBm", inputPowerDb), ...
"Interpreter", "tex");
%% Final figure: Thorlabs gain versus wavelength
figure("Name", "Thorlabs PDFA gain versus wavelength", "Color", "w");
hold on;
colors = lines(numel(pumpLevelsLambda));
gainHandles = gobjects(numel(pumpLevelsLambda), 1);
for pumpIndex = 1:numel(pumpLevelsLambda)
pumpLevel = pumpLevelsLambda(pumpIndex);
signalOutputDb = whThorlabsLambda.getStoValue( ...
"psig_osa", inputPowerDb, wavelengthsNm, pumpLevel);
gainDb = signalOutputDb - correctedInputPowerDb;
gainHandles(pumpIndex) = plot(wavelengthsNm, gainDb, "-", ...
"Color", colors(pumpIndex, :), ...
"LineWidth", 1.6, ...
"DisplayName", sprintf("Pump %g%%: gain", pumpLevel));
end
title(sprintf("Thorlabs PDFA: gain versus wavelength at P_{in} = %g dBm", ...
inputPowerDb), "Interpreter", "tex");
xlabel("Wavelength [nm]", "Interpreter", "none");
ylabel("Gain [dB]", "Interpreter", "none");
grid on;
box on;
legend(gainHandles, "Location", "best", "Interpreter", "none");

View File

@@ -0,0 +1,87 @@
%AVAILABLE_OUTPUT_POWER_VS_WAVELENGTH Plot achieved laser output power.
%
% laser_outputpower_sweep_2_many_averages.mat contains the requested
% wavelength and laser-power values. The corresponding 20-trace
% measurements are stored in the local FIG file and are averaged here.
scriptDir = fileparts(mfilename("fullpath"));
parameterData = load(fullfile(scriptDir, "laser_outputpower_sweep_2_many_averages.mat"));
laserparams = parameterData.laserparams;
averagesPerSetpoint = 20;
measurementFigure = openfig( ...
fullfile(scriptDir, "available_output_power_exfo_t100_oband_20_averages.fig"), ...
"invisible");
measurementAxes = findall(measurementFigure, "Type", "axes");
measurementLines = findall(measurementAxes, "Type", "line");
numberOfSetpoints = numel(laserparams.laserpower);
numberOfWavelengths = numel(laserparams.lambda);
expectedLineCount = numberOfSetpoints * averagesPerSetpoint;
if numel(measurementLines) ~= expectedLineCount
close(measurementFigure);
error("Expected %d measurement traces, found %d.", ...
expectedLineCount, numel(measurementLines));
end
measuredLambda = measurementLines(1).XData(:).';
if numel(measuredLambda) ~= numberOfWavelengths || ...
any(abs(measuredLambda - laserparams.lambda) > 1e-9)
close(measurementFigure);
error("Measurement wavelengths do not match laserparams.lambda.");
end
achievedPowerDbm = zeros(numberOfSetpoints, numberOfWavelengths);
for setpointIndex = 1:numberOfSetpoints
firstLine = (setpointIndex - 1) * averagesPerSetpoint + 1;
traces = zeros(averagesPerSetpoint, numberOfWavelengths);
for averageIndex = 1:averagesPerSetpoint
trace = measurementLines(firstLine + averageIndex - 1).YData(:).';
traces(averageIndex, :) = trace;
end
achievedPowerDbm(setpointIndex, :) = mean(traces, 1, "omitnan");
end
close(measurementFigure);
% The measurement traces are stored in descending output-power order.
% Sort them so they match the ascending laserparams.laserpower order.
[~, sortIndex] = sort(mean(achievedPowerDbm, 2, "omitnan"));
achievedPowerDbm = achievedPowerDbm(sortIndex, :);
figure("Name", "Available laser output power", "Color", "w");
hold on;
colors = lines(numberOfSetpoints);
annotationX = laserparams.lambda(end) - 1;
for setpointIndex = 1:numberOfSetpoints
yline(laserparams.laserpower(setpointIndex), ":", ...
"Color", [0.65 0.65 0.65], "HandleVisibility", "off");
plot(laserparams.lambda, achievedPowerDbm(setpointIndex, :), "-", ...
"Color", colors(setpointIndex, :), "LineWidth", 1.6, ...
"DisplayName", sprintf("Target %g dBm: measured", ...
laserparams.laserpower(setpointIndex)));
text(annotationX, laserparams.laserpower(setpointIndex), ...
sprintf("P_{set} = %g dBm", laserparams.laserpower(setpointIndex)), ...
"Color", colors(setpointIndex, :), "Interpreter", "tex", ...
"HorizontalAlignment", "right", "VerticalAlignment", "bottom", ...
"HandleVisibility", "off");
end
title("Available laser output power versus wavelength", "Interpreter", "none");
xlabel("Wavelength [nm]", "Interpreter", "none");
ylabel("Achieved output power [dBm]", "Interpreter", "none");
grid on;
box on;
legend("Location", "northwest", "Interpreter", "none");
% Export the same figure with the repository-local matlab2tikz version.
repoDir = fileparts(fileparts(scriptDir));
addpath(fullfile(repoDir, "Libs", "mat2tikz", "src"));
tikzOutputFile = fullfile(scriptDir, "available_output_power_vs_wavelength.tikz");
matlab2tikz(char(tikzOutputFile), ...
'showInfo', false, ...
'showHiddenStrings', true, ...
'extraAxisOptions', {'legend style={font=\footnotesize}'});

View File

@@ -0,0 +1,350 @@
% This file was created by matlab2tikz.
%
\definecolor{mycolor1}{rgb}{0.89412,0.10196,0.10980}%
\definecolor{mycolor2}{rgb}{0.21569,0.49412,0.72157}%
\definecolor{mycolor3}{rgb}{0.30196,0.68627,0.29020}%
\definecolor{mycolor4}{rgb}{0.59608,0.30588,0.63922}%
\definecolor{mycolor5}{rgb}{1.00000,0.49804,0.00000}%
\definecolor{mycolor6}{rgb}{0.12941,0.12941,0.12941}%
%
\begin{tikzpicture}
\begin{axis}[%
width=6.458in,
height=5.094in,
at={(1.083in,0.688in)},
scale only axis,
xmin=1260,
xmax=1360,
xlabel style={font=\color{mycolor6}},
xlabel={Wavelength [nm]},
ymin=-2,
ymax=12,
ylabel style={font=\color{mycolor6}},
ylabel={Achieved output power [dBm]},
axis background/.style={fill=white},
title style={font=\bfseries\color{mycolor6}},
title={Available laser output power versus wavelength},
xmajorgrids,
ymajorgrids,
grid style={dashed},
legend style={at={(0.03,0.97)}, anchor=north west, legend cell align=left, align=left},
legend style={font=\footnotesize}
]
\addplot [color=white!65!black, dotted, forget plot]
table[row sep=crcr]{%
1260 0\\
1360 0\\
};
\addplot [color=mycolor1, line width=1.6pt]
table[row sep=crcr]{%
1260 -0.00174427500000025\\
1262 0.0049530015000002\\
1264 0.00258243800000004\\
1266 0.047948272\\
1268 0.035907209\\
1270 0.039917379\\
1272 0.0499743345000001\\
1274 0.05230108\\
1276 0.0532826289999999\\
1278 0.0831857375000001\\
1280 0.0627413579999998\\
1282 0.0488831240000001\\
1284 0.0447253150000001\\
1286 0.0355354289999999\\
1288 -0.00749986499999977\\
1290 0.00176488450000019\\
1292 -0.01833745\\
1294 -0.0421456650000001\\
1296 -0.05284091\\
1298 -0.0608307149999999\\
1300 -0.0629399749999999\\
1302 -0.0542474850000001\\
1304 -0.0943844899999999\\
1306 -0.11916809\\
1308 -0.088293575\\
1310 -0.106838255\\
1312 -0.142186045\\
1314 -0.122909425\\
1316 -0.136120425\\
1318 -0.169754955\\
1320 -0.14376953\\
1322 -0.153371705\\
1324 -0.19083209\\
1326 -0.16585443\\
1328 -0.16738367\\
1330 -0.18946759\\
1332 -0.16708397\\
1334 -0.125356505\\
1336 -0.125601985\\
1338 -0.133365665\\
1340 -0.07764601\\
1342 -0.0349484850000002\\
1344 -0.054868575\\
1346 0.0190894669999999\\
1348 0.0787105715000002\\
1350 0.0655493160000001\\
1352 0.0920785395000001\\
1354 0.139883643\\
1356 0.1389105975\\
1358 0.15040566\\
1360 0.1795823165\\
};
\addlegendentry{Target 0 dBm: measured}
\node[above left, align=right, inner sep=0, font=\color{mycolor1}]
at (axis cs:1359,0) {$\text{P}_{\text{set}}\text{ = 0 dBm}$};
\addplot [color=white!65!black, dotted, forget plot]
table[row sep=crcr]{%
1260 3\\
1360 3\\
};
\addplot [color=mycolor2, line width=1.6pt]
table[row sep=crcr]{%
1260 3.001012505\\
1262 3.006422054\\
1264 3.008953025\\
1266 3.0601149535\\
1268 3.0454856245\\
1270 3.047445995\\
1272 3.060657117\\
1274 3.071385682\\
1276 3.06080489\\
1278 3.079472899\\
1280 3.061335045\\
1282 3.050070503\\
1284 3.055032579\\
1286 3.042449593\\
1288 3.0135694995\\
1290 3.0167014675\\
1292 3.0036652445\\
1294 2.9777426285\\
1296 2.9726872215\\
1298 2.9398711715\\
1300 2.916202023\\
1302 2.927803093\\
1304 2.904134328\\
1306 2.888280249\\
1308 2.8835689165\\
1310 2.882305029\\
1312 2.8897452415\\
1314 2.8942784325\\
1316 2.885521624\\
1318 2.8517907965\\
1320 2.843298358\\
1322 2.855958141\\
1324 2.8360024095\\
1326 2.82446217\\
1328 2.8612967225\\
1330 2.8202706295\\
1332 2.8185087195\\
1334 2.8657736085\\
1336 2.859368219\\
1338 2.8706378305\\
1340 2.9512009455\\
1342 2.9559564765\\
1344 2.9683251885\\
1346 3.0288653\\
1348 3.082960566\\
1350 3.0701560325\\
1352 3.1103352335\\
1354 3.148218037\\
1356 3.137503138\\
1358 3.1635465145\\
1360 3.1890315515\\
};
\addlegendentry{Target 3 dBm: measured}
\node[above left, align=right, inner sep=0, font=\color{mycolor2}]
at (axis cs:1359,3) {$\text{P}_{\text{set}}\text{ = 3 dBm}$};
\addplot [color=white!65!black, dotted, forget plot]
table[row sep=crcr]{%
1260 6\\
1360 6\\
};
\addplot [color=mycolor3, line width=1.6pt]
table[row sep=crcr]{%
1260 5.987181772\\
1262 6.01080855\\
1264 6.0196549575\\
1266 6.0658400605\\
1268 6.0468314485\\
1270 6.0572872865\\
1272 6.0622857945\\
1274 6.0605621715\\
1276 6.056511485\\
1278 6.075664353\\
1280 6.0551906125\\
1282 6.0599489405\\
1284 6.0755984025\\
1286 6.0270936715\\
1288 6.0354817275\\
1290 6.022740339\\
1292 6.002313234\\
1294 5.980794718\\
1296 5.976163088\\
1298 5.950657563\\
1300 5.9223667825\\
1302 5.935413683\\
1304 5.9015896255\\
1306 5.8993098045\\
1308 5.891102095\\
1310 5.881858564\\
1312 5.8940111635\\
1314 5.903941552\\
1316 5.888497728\\
1318 5.857725498\\
1320 5.8567185045\\
1322 5.861671574\\
1324 5.8387118195\\
1326 5.833706753\\
1328 5.868907057\\
1330 5.8273320035\\
1332 5.8389400305\\
1334 5.876151523\\
1336 5.8620020595\\
1338 5.892156925\\
1340 5.9586369025\\
1342 5.9567924075\\
1344 5.982374356\\
1346 6.0375119525\\
1348 6.0890175215\\
1350 6.086149783\\
1352 6.1184705185\\
1354 6.1594823235\\
1356 6.1509530645\\
1358 6.180046173\\
1360 6.201404003\\
};
\addlegendentry{Target 6 dBm: measured}
\node[above left, align=right, inner sep=0, font=\color{mycolor3}]
at (axis cs:1359,6) {$\text{P}_{\text{set}}\text{ = 6 dBm}$};
\addplot [color=white!65!black, dotted, forget plot]
table[row sep=crcr]{%
1260 8\\
1360 8\\
};
\addplot [color=mycolor4, line width=1.6pt]
table[row sep=crcr]{%
1260 6.031553038\\
1262 6.475926562\\
1264 6.9678239935\\
1266 7.3812834605\\
1268 7.6372069015\\
1270 7.9191694625\\
1272 8.0705430165\\
1274 8.075813651\\
1276 8.083409099\\
1278 8.084946346\\
1280 8.070698563\\
1282 8.072232142\\
1284 8.0797347955\\
1286 8.0427192865\\
1288 8.0515738335\\
1290 8.0293537885\\
1292 8.003326676\\
1294 7.9976869805\\
1296 7.99008813\\
1298 7.9492751595\\
1300 7.9346835125\\
1302 7.9539901195\\
1304 7.9017661065\\
1306 7.911127979\\
1308 7.907037724\\
1310 7.896197244\\
1312 7.907979363\\
1314 7.9250423625\\
1316 7.8998658815\\
1318 7.869535715\\
1320 7.8711403195\\
1322 7.8730818405\\
1324 7.845561009\\
1326 7.85929528\\
1328 7.8721726815\\
1330 7.843708335\\
1332 7.872582495\\
1334 7.8930589925\\
1336 7.882136431\\
1338 7.920374461\\
1340 7.9612389245\\
1342 7.968047289\\
1344 8.002771872\\
1346 8.054855834\\
1348 8.0969126465\\
1350 8.1158346705\\
1352 8.1537042675\\
1354 8.181075273\\
1356 8.1743607905\\
1358 8.184566916\\
1360 8.231713996\\
};
\addlegendentry{Target 8 dBm: measured}
\node[above left, align=right, inner sep=0, font=\color{mycolor4}]
at (axis cs:1359,8) {$\text{P}_{\text{set}}\text{ = 8 dBm}$};
\addplot [color=white!65!black, dotted, forget plot]
table[row sep=crcr]{%
1260 10\\
1360 10\\
};
\addplot [color=mycolor5, line width=1.6pt]
table[row sep=crcr]{%
1260 5.7017284215\\
1262 6.1986312315\\
1264 6.7501614225\\
1266 7.1907112785\\
1268 7.478677345\\
1270 7.8060430325\\
1272 8.1067171935\\
1274 8.276471823\\
1276 8.4638177645\\
1278 8.6111332305\\
1280 8.743436757\\
1282 8.832386923\\
1284 8.91278028\\
1286 9.04439012185\\
1288 9.04376373615\\
1290 9.0798988258\\
1292 9.03395511625\\
1294 8.9921718725\\
1296 9.0235644005\\
1298 8.984754639\\
1300 9.14610934125\\
1302 9.13106503035\\
1304 9.2086507995\\
1306 9.19881023815\\
1308 9.265719913\\
1310 9.09915279415\\
1312 9.26554581525\\
1314 9.52831931645\\
1316 9.6669975241\\
1318 9.7107443498\\
1320 9.7898754663\\
1322 9.83403242775\\
1324 9.8726901865\\
1326 9.8822262349\\
1328 9.88426970855\\
1330 9.84195532625\\
1332 9.87982448925\\
1334 9.90844454246\\
1336 9.89706352885\\
1338 9.9274772081\\
1340 9.991583686247\\
1342 9.983836191075\\
1344 10.01432601913\\
1346 10.08652218347\\
1348 10.1116048911\\
1350 10.1206856349\\
1352 10.1690383129\\
1354 10.191050599\\
1356 10.19157208575\\
1358 10.21642409705\\
1360 10.24591036845\\
};
\addlegendentry{Target 10 dBm: measured}
\node[above left, align=right, inner sep=0, font=\color{mycolor5}]
at (axis cs:1359,10) {$\text{P}_{\text{set}}\text{ = 10 dBm}$};
\end{axis}
\end{tikzpicture}%

View File

@@ -1,6 +1,6 @@
wh_aeon = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Lab_analysis\aeon_soa_measurement_lambda_plaser_pump.mat");
wh_aeon = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Lab_oband_laser_and_amp_analysis\aeon_soa_measurement_lambda_plaser_pump.mat");
wh_aeon = wh_aeon.wh;
wh_thor = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Lab_analysis\thorlabs_pdfa_measurement_lambda_plaser_pump.mat");
wh_thor = load("C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\Lab_oband_laser_and_amp_analysis\thorlabs_pdfa_measurement_lambda_plaser_pump.mat");
wh_thor = wh_thor.wh;
% Silas' custom "warehouse" datatype
@@ -36,6 +36,7 @@ for p = 1:numel(plasers)
plot(wavelength_osa,spectrum_osa,'DisplayName',sprintf('P_{in}: %d dB ',plasers(p)),'Color',cols(ccnt+1,:));
ylim([-60, 20]);
beautifyBERplot("logscale",false,"setmarkers",0,"setcolors",0);
end

View File

@@ -0,0 +1,193 @@
if 1
% A) RUN FULL LOOP
M_format = [2,4,6,8];
snr = 10:25;
else
% B) RUN FOR DEBUG AND TEST
M_format = 4;
snr = 20;
end
for m = 1:length(M_format)
% --- Parameters ---
M = M_format(m); % PAM order (e.g., 2,4,8)
Nsym = 1e5; % number of symbols
h = [1, 0.5, 0.2]; % Impulse response to remove
b = log2(M);
if M == 6 b = 5; end
rng(1);
bits_tx = logical(randi([0 1], Nsym, b, 'uint8'));
tx_symbols = pammap(bits_tx,M);
if M == 6
states = unique(tx_symbols);
pam6transitions = combvec(states',states')'; % pam6transitions =
bitmapping = pamdemap(reshape(pam6transitions',1,[])',M);
else
bitmapping = pamdemap(unique(tx_symbols),M);
end
scaling = sqrt(sum(unique(tx_symbols).^2)/numel(unique(tx_symbols)));
tx_symbols = tx_symbols ./ scaling;
% apply impulse response to signal
y_filt = filter(h, 1, tx_symbols);
for s = 1:length(snr)
% apply noise
y = awgn(y_filt,snr(s),"measured",1);
y = Electricalsignal(y);
% apply ml-MLSE
adaptive_mu = 0;
mu_lms = 0.15;
ml_mlse_equalizer = ML_MLSE("epochs_tr",50,"epochs_dd",1,"len_tr",length(y)/2,...
"mu_dd",mu_lms,"mu_tr",mu_lms,"order",11,"sps",1,...
"L",2,"delta",4,"adaptive_mu",adaptive_mu);
[ml_mlse_estimate,~] = ml_mlse_equalizer.process(y,tx_symbols);
rx_symbols = ml_mlse_estimate .* scaling;
bits_rx = pamdemap(rx_symbols,M);
BER_ml(m,s) = nnz(bits_tx ~= bits_rx) / numel(bits_tx);
fprintf('BER = %.2e \n', BER_ml(m,s));
% apply bcjr
BCJR = MLSE("DIR",h,"duobinary_output",0,"M",M,"trellis_states",unique(tx_symbols));
[viterbi_estimate,LLR,GMI(m,s)] = BCJR.process(y,tx_symbols,bits_tx,bitmapping);
% decode LLR's
bits_LLR = LLR > 0;
% demap viterbi symbols sequence
rx_symbols = viterbi_estimate .* scaling;
bits_rx = pamdemap(rx_symbols,M);
% BER calc
BER_vit(m,s) = nnz(bits_tx ~= bits_LLR) / numel(bits_tx);
fprintf('BER LLR = %.2e \n', BER_vit(m,s));
BER_llr(m,s) = nnz(bits_tx ~= bits_rx) / numel(bits_tx);
fprintf('BER = %.2e \n', BER_llr(m,s));
end
end
%%
figure();hold on
for m = 1:length(M_format)
p=plot(snr,BER_llr(m,:),'DisplayName',sprintf('Viterbi: PAM %d',M_format(m)));
plot(snr,BER_ml(m,:),'DisplayName',sprintf('ML-Based: PAM %d',M_format(m)),'LineStyle',':','Color',p.Color);
end
ylabel('BER');
xlabel('SNR')
title('BER vs. SNR');
set(gca, 'XScale', 'linear', ...
'YScale', 'log', ...
'TickLabelInterpreter', 'latex', ...
'FontSize', 11);
%%
figure();hold on
for m = 1:length(M_format)
plot(snr,GMI(m,:),'DisplayName',sprintf('GMI PAM %d',M_format(m)))
end
ylabel('GMI');
xlabel('SNR')
title('GMI vs. SNR');
set(gca, 'XScale', 'linear', ...
'YScale', 'linear', ...
'TickLabelInterpreter', 'latex', ...
'FontSize', 11);
function symbols = pammap(bits,M)
bits = logical(bits);
if M == 2
symbols = bits;
elseif M == 4
symbols= 2*bits(:,1) + (bits(:,1)==bits(:,2));
symbols=2*symbols-3;
elseif M == 6
m = 1;
if size(bits,2)>size(bits,1)
bits = bits'; %vector aufrecht stellen
end
bits = reshape(bits',1,[])';
thres = [-3 5;-1 5;-3 -5;-1 -5;-5 3;-5 1;-5 -3;-5 -1;-1 3;-1 1;-1 -3;-1 -1;-3 3;-3 1;-3 -3;-3 -1;3 5;1 5;3 -5;1 -5;5 3;5 1;5 -3;5 -1;1 3;1 1;1 -3;1 -1;3 3;3 1;3 -3;3 -1];
% LUT based mapping
for k = 1:5:fix(length(bits)/5)*5
symbols(m:m+1,1) = thres(bin2dec(int2str(bits(k:k+4)'))+1,:);
m = m+2;
end
elseif M == 8
x1 = bits(:,1);
x2 = (bits(:,1)==bits(:,3));
x3 = x2~=bits(:,2);
symbols = 4*x1 + 2*x2 + x3;
symbols=2*symbols-7;
end
end
function bits = pamdemap(symbols,M)
if M == 2
thres=0;
elseif M == 4
thres=[-2,0,2];
elseif M == 6
thres = [-3 5;-1 5;-3 -5;-1 -5;-5 3;-5 1;-5 -3;-5 -1;-1 3;-1 1;-1 -3;-1 -1;-3 3;-3 1;-3 -3;-3 -1;3 5;1 5;3 -5;1 -5;5 3;5 1;5 -3;5 -1;1 3;1 1;1 -3;1 -1;3 3;3 1;3 -3;3 -1];
elseif M == 8
thres=-6:2:6;
end
if M ~= 6
symbols = symbols';
a = squeeze(repmat(real(symbols),[1 1 length(thres)])); %Eingangssignal in 3 spalten
b = squeeze(repmat(reshape(thres(:).',[1 1 length(thres)]),[1 length(symbols) 1])); %Threshold in 3 Spalten
comp_real = a > b; %check for each symbol/ sampling if it exeeds the obj.thresholdseshold 1, 2 or 3
comp_real=repmat(real(symbols),[1 1 length(thres)]) > repmat(reshape(thres(:).',[1 1 length(thres)]),[1 length(symbols) 1]);
s1=size(comp_real,1);
s2=size(comp_real,2);
end
if M == 2
data_out=abs(comp_real(:,:,1));
elseif M == 4
data_out=[comp_real(:,:,2); ones(s1,s2) - comp_real(:,:,1) + comp_real(:,:,3)];
elseif M == 6
if size(symbols,2) > 1
symbols = symbols.';
end
if length(symbols)/2 ~= round(length(symbols)/2)
symbols = [symbols;0];
end
m = 1;
for n = 1:2:length(symbols)
dist = sqrt((symbols(n)-thres(:,1)).^2+(symbols(n+1)-thres(:,2)).^2);
[~,dd_idx] = min(dist);
% dec_out(n:n+1) = LUT(dd_idx,:);
data_out(m:m+4) = bitget(dd_idx-1,5:-1:1);
m = m+5;
end
data_out = reshape(data_out',5,[]);
elseif M == 8
data_out=[comp_real(:,:,4);
comp_real(:,:,1)-comp_real(:,:,3)+comp_real(:,:,5)-comp_real(:,:,7);
1-comp_real(:,:,2)+comp_real(:,:,6)];
end
bits = data_out';
end

Binary file not shown.