diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index f88eadb..608510f 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -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); @@ -526,10 +526,10 @@ classdef Signal end y_margin = 0.05 * y_range; - ylim([y_min - y_margin, y_max + y_margin]); - - yticks(-200:10:200); - grid on; + ylim([y_min - y_margin, y_max + y_margin]); + + yticks(-200:10:200); + grid on; % Add legend if not already present if isempty(get(gca, 'Legend')) @@ -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))); diff --git a/Classes/01_transmit/ChannelFreqResp.m b/Classes/01_transmit/ChannelFreqResp.m index 21c068e..c0d4c11 100644 --- a/Classes/01_transmit/ChannelFreqResp.m +++ b/Classes/01_transmit/ChannelFreqResp.m @@ -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(); diff --git a/Classes/04_DSP/Equalizer/ML_MLSE_DUOBINARY.m b/Classes/04_DSP/Equalizer/ML_MLSE_DUOBINARY.m index 2a9c206..d623a71 100644 --- a/Classes/04_DSP/Equalizer/ML_MLSE_DUOBINARY.m +++ b/Classes/04_DSP/Equalizer/ML_MLSE_DUOBINARY.m @@ -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 diff --git a/Datatypes/clr.m b/Datatypes/clr.m index bda20d1..2899191 100644 --- a/Datatypes/clr.m +++ b/Datatypes/clr.m @@ -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], ... diff --git a/Functions/EQ_blocks/duobinary_target.m b/Functions/EQ_blocks/duobinary_target.m index 2e42e75..4edd8b2 100644 --- a/Functions/EQ_blocks/duobinary_target.m +++ b/Functions/EQ_blocks/duobinary_target.m @@ -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); diff --git a/Functions/EQ_blocks/vnle_postfilter_mlse.m b/Functions/EQ_blocks/vnle_postfilter_mlse.m index 5bcf959..6616986 100644 --- a/Functions/EQ_blocks/vnle_postfilter_mlse.m +++ b/Functions/EQ_blocks/vnle_postfilter_mlse.m @@ -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; diff --git a/Functions/EQ_recipes/dsp_400g_recipe.m b/Functions/EQ_recipes/dsp_400g_recipe.m index 751534c..5fe580e 100644 --- a/Functions/EQ_recipes/dsp_400g_recipe.m +++ b/Functions/EQ_recipes/dsp_400g_recipe.m @@ -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) diff --git a/Functions/EQ_visuals/showEQNoisePSD.m b/Functions/EQ_visuals/showEQNoisePSD.m index ba613da..26fda38 100644 --- a/Functions/EQ_visuals/showEQNoisePSD.m +++ b/Functions/EQ_visuals/showEQNoisePSD.m @@ -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 diff --git a/Functions/EQ_visuals/showLevelHistogram.m b/Functions/EQ_visuals/showLevelHistogram.m index c86828d..15395bd 100644 --- a/Functions/EQ_visuals/showLevelHistogram.m +++ b/Functions/EQ_visuals/showLevelHistogram.m @@ -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]); diff --git a/Functions/Job_Processing/loadAndSyncRunSignals.m b/Functions/Job_Processing/loadAndSyncRunSignals.m index d1e03fa..d1fda35 100644 --- a/Functions/Job_Processing/loadAndSyncRunSignals.m +++ b/Functions/Job_Processing/loadAndSyncRunSignals.m @@ -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', ... diff --git a/Functions/Metrics/calc_snr.m b/Functions/Metrics/calc_snr.m index ff25cd9..6151e89 100644 --- a/Functions/Metrics/calc_snr.m +++ b/Functions/Metrics/calc_snr.m @@ -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); - - % 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)); - + snr_per_level = zeros(numel(levels),1); + + 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)); + + 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 diff --git a/Functions/frequency_wavelength/calcWavelengthPlan.m b/Functions/frequency_wavelength/calcWavelengthPlan.m new file mode 100644 index 0000000..dea849f --- /dev/null +++ b/Functions/frequency_wavelength/calcWavelengthPlan.m @@ -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 + + diff --git a/Functions/frequency_wavelength/df2dw.m b/Functions/frequency_wavelength/df2dw.m new file mode 100644 index 0000000..1aefb74 --- /dev/null +++ b/Functions/frequency_wavelength/df2dw.m @@ -0,0 +1,4 @@ +function dw = df2dw(df,f) + %delta frequency to delta wavelength + dw = 299792458 .* df ./ f.^2; +end \ No newline at end of file diff --git a/Functions/frequency_wavelength/dw2df.m b/Functions/frequency_wavelength/dw2df.m new file mode 100644 index 0000000..d9f1fcc --- /dev/null +++ b/Functions/frequency_wavelength/dw2df.m @@ -0,0 +1,4 @@ +function df = dw2df(dw,w) + %delta wavelength to delta frequency + df = 299792458 .* dw ./ w.^2; +end \ No newline at end of file diff --git a/Functions/frequency_wavelength/getSweepWavelengths.m b/Functions/frequency_wavelength/getSweepWavelengths.m new file mode 100644 index 0000000..1c848f4 --- /dev/null +++ b/Functions/frequency_wavelength/getSweepWavelengths.m @@ -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 \ No newline at end of file diff --git a/Functions/frequency_wavelength/hz2nm.m b/Functions/frequency_wavelength/hz2nm.m new file mode 100644 index 0000000..43f8ffc --- /dev/null +++ b/Functions/frequency_wavelength/hz2nm.m @@ -0,0 +1,4 @@ +function nm = hz2nm(hz) + m = (299792458 ./ hz); %wavelength in meter + nm = m .* 10^9; +end diff --git a/Functions/frequency_wavelength/nm2hz.m b/Functions/frequency_wavelength/nm2hz.m new file mode 100644 index 0000000..71c467c --- /dev/null +++ b/Functions/frequency_wavelength/nm2hz.m @@ -0,0 +1,4 @@ +function hz = nm2hz(nm) + wavelen_in_m = nm.* 1e-9; + hz = (299792458 ./ wavelen_in_m); %frequency in Hertz +end \ No newline at end of file diff --git a/Functions/frequency_wavelength/nm2thz.m b/Functions/frequency_wavelength/nm2thz.m new file mode 100644 index 0000000..b56e5e4 --- /dev/null +++ b/Functions/frequency_wavelength/nm2thz.m @@ -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 diff --git a/Functions/frequency_wavelength/thz2nm.m b/Functions/frequency_wavelength/thz2nm.m new file mode 100644 index 0000000..ced382a --- /dev/null +++ b/Functions/frequency_wavelength/thz2nm.m @@ -0,0 +1,4 @@ +function nm = thz2nm(thz) + freq_in_hz = thz.* 1e12; + nm = (299792458 ./ freq_in_hz) .*1e9; %wavelength in nanometer +end \ No newline at end of file diff --git a/Functions/mat2tikz_improved.m b/Functions/mat2tikz_improved.m index bc0b040..2cc5087 100644 --- a/Functions/mat2tikz_improved.m +++ b/Functions/mat2tikz_improved.m @@ -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', ... diff --git a/Functions/saveDirectoryInfo.m b/Functions/saveDirectoryInfo.m new file mode 100644 index 0000000..24f2b0d --- /dev/null +++ b/Functions/saveDirectoryInfo.m @@ -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 \ No newline at end of file diff --git a/Tests/04_DSP/Equalizer/ML_MLSE_DUOBINARY_test.m b/Tests/04_DSP/Equalizer/ML_MLSE_DUOBINARY_test.m index 0c1786c..2f1f7e3 100644 --- a/Tests/04_DSP/Equalizer/ML_MLSE_DUOBINARY_test.m +++ b/Tests/04_DSP/Equalizer/ML_MLSE_DUOBINARY_test.m @@ -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); diff --git a/directory_info_sioe_labor_20260720_114654.mat b/directory_info_sioe_labor_20260720_114654.mat new file mode 100644 index 0000000..ccc495f Binary files /dev/null and b/directory_info_sioe_labor_20260720_114654.mat differ diff --git a/projects/AWG_output_power_simulation/output_power.m b/projects/AWG_output_power_simulation/output_power.m index a31acdc..608839b 100644 --- a/projects/AWG_output_power_simulation/output_power.m +++ b/projects/AWG_output_power_simulation/output_power.m @@ -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); diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_EYES.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_EYES.m index 626862c..bcbdf4d 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_EYES.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_EYES.m @@ -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 = ; - - -% ============================================================ -% 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]; diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_Spectra.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_Spectra.m index 5895cf3..750950f 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_Spectra.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_Spectra.m @@ -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,85 +18,102 @@ 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); - 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); + %%% 2) PLOT PREEMPH. TX SIGNAL + if dbmode == 0 + maxamp = -37; + precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs); - Digi_sig_pre = Digi_sig_pre.resample("fs_out",fdac); + 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.normalize("mode","rms"); + Digi_sig_pre = Digi_sig_pre.resample("fs_out",256e9); - 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= Digi_sig_pre.normalize("mode","rms"); + 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 fp = QueryFilter(); M = 4; - fp.where('Runs', 'pam_level','EQUALS', M); + 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', 'wavelength','EQUALS', 1310); %1327.4 - fp.where('Runs', 'db_mode','EQUALS', dbmode); - fp.where('Runs', 'rop_attenuation','EQUAL', 0); - + 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); + [dataTable,~] = database.queryDB(fp, database.getTableFieldNames('Runs')); - + dataTable = queryRunid(dataTable.run_id, database); fsym = dataTable.symbolrate; M = double(dataTable.pam_level); - + % Load and Sync signal data from DB [Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncRunSignals(dataTable, dsp_options); - + % Preprocess signal Scpe_sig = preprocessSignal(Scpe_cell{1}, Symbols, fsym); 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,48 +145,48 @@ 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; % Do EQ and find alpha's len_tr = 4096*2; - + ffe_order = [50, 5, 5]; dfe_order = [0, 0, 0]; pf_ncoeffs = 1; mu_ffe = [0.0001, 0.0008, 0.001]; mu_dfe = 0.0004; mu_dc = 0.005; - + %%% FULL RESP TARGET 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); pf_1 = Postfilter("ncoeff",1,"useBurg",1); - + [eq_signal_sd, eq_noise] = eq_.process(Scpe_sig, Symbols); % eq_noise.signal = eq_noise.signal - mean(eq_noise.signal); % eq_noise = eq_noise.normalize("mode","rms"); - + [mlse_sig_sd,whitened_noise] = pf_1.process(eq_signal_sd, eq_noise); fig = figure(fignum+dbmode+10); hold on - + [h, w] = freqz(1, pf_1.coefficients, length(eq_noise), "whole", eq_noise.fs); h = h / max(abs(h)); % Normalize the filter response w_ = (w - eq_noise.fs / 2); - + %%% 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); % db_noise = db_noise.normalize("mode","rms"); - + %%% 1-3) Plot EQ Noise EEN figure(fignum+dbmode+10) eq_noise.spectrum("displayname", 'Noise', "fignum", fignum+dbmode+10, "normalizeTo0dB", 0,"color",clr.Set1.green,"normalizeToDC",0,"addDCoffset",0); @@ -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 diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_WAVELENGTH_VS_BAUDRATE.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_WAVELENGTH_VS_BAUDRATE.m deleted file mode 100644 index d987359..0000000 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_WAVELENGTH_VS_BAUDRATE.m +++ /dev/null @@ -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 diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/generate_spectrum_plots.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/generate_spectrum_plots.m index 9e413af..41c936c 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/generate_spectrum_plots.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/generate_spectrum_plots.m @@ -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 diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/power_vs_wavelength.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/power_vs_wavelength.m index b7e38ad..dabbd15 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/power_vs_wavelength.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/power_vs_wavelength.m @@ -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); @@ -93,4 +94,29 @@ end 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); \ No newline at end of file +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); \ No newline at end of file diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/bias_evaluation.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/bias_evaluation.m index c042156..b1a4b91 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/bias_evaluation.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/bias_evaluation.m @@ -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)); diff --git a/projects/Diss/400G_revisit/FIGURE_EQ_NOISE_VS_BAUDRATE.m b/projects/Diss/400G_revisit/FIGURE_EQ_NOISE_VS_BAUDRATE.m new file mode 100644 index 0000000..54d0d88 --- /dev/null +++ b/projects/Diss/400G_revisit/FIGURE_EQ_NOISE_VS_BAUDRATE.m @@ -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 \ No newline at end of file diff --git a/projects/Diss/400G_revisit/PLOT_BER_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m b/projects/Diss/400G_revisit/PLOT_BER_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m index 7522ccf..dc57fc4 100644 --- a/projects/Diss/400G_revisit/PLOT_BER_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m +++ b/projects/Diss/400G_revisit/PLOT_BER_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m @@ -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); diff --git a/projects/Diss/400G_revisit/PLOT_BER_VS_ALGO.m b/projects/Diss/400G_revisit/PLOT_BER_VS_ALGO.m index eed6bc1..34776eb 100644 --- a/projects/Diss/400G_revisit/PLOT_BER_VS_ALGO.m +++ b/projects/Diss/400G_revisit/PLOT_BER_VS_ALGO.m @@ -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"); diff --git a/projects/Diss/400G_revisit/PLOT_DUOBINARY_DETECTION_BER_VS_BITRATE.m b/projects/Diss/400G_revisit/PLOT_DUOBINARY_DETECTION_BER_VS_BITRATE.m index 9c52776..f2a2890 100644 --- a/projects/Diss/400G_revisit/PLOT_DUOBINARY_DETECTION_BER_VS_BITRATE.m +++ b/projects/Diss/400G_revisit/PLOT_DUOBINARY_DETECTION_BER_VS_BITRATE.m @@ -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 diff --git a/projects/Diss/400G_revisit/PLOT_EYES_400G_REVISIT.m b/projects/Diss/400G_revisit/PLOT_EYES_400G_REVISIT.m new file mode 100644 index 0000000..c43d45c --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_EYES_400G_REVISIT.m @@ -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 + 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 diff --git a/projects/Diss/400G_revisit/PLOT_MEMORYLESS_DBTGT_VS_VNLE_DB_MLSE.m b/projects/Diss/400G_revisit/PLOT_MEMORYLESS_DBTGT_VS_VNLE_DB_MLSE.m new file mode 100644 index 0000000..e2d7351 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_MEMORYLESS_DBTGT_VS_VNLE_DB_MLSE.m @@ -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 + 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 + 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 diff --git a/projects/Diss/400G_revisit/PLOT_MLSE_N_TAP_BER_BY_DB_MODE.m b/projects/Diss/400G_revisit/PLOT_MLSE_N_TAP_BER_BY_DB_MODE.m new file mode 100644 index 0000000..671f382 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_MLSE_N_TAP_BER_BY_DB_MODE.m @@ -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 + 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 + 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 + dbModeCol(end+1, 1) = double(runTable.db_mode(rowIdx)); %#ok + pfCol(end+1, 1) = pf; %#ok + bitrateCol(end+1, 1) = double(runTable.bitrate(rowIdx)) .* 1e-9; %#ok + berCol(end+1, 1) = ber; %#ok + berpCol(end+1, 1) = berp; %#ok +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 + dbModeCol(end+1, 1) = double(runTable.db_mode(rowIdx)); %#ok + bitrateCol(end+1, 1) = double(runTable.bitrate(rowIdx)) .* 1e-9; %#ok + berCol(end+1, 1) = ber; %#ok + berpCol(end+1, 1) = berp; %#ok +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 diff --git a/projects/Diss/400G_revisit/PLOT_NGMI_AIR_GMI_RATES_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m b/projects/Diss/400G_revisit/PLOT_NGMI_AIR_GMI_RATES_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m new file mode 100644 index 0000000..e387475 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_NGMI_AIR_GMI_RATES_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m @@ -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 + 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 diff --git a/projects/Diss/400G_revisit/PLOT_REPROCESSED_DSP400G_WAREHOUSE_ANALYSIS.m b/projects/Diss/400G_revisit/PLOT_REPROCESSED_DSP400G_WAREHOUSE_ANALYSIS.m new file mode 100644 index 0000000..b69a898 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_REPROCESSED_DSP400G_WAREHOUSE_ANALYSIS.m @@ -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 diff --git a/projects/Diss/400G_revisit/PLOT_ROP_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m b/projects/Diss/400G_revisit/PLOT_ROP_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m new file mode 100644 index 0000000..4f84fc1 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_ROP_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m @@ -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 + 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 diff --git a/projects/Diss/400G_revisit/PLOT_SNR_BER_DUOBINARY_PARTIAL_RESPONSE.m b/projects/Diss/400G_revisit/PLOT_SNR_BER_DUOBINARY_PARTIAL_RESPONSE.m new file mode 100644 index 0000000..ef8a0a9 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_SNR_BER_DUOBINARY_PARTIAL_RESPONSE.m @@ -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 + pamCol(end+1, 1) = double(runTable.pam_level(runRow)); %#ok + dbModeCol(end+1, 1) = double(runTable.db_mode(runRow)); %#ok + symbolrateCol(end+1, 1) = double(runTable.symbolrate(runRow)) * 1e-9; %#ok + snrCol(end+1, 1) = snr; %#ok + berCol(end+1, 1) = ber; %#ok + 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 diff --git a/projects/Diss/400G_revisit/PLOT_WAREHOUSE_BAUDRATE_ROP_ANALYSIS.m b/projects/Diss/400G_revisit/PLOT_WAREHOUSE_BAUDRATE_ROP_ANALYSIS.m new file mode 100644 index 0000000..f3a76c5 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_WAREHOUSE_BAUDRATE_ROP_ANALYSIS.m @@ -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 + 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 + 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 diff --git a/projects/Diss/400G_revisit/PLOT_WAVELENGTH_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m b/projects/Diss/400G_revisit/PLOT_WAVELENGTH_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m new file mode 100644 index 0000000..a9657e2 --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_WAVELENGTH_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m @@ -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 + 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 diff --git a/projects/Diss/400G_revisit/REPLOT_OPT_FILT_ANALYSIS_BY_PAM.m b/projects/Diss/400G_revisit/REPLOT_OPT_FILT_ANALYSIS_BY_PAM.m new file mode 100644 index 0000000..df437e2 --- /dev/null +++ b/projects/Diss/400G_revisit/REPLOT_OPT_FILT_ANALYSIS_BY_PAM.m @@ -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); diff --git a/projects/Diss/400G_revisit/RUN_REPROCESS_BAUDRATE_FILES_DSP400G.m b/projects/Diss/400G_revisit/RUN_REPROCESS_BAUDRATE_FILES_DSP400G.m new file mode 100644 index 0000000..6ec94b3 --- /dev/null +++ b/projects/Diss/400G_revisit/RUN_REPROCESS_BAUDRATE_FILES_DSP400G.m @@ -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 + "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 diff --git a/projects/Diss/400G_revisit/RX_sprectra.m b/projects/Diss/400G_revisit/RX_sprectra.m new file mode 100644 index 0000000..bb39f77 --- /dev/null +++ b/projects/Diss/400G_revisit/RX_sprectra.m @@ -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); \ No newline at end of file diff --git a/projects/Diss/400G_revisit/TX_spectra.m b/projects/Diss/400G_revisit/TX_spectra.m new file mode 100644 index 0000000..df869bd --- /dev/null +++ b/projects/Diss/400G_revisit/TX_spectra.m @@ -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); diff --git a/projects/Diss/400G_revisit/auswertung_baudrate.mlx b/projects/Diss/400G_revisit/auswertung_baudrate.mlx new file mode 100644 index 0000000..78b0807 Binary files /dev/null and b/projects/Diss/400G_revisit/auswertung_baudrate.mlx differ diff --git a/projects/Diss/400G_revisit/investigate_400g_algorithms.m b/projects/Diss/400G_revisit/investigate_400g_algorithms.m index 56d5856..a8d0dc6 100644 --- a/projects/Diss/400G_revisit/investigate_400g_algorithms.m +++ b/projects/Diss/400G_revisit/investigate_400g_algorithms.m @@ -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 + bitrateCol(end+1, 1) = double(bitrate) * 1e-9; %#ok + berpCol(end+1, 1) = double(berp); %#ok +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); diff --git a/projects/Diss/400G_revisit/mlse_n_tap_pam4.mat b/projects/Diss/400G_revisit/mlse_n_tap_pam4.mat new file mode 100644 index 0000000..56ed143 Binary files /dev/null and b/projects/Diss/400G_revisit/mlse_n_tap_pam4.mat differ diff --git a/projects/Diss/400G_revisit/mlse_n_tap_pam4_2km.mat b/projects/Diss/400G_revisit/mlse_n_tap_pam4_2km.mat new file mode 100644 index 0000000..de2375d Binary files /dev/null and b/projects/Diss/400G_revisit/mlse_n_tap_pam4_2km.mat differ diff --git a/projects/Diss/400G_revisit/opt_filt_analysis.fig b/projects/Diss/400G_revisit/opt_filt_analysis.fig new file mode 100644 index 0000000..c2a9b6b Binary files /dev/null and b/projects/Diss/400G_revisit/opt_filt_analysis.fig differ diff --git a/projects/Diss/400G_revisit/results_duobinary_eq_memoryless_2km_pam468.mat b/projects/Diss/400G_revisit/results_duobinary_eq_memoryless_2km_pam468.mat new file mode 100644 index 0000000..849919a Binary files /dev/null and b/projects/Diss/400G_revisit/results_duobinary_eq_memoryless_2km_pam468.mat differ diff --git a/projects/Diss/400G_revisit/results_snr_duobinary_and_partialresponse.mat b/projects/Diss/400G_revisit/results_snr_duobinary_and_partialresponse.mat new file mode 100644 index 0000000..a4d7793 Binary files /dev/null and b/projects/Diss/400G_revisit/results_snr_duobinary_and_partialresponse.mat differ diff --git a/projects/Diss/400G_revisit/results_snr_duobinary_and_partialresponse_pam468_10km.mat b/projects/Diss/400G_revisit/results_snr_duobinary_and_partialresponse_pam468_10km.mat new file mode 100644 index 0000000..3046bd6 Binary files /dev/null and b/projects/Diss/400G_revisit/results_snr_duobinary_and_partialresponse_pam468_10km.mat differ diff --git a/projects/Diss/400G_revisit/rop_vs_baudrate_sweep.mat b/projects/Diss/400G_revisit/rop_vs_baudrate_sweep.mat new file mode 100644 index 0000000..fc1c9cf Binary files /dev/null and b/projects/Diss/400G_revisit/rop_vs_baudrate_sweep.mat differ diff --git a/projects/Lab_oband_laser_and_amp_analysis/amplifier_gain_analysis.m b/projects/Lab_oband_laser_and_amp_analysis/amplifier_gain_analysis.m new file mode 100644 index 0000000..b3a90f2 --- /dev/null +++ b/projects/Lab_oband_laser_and_amp_analysis/amplifier_gain_analysis.m @@ -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"); diff --git a/projects/Lab_oband_laser_and_amp_analysis/available_output_power_vs_wavelength.m b/projects/Lab_oband_laser_and_amp_analysis/available_output_power_vs_wavelength.m new file mode 100644 index 0000000..3243073 --- /dev/null +++ b/projects/Lab_oband_laser_and_amp_analysis/available_output_power_vs_wavelength.m @@ -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}'}); diff --git a/projects/Lab_oband_laser_and_amp_analysis/available_output_power_vs_wavelength.tikz b/projects/Lab_oband_laser_and_amp_analysis/available_output_power_vs_wavelength.tikz new file mode 100644 index 0000000..414a97e --- /dev/null +++ b/projects/Lab_oband_laser_and_amp_analysis/available_output_power_vs_wavelength.tikz @@ -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}% \ No newline at end of file diff --git a/projects/Lab_oband_laser_and_amp_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m b/projects/Lab_oband_laser_and_amp_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m index c14d43d..d86c9c3 100644 --- a/projects/Lab_oband_laser_and_amp_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m +++ b/projects/Lab_oband_laser_and_amp_analysis/data_analysis_and_plot_scripts/lambda_vs_spectrum.m @@ -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 @@ -34,7 +34,8 @@ for p = 1:numel(plasers) spectrum_osa = wh_aeon.getStoValue('spectrum_osa',plasers(p),lambda,pumps(pmp)); wavelength_osa = wh_aeon.getStoValue('wavelength_osa',plasers(p),lambda,pumps(pmp)); 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); diff --git a/projects/MLSE_ML_based/minimal_example_huawei/minimal_example_ml_mlse.m b/projects/MLSE_ML_based/minimal_example_huawei/minimal_example_ml_mlse.m new file mode 100644 index 0000000..87b9ea9 --- /dev/null +++ b/projects/MLSE_ML_based/minimal_example_huawei/minimal_example_ml_mlse.m @@ -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 \ No newline at end of file diff --git a/workerError.mat b/workerError.mat index f03b730..73a5a0f 100644 Binary files a/workerError.mat and b/workerError.mat differ