diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 608510f..f0852e8 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,20); + p_lin = movmean(p_lin,5); if options.normalizeTo0dB p_lin = p_lin ./ max(p_lin); diff --git a/Classes/04_DSP/Equalizer/ML_MLSE.m b/Classes/04_DSP/Equalizer/ML_MLSE.m index a27712b..ee90526 100644 --- a/Classes/04_DSP/Equalizer/ML_MLSE.m +++ b/Classes/04_DSP/Equalizer/ML_MLSE.m @@ -280,7 +280,7 @@ classdef ML_MLSE < handle if sym_idx>=obj.L if obj.adaptive_mu mu_eff=CE_smooth(symbol); - mu_eff=max(min(mu_eff,0.2),1e-4); + mu_eff=max(min(mu_eff,0.2),1e-5); else mu_eff=mu; end diff --git a/Functions/EQ_blocks/ffe.m b/Functions/EQ_blocks/ffe.m index 07ea93b..5860b15 100644 --- a/Functions/EQ_blocks/ffe.m +++ b/Functions/EQ_blocks/ffe.m @@ -50,10 +50,10 @@ eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); %% Calculate performance metrics [snr, snr_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal); % [gmi] = calc_air(eq_signal_sd, tx_symbols, "skip_front", 10000, "skip_end", 10000); -[gmi] = calc_ngmi(eq_signal_sd,tx_symbols); +[gmi,ngmi] = calc_ngmi(eq_signal_sd,tx_symbols); gmi = max(gmi,0); -air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M)); +air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* ngmi; [evm_total, evm_lvl] = calc_evm(eq_signal_sd, tx_symbols); [std_total, std_lvl] = calc_std(eq_signal_sd, tx_symbols); [std_rxraw_total, std_rxraw_lvl] = calc_std(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols); diff --git a/Functions/EQ_recipes/dsp_400g_recipe.m b/Functions/EQ_recipes/dsp_400g_recipe.m index 5fe580e..b672867 100644 --- a/Functions/EQ_recipes/dsp_400g_recipe.m +++ b/Functions/EQ_recipes/dsp_400g_recipe.m @@ -86,7 +86,7 @@ if options.duob_mode ~= db_mode.db_encoded [vnle_results, equalized_signal] = runFfe(eq_vnle, "VNLE", ... Scpe_sig, Symbols, Tx_bits, options); vnle_results.config.equalizer_structure = equalizer_structure.ffe; - ffe_results.recipe_config = collectRecipeConfig("vnle", eq_ffe, p, options); + ffe_results.recipe_config = collectRecipeConfig("vnle", eq_vnle, p, options); output.vnle_package = vnle_results; if p.plot_output_signals @@ -273,14 +273,14 @@ end function p = defaultRecipeParameters() p = struct(); -p.run_ffe = false; -p.run_vnle = false; -p.run_dfe = false; -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 +p.run_ffe = 0; +p.run_vnle = 0; +p.run_dfe = 0; +p.run_vnle_mlse = 1; +p.run_dbtgt = 0; +p.run_ml_mlse = 0; % non-encoded and precoded branches +p.run_ml_mlse_db = 0; % db_encoded: ML-based MLSE +p.run_mlse_db = 1; % db_encoded: conventional MLSE p.preprocess_mode = "auto"; diff --git a/Theory/Mapping_Coding/partial_response/duobinary_back2back.m b/Theory/Mapping_Coding/partial_response/duobinary_back2back.m index 3fccde2..6e23652 100644 --- a/Theory/Mapping_Coding/partial_response/duobinary_back2back.m +++ b/Theory/Mapping_Coding/partial_response/duobinary_back2back.m @@ -1,4 +1,4 @@ -for M = [4] +for M = [6] bits = Signalgenerator("form", signalform.prms,"M", M,"order", 18).process(); mapper = PAMmapper(M,0); symbols = mapper.map(bits); @@ -13,7 +13,7 @@ for M = [4] symbols_db_ = awgn_channel(symbols_db,"snr_dB",15); figure(1003042); showLevelHistogram(symbols_db_,symbols_db,"ref_symbol_uncoded",symbols); - showLevelHistogram(symbols_db_,symbols_db); + % showLevelHistogram(symbols_db_,symbols_db); symbols_rx = db.decode(symbols_db); bits_rx = PAMmapper(M,0).demap(symbols_rx); diff --git a/Theory/Optical/calcFWM/analytical_calculation_paper.m b/Theory/Optical/calcFWM/analytical_calculation_paper.m index d088424..030b34a 100644 --- a/Theory/Optical/calcFWM/analytical_calculation_paper.m +++ b/Theory/Optical/calcFWM/analytical_calculation_paper.m @@ -16,13 +16,14 @@ for i = 1:length(N_) [Mndg,Mdg] = getProducts(N_(i)); total(i) = sum(Mndg) + sum(Mdg); subplot(1,length(N_),i) + fprintf('%d Channel: Non degenerate: %d; Degenerate: %d \n', N_(i),sum(Mndg),sum(Mdg)); xline(calcWavelengthPlan(N_(i), df_hz, center_nm)); hold on stem(channelplan_nm,(Mndg+Mdg),'filled','LineWidth',1,'Marker','o','MarkerSize',2) stem(channelplan_nm,(Mdg),'filled','LineWidth',1,'Marker','none'); ylim([0,100]) xlim([1260, 1365]); - grid off + % grid off xlabel('O-band wavelength region in nm'); ylabel('Number of FWM products'); title([num2str(N_(i)),' ch.']) diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_NGMI_THESIS_FINAL.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_NGMI_THESIS_FINAL.m new file mode 100644 index 0000000..96bfe9e --- /dev/null +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_NGMI_THESIS_FINAL.m @@ -0,0 +1,458 @@ +%% Final thesis figure: best EQ curves for NGMI, AIR and FEC rates +% Four-panel thesis view based on FIGURE_NGMI.m: +% a) NGMI +% b) AIR +% c) SD+HD FEC NDR +% d) O-FEC and KP4+Hamming NDR +% +% The normal EQ classes use the 2 km data set. DBS + VNLE + MLSE uses the +% post-2026 10 km data set because that is the available DBS measurement set. +% Within each selected EQ/PAM curve, the best valid result is retained for +% every symbol rate. The PAM mapping follows the reference thesis figure: +% VNLE -> PAM-8 +% VNLE + PF + MLSE -> PAM-6 +% VNLE DBt. + MLSE -> PAM-4 +% ML pre-EQ + Viterbi -> PAM-4 +% DBS + VNLE + MLSE -> PAM-4 and PAM-6 + +clear; + +%% 1) Configuration + +normalFiberLengthKm = 2; +duobinaryFiberLengthKm = 2; +selectedWavelengthNm = 1310; +selectedRopAttenuation = 0; +selectedIsMpi = 0; +duobinaryDateCutoff = datetime("2026-01-01 00:00:00"); + +% Thesis colors and curvet /PAM mapping. +curves = struct; +curves(1).name = "VNLE + PF + MLSE"; +curves(1).algorithm_key = "vnle_pf_mlse"; +curves(1).pam = 8; +curves(1).color = clr.Paired.red; +curves(1).marker = "o"; +curves(1).source = "normal"; +curves(1).usePrecodedBer = false; + +curves(2).name = "VNLE + PF + MLSE"; +curves(2).algorithm_key = "vnle_pf_mlse"; +curves(2).pam = 6; +curves(2).color = clr.Paired.green; +curves(2).marker = "square"; +curves(2).source = "normal"; +curves(2).usePrecodedBer = false; + +curves(3).name = "VNLE DBt. + MLSE"; +curves(3).algorithm_key = "vnle_db_mlse"; +curves(3).pam = 4; +curves(3).color = clr.Paired.blue; +curves(3).marker = "diamond"; +curves(3).source = "normal"; +curves(3).usePrecodedBer = true; + +curves(4).name = "ML pre-EQ + Viterbi"; +curves(4).algorithm_key = "ml_mlse"; +curves(4).pam = 4; +curves(4).color = clr.Paired.purple; +curves(4).marker = "^"; +curves(4).source = "normal"; +curves(4).usePrecodedBer = true; + +curves(5).name = "DBS + VNLE + MLSE (PAM-4)"; +curves(5).algorithm_key = "db_encoded"; +curves(5).pam = 4; +curves(5).color = clr.Paired.orange; +curves(5).marker = "v"; +curves(5).source = "duobinary"; +curves(5).usePrecodedBer = false; + +curves(6).name = "DBS + VNLE + MLSE (PAM-6)"; +curves(6).algorithm_key = "db_encoded"; +curves(6).pam = 6; +curves(6).color = clr.Paired.orange; +curves(6).marker = "v"; +curves(6).source = "duobinary"; +curves(6).usePrecodedBer = false; + +%% 2) Query normal 2 km and duobinary 10 km data + +db = DBHandler( ... + "dataBase", "labor_highspeed", ... + "type", "mysql", ... + "server", "192.168.178.192", ... + "user", "silas", ... + "password", "silas"); +db.refresh(); + +selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime'); + +normalRows = queryMeasurementRows(db, selectedFields, normalFiberLengthKm,selectedWavelengthNm); +normalRows = cleanMeasurementRows(normalRows); +normalRows = normalRows(ismember(normalRows.db_mode, ... + [double(db_mode.no_db), double(db_mode.db_precoded)]), :); +normalRows.algorithm_key = lower(string(normalRows.equalizer_structure)); + +duobinaryRows = queryMeasurementRows(db, selectedFields, duobinaryFiberLengthKm,selectedWavelengthNm); +duobinaryRows = cleanMeasurementRows(duobinaryRows); +duobinaryRows = duobinaryRows(duobinaryRows.db_mode == double(db_mode.db_encoded), :); +duobinaryRows = duobinaryRows( ... + equalizerMask(duobinaryRows.equalizer_structure, ... + equalizer_structure.db_encoded), :); + +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)>duobinaryDateCutoff,:); + duobinaryRows.date_of_processing = string(duobinaryRows.date_of_processing); +else + warning("figure_ngmi_thesis:NoProcessingDate", ... + "date_of_processing was not returned; no DBS date filter was applied."); +end +% DBS PAM-8 is intentionally excluded from the final thesis figure. +duobinaryRows = duobinaryRows(ismember(duobinaryRows.pam_level, [4 6]), :); +duobinaryRows.algorithm_key = repmat("db_encoded", height(duobinaryRows), 1); + +fprintf("Normal %g km rows: %d\n", normalFiberLengthKm, height(normalRows)); +fprintf("DBS %g km rows after date filter: %d\n", ... + duobinaryFiberLengthKm, height(duobinaryRows)); + +%% 3) Extract the best curves and calculate FEC rates + +tp = TransmissionPerformance; +results = repmat(emptyResult(), numel(curves), 1); + +for curveIdx = 1:numel(curves) + curve = curves(curveIdx); + + if curve.source == "normal" + curveRows = normalRows; + else + curveRows = duobinaryRows; + end + + curveRows = curveRows(curveRows.pam_level == curve.pam & ... + curveRows.algorithm_key == curve.algorithm_key, :); + + if isempty(curveRows) + warning("figure_ngmi_thesis:NoCurveRows", ... + "No rows found for %s, PAM-%d.", curve.name, curve.pam); + continue + end + + curveRows.BER_plot = curveRows.BER; + curveRows.precode = zeros(height(curveRows), 1); + if curve.usePrecodedBer && ismember("BER_precoded", ... + string(curveRows.Properties.VariableNames)) + usePrecoded = isfinite(curveRows.BER_precoded); + curveRows.BER_plot(usePrecoded) = curveRows.BER_precoded(usePrecoded); + curveRows.precode(usePrecoded) = 1; + end + + curveRows = addDerivedMetrics(curveRows); + + % Each metric is optimized independently. In particular, the BER + % winner does not have to be the NGMI or AIR winner for a given baud + % rate because these quantities are stored/calculated separately. + berSeries = bestSeries(curveRows, "BER_plot"); + ngmiSeries = bestSeries(curveRows, "NGMI"); + airSeries = bestSeries(curveRows, "AIR_Gbps"); + + grossRate = double(curveRows.grossrate); + measuredNgmi = double(curveRows.NGMI); + measuredBer = double(curveRows.BER_plot); + measuredNgmi(~isfinite(measuredNgmi) | measuredNgmi < 0 | ... + measuredNgmi > 1.05) = NaN; + measuredBer(~isfinite(measuredBer) | measuredBer <= 0 | ... + measuredBer > 0.5) = NaN; + + ndr = tp.calculateNetRate(grossRate, ... + "NGMI", measuredNgmi, ... + "BER", measuredBer); + + curveRows.NDR_SDHD = columnVector(ndr.SDHD.NetRate) .* 1e-9; + curveRows.NDR_O_FEC = columnVector(ndr.O_FEC.NetRate) .* 1e-9; + curveRows.NDR_KP4_HAMMING = columnVector(ndr.KP4_hamming.NetRate) .* 1e-9; + + results(curveIdx).name = curve.name; + results(curveIdx).color = curve.color; + results(curveIdx).marker = curve.marker; + results(curveIdx).pam = curve.pam; + results(curveIdx).ber = berSeries; + results(curveIdx).ngmi = ngmiSeries; + results(curveIdx).air = airSeries; + results(curveIdx).sdhd = bestSeries(curveRows, "NDR_SDHD"); + results(curveIdx).ofec = bestSeries(curveRows, "NDR_O_FEC"); + results(curveIdx).kp4Hamming = bestSeries(curveRows, "NDR_KP4_HAMMING"); + + fprintf("%s, PAM-%d: BER=%d, NGMI=%d, AIR=%d, SD+HD=%d, O-FEC=%d, KP4+Hamming=%d points\n", ... + curve.name, curve.pam, ... + numel(results(curveIdx).ber.x), ... + numel(results(curveIdx).ngmi.x), ... + numel(results(curveIdx).air.x), ... + numel(results(curveIdx).sdhd.x), ... + numel(results(curveIdx).ofec.x), ... + numel(results(curveIdx).kp4Hamming.x)); +end + +%% 4) Four-panel thesis figure + +fig = figure(72+normalFiberLengthKm); clf; +t = tiledlayout(fig, 1, 4, ... + "TileSpacing", "compact", ... + "Padding", "compact"); + +axNgmi = nexttile(t, 1); +hold(axNgmi, "on"); +for curveIdx = 1:numel(results) + if isempty(results(curveIdx).ngmi.x) + continue + end + plotSeries(axNgmi, results(curveIdx).ngmi, ... + results(curveIdx), "-", results(curveIdx).marker, results(curveIdx).name); +end +formatAxis(axNgmi, "NGMI", [0.90 1.00], [100 220]); +title(axNgmi, "a) NGMI"); + +axAir = nexttile(t, 2); +hold(axAir, "on"); +for curveIdx = 1:numel(results) + if isempty(results(curveIdx).air.x) + continue + end + plotSeries(axAir, results(curveIdx).air, results(curveIdx), ... + "-", results(curveIdx).marker, results(curveIdx).name); +end +formatAxis(axAir, "AIR [Gb/s]", [280 440], [100 220]); +yline(axAir, 400, "--", "Color", [0.25 0.25 0.25], ... + "HandleVisibility", "off"); +title(axAir, "b) AIR"); + +axSdhd = nexttile(t, 3); +hold(axSdhd, "on"); +for curveIdx = 1:numel(results) + if isempty(results(curveIdx).sdhd.x) + continue + end + plotSeries(axSdhd, results(curveIdx).sdhd, results(curveIdx), ... + "-", results(curveIdx).marker, results(curveIdx).name); +end +formatAxis(axSdhd, "NDR [Gb/s]", [280 430], [100 220]); +yline(axSdhd, 400, "--", "Color", [0.25 0.25 0.25], ... + "HandleVisibility", "off"); +title(axSdhd, "c) SD+HD FEC"); + +axHd = nexttile(t, 4); +hold(axHd, "on"); +hdHandles = gobjects(numel(curves), 1); +for curveIdx = 1:numel(results) + curve = results(curveIdx); + if ~isempty(curve.ofec.x) + hdHandles(curveIdx) = plotSeries(axHd, curve.ofec, curve, ":", curve.marker, ... + curve.name + " — O-FEC"); + end + if ~isempty(curve.kp4Hamming.x) + plotSeries(axHd, curve.kp4Hamming, curve, "--", "diamond", ... + curve.name + " — KP4+Hamming"); + end +end +formatAxis(axHd, "NDR [Gb/s]", [280 430], [100 220]); +yline(axHd, 400, "--", "Color", [0.25 0.25 0.25], ... + "HandleVisibility", "off"); +title(axHd, "d) O-FEC; KP4+Hamming"); + +validLegendHandles = hdHandles(isgraphics(hdHandles)); +if ~isempty(validLegendHandles) + legend(axHd, validLegendHandles, ... + {results(isgraphics(hdHandles)).name}, ... + "Location", "southoutside", ... + "NumColumns", min(3, numel(validLegendHandles)), ... + "Interpreter", "none"); +end + +sgtitle(t, sprintf("Best EQ-class results, lambda=%g nm", selectedWavelengthNm)); +set(fig, "Position", 1e3 .* [0.18 0.55 1.52 0.31]); + +%% Local helpers + +function T = queryMeasurementRows(db, fields, fiberLengthKm,selectedWavelengthNm) +fp = QueryFilter(); +fp.where("Runs", "fiber_length", "EQUALS", fiberLengthKm); +% fp.where("Runs", "wavelength", "EQUALS", selectedWavelengthNm); +fp.where("Runs", "rop_attenuation", "EQUALS", 0); +fp.where("Runs", "is_mpi", "EQUALS", 0); +[T, query] = db.queryDB(fp, fields); +disp(query); +end + +function T = cleanMeasurementRows(T) +numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ... + "symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ... + "rop_attenuation", "numBits", "numBitErr", "BER", ... + "numBitErr_precoded", "BER_precoded", "GMI", "AIR", "NGMI"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(T.Properties.VariableNames)) + T.(char(fieldName)) = numericColumn(T.(char(fieldName))); + end +end +T.equalizer_structure = lower(string(T.equalizer_structure)); +end + +function result = emptyResult() +result = struct( ... + "name", "", ... + "color", [0 0 0], ... + "marker", "o", ... + "pam", NaN, ... + "ber", emptySeries(), ... + "ngmi", emptySeries(), ... + "air", emptySeries(), ... + "sdhd", emptySeries(), ... + "ofec", emptySeries(), ... + "kp4Hamming", emptySeries()); +end + +function series = emptySeries() +series = struct( ... + "x", [], ... + "y", [], ... + "ber", [], ... + "ngmi", [], ... + "air", [], ... + "wavelength", [], ... + "precode", [], ... + "dbMode", []); +end + +function series = bestSeries(T, fieldName) +series = emptySeries(); +if ~ismember(fieldName, string(T.Properties.VariableNames)) + return +end + +values = numericColumn(T.(char(fieldName))); +valid = isfinite(values); +if strcmp(fieldName, "BER_plot") + % BER is minimized; zero is omitted because the axis is logarithmic. + valid = valid & values > 0 & values <= 0.5; +elseif strcmp(fieldName, "NGMI") + valid = valid & values >= 0 & values <= 1.05; +else + valid = valid & values >= 0; +end + +candidate = T(valid, :); +values = values(valid); +if isempty(candidate) + return +end + +[groupId, ~] = findgroups(candidate.symbolrate_GBd); +keepIndex = zeros(max(groupId), 1); +for groupIdx = 1:max(groupId) + rowIndex = find(groupId == groupIdx); + if strcmp(fieldName, "BER_plot") + [~, localIndex] = min(values(rowIndex)); + else + [~, localIndex] = max(values(rowIndex)); + end + keepIndex(groupIdx) = rowIndex(localIndex(1)); +end + +candidate = candidate(keepIndex, :); +series.x = double(candidate.symbolrate_GBd); +series.y = values(keepIndex); +[series.x, order] = sort(series.x); +series.y = series.y(order); + +% Preserve the complete selected database row as data-tip metadata. The +% metadata therefore belongs to the plotted point, even when BER, NGMI, +% AIR, and NDR select different rows at the same baud rate. +series.ber = numericOrNaN(candidate, "BER_plot"); +series.ngmi = numericOrNaN(candidate, "NGMI"); +series.air = numericOrNaN(candidate, "AIR_Gbps"); +series.wavelength = numericOrNaN(candidate, "wavelength"); +series.precode = numericOrNaN(candidate, "precode"); +series.dbMode = numericOrNaN(candidate, "db_mode"); +series.ber = series.ber(order); +series.ngmi = series.ngmi(order); +series.air = series.air(order); +series.wavelength = series.wavelength(order); +series.precode = series.precode(order); +series.dbMode = series.dbMode(order); +end + +function T = addDerivedMetrics(T) +T.symbolrate_GBd = double(T.symbolrate) .* 1e-9; + +T.AIR_Gbps = numericOrNaN(T, "AIR") .* 1e-9; +gmi = numericOrNaN(T, "GMI"); +fallbackAir = gmi .* double(T.symbolrate) .* 1e-9; +grossRate = double(T.grossrate); +useFallback = ~isfinite(T.AIR_Gbps) | T.AIR_Gbps < 0 | ... + (isfinite(grossRate) & T.AIR_Gbps > grossRate .* 1.05e-9); +T.AIR_Gbps(useFallback) = fallbackAir(useFallback); +T.NGMI = numericOrNaN(T, "NGMI"); +end + +function h = plotSeries(ax, series, curve, lineStyle, marker, displayName) +h = plot(ax, series.x, series.y, ... + "LineStyle", lineStyle, ... + "Marker", marker, ... + "MarkerSize", 4, ... + "LineWidth", 1.35, ... + "Color", curve.color, ... + "MarkerFaceColor", curve.color, ... + "MarkerEdgeColor", curve.color, ... + "DisplayName", displayName); +h.DataTipTemplate.DataTipRows = [ ... + dataTipTextRow("X", series.x); ... + dataTipTextRow("Y", series.y); ... + dataTipTextRow("symbol rate [GBd]", series.x); ... + dataTipTextRow("BER", series.ber); ... + dataTipTextRow("NGMI", series.ngmi); ... + dataTipTextRow("AIR [Gb/s]", series.air); ... + dataTipTextRow("wavelength [nm]", series.wavelength); ... + dataTipTextRow("precoding", series.precode); ... + dataTipTextRow("db_mode", series.dbMode)]; +end + +function formatAxis(ax, yLabel, yLimits, xLimits) +set(ax, "FontSize", 8, "TickLabelInterpreter", "none"); +xlabel(ax, "Baud rate [GBd]"); +ylabel(ax, yLabel); +xlim(ax, xLimits); +xticks(ax, 100:15:220); +ylim(ax, yLimits); +grid(ax, "on"); +grid(ax, "minor"); +box(ax, "on"); +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 = numericOrNaN(T, fieldName) +if ismember(fieldName, string(T.Properties.VariableNames)) + values = numericColumn(T.(char(fieldName))); +else + values = NaN(height(T), 1); +end +values = values(:); +end + +function values = columnVector(values) +values = double(values(:)); +end + +function mask = equalizerMask(equalizerColumn, eqValue) +mask = lower(string(equalizerColumn)) == lower(string(eqValue)); +end diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_NGMI_v2.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_NGMI_v2.m index a690207..314f00e 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_NGMI_v2.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_NGMI_v2.m @@ -161,6 +161,7 @@ for ti = 1:3 grid minor; box on; beautifyBERplot; yline(400,'HandleVisibility','off'); + legend end % === FIX FIGURE SIZE FOR TIKZ ========================================== diff --git a/projects/Diss/400G_revisit/CHECK_BEST_PREEMPHASIS_PRECODE_2KM_1310NM.m b/projects/Diss/400G_revisit/CHECK_BEST_PREEMPHASIS_PRECODE_2KM_1310NM.m new file mode 100644 index 0000000..97bedc7 --- /dev/null +++ b/projects/Diss/400G_revisit/CHECK_BEST_PREEMPHASIS_PRECODE_2KM_1310NM.m @@ -0,0 +1,241 @@ +%% Best pre-emphasis and precoding settings at 2 km / 1310 nm +% The four requested technique labels map to the stored equalizer classes: +% VNLE -> vnle +% MLSE -> VNLE + PF + MLSE (vnle_pf_mlse) +% DB tgt. -> VNLE DBt. + MLSE (vnle_db_mlse) +% ML based -> ML pre-EQ + Viterbi (ml_mlse) +% +% For every technique and PAM format, the globally lowest valid BER is +% selected across the available gross rates, db_mode 0/1 variants, and +% BER/BER_precoded variants. The resulting table reports binary settings: +% preemph = 0/1 and precode = 0/1. + +clear; clc; + +%% 1) Query the 2 km / 1310 nm normal-link data + +selectedFiberLengthKm = 2; +selectedWavelengthNm = 1310; +selectedRopAttenuation = 0; +selectedIsMpi = 0; +selectedPamLevels = [4 6 8]; +normalDbModes = [double(db_mode.no_db), double(db_mode.db_precoded)]; + +techniques = table( ... + ["VNLE"; "MLSE"; "DB tgt."; "ML based"], ... + ["vnle"; "vnle_pf_mlse"; "vnle_db_mlse"; "ml_mlse"], ... + 'VariableNames', ["technique", "algorithm_key"]); + +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); + +fields = db.getTableFieldNames('dashboard_ungrouped_alltime'); +fields = appendMissingFields(fields, {'Runs.precomp_amp'; 'Runs.is_mpi'}); +[rawData, query] = db.queryDB(fp, fields); +disp(query); + +data = cleanRows(rawData); +data = data(ismember(data.db_mode, normalDbModes) & ... + ismember(data.pam_level, selectedPamLevels), :); + +if ~ismember("precomp_amp", string(data.Properties.VariableNames)) + warning("check_best_settings:NoPrecompAmp", ... + "Runs.precomp_amp was not returned; using db_mode == 0 as preemph = 1."); + data.pre_emphasis = data.db_mode == double(db_mode.no_db); +else + data.pre_emphasis = derivePreEmphasis(data.precomp_amp, data.db_mode); +end + +%% 2) Build BER and BER_precoded candidates + +baseRows = data(isfinite(data.BER) & data.BER > 0, :); +baseRows.precode = zeros(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) & data.BER_precoded > 0, :); + precodedRows.precode = ones(height(precodedRows), 1); + precodedRows.BER_plot = precodedRows.BER_precoded; + precodedRows.algorithm_key = algorithmKeyFromEqualizer( ... + precodedRows.equalizer_structure); + precodedRows = precodedRows(precodedRows.algorithm_key ~= "", :); + metricRows = [baseRows; precodedRows]; +else + warning("check_best_settings:NoPrecodedBer", ... + "BER_precoded was not returned; only precode = 0 is available."); + metricRows = baseRows; +end + +fprintf("Candidate rows: %d\n", height(metricRows)); + +%% 3) Find the lowest BER and its binary settings + +numResults = height(techniques) * numel(selectedPamLevels); +best = repmat(struct( ... + "technique", "", ... + "algorithm_key", "", ... + "pam", NaN, ... + "BER", NaN, ... + "preemph", NaN, ... + "precode", NaN, ... + "db_mode", NaN, ... + "grossrate_Gbps", NaN, ... + "symbolrate_GBd", NaN), numResults, 1); + +resultIdx = 0; +for techniqueIdx = 1:height(techniques) + for pamIdx = 1:numel(selectedPamLevels) + resultIdx = resultIdx + 1; + algorithmKey = techniques.algorithm_key(techniqueIdx); + pamLevel = selectedPamLevels(pamIdx); + candidates = metricRows( ... + metricRows.algorithm_key == algorithmKey & ... + metricRows.pam_level == pamLevel, :); + + best(resultIdx).technique = techniques.technique(techniqueIdx); + best(resultIdx).algorithm_key = algorithmKey; + best(resultIdx).pam = pamLevel; + if isempty(candidates) + continue + end + + [~, bestRowIdx] = min(candidates.BER_plot); + row = candidates(bestRowIdx, :); + best(resultIdx).BER = row.BER_plot; + best(resultIdx).preemph = double(row.pre_emphasis); + best(resultIdx).precode = double(row.precode); + best(resultIdx).db_mode = double(row.db_mode); + best(resultIdx).grossrate_Gbps = double(row.grossrate) .* 1e-9; + best(resultIdx).symbolrate_GBd = double(row.symbolrate) .* 1e-9; + end +end + +bestRows = struct2table(best); +disp(bestRows); + +%% 4) Compact tables for direct use in the thesis configuration + +binaryTable = table(techniques.technique, ... + NaN(height(techniques), 1), NaN(height(techniques), 1), ... + NaN(height(techniques), 1), NaN(height(techniques), 1), ... + NaN(height(techniques), 1), NaN(height(techniques), 1), ... + 'VariableNames', ["technique", ... + "PAM4_preemph", "PAM4_precode", ... + "PAM6_preemph", "PAM6_precode", ... + "PAM8_preemph", "PAM8_precode"]); + +labelTable = table(techniques.technique, strings(height(techniques), 1), ... + strings(height(techniques), 1), strings(height(techniques), 1), ... + 'VariableNames', ["technique", "PAM4", "PAM6", "PAM8"]); + +for rowIdx = 1:height(bestRows) + techniqueIdx = find(techniques.algorithm_key == bestRows.algorithm_key(rowIdx), 1); + pamField = sprintf("PAM%d", bestRows.pam(rowIdx)); + binaryTable.(sprintf("%s_preemph", pamField))(techniqueIdx) = ... + bestRows.preemph(rowIdx); + binaryTable.(sprintf("%s_precode", pamField))(techniqueIdx) = ... + bestRows.precode(rowIdx); + + if isfinite(bestRows.BER(rowIdx)) + labelTable.(pamField)(techniqueIdx) = sprintf( ... + "preemph=%d, precode=%d (BER=%.3g)", ... + bestRows.preemph(rowIdx), bestRows.precode(rowIdx), ... + bestRows.BER(rowIdx)); + else + labelTable.(pamField)(techniqueIdx) = "no valid BER"; + end +end + +disp("Binary setting table:"); +disp(binaryTable); +disp("Settings with winning BER:"); +disp(labelTable); + +%% 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 T = cleanRows(T) +numericFields = ["result_id", "run_id", "eq_id", "bitrate", ... + "grossrate", "symbolrate", "pam_level", "wavelength", ... + "fiber_length", "db_mode", "rop_attenuation", "precomp_amp", ... + "numBits", "numBitErr", "BER", "numBitErr_precoded", ... + "BER_precoded", "GMI", "AIR", "NGMI"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(T.Properties.VariableNames)) + T.(char(fieldName)) = numericColumn(T.(char(fieldName))); + end +end +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 algorithmKey = algorithmKeyFromEqualizer(equalizerColumn) +eqNumeric = equalizerNumeric(equalizerColumn); +algorithmKey = strings(size(eqNumeric)); +algorithmKey(eqNumeric == double(equalizer_structure.vnle)) = "vnle"; +algorithmKey(eqNumeric == double(equalizer_structure.vnle_pf_mlse)) = ... + "vnle_pf_mlse"; +algorithmKey(eqNumeric == double(equalizer_structure.vnle_db_mlse)) = ... + "vnle_db_mlse"; +algorithmKey(eqNumeric == double(equalizer_structure.ml_mlse)) = "ml_mlse"; +end + +function eqNumeric = equalizerNumeric(equalizerColumn) +if isa(equalizerColumn, "equalizer_structure") || isnumeric(equalizerColumn) + eqNumeric = double(equalizerColumn); + return +end + +equalizerString = string(equalizerColumn); +eqNumeric = str2double(equalizerString); +enumNames = ["vnle", "vnle_pf_mlse", "vnle_db_mlse", "ml_mlse"]; +enumValues = [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(enumNames) + missingNumeric = isnan(eqNumeric); + eqNumeric(missingNumeric & equalizerString == enumNames(idx)) = ... + enumValues(idx); +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 diff --git a/projects/Diss/400G_revisit/FIGURE_EQ_NOISE_VS_BAUDRATE.m b/projects/Diss/400G_revisit/FIGURE_EQ_NOISE_VS_BAUDRATE.m index 54d0d88..9bc827b 100644 --- a/projects/Diss/400G_revisit/FIGURE_EQ_NOISE_VS_BAUDRATE.m +++ b/projects/Diss/400G_revisit/FIGURE_EQ_NOISE_VS_BAUDRATE.m @@ -35,10 +35,6 @@ 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'); @@ -65,6 +61,12 @@ for runIndex = 1:numel(run_ids) runData = queryRunid(run_ids(runIndex), db); fsym = double(runData.symbolrate(1)); bitrate = double(runData.bitrate(1)); + symbolrate = double(runData.symbolrate(1)); + if double(runData.db_mode) == 0 || double(runData.db_mode) == 2 + ispreemph = 1; + elseif double(runData.db_mode) == 1 + ispreemph = 0; + end dsp_options.max_occurences = 20; dsp_options.start_occurence = 2; @@ -94,53 +96,54 @@ for runIndex = 1:numel(run_ids) "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); + % + dbt = 1; + if dbt + Symbols_ = Duobinary().encode(Symbols); + else + Symbols_ = Symbols; + end + [equalized_signal, eq_noise] = eq_vnle.process(Scpe_sig, Symbols_); 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); + if ~dbt + postfilter_order = 3; + pf = Postfilter("ncoeff", postfilter_order, "useBurg", 1); + [mlse_sig_sd,whitened_noise] = pf.process(equalized_signal, eq_noise); + end + + fig = figure(220); + if ~dbt + showEQNoisePSD(eq_noise, ... + "fignum", fig.Number, ... + "displayname", sprintf('EEN: %.0f GBd; preemph: %d', symbolrate*1e-9,ispreemph), ... + "postfilter_taps", pf.coefficients, ... + "colormode", "qualitative","color",clr.Paired.lblue); + whitened_noise = whitened_noise - mean(whitened_noise.signal); + whitened_noise.spectrum("displayname", 'Whitened Noise', "fignum", fig.Number, "normalizeTo0dB", 0,"fft_length",4096,"normalizeToDC",0,"color",clr.Paired.dblue); + else + showEQNoisePSD(eq_noise, ... + "fignum", fig.Number, ... + "displayname", sprintf('EEN: %.0f GBd; preemph: %d', symbolrate*1e-9,ispreemph), ... + "colormode", "qualitative","color",clr.Paired.lblue); + 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'); + + % ylim([-70 -50]); max_fs = max(max_fs, eq_noise.fs); + xlim([-max_fs/2 max_fs/2]*1e-9); 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}; diff --git a/projects/Diss/400G_revisit/PLOT_AIR_BEST_CURVES_OVER_WAVELENGTH_10KM.m b/projects/Diss/400G_revisit/PLOT_AIR_BEST_CURVES_OVER_WAVELENGTH_10KM.m new file mode 100644 index 0000000..984131b --- /dev/null +++ b/projects/Diss/400G_revisit/PLOT_AIR_BEST_CURVES_OVER_WAVELENGTH_10KM.m @@ -0,0 +1,407 @@ +%% Maximum AIR over wavelength at 10 km +% One curve is shown for each selected EQ/PAM contender. For every +% wavelength, the row with the maximum AIR is retained across the available +% symbol rates and normal db_mode 0/1 variants. +% DBS uses db_mode = 2, equalizer_structure = db_encoded, and only results +% processed after 2026-01-01. + +clear; clc; + +%% 1) Configuration + +normalFiberLengthKm = 2; +duobinaryFiberLengthKm = 10; +selectedRopAttenuation = 0; +selectedIsMpi = 0; +duobinaryDateCutoff = datetime("2026-01-01 00:00:00"); + +% Thesis colors and curve/PAM mapping. +curves = struct; + +curves(1).name = "DBS + VNLE + MLSE (PAM-4)"; +curves(1).algorithm_key = "db_encoded"; +curves(1).pam = 4; +curves(1).color = clr.Paired.orange; +curves(1).marker = "square"; +curves(1).source = "duobinary"; +curves(1).usePrecodedBer = false; + +curves(2).name = "VNLE DBt. + MLSE"; +curves(2).algorithm_key = "vnle_db_mlse"; +curves(2).pam = 4; +curves(2).color = clr.Paired.blue; +curves(2).marker = "square"; +curves(2).source = "normal"; +curves(2).usePrecodedBer = true; + +% curves(3).name = "ML pre-EQ + Viterbi"; +% curves(3).algorithm_key = "ml_mlse"; +% curves(3).pam = 4; +% curves(3).color = clr.Paired.purple; +% curves(3).marker = "square"; +% curves(3).source = "normal"; +% curves(3).usePrecodedBer = true; + +curves(3).name = "VNLE + PF + MLSE"; +curves(3).algorithm_key = "vnle_pf_mlse"; +curves(3).pam = 6; +curves(3).color = clr.Paired.green; +curves(3).marker = "v"; +curves(3).source = "normal"; +curves(3).usePrecodedBer = false; + +% curves(5).name = "ML pre-EQ + Viterbi"; +% curves(5).algorithm_key = "ml_mlse"; +% curves(5).pam = 6; +% curves(5).color = clr.Paired.purple; +% curves(5).marker = "v"; +% curves(5).source = "normal"; +% curves(5).usePrecodedBer = true; + +curves(4).name = "DBS + VNLE + MLSE (PAM-6)"; +curves(4).algorithm_key = "db_encoded"; +curves(4).pam = 6; +curves(4).color = clr.Paired.orange; +curves(4).marker = "v"; +curves(4).source = "duobinary"; +curves(4).usePrecodedBer = false; + +curves(5).name = "VNLE + PF + MLSE"; +curves(5).algorithm_key = "vnle_pf_mlse"; +curves(5).pam = 8; +curves(5).color = clr.Paired.red; +curves(5).marker = "o"; +curves(5).source = "normal"; +curves(5).usePrecodedBer = false; + + + +%% 2) Query normal and DBS rows over all wavelengths + +db = DBHandler( ... + "dataBase", "labor_highspeed", ... + "type", "mysql", ... + "server", "192.168.178.192", ... + "user", "silas", ... + "password", "silas"); +db.refresh(); + +selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime'); + +normalRows = queryRows(db, selectedFields, normalFiberLengthKm, ... + selectedRopAttenuation, selectedIsMpi); +normalRows = cleanRows(normalRows); +normalRows = normalRows(ismember(normalRows.db_mode, ... + [double(db_mode.no_db), double(db_mode.db_precoded)]), :); +normalRows.algorithm_key = lower(string(normalRows.equalizer_structure)); + +duobinaryRows = queryRows(db, selectedFields, duobinaryFiberLengthKm, ... + selectedRopAttenuation, selectedIsMpi); +duobinaryRows = cleanRows(duobinaryRows); +duobinaryRows = duobinaryRows( ... + duobinaryRows.db_mode == double(db_mode.db_encoded), :); +duobinaryRows = duobinaryRows( ... + equalizerMask(duobinaryRows.equalizer_structure, ... + equalizer_structure.db_encoded), :); + +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) > duobinaryDateCutoff, :); +else + warning("plot_air_wavelength:NoProcessingDate", ... + "date_of_processing was not returned; no DBS date filter was applied."); +end +duobinaryRows.algorithm_key = repmat("db_encoded", height(duobinaryRows), 1); + +normalRows = addAirMetric(normalRows); +duobinaryRows = addAirMetric(duobinaryRows); + +fprintf("Normal %g km rows: %d\n", normalFiberLengthKm, height(normalRows)); +fprintf("DBS %g km rows after date filter: %d\n", ... + duobinaryFiberLengthKm, height(duobinaryRows)); + +%% 3) Select maximum AIR per wavelength for every configured curve + +tp = TransmissionPerformance; +results = repmat(emptyResult(), numel(curves), 1); +for curveIdx = 1:numel(curves) + curve = curves(curveIdx); + if curve.source == "normal" + sourceRows = normalRows; + else + sourceRows = duobinaryRows; + end + + curveRows = sourceRows(sourceRows.pam_level == curve.pam & ... + sourceRows.algorithm_key == curve.algorithm_key, :); + curveRows.BER_plot = curveRows.BER; + if curve.usePrecodedBer && ismember("BER_precoded", ... + string(curveRows.Properties.VariableNames)) + usePrecoded = isfinite(curveRows.BER_precoded); + curveRows.BER_plot(usePrecoded) = curveRows.BER_precoded(usePrecoded); + end + + measuredNgmi = numericOrNaN(curveRows, "NGMI"); + measuredBer = numericOrNaN(curveRows, "BER_plot"); + grossRate = numericOrNaN(curveRows, "grossrate"); + measuredNgmi(~isfinite(measuredNgmi) | measuredNgmi < 0 | ... + measuredNgmi > 1.05) = NaN; + measuredBer(~isfinite(measuredBer) | measuredBer <= 0 | ... + measuredBer > 0.5) = NaN; + + ndr = tp.calculateNetRate(grossRate, ... + "NGMI", measuredNgmi, ... + "BER", measuredBer); + curveRows.NDR_O_FEC = columnVector(ndr.O_FEC.NetRate) .* 1e-9; + + bestRows = bestMetricByWavelength(curveRows, "AIR_Gbps"); + bestOfecRows = bestMetricByWavelength(curveRows, "NDR_O_FEC"); + + results(curveIdx).name = curve.name; + results(curveIdx).color = curve.color; + results(curveIdx).marker = curve.marker; + results(curveIdx).pam = curve.pam; + results(curveIdx).air = makeMetricSeries( ... + bestRows, "AIR_Gbps", curve.usePrecodedBer); + results(curveIdx).ofec = makeMetricSeries( ... + bestOfecRows, "NDR_O_FEC", curve.usePrecodedBer); + + fprintf("%s, PAM-%d: AIR=%d, O-FEC=%d wavelength points\n", ... + curve.name, curve.pam, ... + numel(results(curveIdx).air.x), ... + numel(results(curveIdx).ofec.x)); +end + +%% 4) Plot maximum AIR and O-FEC NDR over wavelength + +fig = figure(435); clf; +t = tiledlayout(fig, 1, 2, ... + "TileSpacing", "compact", ... + "Padding", "compact"); + +axAir = nexttile(t, 1); +airHandles = plotMetricPanel(axAir, results, "air", ... + "AIR [Gb/s]", "a) Maximum AIR"); + +axOfec = nexttile(t, 2); +ofecHandles = plotMetricPanel(axOfec, results, "ofec", ... + "O-FEC NDR [Gb/s]", "b) O-FEC"); + +allWavelengthCells = cell(numel(results), 1); +for curveIdx = 1:numel(results) + allWavelengthCells{curveIdx} = [ ... + results(curveIdx).air.x(:); ... + results(curveIdx).ofec.x(:)]; +end +allWavelengths = vertcat(allWavelengthCells{:}); +allWavelengths = allWavelengths(isfinite(allWavelengths)); +if ~isempty(allWavelengths) + xLimits = [min(allWavelengths) - 1, max(allWavelengths) + 1]; + for ax = [axAir, axOfec] + xlim(ax, xLimits); + xticks(ax, unique(allWavelengths)); + end +end + +validOfecMask = isgraphics(ofecHandles); +if any(validOfecMask) + legend(axOfec, ofecHandles(validOfecMask), ... + makeLegendNames(results(validOfecMask)), ... + "Location", "southoutside", ... + "NumColumns", min(3, nnz(validOfecMask)), ... + "Interpreter", "none"); +elseif any(isgraphics(airHandles)) + validAirMask = isgraphics(airHandles); + legend(axAir, airHandles(validAirMask), ... + makeLegendNames(results(validAirMask)), ... + "Location", "southoutside", ... + "NumColumns", min(3, nnz(validAirMask)), ... + "Interpreter", "none"); +end + +sgtitle(t, sprintf("Best AIR and O-FEC results over wavelength, %.0f km", ... + normalFiberLengthKm)); +set(fig, "Position", 1e3 .* [0.08 0.42 1.75 0.48]); + +%% Local helpers + +function T = queryRows(db, fields, fiberLengthKm, ropAttenuation, isMpi) +fp = QueryFilter(); +fp.where("Runs", "fiber_length", "EQUALS", fiberLengthKm); +fp.where("Runs", "rop_attenuation", "EQUALS", ropAttenuation); +fp.where("Runs", "is_mpi", "EQUALS", isMpi); +[T, query] = db.queryDB(fp, fields); +disp(query); +end + +function T = cleanRows(T) +numericFields = ["result_id", "run_id", "eq_id", "bitrate", ... + "grossrate", "symbolrate", "pam_level", "wavelength", ... + "fiber_length", "db_mode", "rop_attenuation", "numBits", ... + "numBitErr", "BER", "numBitErr_precoded", "BER_precoded", ... + "GMI", "AIR", "NGMI"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(T.Properties.VariableNames)) + T.(char(fieldName)) = numericColumn(T.(char(fieldName))); + end +end +T.equalizer_structure = lower(string(T.equalizer_structure)); +end + +function T = addAirMetric(T) +T.AIR_Gbps = numericOrNaN(T, "AIR") .* 1e-9; +gmi = numericOrNaN(T, "GMI"); +symbolrate = numericOrNaN(T, "symbolrate"); +grossrate = numericOrNaN(T, "grossrate"); +fallbackAir = gmi .* symbolrate .* 1e-9; +useFallback = ~isfinite(T.AIR_Gbps) | T.AIR_Gbps < 0 | ... + (isfinite(grossrate) & T.AIR_Gbps > grossrate .* 1.05e-9); +T.AIR_Gbps(useFallback) = fallbackAir(useFallback); +end + +function bestRows = bestMetricByWavelength(T, metricField) +if isempty(T) + bestRows = T; + return +end + +values = numericOrNaN(T, metricField); +valid = isfinite(T.wavelength) & isfinite(values) & values >= 0; +candidate = T(valid, :); +values = values(valid); +if isempty(candidate) + bestRows = candidate; + return +end + +[groupId, ~] = findgroups(candidate.wavelength); +keepIdx = zeros(max(groupId), 1); +for groupIdx = 1:max(groupId) + rowIdx = find(groupId == groupIdx); + [~, localIdx] = max(values(rowIdx)); + keepIdx(groupIdx) = rowIdx(localIdx(1)); +end +bestRows = sortrows(candidate(keepIdx, :), "wavelength"); +end + +function series = makeMetricSeries(T, metricField, usePrecodedBer) +series = emptyMetricSeries(); +if isempty(T) + return +end + +series.x = numericOrNaN(T, "wavelength"); +series.y = numericOrNaN(T, metricField); +series.air = numericOrNaN(T, "AIR_Gbps"); +series.ngmi = numericOrNaN(T, "NGMI"); +series.ber = numericOrNaN(T, "BER_plot"); +series.symbolrate = numericOrNaN(T, "symbolrate") .* 1e-9; +series.grossrate = numericOrNaN(T, "grossrate") .* 1e-9; +series.precoding = repmat(double(usePrecodedBer), height(T), 1); +series.dbMode = numericOrNaN(T, "db_mode"); +end + +function series = emptyMetricSeries() +series = struct( ... + "x", [], ... + "y", [], ... + "air", [], ... + "ngmi", [], ... + "ber", [], ... + "symbolrate", [], ... + "grossrate", [], ... + "precoding", [], ... + "dbMode", []); +end + +function handles = plotMetricPanel(ax, results, seriesField, yLabel, panelTitle) +hold(ax, "on"); +handles = gobjects(numel(results), 1); +for curveIdx = 1:numel(results) + result = results(curveIdx); + series = result.(seriesField); + if isempty(series.x) + continue + end + + handles(curveIdx) = plot(ax, series.x, series.y, ... + "LineStyle", "-", ... + "Marker", result.marker, ... + "MarkerSize", 5, ... + "LineWidth", 1.35, ... + "Color", result.color, ... + "MarkerFaceColor", result.color, ... + "MarkerEdgeColor", result.color, ... + "DisplayName", result.name); + + handles(curveIdx).DataTipTemplate.DataTipRows = [ ... + dataTipTextRow(yLabel, series.y); ... + dataTipTextRow("AIR [Gb/s]", series.air); ... + dataTipTextRow("NGMI", series.ngmi); ... + dataTipTextRow("BER", series.ber); ... + dataTipTextRow("Symbol rate [GBd]", series.symbolrate); ... + dataTipTextRow("Gross rate [Gb/s]", series.grossrate); ... + dataTipTextRow("precoding", series.precoding); ... + dataTipTextRow("db_mode", series.dbMode)]; +end + +set(ax, "FontSize", 9, "TickLabelInterpreter", "none"); +xlabel(ax, "Wavelength [nm]"); +ylabel(ax, yLabel); +title(ax, panelTitle); +grid(ax, "on"); +grid(ax, "minor"); +box(ax, "on"); +end + +function names = makeLegendNames(results) +names = strings(numel(results), 1); +for idx = 1:numel(results) + if contains(results(idx).name, "PAM-") + names(idx) = results(idx).name; + else + names(idx) = results(idx).name + ... + " (PAM-" + string(results(idx).pam) + ")"; + end +end +end + +function result = emptyResult() +result = struct( ... + "name", "", ... + "color", [0 0 0], ... + "marker", "o", ... + "pam", NaN, ... + "air", emptyMetricSeries(), ... + "ofec", emptyMetricSeries()); +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 = numericOrNaN(T, fieldName) +if ismember(fieldName, string(T.Properties.VariableNames)) + values = numericColumn(T.(char(fieldName))); +else + values = NaN(height(T), 1); +end +values = values(:); +end + +function values = columnVector(values) +values = double(values(:)); +end + +function mask = equalizerMask(equalizerColumn, eqValue) +mask = lower(string(equalizerColumn)) == lower(string(eqValue)); +end 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 dc57fc4..02e0b13 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 @@ -33,7 +33,7 @@ db.refresh(); fp = QueryFilter(); fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm); -fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm); +% fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm); if ~isempty(selectedRopAttenuation) fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation); end @@ -141,7 +141,7 @@ for pamIdx = 1:numel(selectedPamLevels) algoData = sortrows(bestPlotData(rowMask, :), "grossrate_Gbps"); if showRawEntries - scatter(ax, algoData.grossrate_Gbps, algoData.BER_plot, ... + rawHandle = scatter(ax, algoData.grossrate_Gbps, algoData.BER_plot, ... 9, ... "Marker", ".", ... "MarkerEdgeColor", style.color, ... @@ -149,10 +149,14 @@ for pamIdx = 1:numel(selectedPamLevels) "MarkerEdgeAlpha", 0.25, ... "MarkerFaceAlpha", 0.25, ... "HandleVisibility", "off"); + rawHandle.DataTipTemplate.DataTipRows = [ ... + dataTipTextRow("Gross rate [Gb/s]", algoData.grossrate_Gbps); ... + dataTipTextRow("BER", algoData.BER_plot); ... + dataTipTextRow("db_mode", algoData.db_mode)]; end if showBestLine - plot(ax, algoData.grossrate_Gbps, algoData.BER_plot, ... + bestHandle = plot(ax, algoData.grossrate_Gbps, algoData.BER_plot, ... "LineStyle", style.lineStyle, ... "Marker", style.marker, ... "MarkerSize", 5, ... @@ -161,6 +165,10 @@ for pamIdx = 1:numel(selectedPamLevels) "MarkerFaceColor", style.markerFaceColor, ... "MarkerEdgeColor", style.color, ... "DisplayName", style.name); + bestHandle.DataTipTemplate.DataTipRows = [ ... + dataTipTextRow("Gross rate [Gb/s]", algoData.grossrate_Gbps); ... + dataTipTextRow("BER", algoData.BER_plot); ... + dataTipTextRow("db_mode", algoData.db_mode)]; end end 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 f2a2890..6c20474 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 @@ -81,7 +81,7 @@ if isempty(plotData) return end -plotData.bitrate_Gbps = plotData.bitrate .* 1e-9; +plotData.bitrate_Gbps = plotData.grossrate .* 1e-9; plotData = sortrows(plotData, ... ["pam_level", "wavelength", "detection_type", "bitrate_Gbps", "run_id"]); 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 index e387475..32d7832 100644 --- 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 @@ -4,7 +4,7 @@ % 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. +% BER, NGMI, GMI, AIR, and all requested FEC NDR curves. % NDR values are calculated from the measured BER/NGMI using % TransmissionPerformance. @@ -13,7 +13,7 @@ clear; clc; %% 1) Query data selectedPamLevels = [4, 6, 8]; -selectedFiberLengthKm = 10; +selectedFiberLengthKm = 2; selectedWavelengthNm = 1310; selectedRopAttenuation = 0; selectedIsMpi = 0; @@ -21,7 +21,7 @@ selectedIsMpi = 0; normalDbModes = [double(db_mode.no_db), double(db_mode.db_precoded)]; duobinaryDbMode = double(db_mode.db_encoded); -showRawEntries = false; +showRawEntries = true; % set true to overlay all candidate rows maxBerForPlot = 0.5; algoStyles = defaultAlgorithmStyles(); @@ -76,6 +76,9 @@ end 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"),:); + % Keep the table variable type compatible with the normal rows before + % concatenating both result sets. + duobinaryRows.date_of_processing = string(duobinaryRows.date_of_processing); else warning("plot_best_metrics:NoProcessingDate", ... "date_of_processing was not returned; no duobinary date filtering was applied."); @@ -94,7 +97,8 @@ 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"])); +disp(groupcounts(plotData, ... + ["algorithm_key", "pam_level", "pre_emphasis", "precode"])); %% 3) Calculate BER-/NGMI-dependent net rates @@ -108,30 +112,48 @@ 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; +plotData.NDR_STAIR = columnVector(ndr.STAIR.NetRate) .* 1e-9; +plotData.NDR_HD = columnVector(ndr.HD.NetRate) .* 1e-9; +plotData.NDR_KP4 = columnVector(ndr.KP4.NetRate) .* 1e-9; +plotData.NDR_KP4_HAMMING = columnVector(ndr.KP4_hamming.NetRate) .* 1e-9; +plotData.NDR_O_FEC = columnVector(ndr.O_FEC.NetRate) .* 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"}); +metricDefinitions = struct; +metricDefinitions(1).fields = "BER_plot"; +metricDefinitions(1).label = "BER"; +metricDefinitions(1).seriesLabels = "BER"; +metricDefinitions(1).seriesLineStyles = "-"; +metricDefinitions(2).fields = "NGMI"; +metricDefinitions(2).label = "NGMI"; +metricDefinitions(2).seriesLabels = "NGMI"; +metricDefinitions(2).seriesLineStyles = "-"; +metricDefinitions(3).fields = "GMI"; +metricDefinitions(3).label = "GMI [bit/sym]"; +metricDefinitions(3).seriesLabels = "GMI"; +metricDefinitions(3).seriesLineStyles = "-"; +metricDefinitions(4).fields = "AIR_Gbps"; +metricDefinitions(4).label = "AIR [Gb/s]"; +metricDefinitions(4).seriesLabels = "AIR"; +metricDefinitions(4).seriesLineStyles = "-"; +metricDefinitions(5).fields = ["NDR_SDHD", "NDR_STAIR", "NDR_HD", ... + "NDR_KP4", "NDR_KP4_HAMMING", "NDR_O_FEC"]; +metricDefinitions(5).label = "NDR [Gb/s]"; +metricDefinitions(5).seriesLabels = ["SD+HD", "Staircase", "HD-FEC", ... + "KP4", "KP4 + Hamming", "O-FEC"]; +metricDefinitions(5).seriesLineStyles = ["-"; "--"; ":"; "-."; "-"; "--"]; bestMetricData = cell(numel(metricDefinitions), 1); for metricIdx = 1:numel(metricDefinitions) - bestMetricData{metricIdx} = bestMetricRows(plotData, ... - metricDefinitions(metricIdx).field); + bestMetricData{metricIdx} = cell(numel(metricDefinitions(metricIdx).fields), 1); + for seriesIdx = 1:numel(metricDefinitions(metricIdx).fields) + bestMetricData{metricIdx}{seriesIdx} = bestMetricRows(plotData, ... + metricDefinitions(metricIdx).fields(seriesIdx)); + end end -availableStyles = algoStyles(hasAlgorithmRows(plotData, algoStyles), :); +availableStyles = algoStyles(hasPlotStyleRows(plotData, algoStyles), :); if isempty(availableStyles) warning("plot_best_metrics:NoSelectedAlgorithms", ... "None of the configured algorithm styles match the queried rows."); @@ -147,53 +169,60 @@ t = tiledlayout(fig, numel(metricDefinitions), numel(selectedPamLevels), ... 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 seriesIdx = 1:numel(metric.fields) + metricData = bestMetricData{metricIdx}{seriesIdx}; + metricField = metric.fields(seriesIdx); + 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 + for styleIdx = 1:height(availableStyles) + style = availableStyles(styleIdx, :); + rowMask = pamMask & ... + metricData.algorithm_key == style.algorithm_key & ... + metricData.pre_emphasis == style.pre_emphasis; + if ~any(rowMask) + continue + end - algoData = sortrows(metricData(rowMask, :), "symbolrate_GBd"); - y = algoData.(metric.field); - valid = isfinite(y); - if ~any(valid) - continue - end + algoData = sortrows(metricData(rowMask, :), "symbolrate_GBd"); + y = algoData.(char(metricField)); + valid = isfinite(y); + if ~any(valid) + continue + end - if showRawEntries - scatter(ax, algoData.symbolrate_GBd(valid), y(valid), ... - 9, ... - "Marker", ".", ... + if isequal(showRawEntries, true) + 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", metric.seriesLineStyles(seriesIdx), ... + "Marker", style.marker, ... + "MarkerSize", 4, ... + "LineWidth", 1.35, ... + "Color", style.color, ... + "MarkerFaceColor", style.markerFaceColor, ... "MarkerEdgeColor", style.color, ... - "MarkerFaceColor", style.color, ... - "MarkerEdgeAlpha", 0.25, ... - "MarkerFaceAlpha", 0.25, ... - "HandleVisibility", "off"); + "DisplayName", style.displayName); 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"); + elseif metricIdx == numel(metricDefinitions) && pamIdx == 1 + addNdrLegend(ax, metric.seriesLabels, metric.seriesLineStyles); end end end @@ -229,17 +258,27 @@ 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"]); +algorithmKey = repelem(["vnle"; "vnle_pf_mlse"; "vnle_db_mlse"; ... + "ml_mlse"; "db_encoded"], 2, 1); +name = repelem(["VNLE"; "VNLE + PF + MLSE"; "VNLE DBt. + MLSE"; ... + "ML pre-EQ + Viterbi"; "DBS + VNLE + MLSE"], 2, 1); +preEmphasis = repmat([false; true], 5, 1); +preEmphasisLabel = repmat(["w/o pre-emph."; "w/ pre-emph."], 5, 1); +preEmphasisLabel(algorithmKey == "db_encoded") = "duobinary"; + +styles = table(algorithmKey, name, preEmphasis, preEmphasisLabel, ... + repelem(["o"; "square"; "diamond"; "^"; "v"], 2, 1), ... + repmat("-", 10, 1), ... + repmat("w", 10, 1), ... + [clr.Paired.lred; clr.Paired.dred; ... + clr.Paired.lgreen; clr.Paired.dgreen; ... + clr.Paired.lblue; clr.Paired.dblue; ... + clr.Paired.llila; clr.Paired.dlila; ... + clr.Paired.lorange; clr.Paired.dorange], ... + 'VariableNames', ["algorithm_key", "name", "pre_emphasis", ... + "preEmphasisLabel", "marker", "lineStyle", ... + "markerFaceColor", "color"]); +styles.displayName = styles.name + "; " + styles.preEmphasisLabel; end function plotData = buildNormalMetricRows(data) @@ -247,6 +286,7 @@ baseRows = data(isfinite(data.BER), :); baseRows.precode = false(height(baseRows), 1); baseRows.BER_plot = baseRows.BER; baseRows.algorithm_key = algorithmKeyFromEqualizer(baseRows.equalizer_structure); +baseRows.pre_emphasis = baseRows.db_mode == double(db_mode.db_precoded); baseRows = baseRows(baseRows.algorithm_key ~= "", :); if ismember("BER_precoded", string(data.Properties.VariableNames)) @@ -255,6 +295,7 @@ if ismember("BER_precoded", string(data.Properties.VariableNames)) precodedRows.BER_plot = precodedRows.BER_precoded; precodedRows.algorithm_key = algorithmKeyFromEqualizer( ... precodedRows.equalizer_structure); + precodedRows.pre_emphasis = precodedRows.db_mode == double(db_mode.db_precoded); precodedRows = precodedRows(precodedRows.algorithm_key ~= "", :); plotData = [baseRows; precodedRows]; else @@ -269,6 +310,7 @@ 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); +plotData.pre_emphasis = false(height(plotData), 1); end function plotData = addDerivedMetrics(plotData) @@ -343,7 +385,7 @@ if isempty(candidateData) return end -groupVars = ["pam_level", "algorithm_key", "symbolrate_GBd"]; +groupVars = ["pam_level", "algorithm_key", "pre_emphasis", "symbolrate_GBd"]; groupId = findgroups(candidateData(:, groupVars)); keepIdx = zeros(max(groupId), 1); for curGroup = 1:max(groupId) @@ -359,10 +401,11 @@ end bestData = sortrows(candidateData(keepIdx, :), groupVars); end -function keep = hasAlgorithmRows(data, algoStyles) +function keep = hasPlotStyleRows(data, algoStyles) keep = false(height(algoStyles), 1); for idx = 1:height(algoStyles) - keep(idx) = any(data.algorithm_key == algoStyles.algorithm_key(idx)); + keep(idx) = any(data.algorithm_key == algoStyles.algorithm_key(idx) & ... + data.pre_emphasis == algoStyles.pre_emphasis(idx)); end end @@ -376,7 +419,7 @@ grid(ax, "on"); grid(ax, "minor"); box(ax, "on"); -switch metric.field +switch metric.fields(1) case "BER_plot" set(ax, "YScale", "log"); ylim(ax, [1e-5 0.2]); @@ -405,3 +448,15 @@ else title(ax, sprintf("PAM-%d", pamLevel)); end end + +function addNdrLegend(ax, seriesLabels, seriesLineStyles) +legendHandles = gobjects(numel(seriesLabels), 1); +for idx = 1:numel(seriesLabels) + legendHandles(idx) = plot(ax, NaN, NaN, ... + "Color", [0.2 0.2 0.2], ... + "LineStyle", seriesLineStyles(idx), ... + "LineWidth", 1.35, ... + "DisplayName", seriesLabels(idx)); +end +legend(ax, legendHandles, "Location", "southwest", "Interpreter", "none"); +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 index 4f84fc1..58d743b 100644 --- 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 @@ -9,9 +9,9 @@ clear; clc; %% 1) Query data selectedPamLevels = [4, 6, 8]; -selectedFiberLengthKm = 1; +selectedFiberLengthKm = 10; selectedWavelengthNm = 1310; -selectedBitrateGbps = 360; +selectedBitrateGbps = 420; selectedIsMpi = 0; % set [] to use all entries maxPowerPdIn = []; % set a numeric limit to enable 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 index a9657e2..78d05aa 100644 --- 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 @@ -9,7 +9,7 @@ clear; clc; %% 1) Query data selectedPamLevels = [4, 6, 8]; -selectedFiberLengthKm = 5; +selectedFiberLengthKm = 10; 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 diff --git a/projects/Diss/400G_revisit/RX_sprectra.m b/projects/Diss/400G_revisit/RX_sprectra.m index bb39f77..0b42458 100644 --- a/projects/Diss/400G_revisit/RX_sprectra.m +++ b/projects/Diss/400G_revisit/RX_sprectra.m @@ -8,7 +8,7 @@ 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.max_occurences = 16; dsp_options.debug_plots = false; @@ -33,65 +33,148 @@ db = DBHandler("dataBase", [dsp_options.dataBase], ... "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); +cols = cbrewer2('Spectral',6); +cols = linspecer(2); +len = [1 2 3 5 6 8 10]; +idx = 1; +for lambda = 1302%[1297.5] + + fp = QueryFilter(); + len = 10; + fp.where('Runs','fiber_length','EQUALS', 10); % 1 2 3 5 6 8 10 + fp.where('Runs','wavelength','EQUALS', lambda); % + fp.where('Runs','bitrate','EQUALS', 450e9); + 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, S, ~] = loadAndSyncRunSignals(dataTable(1,:), dsp_options); + + average_signals = 1; + if average_signals + Scpe_sig_avg = S{1}; + scope_mean = zeros(size(S{1}.signal)); + for n=1:numel(S) + scope_mean = scope_mean + S{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig_avg.signal = scope_mean; + ScopeSignal_preemph = preprocessSignal(Scpe_sig_avg, Symbols_preemph, Symbols_preemph.fs,"gaussian_cutoff_factor",0.8); + else + ScopeSignal = S{1}; + ScopeSignal_preemph = preprocessSignal(ScopeSignal, Symbols_preemph, Symbols_preemph.fs); + end + + fignum = 2; + dn = sprintf("Rx Spectrum; %d nm",floor(lambda)); + ScopeSignal_preemph.spectrum("displayname",dn,'fignum',fignum,'normalizeTo0dB',1,'fft_length',2^15,'show_onesided',1,'color',cols(idx,:)); -[~, Symbols_preemph, Scpe_cell_preemph, ~] = loadAndSyncRunSignals(dataTable(1,:), dsp_options); -ScopeSignal = Scpe_cell_preemph{1}; -ScopeSignal_preemph = preprocessSignal(ScopeSignal, Symbols_preemph, Symbols_preemph.fs); + if len == 10 + + fp.where('Runs', 'db_mode','EQUALS', 2); + fields = db.getTableFieldNames('Runs'); + [dataTable, query] = db.queryDB(fp, fields); + + [~, Symbols_db, S, found_sync] = loadAndSyncRunSignals(dataTable(1,:), dsp_options); + ScopeSignal = S{1}; + average_signals = 1; + if average_signals + Scpe_sig_avg = S{1}; + scope_mean = zeros(size(S{1}.signal)); + for n=1:numel(S) + scope_mean = scope_mean + S{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig_avg.signal = scope_mean; + ScopeSignal = preprocessSignal(Scpe_sig_avg, Symbols_preemph, Symbols_preemph.fs,"gaussian_cutoff_factor",0.8); + else + ScopeSignal = S{1}; + ScopeSignal = preprocessSignal(ScopeSignal, Symbols_preemph, Symbols_preemph.fs); + end + + ScopeSignal_DB = preprocessSignal(ScopeSignal, Symbols_db, Symbols_db.fs); + end + + if len == 10 + fignum = 2; + dn = sprintf("DBS Rx Spectrum; %d nm",floor(lambda)); + ScopeSignal_DB.spectrum("displayname",dn,'fignum',fignum,'normalizeTo0dB',1,'fft_length',2^15,'show_onesided',1,'color',cols(idx+1,:)); + % Symbols_db.spectrum("displayname",dn,'fignum',fignum,'normalizeTo0dB',1,'fft_length',1024,'show_onesided',1,'color',cols(idx+1,:)); + end + + + + % Fiber and system parameters + lambda0 = 1314e-9; % zero-dispersion wavelength [m] + lambda = lambda*1e-9; % operating wavelength [m] + S0 = 0.092; % dispersion slope [ps/(nm²·km)] + L = 10e3; % fiber length [m] + c = physconst('lightspeed'); + + % Derived quantities + S0_si = S0 * 1e3; % → s/m³ + D_lambda = (S0/4) * (lambda*1e9 - (lambda0*1e9)^4/(lambda*1e9)^3); % ps/(nm·km) + D_si = D_lambda * 1e-6; % → s/m² + b2 = -D_si * lambda^2 / (2*pi*c); % s²/m + + Dacc = D_lambda * L; + fprintf('Accumulated Dispersion: %.2f ps/nm \n', Dacc / 1e3); + + % Frequency grid + f_max = 200e9; + f = linspace(0, f_max, 5000); % [Hz] + + % IM/DD transfer function (power fading) + phi = 2*pi^2 * b2 * f.^2 * L; + H = abs(cos(phi)); + + % Plot + plot(f/1e9, 10*log10(H), 'LineWidth', 1,'Color','black','DisplayName','CD Transfer function','HandleVisibility','on'); + grid on; box on; + xlabel('Frequency [GHz]'); + ylabel('Magnitude [dB]'); + % title(sprintf('IM/DD Power Fading: 10 km; 1275nm', lambda*1e9, L/1000),"Interpreter","latex"); + % ylim([-20 0]); + + idx = idx +1; + +end %% 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); +fp.where('Runs', 'db_mode','EQUALS', 0); 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); +fignum = len; +ScopeSignal_preemph.spectrum("displayname",'Full Response w/ preemphasis','fignum',fignum,'normalizeTo0dB',0,'color',clr.Paired.dgreen,'fft_length',4096*6); %% 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); +if len == 10 + 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); +end -fields = db.getTableFieldNames('Runs'); -[dataTable, query] = db.queryDB(fp, fields); +if len == 10 + fignum = 2; + dn = sprintf("DBS Rx Spectrum; %d nm",floor(lambda)); + ScopeSignal_DB.spectrum("displayname",dn,'fignum',fignum,'normalizeTo0dB',1,'fft_length',2^15,'show_onesided',1,'color',cols(idx,:)); +end -[~, 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 +%% \ No newline at end of file diff --git a/projects/Diss/400G_revisit/TABLE_BEST_DATARATES_2KM_10KM_1310NM.m b/projects/Diss/400G_revisit/TABLE_BEST_DATARATES_2KM_10KM_1310NM.m new file mode 100644 index 0000000..94328c3 --- /dev/null +++ b/projects/Diss/400G_revisit/TABLE_BEST_DATARATES_2KM_10KM_1310NM.m @@ -0,0 +1,261 @@ +%% Best AIR and FEC data rates over all wavelengths for 2 km and 10 km +% The table follows the contender mapping used in the thesis figures: +% PAM-4: DB target + MLSE, with DBS + VNLE + MLSE at 10 km +% PAM-6: VNLE + PF + MLSE, with DBS + VNLE + MLSE at 10 km +% PAM-8: VNLE + PF + MLSE +% +% For every configured contender, the maximum valid rate is selected across +% the available gross rates and symbol rates. BER_precoded is used for the +% DB-target curve, as in FIGURE_NGMI_THESIS_FINAL.m. + +clear; + +%% Configuration + +normalFiberLengthsKm = [2 10]; +duobinaryFiberLengthKm = 10; +selectedRopAttenuation = 0; +selectedIsMpi = 0; +duobinaryDateCutoff = datetime("2026-01-01 00:00:00"); + +curves = struct( ... + "name", {"DBt. + MLSE"; "DBS + VNLE + MLSE"; ... + "VNLE + PF + MLSE"; "DBS + VNLE + MLSE"; ... + "VNLE + PF + MLSE"}, ... + "algorithm_key", {"vnle_db_mlse"; "db_encoded"; ... + "vnle_pf_mlse"; "db_encoded"; ... + "vnle_pf_mlse"}, ... + "pam", {4; 4; 6; 6; 8}, ... + "source", {"normal"; "duobinary"; "normal"; "duobinary"; "normal"}, ... + "usePrecodedBer", {true; false; false; false; false}); + +%% Query the normal and DBS data + +db = DBHandler( ... + "dataBase", "labor_highspeed", ... + "type", "mysql", ... + "server", "192.168.178.192", ... + "user", "silas", ... + "password", "silas"); +db.refresh(); +selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime'); + +normalRows = cell(size(normalFiberLengthsKm)); +for lengthIdx = 1:numel(normalFiberLengthsKm) + normalRows{lengthIdx} = queryRows(db, selectedFields, ... + normalFiberLengthsKm(lengthIdx), ... + selectedRopAttenuation, selectedIsMpi); + normalRows{lengthIdx} = cleanRows(normalRows{lengthIdx}); + normalRows{lengthIdx} = normalRows{lengthIdx}(ismember( ... + normalRows{lengthIdx}.db_mode, ... + [double(db_mode.no_db), double(db_mode.db_precoded)]), :); + normalRows{lengthIdx}.algorithm_key = ... + lower(string(normalRows{lengthIdx}.equalizer_structure)); + fprintf("Normal %g km rows: %d\n", normalFiberLengthsKm(lengthIdx), ... + height(normalRows{lengthIdx})); +end + +duobinaryRows = queryRows(db, selectedFields, duobinaryFiberLengthKm, ... + selectedRopAttenuation, selectedIsMpi); +duobinaryRows = cleanRows(duobinaryRows); +duobinaryRows = duobinaryRows(duobinaryRows.db_mode == ... + double(db_mode.db_encoded), :); +duobinaryRows = duobinaryRows( ... + duobinaryRows.equalizer_structure == "db_encoded", :); +if ismember("date_of_processing", string(duobinaryRows.Properties.VariableNames)) + processedDate = datetime(string(duobinaryRows.date_of_processing)); + duobinaryRows = duobinaryRows(processedDate > duobinaryDateCutoff, :); +else + warning("table_best_rates:NoProcessingDate", ... + "date_of_processing was not returned; no DBS date filter was applied."); +end +duobinaryRows.algorithm_key = repmat("db_encoded", height(duobinaryRows), 1); +fprintf("DBS %g km rows after date filter: %d\n", ... + duobinaryFiberLengthKm, height(duobinaryRows)); + +%% Calculate maximum rates + +tp = TransmissionPerformance; +resultRows = repmat(emptyResultRow(), 0, 1); + +for lengthIdx = 1:numel(normalFiberLengthsKm) + for curveIdx = 1:numel(curves) + curve = curves(curveIdx); + if curve.source ~= "normal" + continue + end + resultRows(end+1) = evaluateCurve(normalRows{lengthIdx}, curve, ... + normalFiberLengthsKm(lengthIdx), tp); %#ok + end +end + +for curveIdx = 1:numel(curves) + curve = curves(curveIdx); + if curve.source ~= "duobinary" + continue + end + resultRows(end+1) = evaluateCurve(duobinaryRows, curve, ... + duobinaryFiberLengthKm, tp); %#ok +end + +ratesTable = struct2table(resultRows); +ratesTable = ratesTable(:, ... + ["pam", "fiber_length_km", "technique", "AIR_Gbps", ... + "SDHD_Gbps", "HD_Gbps", "AIR_grossrate_Gbps", ... + "SDHD_grossrate_Gbps", "HD_grossrate_Gbps", ... + "AIR_symbolrate_GBd", "SDHD_symbolrate_GBd", ... + "HD_symbolrate_GBd"]); + +disp("Best rates (Gb/s):"); +disp(ratesTable(:, ["pam", "fiber_length_km", "technique", ... + "AIR_Gbps", "SDHD_Gbps", "HD_Gbps"])); + +displayTable = table( ... + ratesTable.pam(:), ... + ratesTable.fiber_length_km(:), ... + string(ratesTable.technique(:)), ... + formatRateWithSymbol(ratesTable.AIR_Gbps(:), ratesTable.AIR_symbolrate_GBd(:)), ... + formatRateWithSymbol(ratesTable.SDHD_Gbps(:), ratesTable.SDHD_symbolrate_GBd(:)), ... + formatRateWithSymbol(ratesTable.HD_Gbps(:), ratesTable.HD_symbolrate_GBd(:)), ... + 'VariableNames', {'PAM', 'Distance_km', 'EQ_contender', ... + 'AIR', 'SDHD_FEC', 'HD_FEC'}); +disp("Best rates with the symbol rate used for each metric:"); +disp(displayTable); + +%% Local functions + +function T = queryRows(db, fields, fiberLengthKm, ropAttenuation, isMpi) +fp = QueryFilter(); +fp.where("Runs", "fiber_length", "EQUALS", fiberLengthKm); +fp.where("Runs", "rop_attenuation", "EQUALS", ropAttenuation); +fp.where("Runs", "is_mpi", "EQUALS", isMpi); +[T, query] = db.queryDB(fp, fields); +disp(query); +end + +function T = cleanRows(T) +numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ... + "symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ... + "rop_attenuation", "numBits", "numBitErr", "BER", ... + "numBitErr_precoded", "BER_precoded", "GMI", "AIR", "NGMI"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(T.Properties.VariableNames)) + T.(char(fieldName)) = numericColumn(T.(char(fieldName))); + end +end +T.equalizer_structure = lower(string(T.equalizer_structure)); +end + +function result = evaluateCurve(sourceRows, curve, fiberLengthKm, tp) +result = emptyResultRow(); +result.pam = curve.pam; +result.fiber_length_km = fiberLengthKm; +result.technique = curve.name; + +curveRows = sourceRows(sourceRows.pam_level == curve.pam & ... + sourceRows.algorithm_key == curve.algorithm_key, :); +if isempty(curveRows) + warning("table_best_rates:NoRows", ... + "No rows found for %s, PAM-%d, %g km.", ... + curve.name, curve.pam, fiberLengthKm); + return +end + +curveRows.BER_plot = numericOrNaN(curveRows, "BER"); +if curve.usePrecodedBer && ismember("BER_precoded", ... + string(curveRows.Properties.VariableNames)) + precodedBer = numericOrNaN(curveRows, "BER_precoded"); + usePrecoded = isfinite(precodedBer); + curveRows.BER_plot(usePrecoded) = precodedBer(usePrecoded); +end + +curveRows = addAirMetric(curveRows); +grossRate = numericOrNaN(curveRows, "grossrate"); +ngmi = numericOrNaN(curveRows, "NGMI"); +ber = numericOrNaN(curveRows, "BER_plot"); +ngmi(~isfinite(ngmi) | ngmi < 0 | ngmi > 1.05) = NaN; +ber(~isfinite(ber) | ber <= 0 | ber > 0.5) = NaN; + +ndr = tp.calculateNetRate(grossRate, "NGMI", ngmi, "BER", ber); +curveRows.SDHD_Gbps = columnVector(ndr.SDHD.NetRate) .* 1e-9; +curveRows.HD_Gbps = columnVector(ndr.O_FEC.NetRate) .* 1e-9; + +result = assignMaximum(result, curveRows, "AIR_Gbps"); +result = assignMaximum(result, curveRows, "SDHD_Gbps"); +result = assignMaximum(result, curveRows, "HD_Gbps"); +end + +function result = assignMaximum(result, T, fieldName) +values = numericOrNaN(T, fieldName); +valid = isfinite(values) & values >= 0; +if ~any(valid) + return +end +[maximum, localIndex] = max(values(valid)); +validRows = find(valid); +rowIndex = validRows(localIndex); +prefix = erase(fieldName, "_Gbps"); +result.(fieldName) = maximum; +result.(char(prefix + "_grossrate_Gbps")) = ... + numericOrNaN(T(rowIndex, :), "grossrate") .* 1e-9; +result.(char(prefix + "_symbolrate_GBd")) = ... + numericOrNaN(T(rowIndex, :), "symbolrate") .* 1e-9; +end + +function T = addAirMetric(T) +T.AIR_Gbps = numericOrNaN(T, "AIR") .* 1e-9; +gmi = numericOrNaN(T, "GMI"); +fallbackAir = gmi .* numericOrNaN(T, "symbolrate") .* 1e-9; +grossRate = numericOrNaN(T, "grossrate"); +useFallback = ~isfinite(T.AIR_Gbps) | T.AIR_Gbps < 0 | ... + (isfinite(grossRate) & T.AIR_Gbps > grossRate .* 1.05e-9); +T.AIR_Gbps(useFallback) = fallbackAir(useFallback); +end + +function result = emptyResultRow() +result = struct( ... + "pam", NaN, ... + "fiber_length_km", NaN, ... + "technique", "", ... + "AIR_Gbps", NaN, ... + "SDHD_Gbps", NaN, ... + "HD_Gbps", NaN, ... + "AIR_grossrate_Gbps", NaN, ... + "SDHD_grossrate_Gbps", NaN, ... + "HD_grossrate_Gbps", NaN, ... + "AIR_symbolrate_GBd", NaN, ... + "SDHD_symbolrate_GBd", NaN, ... + "HD_symbolrate_GBd", NaN); +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 = numericOrNaN(T, fieldName) +if ismember(fieldName, string(T.Properties.VariableNames)) + values = numericColumn(T.(char(fieldName))); +else + values = NaN(height(T), 1); +end +values = values(:); +end + +function values = columnVector(values) +values = double(values(:)); +end + +function text = formatRateWithSymbol(rateGbps, symbolrateGBd) +text = strings(size(rateGbps)); +valid = isfinite(rateGbps) & isfinite(symbolrateGBd); +text(valid) = compose("%.0f (%.0f GBd)", ... + rateGbps(valid), symbolrateGBd(valid)); +text(~valid) = "--"; +end diff --git a/projects/Diss/400G_revisit/TX_spectra.m b/projects/Diss/400G_revisit/TX_spectra.m index df869bd..70779ef 100644 --- a/projects/Diss/400G_revisit/TX_spectra.m +++ b/projects/Diss/400G_revisit/TX_spectra.m @@ -1,5 +1,5 @@ -rates = [420e9]; +rates = [400e9]; rcalpha = 0.05; fsym = rates/2; apply_pulsef = 1; @@ -56,8 +56,6 @@ Digi_sig_rx.spectrum(... - - Digi_sig_pre_rx.spectrum(... "displayname",'Rx w/ pre-emphasis',... "fignum",2,"normalizeTo0dB",0,"color",clr.Paired.dblue); @@ -80,6 +78,9 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"al "mrds_code",0,"mrds_blocklength",512).process(); Digi_sig_DB= Digi_sig_DB.normalize("mode","oneone"); +Digi_sig_DB_rx = precomp_est.apply(... + Digi_sig_DB,'maxampdb',0,'loadPath',precomp_path,'fileName',precomp_fn); + 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"; @@ -89,27 +90,33 @@ 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"); +Digi_sig_DB_pre_rx = precomp_est.apply(... + Digi_sig_DB_pre,'maxampdb',0,'loadPath',precomp_path,'fileName',precomp_fn); + %% 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_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_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_rx = AWG_.process(Digi_sig_DB_pre_rx); +Digi_sig_DB.spectrum("displayname",'DB Response w/o preemphasis','fignum',1,'normalizeTo0dB',1,'color',clr.Paired.lorange); +Digi_sig_DB_pre.spectrum("displayname",'DB Response w/ preemphasis','fignum',1,'normalizeTo0dB',1,'color',clr.Paired.dorange); -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); +Digi_sig_DB_rx.spectrum("displayname",'DB Response w/o preemphasis at Rx','fignum',1,'normalizeTo0dB',1,'color',clr.Paired.lred); +Digi_sig_DB_pre_rx.spectrum("displayname",'DB Response w/ preemphasis at Rx','fignum',1,'normalizeTo0dB',1,'color',clr.Paired.dred); \ No newline at end of file diff --git a/projects/Diss/400G_revisit/investigate_400g_algorithms.m b/projects/Diss/400G_revisit/investigate_400g_algorithms.m index a8d0dc6..882ba74 100644 --- a/projects/Diss/400G_revisit/investigate_400g_algorithms.m +++ b/projects/Diss/400G_revisit/investigate_400g_algorithms.m @@ -5,7 +5,7 @@ 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.max_occurences = 1; dsp_options.debug_plots = false; @@ -36,11 +36,12 @@ 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','LESS_THAN', 480e9); -% fp.where('Runs','pam_level','EQUALS', 4); +% fp.where('Runs','bitrate','LESS_THAN', 330e9); +fp.where('Runs','symbolrate','EQUALS', 174e9); +fp.where('Runs','pam_level','EQUALS', 6); 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', 1); fields = db.getTableFieldNames('Runs'); [dataTable, query] = db.queryDB(fp, fields); @@ -67,7 +68,7 @@ dsp_options.userParameters = struct(); % dsp_options.userParameters.len_tr = 4096*2; % dsp_options.userParameters.pf_ncoeffs = [1,2,3]; -dsp_options.userParameters.decoding_mode = [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]; @@ -86,7 +87,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.parallel, ... +[results, wh] = submitJobs(run_ids, dsp_options, processingMode.serial, ... "wh", wh, ... "waitbar", true); diff --git a/projects/MLSE_ML_based/minimal_example_huawei/ml_mlse_pam.m b/projects/MLSE_ML_based/minimal_example_huawei/ml_mlse_pam.m index f33d40d..3de3cab 100644 --- a/projects/MLSE_ML_based/minimal_example_huawei/ml_mlse_pam.m +++ b/projects/MLSE_ML_based/minimal_example_huawei/ml_mlse_pam.m @@ -210,8 +210,8 @@ classdef ml_mlse_pam < handle % ============================================================== % ML-Based Branch Metric Estimation + Viterbi % ============================================================== - % debug = 1; - % showPlots = 1; + debug = 0; + showPlots = 0; nSymbols = ceil(N/obj.sps); diff --git a/projects/MLSE_ML_based/theoretic_channel_evaluation.m b/projects/MLSE_ML_based/theoretic_channel_evaluation.m index cdcdc81..a0ad445 100644 --- a/projects/MLSE_ML_based/theoretic_channel_evaluation.m +++ b/projects/MLSE_ML_based/theoretic_channel_evaluation.m @@ -47,7 +47,7 @@ ber_ml_mlse_l4 = zeros(size(SNR_dB)); epochs_training = 100; -parfor i = 1:numel(SNR_dB) +for i = 1:numel(SNR_dB) symbols_noi = symbols_filt; symbols_noi.signal = awgn(symbols_filt.signal, SNR_dB(i), 'measured'); % AWGN with given SNR @@ -73,46 +73,46 @@ parfor i = 1:numel(SNR_dB) [~, ~, ber_ffe(i), ~] = calc_ber(Eq_bits.signal, Bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); fprintf('FFE: %.2e \n',ber_ffe(i)); - % Postfilter - [y_white,~] = pf_.process(y_ffe, ffe_noise); - - % Sequence Est - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',pf_.coefficients); - [y_mlse] = mlse_.process(y_white,Symbols); - mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); - [~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); - fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i)); - - % ML-base MLSE L=2 - adaptive_mu = 0; - mu_lms = 0.15; - ml_mlse_equalizer = ML_MLSE("epochs_tr",epochs_training,"epochs_dd",1,"len_tr",2^15,... - "mu_dd",mu_lms,"mu_tr",mu_lms,"order",11,"sps",1,... - "traceback_depth",128,"L",2,"delta",4,"adaptive_mu",adaptive_mu); - [y_ml_mlse,~] = ml_mlse_equalizer.process(symbols_noi,Symbols); - ml_mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ml_mlse); - [~, errors, ber_ml_mlse_l2(i), errpos] = calc_ber(ml_mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); - fprintf('ML MLSE BER: %.2e \n',ber_ml_mlse_l2(i)); - - % ML-base MLSE L=3 - mu_lms = 0.15; - ml_mlse_equalizer = ML_MLSE("epochs_tr",epochs_training,"epochs_dd",1,"len_tr",2^16,... - "mu_dd",mu_lms,"mu_tr",mu_lms,"order",11,"sps",1,... - "traceback_depth",128,"L",3,"delta",4,"adaptive_mu",adaptive_mu); - [y_ml_mlse,~] = ml_mlse_equalizer.process(symbols_noi,Symbols); - ml_mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ml_mlse); - [~, errors, ber_ml_mlse_l3(i), errpos] = calc_ber(ml_mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); - fprintf('ML MLSE BER: %.2e \n',ber_ml_mlse_l3(i)); - - % % ML-base MLSE L=5 + % % Postfilter + % [y_white,~] = pf_.process(y_ffe, ffe_noise); + % + % % Sequence Est + % mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',pf_.coefficients); + % [y_mlse] = mlse_.process(y_white,Symbols); + % mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse); + % [~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + % fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i)); + % + % % ML-base MLSE L=2 + % adaptive_mu = 0; % mu_lms = 0.15; % ml_mlse_equalizer = ML_MLSE("epochs_tr",epochs_training,"epochs_dd",1,"len_tr",2^15,... % "mu_dd",mu_lms,"mu_tr",mu_lms,"order",11,"sps",1,... - % "traceback_depth",128,"L",5,"delta",4); + % "traceback_depth",128,"L",2,"delta",4,"adaptive_mu",adaptive_mu); % [y_ml_mlse,~] = ml_mlse_equalizer.process(symbols_noi,Symbols); % ml_mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ml_mlse); - % [~, errors, ber_ml_mlse_l5(i), errpos] = calc_ber(ml_mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); - % fprintf('ML MLSE BER: %.2e \n',ber_ml_mlse_l5(i)); + % [~, errors, ber_ml_mlse_l2(i), errpos] = calc_ber(ml_mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + % fprintf('ML MLSE BER: %.2e \n',ber_ml_mlse_l2(i)); + % + % % ML-base MLSE L=3 + % mu_lms = 0.15; + % ml_mlse_equalizer = ML_MLSE("epochs_tr",epochs_training,"epochs_dd",1,"len_tr",2^16,... + % "mu_dd",mu_lms,"mu_tr",mu_lms,"order",11,"sps",1,... + % "traceback_depth",128,"L",3,"delta",4,"adaptive_mu",adaptive_mu); + % [y_ml_mlse,~] = ml_mlse_equalizer.process(symbols_noi,Symbols); + % ml_mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ml_mlse); + % [~, errors, ber_ml_mlse_l3(i), errpos] = calc_ber(ml_mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + % fprintf('ML MLSE BER: %.2e \n',ber_ml_mlse_l3(i)); + + % ML-base MLSE L=5 + mu_lms = 0.15; + ml_mlse_equalizer = ML_MLSE("epochs_tr",epochs_training,"epochs_dd",1,"len_tr",2^15,... + "mu_dd",mu_lms,"mu_tr",mu_lms,"order",11,"sps",1,... + "traceback_depth",128,"L",5,"delta",4); + [y_ml_mlse,~] = ml_mlse_equalizer.process(symbols_noi,Symbols); + ml_mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ml_mlse); + [~, errors, ber_ml_mlse_l5(i), errpos] = calc_ber(ml_mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); + fprintf('ML MLSE BER: %.2e \n',ber_ml_mlse_l5(i)); end diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m index 1435d1f..6456815 100644 --- a/projects/WDM/WDM_auswertung.m +++ b/projects/WDM/WDM_auswertung.m @@ -1,5 +1,5 @@ -base = "C:\Users\Silas\Nextcloud4\Cluster"; +base = "C:\Users\Silas\Nextcloud2\CAU\Cluster"; all_files = dir(fullfile(base, "**/*.mat")); schemes = ["co","pair","alt","seg"]; @@ -77,7 +77,7 @@ res_all = drop_empty_realizations(res_all); S = plot_BER_vs_ROP(res_all, 'fec', 3.8e-3); %% Routine B: violin plot (independent) -plot_FEC_violin(res_all, 'tech','VNLE', 'fec',3.8e-3, 'ylim',[-10 0],'eval_ptr',5); +plot_FEC_violin(res_all, 'tech','VNLE', 'fec',3.8e-3, 'ylim',[-10 0],'eval_ptr',10); %% diff --git a/projects/WDM/WDM_model.m b/projects/WDM/WDM_model.m index aba8593..193b0e9 100644 --- a/projects/WDM/WDM_model.m +++ b/projects/WDM/WDM_model.m @@ -6,10 +6,10 @@ function WDM_model(options) arguments options.num_channels = 16; options.channel_spacing = 400e9; - options.fiber_length_km = 0; + options.fiber_length_km = 1; options.rand_key = 1; options.num_realiz = 1; - options.fwm_mitigation_technique = "co"; + options.fwm_mitigation_technique = "pair"; end %% @@ -72,7 +72,6 @@ host = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end fname = sprintf('WDM_%s_%s_%s_%dkm_%dch_%dghz_%s.mat', char(t), host, jobid, options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique); -%% s.num_realiz = options.num_realiz; % s.wavelengthplan = calcWavelengthPlan(16,400e9,1310); s.wavelengthplan = calcWavelengthPlan(options.num_channels,options.channel_spacing,1310); @@ -162,7 +161,7 @@ end for realiz = 1:s.num_realiz - parfor l = 1:N + for l = 1:N [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource(... "fsym",fsym,"M",s.M,"order",18,"useprbs",0,... @@ -199,9 +198,9 @@ for realiz = 1:s.num_realiz Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm); - % Opt_sig_wdm.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0); + Opt_sig_wdm.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',1); - % Opt_sig_wdm.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',1,'max_num_lines',2); + %% Opt_sig_wdm.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',1,'max_num_lines',2); %%%%%% Fiber %%%%%% Opt_sig_wdm_fib=Opt_sig_wdm; @@ -219,20 +218,17 @@ for realiz = 1:s.num_realiz "beat_len",10,"corr_len",100,"dz",1,"manakov",0,... "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01,... "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib); - propdist = segment_length; - - end - %%%%%% Demux after 2 km %%%%%% + %% %%%% Demux after 2 km %%%%%% Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib); for ri = 1:length(s.rop) - parfor l = 1:N + for l = 1:N %%%%%% ROP %%%%%% Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",s.rop(ri)).process(Opt_sig_wdm_demux{l}); % rop+10*log10(N) @@ -276,14 +272,13 @@ for realiz = 1:s.num_realiz output_ffe{l,ri,realiz} = ffe_results; - - %VNLE pf_ncoeffs = 1; ffe_order = [50, 5, 5]; eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + useviterbi = 0; if useviterbi mlse_ = MLSE_viterbi("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels); diff --git a/tmp/pdfs/paper-2.png b/tmp/pdfs/paper-2.png new file mode 100644 index 0000000..c235fd3 Binary files /dev/null and b/tmp/pdfs/paper-2.png differ diff --git a/tmp/pdfs/paper-3.png b/tmp/pdfs/paper-3.png new file mode 100644 index 0000000..6d4c33f Binary files /dev/null and b/tmp/pdfs/paper-3.png differ diff --git a/workerError.mat b/workerError.mat index 73a5a0f..0325c3a 100644 Binary files a/workerError.mat and b/workerError.mat differ