From 4099f6820f7a98899aedf6599e37d97c1be9ac45 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Wed, 17 Sep 2025 13:58:58 +0200 Subject: [PATCH] BCJR implementation WDM code added (Pol Cont., Opt MUX/DEMUX, Opt Atten, DP_Fiber) -> the codebase is not optimized to always work with dp signals! --- Classes/00_signals/Opticalsignal.m | 2 + Classes/00_signals/Signal.m | 130 ++-- Classes/02_etc/Amplifier.m | 2 +- Classes/02_optical/DP_Fiber.m | 240 ++++++++ Classes/02_optical/EML.m | 2 +- Classes/02_optical/Optical_Demultiplex.m | 122 ++++ Classes/02_optical/Optical_Multiplex.m | 171 ++++++ Classes/02_optical/Polarization_Controller.m | 92 +++ Classes/04_DSP/Coding/Duobinary.m | 1 + Classes/04_DSP/Sequence Detection/MLSE.m | 261 +++++--- Classes/04_DSP/TransmissionPerformance.m | 52 +- Classes/DataBaseHandler/DBHandler.m | 184 ++++-- Classes/Warehouse_class/classes/DataStorage.m | 5 +- .../Warehouse_class/classes/minimal_example.m | 40 ++ Datatypes/polarization_control_mode.m | 10 + Functions/EQ_structures/dsp_runid.m | 81 ++- Functions/EQ_structures/duobinary_target.m | 40 +- Functions/EQ_structures/ffe.m | 5 +- .../EQ_structures/vnle_postfilter_mlse.m | 42 +- Functions/EQ_visuals/show2Dconstellation.m | 85 +++ Functions/EQ_visuals/showLevelScatter.m | 10 +- .../loadAndSyncSignalDataFromDb.m | 4 +- Functions/Metrics/calc_snr.m | 1 - Functions/Metrics/count_error_bursts.m | 4 +- .../ECOC_2025/dsp_test/hyperparam_tuning.m | 18 +- .../plots_from_database/plot_mpi_trial.m | 45 +- .../theory/analytic_mpi_evaluation.m | 82 +++ .../ECOC_2025/theory/coherence_length_plot.m | 32 + .../Auswertung_JLT/ber_vs_rate.m | 301 ++++++++-- .../Auswertung_JLT/run_dsp_from_db.m | 437 +++++++++----- projects/IMDD_base_system/simulation_bwl.m | 374 +++++++++--- projects/WDM/WDM_model.m | 268 +++++++++ test/bcjr_pam.m | 555 ++++++++++++++++++ test/duobinary_emulation_minimal_example.m | 76 +-- test/minimal_example_bcjr.m | 171 ++++++ test/pam_6_differential_code_understand.m | 79 +++ test/pam_6_states_analysis.m | 212 +++++++ 37 files changed, 3643 insertions(+), 593 deletions(-) create mode 100644 Classes/02_optical/DP_Fiber.m create mode 100644 Classes/02_optical/Optical_Demultiplex.m create mode 100644 Classes/02_optical/Optical_Multiplex.m create mode 100644 Classes/02_optical/Polarization_Controller.m create mode 100644 Classes/Warehouse_class/classes/minimal_example.m create mode 100644 Datatypes/polarization_control_mode.m create mode 100644 Functions/EQ_visuals/show2Dconstellation.m create mode 100644 projects/ECOC_2025/theory/analytic_mpi_evaluation.m create mode 100644 projects/ECOC_2025/theory/coherence_length_plot.m create mode 100644 projects/WDM/WDM_model.m create mode 100644 test/bcjr_pam.m create mode 100644 test/minimal_example_bcjr.m create mode 100644 test/pam_6_differential_code_understand.m create mode 100644 test/pam_6_states_analysis.m diff --git a/Classes/00_signals/Opticalsignal.m b/Classes/00_signals/Opticalsignal.m index 25900a4..b39511e 100644 --- a/Classes/00_signals/Opticalsignal.m +++ b/Classes/00_signals/Opticalsignal.m @@ -5,6 +5,7 @@ classdef Opticalsignal < Signal properties nase lambda + polrot end @@ -18,6 +19,7 @@ classdef Opticalsignal < Signal options.logbook options.lambda options.nase + options.polrot end obj = obj@Signal(signal); diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 9b0c3e5..c84ebf6 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -108,6 +108,7 @@ classdef Signal options.logbook options.nase options.lambda + options.polrot end fn = fieldnames(options); @@ -122,7 +123,7 @@ classdef Signal %convert to optical - o_sig = Opticalsignal(obj.signal,"fs",obj.fs,"lambda",options.lambda,"logbook",obj.logbook,"nase",options.nase); + o_sig = Opticalsignal(obj.signal,"fs",obj.fs,"lambda",options.lambda,"logbook",obj.logbook,"nase",options.nase,"polrot",options.polrot); elseif isa(obj,'Informationsignal') @@ -141,7 +142,7 @@ classdef Signal arguments obj - options.fignum = [] + options.fignum = randi(1000) options.displayname = []; options.timeframe = 0; options.clear = 0; @@ -187,7 +188,7 @@ classdef Signal ylabel('Amplitude'); end - + % Add legend if not already present if isempty(get(gca, 'Legend')) legend; @@ -350,21 +351,30 @@ classdef Signal options.normalizeTo0dB = 0; options.max_num_lines = []; % Leave empty or omit to disable line rotation options.fft_length = []; + % --- NEW options --- + options.useWavelengthAxis (1,1) logical = false % plot x-axis in wavelength + options.lambda0_nm (1,1) double = 1310 % center wavelength [nm] end if isempty(options.fft_length) options.fft_length = 2^(nextpow2(length(obj.signal))-9); end - + if options.normalizeToNyquist == 0 - [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,obj.fs,"centered","power","mean"); - w = w.*1e-9; + [p_lin,f_Hz] = pwelch(obj.signal, hanning(options.fft_length), ... + options.fft_length/2, options.fft_length, ... + obj.fs, "centered", "power", "mean"); + f_GHz = f_Hz*1e-9; % keep frequency vector for frequency axis else - [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,"centered","power","mean"); + [p_lin,f_rad] = pwelch(obj.signal, hanning(options.fft_length), ... + options.fft_length/2, options.fft_length, ... + "centered", "power", "mean"); + % In normalized mode, pwelch returns rad/sample centered on 0. + % We'll keep f_rad for the x-axis in that mode. end if options.normalizeTo0dB - p_lin = p_lin./ max(p_lin); + p_lin = p_lin ./ max(p_lin); p_dbm = 10*log10(p_lin); % normalized to 0 dB ylab = "Normalized PSD"; else @@ -372,36 +382,69 @@ classdef Signal ylab = "Power (dB/Hz)"; end + % --- If requested, build wavelength axis from frequency offset --- + if options.useWavelengthAxis && options.normalizeToNyquist == 0 + c = physconst('LightSpeed'); % [m/s] + lambda0_m = options.lambda0_nm*1e-9; % center wavelength [m] + f_c = c / lambda0_m; % carrier frequency [Hz] + + % exact mapping + f_abs = f_c + f_Hz; % absolute frequency [Hz] + lambda_m = c ./ f_abs; % wavelength [m] + lambda_nm = lambda_m * 1e9; % wavelength [nm] + + % assign axis + x_vec = lambda_nm(:); + x_label = "Wavelength [nm]"; + + % Sort to ensure axis is ascending + [x_vec, sortIdx] = sort(x_vec, 'ascend'); + p_dbm = p_dbm(sortIdx, :); + else + % Frequency or normalized axes + if options.normalizeToNyquist == 0 + x_vec = f_GHz; + x_label = "Frequency in GHz"; + else + x_vec = f_rad; % normalized frequency in rad/sample + x_label = "Normalized Frequency"; + end + end + figure(options.fignum); ax = gca; hold on - if isempty(options.color) - hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1); - else - hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color); + for s = 1:min(size(p_dbm)) + if isempty(options.color) + plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1); + else + plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1, 'Color', options.color); + end end - % If user wants to limit the number of lines, check and remove old lines + % Limit number of lines if requested if ~isempty(options.max_num_lines) && options.max_num_lines > 0 allLines = findall(ax, 'Type', 'Line'); if length(allLines) > options.max_num_lines - % Sort lines by creation order. Usually, the oldest lines appear first in allLines. - % If needed, you can sort by UserData or other criteria. numToRemove = length(allLines) - options.max_num_lines; delete(allLines(1:numToRemove)); end end - if options.normalizeToNyquist == 0 - xlabel("Frequency in GHz"); - edgetick = 2^(nextpow2(obj.fs*1e-9)); - xticks(-edgetick:32:edgetick); - xlim([100*round(min(w)/100,1)-10, 100*round(max(w)/100,1)+10]) - xlim([-128 128]);%256GSa/s + % Axis labels and limits + xlabel(x_label); + + if options.useWavelengthAxis && options.normalizeToNyquist == 0 + xlim([min(x_vec) max(x_vec)]); else - xlabel("Normalized Frequency"); - xlim([-pi, pi]); + if options.normalizeToNyquist == 0 + % Keep your existing freq handling (you can fine-tune as needed) + % xlim([-128 128]); % example for 256 GSa/s if desired + xlim([min(x_vec) max(x_vec)]); + else + xlim([-pi, pi]); + end end ylabel(ylab); @@ -409,24 +452,23 @@ classdef Signal try ylim([max(min(floor(min(p_dbm))-3, ax.YLim(1)),-40), min(max(ceil(max(p_dbm))+3, ax.YLim(2)),10)]); catch - ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]); + ylim([floor(min(p_dbm,[],'all'))-3, ceil(max(p_dbm,[],'all'))+3]); end if options.normalizeTo0dB - % ylim([-40, 3]); - ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]); + ylim([floor(min(p_dbm,[],'all'))-3, ceil(max(p_dbm,[],'all'))+3]); else - ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]); + ylim([floor(min(p_dbm,[],'all'))-3, ceil(max(p_dbm,[],'all'))+3]); end yticks(-200:10:10); - grid on;% grid minor; + grid on; legend - % legend('Interpreter','none'); end + function move_it_spectrum(obj,options) arguments obj @@ -554,7 +596,7 @@ classdef Signal options.unit power_notation = power_notation.dBm end - pow = mean(abs(obj.signal).^2); + pow = sum(mean(abs(obj.signal).^2)); switch options.unit case power_notation.dBm @@ -669,9 +711,9 @@ classdef Signal options.fs_ref = 0; options.debug_plots = 0; end - - - + + + S = {}; inverted = -1; sequenceFound = 0; @@ -715,25 +757,25 @@ classdef Signal findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend') end - + shifts = lags(pkpos); sequenceStarts = shifts; shifts = shifts(shifts>=0); if numel(shifts) > 0 - + %Cut occurences of ref signal from signal (only positive shifts) if all(sign(co(pkpos))==-1) inverted = 1; end - + for c = shifts sig = obj.delay(-c,'mode','samples'); sig.signal = sig.signal(1:length(b));% .* -inverted; S{end+1,1} = sig; end - + % %return/keep the sinal with the highest correlation (only within positive shifts) % [~,idx]=max(pks); % obj.signal = S{idx}.signal; @@ -741,7 +783,7 @@ classdef Signal % swap = S{1}; % S{1} = S{idx}; % S{idx} = swap; - + for c = 1:numel(shifts) S{c}.logbook = []; end @@ -770,10 +812,6 @@ classdef Signal end - - - - %% function obj = filter(obj,a,b) @@ -934,8 +972,8 @@ classdef Signal nn=histcounts(data_ind_y(n,:),1:histpoints+1); hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-( end - - + + plot_data = 20*log10(hist_data); plot_data(plot_data==-Inf) = 0; @@ -970,7 +1008,7 @@ classdef Signal end xlabel('Time in ps') - + % add information @@ -1091,7 +1129,7 @@ classdef Signal x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12); xticklabels(x_tickstring); - grid off + grid off end diff --git a/Classes/02_etc/Amplifier.m b/Classes/02_etc/Amplifier.m index 31324e2..71db478 100644 --- a/Classes/02_etc/Amplifier.m +++ b/Classes/02_etc/Amplifier.m @@ -108,7 +108,7 @@ classdef Amplifier if obj.gain_mode == gain_mode.output_power %get linear gain for output power mode - pow_in = mean(abs(xin.^2)) ; % lin input power + pow_in = sum(mean(abs(xin.^2)),2) ; % lin input power pow_out = 10^(obj.amplification_db/10 - 3) ; % dBm to lin a_lin = sqrt(pow_out/pow_in) ; diff --git a/Classes/02_optical/DP_Fiber.m b/Classes/02_optical/DP_Fiber.m new file mode 100644 index 0000000..a5f3dce --- /dev/null +++ b/Classes/02_optical/DP_Fiber.m @@ -0,0 +1,240 @@ +classdef DP_Fiber + % Dual-Polarization fiber propagation (CNLSE / Manakov) — class version + % Runs full setup + propagation inside process_ (no external loop). + + properties (Access=public) + % ---- User options (public) ---- + L % [km] fiber length + dz % [m] step size target (kept for compatibility; adaptive dz uses SS_* below) + lambda % [nm] reference wavelength + rng % RNG seed + gamma % [1/W/m] nonlinear coefficient + + % Optional / advanced options (match legacy names where possible) + fa % [Hz] sampling frequency + X_alpha % [dB/100km] attenuation per 100 km (per-pol, used for X and Y) + X_beta % [1..4] dispersion coefficients vector for X (Y mirrors X) + D % [ps/(nm*km)] + Ds % [ps/(nm^2*km)] dispersion slope + Dpmd % [ps/sqrt(km)] PMD coefficient + beat_len % [m] beat length (for beta(1) if X_beta is zero) + corr_len % [m] correlation length (not directly used; kept for compatibility) + manakov % 1=Manakov (legacy behavior: manakov=eq-1) + SS_dphimax % [rad] max nonlinear phase per step (adaptive SSFM) + SS_dzmax % [m] max dz (adaptive SSFM) + SS_dzmin % [m] min dz (adaptive SSFM) + n_waveplates % number of PMD waveplates + + % ---- Internal state (persistent between calls) ---- + state % struct mirroring legacy 'state' + end + + methods (Access=public) + function obj = DP_Fiber(options) + % Constructor — copies fields from 'options' and sets defaults. + + arguments + options.L + options.dz + options.lambda + options.rng = 0 + options.gamma + + % optional but recommended + options.fa + + % legacy-compatible optional params + options.X_alpha = 4.605170185988092e-05 % dB/100km + options.X_beta = [0,0,-2.16826193914149e-26,3.56839456298263e-41] + options.D = 17 % ps/(nm*km) + options.Ds = 0.06 % ps/(nm^2*km) + options.Dpmd = 3 % ps/sqrt(km) + options.beat_len = 50 % m + options.corr_len = 50 % m (kept) + options.manakov = 0 % 1=Manakov + options.SS_dphimax = 5e-3 % rad + options.SS_dzmax = 2e4 % m + options.SS_dzmin = 100 % m + options.n_waveplates = 100 + end + + % Copy provided options into properties + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + % Initialize empty state; will be built on first process_ call + obj.state = struct(); + end + + function signalclass_out = process(obj, signalclass_in) + % Public entry point: takes a signal class with .signal (2xN) + % and writes back the propagated signal. + + signalclass_in.signal = obj.process_(signalclass_in.signal, signalclass_in.fs); + + % logbook (kept as in new framework skeleton) + lbdesc = 'DP_Fiber propagation (CNLSE_plain)'; + if ismethod(signalclass_in, 'logbookentry') + signalclass_in = signalclass_in.logbookentry(lbdesc); + end + + signalclass_out = signalclass_in; + end + + function signal_out = process_(obj, signal_in,fs) + % Core processing — builds legacy 'state' and calls CNLSE_plain + % data_in: [2 x N] complex, dual-pol envelope + arguments (Input) + obj + signal_in + fs + end + + % ---- Basic checks + if isempty(signal_in) || size(signal_in,2) ~= 2 + error('DP_Fiber:Input','Expected data_in of size [2 x N].'); + end + if isempty(fs) + error('DP_Fiber:Config','Sampling frequency options.fa is required.'); + end + + obj.fa = fs; + + % ---- RNG (legacy behavior) + R = RandStream("twister","Seed",obj.rng); + + % ---- (Re)build state if empty or size-dependent fields changed + need_rebuild = ~isfield(obj.state,'nt') || (obj.state.nt ~= size(signal_in,2)); + + if need_rebuild + % Constants + c0 = 299792458; % [m/s] + + % Legacy state mapping + st = struct(); + + % High-level + st.L = obj.L * 1000; % [m] legacy expects meters + st.polNames = {'X','Y'}; + st.dt = 1/obj.fa; + st.nt = max(size(signal_in)); + st.omega = 2*pi*[(0:st.nt/2-1),(-st.nt/2:-1)]/(st.dt*st.nt); + st.lambda = obj.lambda * 1e-9; % [m] + st.D = obj.D * 1e-6; % ps/(nm*km) -> s/(nm*m) + st.Dpmd = obj.Dpmd * 1e-6; % ps/sqrt(km) -> s/sqrt(km) + st.Ds = obj.Ds * 1e3; % ps/(nm^2*km) -> s/(nm^2*m) + st.beat_len = obj.beat_len; + st.SS_dzmax = obj.SS_dzmax; + st.SS_dzmin = obj.SS_dzmin; + st.SS_dphimax= obj.SS_dphimax; + + st.chi = 0; % legacy placeholders + st.psi = 0; + + st.manakov = obj.manakov; % 1 if eq==2 (Manakov), 0 if eq==1 (CNLSE) + st.wave_plates = obj.n_waveplates; + + % Alpha (same X/Y) + st.alpha.X = obj.X_alpha; + st.alpha_lin.X = st.alpha.X/10*log(10)/1000; + st.alpha.Y = st.alpha.X; + st.alpha_lin.Y = st.alpha_lin.X; + + % Beta (dispersion) for X (Y mirrors X) + if ~any(obj.X_beta) + % Populate from (beat_len, D, Ds, lambda) like legacy + if obj.beat_len ~= 0 + b1 = pi/obj.beat_len; + else + b1 = 0; + end + b2 = 0; % PMD freq term handled elsewhere in legacy + b3 = -(st.lambda.^2/(2*pi*c0))*st.D; + b4 = (st.lambda^2/(2*pi*c0))^2*st.Ds + (2/st.lambda)*(st.lambda.^2/(2*pi*c0)).^2*st.D; + + st.beta.X = [b1, b2, b3, b4]; % keep 4 terms, legacy had 0 for beta0; b2 unused + % Note: legacy stored [b0,b1,b2,b3]? Here we mirror their usage. + % We follow their X_beta layout length=4. + else + st.beta.X = obj.X_beta; + end + st.beta.Y = st.beta.X; + + % Gamma + st.gamma = obj.gamma; + + % PMD / birefringence (legacy waveplate model) + st.corr_length = st.L / st.wave_plates; + + % DGD formula (Agrawal 1.1.18) + st.dgd = st.Dpmd * sqrt(st.L/1000); % Dpmd in s/sqrt(km), L in m -> convert: sqrt(m/1000) + % For exact legacy match, they used: state.dgd = Dpmd * sqrt(L) wisth L in meters and Dpmd already scaled. + % Using their pattern: + st.dgd = obj.Dpmd*1e-6 * sqrt(st.L); % match old: state.Dpmd already 1e-6*ps/sqrt(km); they used sqrt(L) with L [m] + + brf_multiplier = 1; + if st.Dpmd == 0 + brf_multiplier = 0; + end + + % Waveplate random parameters + st.pauli_mats.s0 = eye(2); + st.pauli_mats.s2 = [0 1; 1 0]; + st.pauli_mats.s3i = [0 1; -1 0]; + + st.brf.theta = (R.rand(st.wave_plates,1)*pi - 0.5*pi) * brf_multiplier; + st.brf.epsilon = 0.5*asin(R.rand(st.wave_plates,1)*2-1) * brf_multiplier; + + st.brf.stokes = NaN(st.wave_plates,3); + st.Ttest = zeros(1, st.wave_plates); + st.brf.matR = cell(st.wave_plates,1); + + for n=1:st.wave_plates + matRth = cos(st.brf.theta(n)) * st.pauli_mats.s0 - sin(st.brf.theta(n)) * st.pauli_mats.s3i; + matRepsilon = complex(cos(st.brf.epsilon(n))*st.pauli_mats.s0, sin(st.brf.epsilon(n))*st.pauli_mats.s2); + matR = matRth * matRepsilon; + + st.brf.matR{n} = matR; + + u1 = matR(1,1); + u2 = matR(1,2); + + st.Ttest(n) = abs(u1).^2 + abs(u2).^2; + st.brf.stokes(n,:) = [abs(u1).^2 - abs(u2).^2 , (u1) * conj(u2) + conj(u1) .* u2 , 1i*(u1) * conj(u2) - conj(u1) .* u2]; + end + + % Frequency-dependent PMD phase term (legacy form) + st.brf.db0 = (R.rand(st.wave_plates,1)*2*pi - pi) * brf_multiplier; + st.brf.db1 = sqrt(3*pi/8)*(st.dgd/obj.fa)/st.wave_plates .* st.omega; + st.brf.simdgd = 0; + % cumsum used in legacy only for debug; keep compatibility variable: + ~cumsum(st.brf.db0); % no-op to mirror legacy path + + % Bookkeeping + st.missing_dz = 0; + st.n_plates_done = 0; + st.test_plates = []; + st.test_plate_numbers = []; + st.lin_z_test = 0; + st.propagated_length = 0; + + % Cache + obj.state = st; + end + + % ---- Call the exact same CNLSE_plain as in the old framework + x_in = signal_in(:,1).'; + y_in = signal_in(:,2).'; + + [x_out, y_out, obj.state] = CNLSE_plain(x_in, y_in, obj.state); + + obj.state.propagated_length = obj.state.propagated_length + obj.state.L; + + signal_out = [x_out; y_out].'; + + + end + end +end diff --git a/Classes/02_optical/EML.m b/Classes/02_optical/EML.m index 96765da..32a902c 100644 --- a/Classes/02_optical/EML.m +++ b/Classes/02_optical/EML.m @@ -70,7 +70,7 @@ classdef EML [signalclass_in.signal,obj] = obj.process_(signalclass_in.signal); % cast the inform. signal to electrical signal - signalclass_in = Opticalsignal(signalclass_in,"fs",obj.fsimu,"logbook",signalclass_in.logbook,"lambda",obj.lambda*1e-9,"nase",0); + signalclass_in = Opticalsignal(signalclass_in,"fs",obj.fsimu,"logbook",signalclass_in.logbook,"lambda",obj.lambda*1e-9,"nase",0,"polrot",0); % append to logbook lbdesc = [num2str(obj.lambda),' nm Laser with ',num2str(obj.power),' dBm P_out. Linew.=',num2str(obj.linewidth*1e-6),' MHz. Modulation mode: ',char(obj.mode) ]; diff --git a/Classes/02_optical/Optical_Demultiplex.m b/Classes/02_optical/Optical_Demultiplex.m new file mode 100644 index 0000000..fc2bc36 --- /dev/null +++ b/Classes/02_optical/Optical_Demultiplex.m @@ -0,0 +1,122 @@ +classdef Optical_Demultiplex < handle + % Dual-Polarization optical demultiplexer + % - Input: total-field signal + % - Output: single-channel dual-pol signal objects in cell array + % + % Notes: + % Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",Opt_sig_wdm_rx.fs/4,"fs_in",Opt_sig_wdm_rx.fs,"lambda_center",1310).process(Opt_sig_wdm_rx); + % Opt_sig_wdm_demux{1}.spectrum("fignum",1100,"displayname",'bla','normalizeTo0dB',0,'max_num_lines',4); + % Opt_sig_wdm_demux{2}.spectrum("fignum",1100,"displayname",'bla','normalizeTo0dB',0,'max_num_lines',4); + + properties (Access=public) + fs_in % [Hz] (optional; inferred from data_in.fs if omitted) + fs_out % [Hz] + lambda_center % [nm] center wavelength of the WDM grid + wavelengthplan % [nm] + attenuation = 0 % [dB] insertion loss + filtype = 1 % 1=Gaussian, 2=Rectangle, 3=No filter + B = 200e9 % [Hz] 3 dB bandwidth (Gaussian) or width (Rect) + mgauss = 3 % Gaussian order (multiple of 1/2) + + % Derived/utility + c = physconst('lightspeed') % [m/s] + end + + methods (Access=public) + function obj = Optical_Demultiplex(options) + arguments + options.fs_in = [] + options.fs_out + options.lambda_center + options.wavelengthplan + options.attenuation = 0 + options.filtype = 1 + options.B = 2.5e10 + options.mgauss = 3 + end + fn = fieldnames(options); + for n = 1:numel(fn) + try obj.(fn{n}) = options.(fn{n}); end + end + end + + function signalclasses_out = process(obj, signalclass_in) + + % ---- Infer wavelength: either given or from input total signal + if isempty(obj.wavelengthplan) + obj.wavelengthplan = signalclass_in.lambda; %meter + else + if all(500e-9 < obj.wavelengthplan) && all(obj.wavelengthplan < 1500e-9) %check if given in nm + obj.wavelengthplan = obj.wavelengthplan.*1e-9; + end + end + + % ---- Infer input sampling rates + if isempty(obj.fs_in) + assert(isprop(signalclass_in,'fs') && ~isempty(signalclass_in.fs), ... + 'Dual_Pol_Demultiplexer: data_in.fs missing and options.fs_in not provided.'); + obj.fs_in = signalclass_in.fs; + end + + % Runs demultiplexing in one go and appends a logbook entry. + [x_envelopes,y_envelopes] = obj.process_(signalclass_in.signal); + + for n = 1:min(size(x_envelopes)) + signalclasses_out{n} = signalclass_in; + signalclasses_out{n}.signal = [x_envelopes(:,n), y_envelopes(:,n)]; + signalclasses_out{n} = signalclasses_out{n}.resample("fs_in",obj.fs_in,"fs_out",obj.fs_out); + signalclasses_out{n}.lambda = obj.wavelengthplan(n); + lbdesc = ['Opt. Demux ', num2str( obj.wavelengthplan(n)),' nm']; + signalclasses_out{n} = signalclasses_out{n}.logbookentry(lbdesc); + end + + end + + function [x_envelopes,y_envelopes] = process_(obj, signal_in) + % Core demux: + % - frequency translate target channel to baseband + % - apply optical filter H + % - resample to fs_out + arguments (Input) + obj + signal_in + end + + w = obj.fs_out ./ obj.fs_in ; + blocklen_in = length(signal_in); + blocklen_out = w*blocklen_in; + + att = 1/10^(obj.attenuation/10); + faxis=linspace( -obj.fs_in/2 , obj.fs_in/2 , blocklen_in+1 ); + faxis=ifftshift(faxis(1:end-1)); + + switch obj.filtype + case 1 + H=exp(-(faxis/obj.B).^(2*obj.mgauss)*log(2)*2^(2*obj.mgauss-1)).'; + case 2 + %all zero filter + H=zeros(1,length(faxis)).'; + %set filter = 1 inside bandwidth -B/2 <-> B/2 + H(abs(faxis)<=obj.B/2)=1; + case 3 + H = 1; + end + + f_mid = obj.c/(obj.lambda_center*1e-9); % center frequency of WDM grid [Hz] + + f_channels = obj.c./(obj.wavelengthplan) ; + + N = numel(f_channels); + + df_T = f_mid - f_channels; + + pha = mod(-2*pi*(0:blocklen_in-1).'.*df_T/obj.fs_in, 2*pi); + lo = cos(pha)+1i*sin(pha); + + x_envelopes = ifft(fft(att.*signal_in(:,1).*lo).*H); + y_envelopes = ifft(fft(att.*signal_in(:,2).*lo).*H); + + end + end +end + diff --git a/Classes/02_optical/Optical_Multiplex.m b/Classes/02_optical/Optical_Multiplex.m new file mode 100644 index 0000000..1d722b0 --- /dev/null +++ b/Classes/02_optical/Optical_Multiplex.m @@ -0,0 +1,171 @@ +classdef Optical_Multiplex < handle + % Takes a cell array of signals + % returns a total field signal + % WDM spacing is given in wavelength plan OR via delta_F + + % The grid is stored in the output signal -> the demux will ideally + % look this up and use this as the demux frequencies... + + % signal_cell = {Opt_sig_1, Opt_sig_2}; + % Opt_sig_wdm = Optical_Multiplex("fs_in",Opt_sig.fs,"fs_out",4*Opt_sig.fs,... + % "lambda_center",1310,"random_key",0,"filtype",1,"B",200e9,"delta_f",400e9).process(signal_cell); + + properties(Access=public) + fs_in + fs_out + lambda_center + delta_f + random_key + attenuation + B + mgauss + filtype + + c = physconst('lightspeed') + f_center + f_T + lambda_T + df_T + end + + methods (Access=public) + function obj = Optical_Multiplex(options) + %NAME Construct an instance of this class + % Detailed explanation goes here + + arguments + options.fs_in + options.fs_out + options.lambda_center + options.B = 200e9 + options.mgauss = 3 + options.filtype = 2 + + options.delta_f = 0 + options.random_key + options.attenuation = 0; + end + + % + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end + + + end + + function signalclass_out = process(obj,signalclasses_in) + + % actual processing of the signal (steps 1. - 3.) + signalclass_out = obj.process_(signalclasses_in); + + % append to logbook + lbdesc = ['Opt. Mux. ']; + signalclass_out = signalclass_out.logbookentry(lbdesc); + + end + + function data_out = process_(obj,data_in) + %METHOD1 Summary of this method goes here + % Detailed explanation goes here + arguments(Input) + obj + data_in cell + end + + % assert(data_in{1}.fs == obj.fs_in,'Sampling rate'); + att = 1/10^(obj.attenuation/10); + N = numel(data_in); + w = obj.fs_out/data_in{1}.fs; + blocklen_in = length(data_in{1}); + blocklen_out = w*blocklen_in; + freqaxis = linspace(-obj.fs_out/2, obj.fs_out/2, blocklen_out+1); + obj.f_center = obj.c/(obj.lambda_center.*1e-9); + + if obj.random_key ~= 0 + res = freqaxis(2)-freqaxis(1); + R = RandStream("twister","Seed",obj.random_key); + laser_frequency_imperfection = res .* round(R.randn(N,1)*10); %in mutliples of the fft resolution, i.e. the distance between two freq. bins + else + laser_frequency_imperfection = zeros(blocklen_in,1); + end + + obj.f_T = []; + obj.df_T = []; + polrots = []; + for o = 1:N + + if obj.delta_f ~= 0 + % user defined a channel spacing in GHz. Build plan + % left and right from zero + obj.df_T(o) = (-length(data_in)/2-0.5+o) .* obj.delta_f; + obj.df_T(o) = obj.df_T(o)+ laser_frequency_imperfection(o); + + obj.f_T = [obj.f_T obj.f_center+obj.df_T(o)]; + else + %center frequencies of channels + obj.f_T = [obj.f_T obj.c/(data_in{o}.lambda)]; + obj.lambda_T = [obj.lambda_T data_in{o}.lambda]; + + %difference between mid frequency of MUX and channels + obj.df_T = [obj.df_T obj.f_center - obj.f_T(o)]; + end + + % adapt frequency shifts to match the FFT grid! Find nearest grid point + [glitch(o),pos] = min(abs( freqaxis-obj.df_T(o) )); + obj.df_T(o) = freqaxis(pos); + + polrots = [polrots, data_in{o}.polrot]; + end + + obj.lambda_T = obj.c ./ (obj.f_center-obj.df_T); + + obj.B = 200e9; %200GHz + faxis = linspace(-obj.fs_out/2,obj.fs_out/2, blocklen_out+1);%generates arow vector faxis of blocklen+1 points linearly spaced between and including -para.fs/2 and para.fs/2 + faxis = ifftshift(faxis(1:end-1)); + + switch obj.filtype + case 1 + H =exp(-(faxis/obj.B).^(2*obj.mgauss)*log(2)*2^(2*obj.mgauss-1)).'; + case 2 + H=zeros(length(faxis),1); + H(find(abs(faxis)<=obj.B/2))=1; + case 3 + H = 1; + end + + x_envelopes = NaN([blocklen_out N]); + y_envelopes = x_envelopes; + + for o = 1:N + + pha = mod(2*pi*(0:blocklen_out-1)*obj.df_T(o)/obj.fs_out,2*pi).'; + lo = cos(pha)+1i*sin(pha); + data_in_resampled = data_in{o}.resample("fs_out",obj.fs_out); + + res_env = ifft(fft(data_in_resampled.signal(:,1)).*H); + x_envelopes(:,o) = att.*res_env.*lo; + + res_env = ifft(fft(data_in_resampled.signal(:,2)).*H); + y_envelopes(:,o) = att.*res_env.*lo; + + end + + data_out = data_in_resampled; + data_out.signal = [sum(x_envelopes,2), sum(y_envelopes,2)]; + data_out.lambda = obj.lambda_T; + data_out.polrot = polrots; + end + + end + + methods (Access=private) + % Cant be seen from outside! So put all your functions here that can/ + % shall not be called from outside + + + end +end diff --git a/Classes/02_optical/Polarization_Controller.m b/Classes/02_optical/Polarization_Controller.m new file mode 100644 index 0000000..62792b9 --- /dev/null +++ b/Classes/02_optical/Polarization_Controller.m @@ -0,0 +1,92 @@ +classdef Polarization_Controller + %Input can be "normal" - output will be DP! + + properties(Access=public) + + mode + desired_angle + desired_power + + rotation_angle + rotation_matrix + + end + + methods (Access=public) + function obj = Polarization_Controller(options) + %NAME Construct an instance of this class + % Detailed explanation goes here + + arguments + options.mode polarization_control_mode = polarization_control_mode.rot_power + options.desired_angle + options.desired_power + + end + + % + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end + + % do more stuff + + end + + function signalclass_out = process(obj,signalclass_in) + + % actual processing of the signal (steps 1. - 3.) + [signalclass_in.signal,signalclass_in.polrot] = obj.process_(signalclass_in.signal,signalclass_in.polrot); + + % append to logbook + lbdesc = ['Logbookentry']; + signalclass_in = signalclass_in.logbookentry(lbdesc); + + % write to output + signalclass_out = signalclass_in; + + end + + function [data_out, polrot_out] = process_(obj,data_in,polrot_in) + % Rotate polarization of am opt signal + + arguments(Input) + obj + data_in double + polrot_in double + end + + + if obj.mode ~= polarization_control_mode.deactivate + + switch obj.mode + case polarization_control_mode.random + obj.rotation_angle = 2*pi*rand ; + + case polarization_control_mode.rot_angle + obj.rotation_angle = obj.desired_angle*pi/180 ; + + case polarization_control_mode.rot_power + obj.rotation_angle = -polrot_in + acos(sqrt(obj.desired_power/100)) ; + end + + obj.rotation_matrix = [cos(obj.rotation_angle) -sin(obj.rotation_angle) ; sin(obj.rotation_angle) cos(obj.rotation_angle)].' ; + + if min(size(data_in)) == 1 + data_in = reshape(data_in,[],1); + data_in = [data_in, zeros(length(data_in),1)]; + end + + data_out = data_in * obj.rotation_matrix; + + polrot_out = polrot_in + obj.rotation_angle ; + + end + + + end + end +end diff --git a/Classes/04_DSP/Coding/Duobinary.m b/Classes/04_DSP/Coding/Duobinary.m index 6f3460b..8590596 100644 --- a/Classes/04_DSP/Coding/Duobinary.m +++ b/Classes/04_DSP/Coding/Duobinary.m @@ -53,6 +53,7 @@ classdef Duobinary % bk(k+1) =mod(data(k)-bk(k),M); % end + % THIS WAS USED! for k = 2:numel(data) bk(k) = mod(data(k)-bk(k-1),M); end diff --git a/Classes/04_DSP/Sequence Detection/MLSE.m b/Classes/04_DSP/Sequence Detection/MLSE.m index 3883877..e51d61b 100644 --- a/Classes/04_DSP/Sequence Detection/MLSE.m +++ b/Classes/04_DSP/Sequence Detection/MLSE.m @@ -29,9 +29,6 @@ classdef MLSE < handle obj.(fn{n}) = options.(fn{n}); end end - - % do more stuff - end function [signalclass_hd,LLR,GMI] = process(obj,signalclass,ref_symbolclass) @@ -57,6 +54,26 @@ classdef MLSE < handle function [VITERBI_ESTIMATION_SYMBOLS,LLR_maxlogmap,GMI] = process_(obj,data_in,data_ref) + debug = 0; + + trellis_state_mode = 2; % General: States should match the target states of the prev. EQ (EQ's job was to reduce the error between signal and the target) + % 0 = use provided states (MUST provide the correct states); + % 1 = normalize to = 1 rms; + % 2 = use target symbols; + % 3 = use statistical levels + % 3 analyzes avg of rx signal levels - can help with nonlinear impairments + + trellis_exclusion = 0; % PAM-6 only (only if data is NOT precoded!) + + scale_mode = 2; % scale_mode: + % 0 = no scaling, + % 1 = RMS→scale MODEL, + % 2 = MMSE/time-corr→scale MODEL, + % 3 = RMS→scale DATA, + % 4 = MMSE/time-corr→scale DATA + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%% PREPARATIONS %%%%%%%% % remove unnecessary zeros at start of impulse response to keep % number of trellis states minimal @@ -69,22 +86,42 @@ classdef MLSE < handle obj.DIR = [0 obj.DIR]; end - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%%%%% PREPARATIONS %%%%%%%% - - %%%% Separate the equalized signal into the respective levels based on the actually transmitted level - constellation = unique(data_ref); - decisionLevels = (constellation(1:end-1) + constellation(2:end)) / 2; - N = length(data_in); - - tx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(data_ref); - % impulse respnse to remove from signal obj.DIR = flip(obj.DIR); %i.e. -0.2676 -0.0478 1.0000 - % Normalize the Trellis states to =1 RMS - obj.trellis_states = obj.trellis_states ./ rms(obj.trellis_states); + tx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(data_ref); + + % Trellis States + if trellis_state_mode == 1 % Normalize the Trellis states to =1 RMS + + obj.trellis_states = obj.trellis_states ./ rms(obj.trellis_states); + + elseif trellis_state_mode == 2 %simply use the states from the ref signal (should be a robust option) + + obj.trellis_states = reshape(unique(data_ref),size(obj.trellis_states)); + + elseif trellis_state_mode == 3 %use_statistical_levels + + %%%% Separate the equalized signal into the respective levels based on the actually transmitted level + constellation = unique(data_ref); + + % find actual levels from rx signal + symbols_for_lvl = NaN(numel(constellation),length(data_ref)); + for l = 1:numel(constellation) + level_amplitude = constellation(l); + symbols_for_lvl(l,data_ref==level_amplitude) = data_in(data_ref==level_amplitude); + end + + %replace the trellis states + avg_levels = mean(symbols_for_lvl,2,'omitnan'); + obj.trellis_states = sort(avg_levels)'; + + %also replace the whole ref signal (PAM-M) levels + [~, idx] = ismember(data_ref, unique(data_ref)); + data_ref = avg_levels(idx); + + end + % seems to be the only way to use combvec for a flexible amount % of vectors. 'combs' contains all trellis states @@ -97,74 +134,84 @@ classdef MLSE < handle states = sum(combs,2); nStates = length(last_sym); - % Calculate all possible input symbols for the desired impulse - % response. Row number is the index of the previous state, - % column number is the index of the next state - % noise free received == branch metrics - noise_free_received = zeros(nStates,nStates); - count_row = 1; - count_col = 1; - for l1 = 1:nStates - for l2 = 1:nStates - if sum(combs(l2,2:end) == combs(l1,1:end-1)) == size(combs,2)-1 - noise_free_received(count_row,count_col) = sum(combs(l2,:).*obj.DIR(end:-1:2)) + last_sym(l1)*obj.DIR(1); - else - noise_free_received(count_row,count_col) = inf; + % % Calculate all possible input symbols for the desired impulse + % % response. Row number is the index of the previous state, + % % column number is the index of the next state + % % noise free received == branch metrics + % assumes: last_sym = combs(:,end); % already defined earlier + levels = sort(unique(obj.trellis_states(:)).'); + edges = [levels(1) levels(end)]; % edge levels (0 and 5 in PAM6) + + noise_free_received = inf(nStates,nStates); % rows: to, cols: from + edge_edge_mask = false(nStates,nStates); % rows: to, cols: from + + for from = 1:nStates + for to = 1:nStates + % valid transition if shift-register overlap holds + if all(combs(to,2:end) == combs(from,1:end-1)) + % noiseless sample for the 'to' state reached from 'from' + noise_free_received(to,from) = ... + dot(combs(to,:), obj.DIR(end:-1:2)) + last_sym(from)*obj.DIR(1); + + % mark edge→edge candidate (to be excluded only on even→odd steps) + edge_edge_mask(to,from) = ... + (last_sym(from)==edges(1) || last_sym(from)==edges(2)) && ... + (last_sym(to) ==edges(1) || last_sym(to) ==edges(2)); end - count_row = count_row + 1; end - count_col = count_col + 1; - count_row = 1; end - % Was used earlier from Tom Wettlin etc. Sufficient for Viterbi - % but the BCJR is sensitive to the scaling, so I use another - % apporach - % % first: RMS normalization of input data (rms==1) - % data_in = data_in ./ rms(data_in); - % - % % then, match amplitude levels of input signal to those of the calculated ideal symbols - % % i.e. match the rms values of data_in to noise_free_received (rms=1.xx) - % if isreal(data_in) - % if obj.M == round(obj.M) - % data_in = data_in * rms(noise_free_received(noise_free_received ~= inf),'all','omitnan'); - % end - % end - - % quick and dirty alignment with xcorr - y = data_in; - x = data_ref; h = flip(obj.DIR(:)).'; - y_ideal = conv(x, h, "same"); - [c,lags] = xcorr(y, y_ideal, 64); % small max lag is enough - [~,ix] = max(abs(c)); - lag = lags(ix); - y_ideal = circshift(y_ideal, lag); + data_in = data_in(:); + y_ideal = conv(data_ref(:), h, "same"); - % center, skip transients, rm mean - skip = max(100, numel(obj.DIR)+50); - y_prep = y(1+skip:end-skip) - mean(y(1+skip:end-skip)); - y_ideal_prep = y_ideal(1+skip:end-skip) - mean(y_ideal(1+skip:end-skip)); - - % one-tap LS/MMSE gain and optional bias - g = (y_ideal_prep' * y_prep) / (y_ideal_prep' * y_ideal_prep); % complex allowed - b = mean(data_in) - g*mean(y_ideal); % intercept (if you keep means) - - dp = dot(y_ideal_prep, y_prep - g*y_ideal_prep); %shall be close to zero - - sigma2 = mean(abs(y_prep - g*y_ideal_prep).^2); - inv2s2 = 1/(2*sigma2); - - debug = 0; - if debug - figure(100); hold on - plot(1:length(y),y,'DisplayName','Y: Received Signal','LineStyle','none','Marker','.','MarkerSize',1); - plot(1:length(y_ideal),y_ideal,'DisplayName','Y ideal: x * DIR','LineStyle','none','Marker','.','MarkerSize',1); - yline(noise_free_received(:),'DisplayName','All Transition States','HandleVisibility','off','Color','red'); - yline(g*noise_free_received(:)+ b,'DisplayName','Scaled Transition States','HandleVisibility','off'); + switch scale_mode + case 0 + g = 1; b = 0; + case 1 % RMS: scale model to data + g = rms(data_in)/rms(y_ideal); b = mean(data_in) - g*mean(y_ideal); + case 2 % MMSE/time-corr: scale states to data + [c,lags] = xcorr(data_in(:), y_ideal, 64); + [~,ix] = max(abs(c)); + lag = lags(ix); + y_ideal = circshift(y_ideal, lag); + mu_y = mean(data_in(:)); + mu_i = mean(y_ideal); + y_c = data_in(:)-mu_y; + yi_c = y_ideal-mu_i; + g = (yi_c'*y_c)/(yi_c'*yi_c); + b = mu_y - g*mu_i; + case 3 % RMS flipped: scale data to model + gd = rms(y_ideal)/rms(data_in); bd = mean(y_ideal) - gd*mean(data_in); + data_in = gd*data_in + bd; + g = 1; b = 0; + case 4 % MMSE/time-corr flipped: scale data to states + [c,lags] = xcorr(data_in(:), y_ideal(:), 64); + [~,ix] = max(abs(c)); + lag = lags(ix); + y_ideal = circshift(y_ideal(:), lag); + mu_y = mean(data_in(:)); + mu_i = mean(y_ideal); + y_c = data_in(:) - mu_y; % data_in centered + yi_c = y_ideal - mu_i; % ideal centered + g = (y_c' * yi_c) / (y_c' * y_c); + b = mu_i - g * mu_y; + data_in = g * data_in(:) + b; + g = 1; b = 0; end - noise_free_received = (g*noise_free_received + b); + % apply (g,b) to model and compute common sigma + noise_free_received = g*noise_free_received + b; + last_sym = g*last_sym + b; + sigma2 = mean(abs(data_in - (g*y_ideal + b)).^2); + inv2s2 = 1/(2*sigma2); + + if debug + figure(100); clf; hold on + showLevelScatter(data_in, data_ref, "fignum", 100); + yline(noise_free_received(:), 'DisplayName','Transition States','Color','red','HandleVisibility','off'); + yline(obj.trellis_states(:), 'DisplayName','Transition States','Color','green','LineWidth',2,'HandleVisibility','off') + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% FORWARD PASS (VITERBI -Alpha's) %%%%% @@ -185,8 +232,14 @@ classdef MLSE < handle for n = 2:length(data_in) bm = -(data_in(n) - noise_free_received).^2 * inv2s2; + + % exclude edge→edge transitions only for even->odd steps && PAM-6 + if mod(n,2) == 0 && obj.M == 6 && trellis_exclusion + bm(edge_edge_mask) = -Inf; + end + pm = pm + bm; - [alpha(:,n),pm_survivor_fw_idx(:,n)] = max(pm,[],2); % choose lowest path metric as new state + [alpha(:,n),pm_survivor_fw_idx(:,n)] = max(pm,[],2); % choose lowest path metric as new state (get min distance for all state transitions towards a new state) pm = repmat(alpha(:,n).',nStates,1); % update pm (chosen state to 2nd dimension -> FROM state) bm_fw(:,:,n) = bm; @@ -202,6 +255,18 @@ classdef MLSE < handle viterbi_path(n-1) = pm_survivor_fw_idx(viterbi_path(n),n); end + + if debug + alpha_ = alpha - min(alpha) + eps; + figure();hold on; + n = 10; + scatter(1:n,obj.trellis_states(repmat([1:numel(obj.trellis_states)]',1,n)),abs(alpha_(:,end-n+1:end)),'Marker','o','LineWidth',1); + scatter(1:n,obj.trellis_states(viterbi_path(end-n+1:end)),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','green'); + % scatter(1:n,data_ref(end-n+1:end),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','red'); + yticks(obj.trellis_states); + ylim([min(obj.trellis_states)-1 max(obj.trellis_states)+1]); + end + VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = first_sym(viterbi_path); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -219,6 +284,12 @@ classdef MLSE < handle for h = length(data_in)-1:-1:1 bm = -(data_in(h+1) - noise_free_received).^2 * inv2s2; + + % exclude edge→edge transitions only for even->odd steps && PAM-6 + if mod(h+1, 2) == 0 && obj.M == 6 && trellis_exclusion + bm(edge_edge_mask) = -Inf; + end + pm = pm + bm.'; [beta(:,h),pm_survivor_bw_idx(:,h)] = max(pm,[],2); % choose lowest path metric as new state pm = repmat(beta(:,h).',nStates,1); % update pm (chosen state to 2nd dimension -> FROM state) @@ -255,6 +326,7 @@ classdef MLSE < handle %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% FORWARD PASS PAM2,4,8 %%%%% + % These are interchangeable... second is chatgpt: nml_LLP = LLP - max(LLP); %subtract highest value for better numerical stability, LLP's are not always close to zero expLLP = exp(nml_LLP); state_prob = expLLP ./ sum(expLLP); % sums to one (or numerically close to one) @@ -266,9 +338,6 @@ classdef MLSE < handle state_prob = exp(logPstate); % exact, sums to 1 - - - if obj.M == 6 num_bits = 5; @@ -357,17 +426,19 @@ classdef MLSE < handle for bit_idx = 1:num_bits % Find indices where bit is 0 and where it is 1 - idx_sym_1 = bit_mapping(:,bit_idx) == 0; + idx_bit_0 = bit_mapping(:,bit_idx) == 0; idx_bit_1 = bit_mapping(:,bit_idx) == 1; - % Sum over log-probabilities (Max-Log approximation: using max instead of sum) - LLR_maxlogmap(:,bit_idx) = max(LLP(idx_bit_1,:), [], 1) - max(LLP(idx_sym_1,:), [], 1); + % Sum over log-probabilities + % Max-Log approximation: using max instead of sum) + LLR_maxlogmap(:,bit_idx) = max(LLP(idx_bit_1,:), [], 1) - max(LLP(idx_bit_0,:), [], 1); % Sum probabilities over states for which the bit is 1 and 0, respectively. - P0 = sum(state_prob(idx_sym_1, :),1); + P0 = sum(state_prob(idx_bit_0, :),1); P1 = sum(state_prob(idx_bit_1, :),1); LLR_exact(:,bit_idx) = log(P1./P0); % N x num_bits + end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -377,12 +448,12 @@ classdef MLSE < handle LLR_exact = LLR_exact; for k = 1:num_bits - idx_bit_1 = (tx_bits(:,k) == 0); %wo sind die 1en - idx_sym_1 = (tx_bits(:,k) == 1); %wo sind die 0en + idx_bit_0 = (tx_bits(:,k) == 0); %wo sind die 1en + idx_bit_1 = (tx_bits(:,k) == 1); %wo sind die 0en %LLR's for all actually transmitted ones or zeros - llr0 = LLR_exact(idx_bit_1,k); - llr1 = LLR_exact(idx_sym_1,k); + llr0 = LLR_exact(idx_bit_0,k); + llr1 = LLR_exact(idx_bit_1,k); % Calculate mutual information for bit position k I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1 @@ -396,7 +467,6 @@ classdef MLSE < handle VITERBI_ESTIMATION_SYMBOLS = VITERBI_ESTIMATION_SYMBOLS./rms(VITERBI_ESTIMATION_SYMBOLS); - if debug %%% DEBUG PLOT LIKELIHOOD RATIOS %%% figure(115);clf @@ -415,6 +485,8 @@ classdef MLSE < handle legend end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% CHECK BER's %%%%% @@ -430,6 +502,13 @@ classdef MLSE < handle fprintf('Viterbi BER: %.2e \n',ber_viterbi); fprintf('Viterbi Errors = %d\n', numErr); + pairs = reshape(VITERBI_ESTIMATION_SYMBOLS,2,[]).'; + levels = sort(unique(VITERBI_ESTIMATION_SYMBOLS)); + isedge = ismember(pairs, [levels(1) levels(end)]); + isforbidden = sum(isedge,2)==2; + fprintf('Found %d forbidden transitions (even→odd edges).\n', nnz(isforbidden)); + + % Convert LLR values to a hard-decision bit stream bit_stream = LLR_maxlogmap > 0; %ratio separates lower or higher than =0 -> simply decode for the negative values bit_stream = reshape(bit_stream',[],1); @@ -458,7 +537,7 @@ classdef MLSE < handle rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(FW_EST'); rx_bits = reshape(rx_bits',[],1); [~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('FW BER: %.2e \n',ber_fw); + fprintf('FW BER : %.2e \n',ber_fw); disp('Stop DEBUG MLSE:') fprintf('\n') end @@ -471,7 +550,5 @@ classdef MLSE < handle s = amax + log(sum(exp(a - amax), dim)); end - end - end diff --git a/Classes/04_DSP/TransmissionPerformance.m b/Classes/04_DSP/TransmissionPerformance.m index 1550e40..d87d691 100644 --- a/Classes/04_DSP/TransmissionPerformance.m +++ b/Classes/04_DSP/TransmissionPerformance.m @@ -42,6 +42,7 @@ classdef TransmissionPerformance 0.8733, 0.8790, 0.8848, 0.8905, 0.8962, ... 0.9019, 0.9077, 0.9134, 0.9191, 0.9248, ... 0.9306, 0.9363, 0.9420, 0.9477, 0.9535]; + NGMITHRESHOLDS_SDHD = [0.8116, 0.8167, 0.8241, 0.8317, 0.8401, ... 0.8459, 0.8512, 0.8574, 0.8685, 0.8746, ... 0.8829, 0.8892, 0.8958, 0.9022, 0.9090,... @@ -60,7 +61,7 @@ classdef TransmissionPerformance %% LUT for KP4-FEC and Inner Code https://grouper.ieee.org/groups/802/3/dj/public/23_03/patra_3dj_01b_2303.pdf - CODE_RATE_KP4_AND_INNER = [0.885799]; + CODE_RATE_KP4_AND_INNER = [0.885799]; %https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9979198 -> 199.8 / 224 = 0.891 als code rate BERTHRESHOLDS_KP4_AND_INNER = 4.85e-3; % https://www.ieee802.org/3/bs/public/14_11/parthasarathy_3bs_01a_1114.pdf @@ -153,11 +154,21 @@ classdef TransmissionPerformance netrates.SDHD.Threshold = NaN(1, numMeasurements); end if ~isempty(ber) + netrates.STAIR.GrossRate = NaN(1, numMeasurements); + netrates.STAIR.NetRate = NaN(1, numMeasurements); + netrates.STAIR.CodeRate = NaN(1, numMeasurements); + netrates.STAIR.Threshold = NaN(1, numMeasurements); + netrates.HD.GrossRate = NaN(1, numMeasurements); netrates.HD.NetRate = NaN(1, numMeasurements); netrates.HD.CodeRate = NaN(1, numMeasurements); netrates.HD.Threshold = NaN(1, numMeasurements); + netrates.KP4.GrossRate = NaN(1, numMeasurements); + netrates.KP4.NetRate = NaN(1, numMeasurements); + netrates.KP4.CodeRate = NaN(1, numMeasurements); + netrates.KP4.Threshold = NaN(1, numMeasurements); + netrates.KP4_hamming.GrossRate = NaN(1, numMeasurements); netrates.KP4_hamming.NetRate = NaN(1, numMeasurements); netrates.KP4_hamming.CodeRate = NaN(1, numMeasurements); @@ -200,11 +211,47 @@ classdef TransmissionPerformance end if ~isempty(idxBER) codeRate = obj.CODE_RATES_HD(idxBER); + netrates.STAIR.NetRate(i) = grossRate(i) * codeRate; + netrates.STAIR.GrossRate(i) = grossRate(i) ; + netrates.STAIR.CodeRate(i) = codeRate; + netrates.STAIR.Threshold(i) = obj.BERTHRESHOLDS_HD(idxBER); + end + + % HD FEC 3,8e-3 + idxBER = []; + for j = length(obj.BERTHRESHOLDS_HDFEC):-1:1 + if ber(i) <= obj.BERTHRESHOLDS_HDFEC(j) + idxBER = j; + break; + end + end + if ~isempty(idxBER) + codeRate = obj.CODE_RATE_HDFEC(idxBER); netrates.HD.NetRate(i) = grossRate(i) * codeRate; netrates.HD.GrossRate(i) = grossRate(i) ; netrates.HD.CodeRate(i) = codeRate; - netrates.HD.Threshold(i) = obj.BERTHRESHOLDS_HD(idxBER); + netrates.HD.Threshold(i) = obj.CODE_RATE_HDFEC(idxBER); end + + + % KP4 + idxBER = []; + for j = length(obj.BERTHRESHOLDS_KP4):-1:1 + if ber(i) <= obj.BERTHRESHOLDS_KP4(j) + idxBER = j; + break; + end + end + if ~isempty(idxBER) + codeRate = obj.CODE_RATE_KP4(idxBER); + netrates.KP4.NetRate(i) = grossRate(i) * codeRate; + netrates.KP4.GrossRate(i) = grossRate(i) ; + netrates.KP4.CodeRate(i) = codeRate; + netrates.KP4.Threshold(i) = obj.CODE_RATE_KP4(idxBER); + end + + + % KP4 + inner Hamming idxBER = []; for j = length(obj.BERTHRESHOLDS_KP4_AND_INNER):-1:1 if ber(i) <= obj.BERTHRESHOLDS_KP4_AND_INNER(j) @@ -220,6 +267,7 @@ classdef TransmissionPerformance netrates.KP4_hamming.Threshold(i) = obj.BERTHRESHOLDS_KP4_AND_INNER(idxBER); end + % O FEC idxBER = []; for j = length(obj.BERTHRESHOLDS_O_FEC):-1:1 if ber(i) <= obj.BERTHRESHOLDS_O_FEC(j) diff --git a/Classes/DataBaseHandler/DBHandler.m b/Classes/DataBaseHandler/DBHandler.m index d21d3be..f1ca74e 100644 --- a/Classes/DataBaseHandler/DBHandler.m +++ b/Classes/DataBaseHandler/DBHandler.m @@ -63,7 +63,9 @@ classdef DBHandler < handle obj.refresh(); else + error('DB seems to be corrupt') + end end @@ -183,6 +185,7 @@ classdef DBHandler < handle fprintf('Raw Rx Paths: Found %d duplictaes of %s \n',duplictae_raw.occurrences(i),duplictae_raw.rx_raw_path(i)); end end + healthyDB = true; end @@ -713,68 +716,167 @@ classdef DBHandler < handle selectedFields end - % Step 1: Handle selectedFields conversion + % -------- Step 1: Normalize selectedFields to {'Table.field', ...} -------- if isempty(selectedFields) || (ischar(selectedFields) && strcmpi(selectedFields, 'all')) - selectedFields = obj.getTableFieldNames('Runs'); % Default to Runs table + selectedFields = obj.getTableFieldNames('Runs'); % default elseif isstruct(selectedFields) - % Convert struct to cell array of 'Table.field' format newFields = {}; tableNames = fieldnames(selectedFields); for t = 1:numel(tableNames) tableStruct = selectedFields.(tableNames{t}); - fieldNames = fieldnames(tableStruct); - for f = 1:numel(fieldNames) - if isequal(tableStruct.(fieldNames{f}), 1) - newFields{end+1} = sprintf('%s.%s', tableNames{t}, fieldNames{f}); + fns = fieldnames(tableStruct); + for f = 1:numel(fns) + if isequal(tableStruct.(fns{f}), 1) + newFields{end+1} = sprintf('%s.%s', tableNames{t}, fns{f}); %#ok end end end selectedFields = newFields; end - % Step 2: Generate COALESCE string for SELECT clause + % Parse the table names actually referenced by the SELECT + reqTables = unique(cellfun(@(s) extractBefore(s, '.'), selectedFields, ... + 'UniformOutput', false)); + + % -------- Step 2: Build SELECT with COALESCE wrapper as you already do ---- selectClause = obj.generateCoalesceString(selectedFields); - % Step 3: Generate FROM clause with joins - mainTable = 'Runs'; - fromClause = ['FROM ', mainTable, ' ']; + % -------- Step 3: FROM and minimal JOIN plan ------------------------------ + % Decide main table: prefer the first explicitly referenced table, else 'Runs' + if ~isempty(reqTables) + mainTable = reqTables{1}; + else + mainTable = 'Runs'; + end - % Prepare join buffers - normalJoins = ''; - equalizerJoin = ''; - - tableNamesAll = fieldnames(obj.tables); - for t = 1:numel(tableNamesAll) - tableName = tableNamesAll{t}; - - if strcmpi(tableName, mainTable) || strcmpi(tableName, 'sqlite_sequence') - continue; - end - - if isfield(obj.tables.(tableName), 'run_id') - % Direct join to Runs - normalJoins = [normalJoins, 'LEFT JOIN ', tableName, ... - ' ON ', mainTable, '.run_id = ', tableName, '.run_id ']; - elseif isfield(obj.tables.(tableName), 'eq_id') - % Equalizer depends on Results, collect this join separately - equalizerJoin = ['LEFT JOIN ', tableName, ... - ' ON Results.eq_id = ', tableName, '.eq_id ']; + % If WHERE references a table not in reqTables (e.g., Runs.*), make sure it’s present. + whereClause = ''; + if ~isempty(filterParams) + whereClause = obj.generateWhereClause(filterParams); + % Heuristic: add 'Runs' if WHERE clause mentions 'Runs.' + if contains(whereClause, 'Runs.') + reqTables = unique([reqTables; {'Runs'}]); %#ok end end - % Combine joins: normal joins first, Equalizer last - fromClause = [fromClause, normalJoins, equalizerJoin]; + % Ensure main table is included + if ~ismember(mainTable, reqTables) + reqTables = unique([mainTable; reqTables]); %#ok + end - % Step 4: Generate WHERE clause if filter parameters exist - if ~isempty(filterParams) - whereClause = obj.generateWhereClause(filterParams); - if ~isempty(whereClause) - query = [selectClause, ' ', fromClause, 'WHERE ', whereClause]; - else - query = [selectClause, ' ', fromClause]; + % We’ll build joins only for the required tables (minus the main) + otherTables = setdiff(reqTables, {mainTable}); + + % Keep track of what’s already in the FROM graph (start with main) + present = string(mainTable); + joins = strings(0,1); + + % Helper lambdas + hasField = @(tbl, fld) isfield(obj.tables.(char(tbl)), char(fld)); + canJoinBy = @(left, right, key) hasField(left, key) && hasField(right, key); + + + % A small helper that adds a LEFT JOIN if the right table isn't present yet + function addJoinByKey(rightTbl, key) + if any(present == string(rightTbl)) + return; % already joined end + % Prefer to join against an already-present table that has the key + anchor = ''; + for k = 1:numel(present) + if canJoinBy(char(present(k)), rightTbl, key) + anchor = char(present(k)); + break; + end + end + if isempty(anchor) + % No anchor in current graph; if the right table is 'Equalizer' and key is eq_id, + % try to ensure a bridge table with eq_id exists (Results or a dashboard view). + if strcmpi(rightTbl,'Equalizer') && strcmpi(key,'eq_id') + % Bring in one eq_id-capable table if it is requested + bridgeOrder = {'Results','dashboard_old','dashboard_new','dashboard_ungrouped'}; + for b = 1:numel(bridgeOrder) + br = bridgeOrder{b}; + if ismember(br, reqTables) && ~any(present == string(br)) && hasField(obj.tables.(br),'eq_id') + % Attach bridge by run_id if possible, otherwise leave for eq_id + if any(present == "Runs") && hasField(obj.tables.(br),'run_id') && hasField(obj.tables.('Runs'),'run_id') + joins(end+1,1) = "LEFT JOIN " + br + " ON Runs.run_id = " + br + ".run_id"; + present(end+1,1) = string(br); + anchor = br; % we can now anchor Equalizer on eq_id to this + break; + else + % Fallback: anchor to main if it shares eq_id + for k = 1:numel(present) + pk = char(present(k)); + if canJoinBy(pk, br, 'eq_id') + joins(end+1,1) = "LEFT JOIN " + br + " ON " + pk + ".eq_id = " + br + ".eq_id"; + present(end+1,1) = string(br); + anchor = br; + break; + end + end + if ~isempty(anchor), break; end + end + end + end + end + end + % Re-scan for an anchor (maybe the bridge helped) + if isempty(anchor) + for k = 1:numel(present) + if canJoinBy(char(present(k)), rightTbl, key) + anchor = char(present(k)); + break; + end + end + end + if isempty(anchor) + % As a final fallback, if the main table is Runs and right has run_id, join by run_id + if strcmpi(mainTable,'Runs') && hasField(obj.tables.(rightTbl),'run_id') && hasField(obj.tables.('Runs'),'run_id') + anchor = 'Runs'; + key = 'run_id'; + end + end + if isempty(anchor) + % Could not find a path; skip join silently (or throw if you prefer strict) + return; + end + joins(end+1,1) = "LEFT JOIN " + rightTbl + " ON " + anchor + "." + key + " = " + rightTbl + "." + key; + present(end+1,1) = string(rightTbl); + end + + % First pass: if WHERE uses Runs.* and mainTable isn’t Runs, ensure Runs is in the graph + if contains(string(whereClause), "Runs.") && ~any(present == "Runs") + % Try to join Runs to whatever has run_id (mainTable ideally) + if hasField(mainTable, 'run_id') && hasField('Runs', 'run_id') + joins(end+1,1) = "LEFT JOIN Runs ON " + string(mainTable) + ".run_id = Runs.run_id"; + present(end+1,1) = "Runs"; + end + end + + % Join the required tables with minimal edges + for i = 1:numel(otherTables) + tbl = otherTables{i}; + % Prefer run_id join if possible, else eq_id, else skip + if any(present == "Runs") && hasField(tbl,'run_id') + addJoinByKey(tbl, 'run_id'); + elseif hasField(tbl,'run_id') && hasField(mainTable,'run_id') + addJoinByKey(tbl, 'run_id'); + elseif hasField(tbl,'eq_id') + addJoinByKey(tbl, 'eq_id'); + else + % no obvious key; skip + end + end + + % Build the FROM clause + fromClause = "FROM " + string(mainTable) + " " + strjoin(joins, " "); + + % -------- Step 4: WHERE (unchanged logic) ------------------------------ + if ~isempty(whereClause) + query = char(strjoin([selectClause, fromClause, "WHERE " + string(whereClause)], " ")); else - query = [selectClause, ' ', fromClause]; + query = char(strjoin([selectClause, fromClause], " ")); end end diff --git a/Classes/Warehouse_class/classes/DataStorage.m b/Classes/Warehouse_class/classes/DataStorage.m index f057c54..5fd288d 100644 --- a/Classes/Warehouse_class/classes/DataStorage.m +++ b/Classes/Warehouse_class/classes/DataStorage.m @@ -51,8 +51,9 @@ classdef DataStorage < handle function save(obj,path) try save(path,"obj"); - catch - + catch e + disp(e.message) + disp('Provide save path') end end diff --git a/Classes/Warehouse_class/classes/minimal_example.m b/Classes/Warehouse_class/classes/minimal_example.m new file mode 100644 index 0000000..b106929 --- /dev/null +++ b/Classes/Warehouse_class/classes/minimal_example.m @@ -0,0 +1,40 @@ + +params = struct; +params.sir = [20:2:36]; +params.laser_linewidth = [1e5 1e6 10e6]; +params.pn_key = [1:3]; +params.vp = [0.25,0.5,0.75,1]; +params.vb = [1:0.1:1.8]; +params.rop = -5:0; + +wh = DataStorage(params); +wh.addStorage("ber"); +wh.addStorage("rop_save"); +wh.addStorage("cspr"); +wh.addStorage("mod_out_pow"); + + +cnt = 1; +for sir = params.sir + for lw = params.laser_linewidth + for pnk = params.pn_key + for vp = params.vp + for vb = params.vb + for rop = params.rop + + current_ber = randn; + wh.addValueToStorage(current_ber,'ber',sir,lw,pnk,vp,vb,rop); + wh.addValueToStorage(-10,'rop_save',sir,lw,pnk,vp,vb,rop); + wh.addValueToStorage(10,'cspr',sir,lw,pnk,vp,vb,rop); + wh.addValueToStorage(-6,'mod_out_pow',sir,lw,pnk,vp,vb,rop); + + end + end + end + end + end +end + +wh.getStoValue('ber',20,1e5,2,0.25,1.1,-3) + +wh.showInfo \ No newline at end of file diff --git a/Datatypes/polarization_control_mode.m b/Datatypes/polarization_control_mode.m new file mode 100644 index 0000000..c968fcd --- /dev/null +++ b/Datatypes/polarization_control_mode.m @@ -0,0 +1,10 @@ +classdef polarization_control_mode < int32 + + enumeration + random (1) + rot_angle (2) + rot_power (3) + deactivate (4) + end + +end \ No newline at end of file diff --git a/Functions/EQ_structures/dsp_runid.m b/Functions/EQ_structures/dsp_runid.m index a698314..09e5e27 100644 --- a/Functions/EQ_structures/dsp_runid.m +++ b/Functions/EQ_structures/dsp_runid.m @@ -96,11 +96,11 @@ try adaption= 1; use_dd_mode = 1; - use_ffe = 0; - use_dfe = 0; + use_ffe = 1; + use_dfe = 1; use_vnle_mlse = 1; use_dbtgt = 1; - use_dbenc = 1; + use_dbenc = 0; addProcessingResultToDatabase = 0; @@ -125,13 +125,12 @@ try % - 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); % mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0,"adaption_technique","lms"); - + eq_post =FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption",adaption_method(adaption),"dd_mode",use_dd_mode); % Duobinary signaling (db encoded) mlse_db_enc = MLSE("DIR", [1,1], "duobinary_output", 0, "M", M, "trellis_states", PAMmapper(M,0).levels); eq_db_enc = EQ("Ne", ffe_order, "Nb", dfe_order, "training_length", len_tr, ... @@ -194,60 +193,78 @@ try if use_vnle_mlse 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); - % mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + + useviterbi = 0; + if useviterbi + mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + else + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + end [ffe_results, mlse_results] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... "precode_mode", duob_mode,... - 'showAnalysis', 1, ... + 'showAnalysis', 0, ... "postFFE", [],... "eth_style_symbol_mapping", 0); ffe_results.metrics.print; ffe_results.config.equalizer_structure = "vnle"; mlse_results.metrics.print; - + output.mlse_package{r} = mlse_results; output.vnle_package{r} = ffe_results; - + if options.append_to_db database.addProcessingResult(run_id, mlse_results.metrics, mlse_results.config); database.addProcessingResult(run_id, ffe_results.metrics, ffe_results.config); end - pf_ncoeffs = 2; - 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); - % mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + % pf_ncoeffs = 2; + % 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); + % % mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + % mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + % + % [ffe_results, mlse_results] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... + % "precode_mode", duob_mode,... + % 'showAnalysis', 0, ... + % "postFFE", [],... + % "eth_style_symbol_mapping", 0); + % + % ffe_results.metrics.print; + % mlse_results.metrics.print; + % + % output.mlse_package{r} = mlse_results; + % output.vnle_package{r} = ffe_results; + % + % if options.append_to_db + % database.addProcessingResult(run_id, mlse_results.metrics, mlse_results.config); + % database.addProcessingResult(run_id, ffe_results.metrics, ffe_results.config); + % end - [ffe_results, mlse_results] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... - "precode_mode", duob_mode,... - 'showAnalysis', 1, ... - "postFFE", [],... - "eth_style_symbol_mapping", 0); - - ffe_results.metrics.print; - mlse_results.metrics.print; - - output.mlse_package{r} = mlse_results; - output.vnle_package{r} = ffe_results; - - if options.append_to_db - database.addProcessingResult(run_id, mlse_results.metrics, mlse_results.config); - database.addProcessingResult(run_id, ffe_results.metrics, ffe_results.config); - end end if use_dbtgt + useviterbi = 0; + if useviterbi + mlse_db_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + else + mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + end + ffe_order = [50, 5, 5]; + eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + dbt_results = duobinary_target(eq_, mlse_db_, M, Scpe_sig, Symbols, Tx_bits, ... "precode_mode", duob_mode, ... 'showAnalysis', 0,... "postFFE", []); + dbt_results.metrics.print; + output.dbtgt_package{r} = dbt_results; if options.append_to_db diff --git a/Functions/EQ_structures/duobinary_target.m b/Functions/EQ_structures/duobinary_target.m index 687d6ca..3a1389f 100644 --- a/Functions/EQ_structures/duobinary_target.m +++ b/Functions/EQ_structures/duobinary_target.m @@ -23,8 +23,14 @@ if ~isempty(options.postFFE) end mlse_.DIR = [1,1]; -% mlse_sig_sd = mlse_.process(eq_signal); -[mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols); +% + +if isa(mlse_,'MLSE_viterbi') + mlse_sig_sd = mlse_.process(eq_signal); +else + [mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols); +end + mlse_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd); % precoding to mitigate error propagation, most prominently used in @@ -43,8 +49,8 @@ switch options.precode_mode tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); tx_bits_precoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded); - rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_precoded); + [~,errors_db_diff_precoded,ber_db_diff_precoded,~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); %B) Just determine BER @@ -59,24 +65,42 @@ switch options.precode_mode mlse_sig_hd_decoded = Duobinary().encode(mlse_sig_hd,"M",M); mlse_sig_hd_decoded = Duobinary().decode(mlse_sig_hd_decoded,"M",M); rx_bits_mlse_decoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_decoded); - [~,errors_db_diff_precoded,ber_db_diff_precoded,~] = calc_ber(rx_bits_mlse_decoded.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - + [~,errors_db_diff_precoded,ber_db_diff_precoded,a] = calc_ber(rx_bits_mlse_decoded.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + burst_db_precoded = count_error_bursts(a, 40); % B) Omit the Coding by comparing with demapped TX symbol sequence tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols); rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); - [bits_db,errors_db,ber_db,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + [bits_db,errors_db,ber_db,a] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + burst_db = count_error_bursts(a, 40); + cols = linspecer(8); + figure();hold on; + stem(1:40,burst_db,'LineWidth',1,'Color',cols(4,:),'Marker','_','DisplayName','w/o diff. precoder'); + stem(1:40,burst_db_precoded,'LineWidth',1,'Color',cols(3,:),'Marker','.','LineStyle','-','DisplayName','w diff. precoder'); + xlabel('Bit Error Burst Length') + ylabel('Occurence') + set(gca, 'yscale', 'log'); end % M = numel(unique(tx_symbols.signal)); rx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); [bits_db,errors_db,ber_db,errorIndice_db] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + + + + alpha = arburg(eq_noise.signal,1);%pf_.coefficients(2); alpha = alpha(2); -gmi_mlse = GMI_MLSE; -air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_mlse ./ log2(double(M)); +if isa(mlse_,'MLSE_viterbi') + gmi_mlse = NaN; + air_mlse = NaN; +else + gmi_mlse = GMI_MLSE; + air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_mlse ./ log2(double(M)); +end db_results = struct(); db_results.metrics = Metricstruct; diff --git a/Functions/EQ_structures/ffe.m b/Functions/EQ_structures/ffe.m index e3ecc04..1f927b1 100644 --- a/Functions/EQ_structures/ffe.m +++ b/Functions/EQ_structures/ffe.m @@ -49,7 +49,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_air(eq_signal_sd, tx_symbols, "skip_front", 10000, "skip_end", 10000); +[gmi] = 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)); [evm_total, evm_lvl] = calc_evm(eq_signal_sd, tx_symbols); [std_total, std_lvl] = calc_std(eq_signal_sd, tx_symbols); diff --git a/Functions/EQ_structures/vnle_postfilter_mlse.m b/Functions/EQ_structures/vnle_postfilter_mlse.m index 2ad8ea2..b86dc93 100644 --- a/Functions/EQ_structures/vnle_postfilter_mlse.m +++ b/Functions/EQ_structures/vnle_postfilter_mlse.m @@ -34,7 +34,7 @@ end % FFE or VNLE [eq_signal_sd, eq_noise] = eq_.process(rx_signal, tx_symbols); -% Apply post-FFE if provided +% Apply post-FFE if provided (does not work properly at the moment...very sad) if ~isempty(options.postFFE) tic [eq_signal_sd, eq_noise] = options.postFFE.process(eq_signal_sd, tx_symbols); @@ -43,13 +43,18 @@ end % Hard decision on VNLE output eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); -eq_signal_hd.spectrum("displayname",'after full response FFE','fignum',2025,'normalizeTo0dB',1); % Process through postfilter and MLSE [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); mlse_.DIR = pf_.coefficients; -[mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(mlse_sig_sd,tx_symbols); +GMI_MLSE = NaN; +if isa(mlse_,'MLSE_viterbi') + [mlse_sig_sd] = mlse_.process(mlse_sig_sd); +else + [mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(mlse_sig_sd,tx_symbols); +end + mlse_sig_hd = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).quantize(mlse_sig_sd); %% Calculate BER based on precoding mode @@ -58,7 +63,12 @@ mlse_sig_hd = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).qua %% Calculate performance metrics % VNLE metrics [snr_vnle, snr_vnle_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal); -[gmi_vnle] = calc_air(eq_signal_sd, tx_symbols, "skip_front", 10000, "skip_end", 10000); +% [gmi_vnle] = calc_air(eq_signal_sd, tx_symbols, "skip_front", 10000, "skip_end", 10000); + +%calculate bitwise GMI +[gmi_vnle] = calc_ngmi(eq_signal_sd,tx_symbols); +gmi_vnle = max(gmi_vnle,0); + air_vnle = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_vnle ./ log2(double(M)); [evm_vnle_total, evm_vnle_lvl] = calc_evm(eq_signal_sd, tx_symbols); [std_vnle_total, std_vnle_lvl] = calc_std(eq_signal_sd, tx_symbols); @@ -67,7 +77,7 @@ air_vnle = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_vnle ./ log2(dou % MLSE metrics alpha = arburg(eq_noise.signal,1);%pf_.coefficients(2); alpha = alpha(2); -gmi_mlse = GMI_MLSE; +gmi_mlse = max(GMI_MLSE,0); air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_mlse ./ log2(double(M)); %% Display analysis if requested @@ -77,6 +87,10 @@ end + + + + %% Prepare output structures % Determine postFFE order if ~isempty(options.postFFE) @@ -119,8 +133,6 @@ ffe_results.config.eq = jsonencode(eq_); ffe_results.config.equalizer_structure = int32(equalizer_structure.vnle); ffe_results.config.comment = 'function: vnle_postfilter_mlse - FFE part'; - - mlse_results = struct(); mlse_results.metrics = Metricstruct; mlse_results.metrics.result_id = NaN; @@ -141,7 +153,6 @@ mlse_results.metrics.Alpha = alpha; mlse_results.metrics.MLSE_dir = mlse_.DIR; % Create MLSE results structure - mlse_results.config = Equalizerstruct(); mlse_results.config.eq = jsonencode(eq_); mlse_.DIR = length(mlse_.DIR)-1; @@ -203,7 +214,7 @@ switch precode_mode mlse_sig_hd_decoded = Duobinary().encode(mlse_sig_hd, "M", M); mlse_sig_hd_decoded = Duobinary().decode(mlse_sig_hd_decoded, "M", M); rx_bits_mlse_decoded = mapper.demap(mlse_sig_hd_decoded); - [~, errors.mlse_precoded, ber.mlse_precoded, ~] = calc_ber(rx_bits_mlse_decoded.signal, tx_bits.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1); + [~, errors.mlse_precoded, ber.mlse_precoded, err_loc_precoded] = calc_ber(rx_bits_mlse_decoded.signal, tx_bits.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1); % B) Omit the Coding by comparing with demapped TX symbol sequence tx_bits_demapped = mapper.demap(tx_symbols); @@ -213,6 +224,19 @@ switch precode_mode rx_bits_mlse = mapper.demap(mlse_sig_hd); [numbits.mlse, errors.mlse, ber.mlse, error_locations.mlse] = calc_ber(rx_bits_mlse.signal, tx_bits_demapped.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1); end + +% cols = linspecer(8); +% burst_precoded = count_error_bursts(err_loc_precoded, 40); +% burst_normal = count_error_bursts(error_locations.mlse, 40); +% figure();hold on; +% stem(1:40,burst_normal,'LineWidth',1,'Color',cols(1,:),'Marker','_','DisplayName','w/o diff. precoder'); +% stem(1:40,burst_precoded,'LineWidth',1,'Color',cols(2,:),'Marker','.','LineStyle','-','DisplayName','w diff. precoder'); +% xlabel('Bit Error Burst Length') +% ylabel('Occurence') +% set(gca, 'yscale', 'log'); + + + end diff --git a/Functions/EQ_visuals/show2Dconstellation.m b/Functions/EQ_visuals/show2Dconstellation.m new file mode 100644 index 0000000..8596676 --- /dev/null +++ b/Functions/EQ_visuals/show2Dconstellation.m @@ -0,0 +1,85 @@ +function [symbols_for_lvl,avg_for_lvl] = show2Dconstellation(eq_signal,ref_symbols,options) +arguments + eq_signal + ref_symbols + options.fignum (1,1) double = NaN % Default to NaN if not provided + options.displayname (1,:) char = '' % Default to an empty string if not provided +end + +plot_shit = 1; + +if isa(eq_signal,'Signal') + eq_signal = eq_signal.signal; +end +if isa(ref_symbols,'Signal') + ref_symbols = ref_symbols.signal; +end + + +if plot_shit + % Determine the figure number to use or create a new figure + if isnan(options.fignum) + fig = figure; % Create a new figure and get its handle + else + fig = figure(options.fignum); % Use the specified figure number + clf; + end +end + + +rx_symbols = eq_signal; %./ rms(eq_signal); +correct_symbols = ref_symbols; + +col = cbrewer2('Paired',numel(unique(correct_symbols))*2); +ccnt = -1; + +levels = unique(correct_symbols); +M = numel(levels); +symbols_for_lvl = NaN(numel(levels),length(correct_symbols)); +start = 1; +ende = length(correct_symbols); + + + + +for l = 1:numel(levels) + ccnt = ccnt+2; + + level_amplitude = levels(l); + + symbols_for_lvl(l,correct_symbols==level_amplitude) = rx_symbols(correct_symbols==level_amplitude); + + statistical_mean(l) = mean(symbols_for_lvl(l,:),'omitnan'); + +end + +tx_even = correct_symbols(1:2:end); +tx_odd = correct_symbols(2:2:end); + +rx_even = rx_symbols(1:2:end); +rx_odd = rx_symbols(2:2:end); + +D_even = tx_even - rx_even; +D_odd = tx_odd - rx_odd; + +D = sqrt(D_even.^2 + D_odd.^2); + +[X,Y] = meshgrid(levels, levels); +[X_,Y_] = meshgrid(statistical_mean, statistical_mean); + +hold on; +scatter(rx_even,rx_odd,5*ones(1,length(D)),D,'.','DisplayName',['Even/Odd PAM-',num2str(M)],'MarkerEdgeColor',col(2,:)); +% colormap(gca,flip(cbrewer2('Spectral',100))) +colormap(gca,'hsv'); +scatter(X_(:), Y_(:), 2, 'x', 'LineWidth', 10, 'MarkerEdgeColor', col(6,:),'DisplayName','Statistical Rx Levels'); +scatter(X(:), Y(:), 2, 'x', 'LineWidth', 10, 'MarkerEdgeColor', col(4,:),'DisplayName','Tx Levels'); +xlim([floor(min(X_(:)))-1, ceil(max(X_(:)))+1]); +ylim([floor(min(Y_(:)))-1, ceil(max(Y_(:)))+1]); + +yticks((levels(1:end-1) + levels(2:end)) / 2); +xticks((levels(1:end-1) + levels(2:end)) / 2); +legend +xlabel('even symbols'); +ylabel('odd symbols'); +axis equal; grid on; +end diff --git a/Functions/EQ_visuals/showLevelScatter.m b/Functions/EQ_visuals/showLevelScatter.m index 70ddb39..f39dcb8 100644 --- a/Functions/EQ_visuals/showLevelScatter.m +++ b/Functions/EQ_visuals/showLevelScatter.m @@ -4,7 +4,7 @@ arguments ref_symbols options.fignum (1,1) double = NaN % Default to NaN if not provided options.displayname (1,:) char = '' % Default to an empty string if not provided - options.f_sym = []; + options.f_sym =1e6; end plot_shit = 1; @@ -30,7 +30,7 @@ if plot_shit end -rx_symbols = eq_signal ./ rms(eq_signal); +rx_symbols = eq_signal; %./ rms(eq_signal); correct_symbols = ref_symbols; f_sym = options.f_sym; @@ -47,7 +47,6 @@ for l = 1:numel(levels) level_amplitude = levels(l); - symbols_for_lvl(l,correct_symbols==level_amplitude) = rx_symbols(correct_symbols==level_amplitude); std_lvl(l) = std(symbols_for_lvl(l,:),'omitnan'); xax_in_sec = ((1:length(correct_symbols)) / f_sym) * 1e6; @@ -68,7 +67,8 @@ for l = 1:numel(levels) ccnt = ccnt+2; level_amplitude = levels(l); - movmean = 1/250 .* movsum(rx_symbols(correct_symbols==level_amplitude),[250/2,250/2], 'Endpoints', 'fill'); + L = 500; + movmean = 1/L .* movsum(rx_symbols(correct_symbols==level_amplitude),[L/2,L/2], 'Endpoints', 'fill'); avg_for_lvl(l,correct_symbols==level_amplitude) = movmean; @@ -125,7 +125,7 @@ if 0 end if plot_shit -yline(levels); +% yline(levels); xlabel('Time in $\mu$s'); ylabel('Normalized Amplitude'); ylim([-3 3]); diff --git a/Functions/Job_Processing/loadAndSyncSignalDataFromDb.m b/Functions/Job_Processing/loadAndSyncSignalDataFromDb.m index 1b76fd4..c10f1d7 100644 --- a/Functions/Job_Processing/loadAndSyncSignalDataFromDb.m +++ b/Functions/Job_Processing/loadAndSyncSignalDataFromDb.m @@ -12,12 +12,12 @@ function [Bits, Symbols, Scpe_cell, found_sync] = loadAndSyncSignalDataFromDb(da % found_sync - Boolean indicating if synchronization was successful found_sync = 0; -tempLocalStorage = 0; +tempLocalStorage = 1; % Define the fixed storage directory relative to the user's MATLAB preferences directory storage_dir = fullfile(prefdir, 'temp_sync_data'); -% Part A: Check and load from local storage if available +% Part A: Check and load from local storage if available- if tempLocalStorage == 1 local_filename = fullfile(storage_dir, sprintf('sync_data_run_%s.mat', num2str(dataTable.run_id))); if exist(local_filename, 'file') diff --git a/Functions/Metrics/calc_snr.m b/Functions/Metrics/calc_snr.m index 2a06f49..ff25cd9 100644 --- a/Functions/Metrics/calc_snr.m +++ b/Functions/Metrics/calc_snr.m @@ -39,6 +39,5 @@ function [snr_all, snr_per_level] = calc_snr(tx_signal, eq_noise) % Compute the SNR for these indices snr_per_level(i) = snr(tx_signal(idx), eq_noise(idx)); - histogram(eq_noise(idx)); end end diff --git a/Functions/Metrics/count_error_bursts.m b/Functions/Metrics/count_error_bursts.m index fdc5fc1..47d5439 100644 --- a/Functions/Metrics/count_error_bursts.m +++ b/Functions/Metrics/count_error_bursts.m @@ -21,8 +21,10 @@ function burst_count = count_error_bursts(err_pos, max_burst_length) % Check if the burst length exceeds any of the thresholds for i = 1:max_burst_length - if burst_length > i + if burst_length == i burst_count(i) = burst_count(i) + 1; + else + end end end diff --git a/projects/ECOC_2025/dsp_test/hyperparam_tuning.m b/projects/ECOC_2025/dsp_test/hyperparam_tuning.m index 75f415e..aec4165 100644 --- a/projects/ECOC_2025/dsp_test/hyperparam_tuning.m +++ b/projects/ECOC_2025/dsp_test/hyperparam_tuning.m @@ -6,10 +6,10 @@ Symbols = load("imdd_simulation\projects\ECOC_2025\dsp_test\pam4_symbols.mat");S eq_ = FFE_DCremoval("epochs_tr",5,"epochs_dd",3,"len_tr",4096*2,"mu_dd",1e-5,"mu_tr",0,"order",50,"sps",2,"decide",0,"mu_dc",0.005,"dc_buffer_len",1); -savePath = 'Z:\2025\ECOC Silas\ecoc_2025\'; -databasePath = 'C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; -database_name = 'ecoc2025_loops.db'; -db = DBHandler("type","mysql"); +% savePath = 'Z:\2025\ECOC Silas\ecoc_2025\'; +% databasePath = 'C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; +% database_name = 'ecoc2025_loops.db'; +% db = DBHandler("type","mysql"); params = (logspace(-6,-2,20)); params = floor((logspace(2,3,20))); @@ -19,16 +19,16 @@ for i = 1:length(params) eq_ = FFE_DCremoval_adaptive_mu("epochs_tr",5,"epochs_dd",3,"len_tr",4096*2,"mu_dd",... 0.0002,"mu_tr",0,"order",25,"sps",2,"decide",0,... - "mu_dc",0.005,"dc_buffer_len",224, ... + "mu_dc",0.005,"dc_buffer_len",1, ... "ffe_buffer_len",1,... "smoothing_buffer_length",0,... - "smoothing_buffer_update",224); + "smoothing_buffer_update",1); - result = vnle(eq_,4,Scpe_sig,Symbols,Tx_bits,"precode_mode",db_mode.no_db,'showAnalysis',1,"postFFE",[],"eth_style_symbol_mapping",0); - ber_ffe(i) = result.ber_vnle; + result = ffe(eq_,4,Scpe_sig,Symbols,Tx_bits,"precode_mode",db_mode.no_db,'showAnalysis',1,"postFFE",[],"eth_style_symbol_mapping",0); + ber_ffe(i) = result.metrics; fprintf(" FFE Results: %.2e\n", ber_ffe(i)); - db.addProcessingResult(run_id, result.resultsVNLE, result.equalizerConfigVNLE); + % db.addProcessingResult(run_id, result.resultsVNLE, result.equalizerConfigVNLE); % eq_ = FFE_adaptive_decision("epochs_tr",5,"epochs_dd",3,"len_tr",4096*2,"mu_dd",... % 0.0003,"mu_tr",0,"order",50,"sps",2,"decide",1,"buffer_length",params(i)); diff --git a/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m index 3b2696c..297f17b 100644 --- a/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m +++ b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m @@ -1,25 +1,24 @@ -basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +% basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; % database = DBHandler("pathToDB",[basePath,'silas_labor.db']); -database = DBHandler("type",'mysql'); +database = DBHandler("type",'mysql','dataBase','labor'); filterParams = database.tables; -%filterParams.Runs.loop_id = 209; -filterParams.Configurations = struct( ... - 'symbolrate', 112e9, ... %[224,336,360,390,420,448] - 'fiber_length', 0, ... - 'db_mode', '"no_db"', ... - 'interference_attenuation', [], ... - 'interference_path_length', [], ... - 'is_mpi', 1, ... - 'pam_level', 4, ... - 'wavelength', 1310, ... - 'precomp_amp', [], ... - 'signal_attenuation', [], ... - 'v_awg', [], ... - 'v_bias', [] ... - ); - +filterParams.Runs.loop_id = 209; +% filterParams.Configurations = struct( ... +% 'symbolrate', 112e9, ... %[224,336,360,390,420,448] +% 'fiber_length', 0, ... +% 'db_mode', '"no_db"', ... +% 'interference_attenuation', [], ... +% 'interference_path_length', 1000, ... +% 'is_mpi', 1, ... +% 'pam_level', 4, ... +% 'wavelength', 1310, ... +% 'precomp_amp', [], ... +% 'signal_attenuation', [], ... +% 'v_awg', [], ... +% 'v_bias', [] ... +% ); % if 1 % % filterParams.EqualizerParameters.dc_buffer_len = 1; @@ -28,9 +27,15 @@ filterParams.Configurations = struct( ... % filterParams.EqualizerParameters.smoothing_buffer_update = 224; % filterParams.EqualizerParameters.DCmu = 0; % end +a = database.getTableFieldNames('Runs'); +b = database.getTableFieldNames('Results'); +c = database.getTableFieldNames('EqualizerParameters'); +d = [a;b;c]; -selectedFields = {'Configurations.run_id' 'Runs.loop_id' 'Runs.date_of_run' 'Runs.rx_raw_path' 'Configurations.bitrate' 'Configurations.v_bias' 'Configurations.v_awg' 'Configurations.precomp_amp' 'Configurations.symbolrate' 'Configurations.pam_level'... - 'Configurations.db_mode' 'Configurations.rop_attenuation' 'Configurations.is_mpi' 'Configurations.interference_attenuation' 'Configurations.interference_path_length' 'Configurations.signal_attenuation' ... +[dataTable,~] = database.queryDB(filterParams, d); + +selectedFields = {'Configurations.run_id' 'Runs.loop_id' 'Runs.date_of_run' 'Runs.rx_raw_path' 'Runs.bitrate' 'Runs.v_bias' 'Runs.v_awg' 'Runs.precomp_amp' 'Runs.symbolrate' 'Runs.pam_level'... + 'Runs.db_mode' 'Runs.rop_attenuation' 'Runs.is_mpi' 'Runs.interference_attenuation' 'Runs.interference_path_length' 'Runs.signal_attenuation' ... 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'EqualizerParameters.dc_buffer_len' 'EqualizerParameters.ffe_buffer_len' 'EqualizerParameters.smoothing_buffer_len' 'EqualizerParameters.smoothing_buffer_update' 'EqualizerParameters.DCmu' 'Measurements.power_pd_in' ... 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.EVM' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'}; diff --git a/projects/ECOC_2025/theory/analytic_mpi_evaluation.m b/projects/ECOC_2025/theory/analytic_mpi_evaluation.m new file mode 100644 index 0000000..9df4f5b --- /dev/null +++ b/projects/ECOC_2025/theory/analytic_mpi_evaluation.m @@ -0,0 +1,82 @@ +% This script is used to evaluate Fig. 1b) in the paper "Adaptive Removal of Multipath Interference in Short Reach 112 GBd PAM-4 IM/DD Systems" + +%% Parameters +df = 1e6;%150e3; % Laser linewidth [Hz] +SIR_dB = 20; % Interference attenuation [dB] +alpha = 10^(-SIR_dB/20); % Interference attenuation [linear] +n_fiber = 1.467; % Refractive index +c = physconst('lightspeed'); % [m/s] + +L = linspace(0,400,40); % Interference delay [m] +tau = n_fiber./c.*L; % Interference time (= tau) [s] + +tau_c = 1/(pi*df); % laser coherence time [s] +L_c = (c/n_fiber)*tau_c; % laser coherence length [m] + +var_sat = 2*alpha^2; % Analytical saturation of variance + +%% Monte–Carlo Simulation +fs = 100e9; % sampling rate [Hz] +Tsim = 50e-6; % sim duration [s] +N = round(Tsim*fs); % number of samples for each realization +max_delay_samples = round(max(tau)*fs); % largest delay that is evaluated (based on max. Interference delay) +phase_noise_std = sqrt(2*pi*df/fs); % standard dev. phase noise + +num_realizations = 10; % number of parallel runs +monte_carlo_variance = zeros(num_realizations, length(L)); +parfor r = 1:num_realizations + + % generate a realization of phase noise random walk + dphi = phase_noise_std * randn(1, N + max_delay_samples); % matlab randn process has std = 1 + phi = cumsum(dphi); + phi_direct = phi(max_delay_samples+1 : max_delay_samples+N); + var_k = zeros(1, length(L)); + for t = 1:length(tau) + + nd = round( tau(t)*fs ); % delay in samples for current interference time + phi_delayed = phi(max_delay_samples+1-nd : max_delay_samples+N-nd); %cut out interfering signal part (was earlier) + + E = exp(1j*phi_direct) + alpha*exp(1j*phi_delayed); % E-fields combined + I = abs(E).^2; % photo current as magnitude square of E-field + var_k(t) = var(I); + + end + monte_carlo_variance(r, :) = var_k; +end + +avg_of_mc_variances = mean(monte_carlo_variance, 1); +std_of_mc_variances = std(monte_carlo_variance, 0, 1); + +%% Analytic variance +analytic_variance = 2*alpha^2 * (1 - exp(-2*pi*df.*tau)).^2; + +%% Plot +cols = [0.3467 0.5360 0.6907 + 0.9153 0.2816 0.2878 + 0.4416 0.7490 0.4322]; + +coherence_length_multiples = 0.5:0.5:ceil(L(end)/L_c); + +figure(); +hold on; +plot(L, avg_of_mc_variances, 'LineWidth',2, 'DisplayName','Simulation','Color',cols(1,:),'LineStyle','-'); +errorbar(L, avg_of_mc_variances,std_of_mc_variances, 'LineWidth',0.7,'LineStyle','none', 'DisplayName','Simulation','Color',cols(1,:),'HandleVisibility','off'); + +plot(L, analytic_variance, 'LineWidth',2, 'DisplayName','Analytic','Color',cols(2,:),'LineStyle','-'); +xticks(coherence_length_multiples.*L_c); +xticklabels(round(coherence_length_multiples.*L_c)); + +norm_to_coherence_len = 1; +if norm_to_coherence_len + xticklabels(coherence_length_multiples); + xlabel('$\tau_c$', 'FontSize',12); +end + +xline(L_c.*coherence_length_multiples, 'LineWidth',1.5, 'DisplayName','Coh. Length','HandleVisibility','off','Color',[0.7,0.7,0.7],'LineStyle','-'); +xlim([0,L(end)]); +yline(var_sat, '-.k','LineWidth',1.5, 'DisplayName','Saturation: 2$\alpha ^2$'); +xlabel('Interference Delay [m]', 'FontSize',12); +grid on; +ylabel('Intensity Variance', 'FontSize',12); +title(sprintf('MPI Variance; %d MHz; SIR: %d dB',df.*1e-6,SIR_dB), 'FontSize',14); +legend('Location','southeast'); diff --git a/projects/ECOC_2025/theory/coherence_length_plot.m b/projects/ECOC_2025/theory/coherence_length_plot.m new file mode 100644 index 0000000..e6d040d --- /dev/null +++ b/projects/ECOC_2025/theory/coherence_length_plot.m @@ -0,0 +1,32 @@ +%% Parameters +df = linspace(100e3,50e6,10000); % Laser FWHM linewidth [Hz] +n_fiber = 1.467; % Fiber group index +c = 3e8; % Speed of light [m/s] + +% Compute coherence length (1/e of mean-fringe decay) +tau_c = 1./(pi*df); +L_c = (c/n_fiber) .* tau_c; % Coherence length [m] + +%% Plot +figure('Color','w'); +loglog(df/1e6, L_c, 'LineWidth',2,'LineStyle','-'); % linewidth in MHz +xticks([0.1, 1, 10, 50]); +yticks([1, 10, 100, 1000]); +yticklabels({'1','10','100','1000'}) +grid on; box on; +xlabel('Laser linewidth [MHz]','FontSize',12,'Interpreter','none'); +ylabel('Coherence length [m]','FontSize',12,'Interpreter','none'); +title('Coherence Length vs. Laser Linewidth','FontSize',14,'Interpreter','none'); + +%% Annotate some key points +hold on; +freqs = [150e3, 1e6, 10e6, 50e6]; % [Hz] +for f = freqs + x = f/1e6; + y = (c/n_fiber) * (1/(pi*f)); + scatter(x,y,'Marker','x','LineWidth',1,'MarkerEdgeColor','black'); + + text(x*1.1,y, sprintf('%.2f MHz', f/1e6), ... + 'FontSize',10,'HorizontalAlignment','left'); + +end diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/ber_vs_rate.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/ber_vs_rate.m index 7039713..d39a24c 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/ber_vs_rate.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/ber_vs_rate.m @@ -4,30 +4,39 @@ dataBase = 'labor_highspeed';%'C:\Users\Silas\Documents\MATLAB\Datensätze\sio db = DBHandler("dataBase", [dataBase], "type", database_type); fp = QueryFilter(); -% fp.where('Runs', 'run_id','EQUALS', 987); +% fp.where('Runs', 'run_id','EQUALS', 987); M = 8; -fp.where('Runs', 'pam_level','EQUALS', M); -% fp.where('Runs', 'bitrate','LESS_THAN', 310e9); -fp.where('Runs', 'fiber_length','EQUALS', 2); -fp.where('Runs', 'is_mpi','EQUALS', 0); -% fp.where('Runs', 'interference_path_length','EQUALS', 1000); -% fp.where('Runs', 'loop_id','GREATER_THAN', 11); -% fp.where('Runs', 'sir','EQUALS',18); -fp.where('Runs', 'wavelength','EQUALS', 1310); -% fp.where('Runs', 'db_mode','EQUALS', 1); -fp.where('Runs', 'rop_attenuation','EQUALS', 0); +fp.where('Runs', 'pam_level','EQUALS', M); +% fp.where('Runs', 'bitrate','LESS_THAN', 310e9); +fp.where('Runs', 'fiber_length','EQUALS', 2); +fp.where('Runs', 'is_mpi','EQUALS', 0); +% fp.where('Runs', 'interference_path_length','EQUALS', 1000); +% fp.where('Runs', 'loop_id','GREATER_THAN', 11); +% fp.where('Runs', 'sir','EQUALS',18); +fp.where('Runs', 'wavelength','EQUALS', 1293); +% fp.where('Runs', 'db_mode','EQUALS', 0); +fp.where('Runs', 'rop_attenuation','EQUALS', 0); -% [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('dashboard')); +fields = db.getTableFieldNames('power_state_info'); +fields = [fields; db.getTableFieldNames('dashboard_ungrouped_new')]; +[dataTable,~] = db.queryDB(fp, fields); eqstructures = unique(dataTable.equalizer_structure); % Create the figure -figure(10); +showFiltered = true; +showPrecoded = false; +show_bitrate = true; +figure(5); hold on -for pre_emph = [1] - - dbmode_filtered = dataTable(dataTable.db_mode == ~pre_emph,:); - for eqs = [equalizer_structure.vnle, equalizer_structure.vnle_pf_mlse , equalizer_structure.vnle_db_mlse] +for eqs = [equalizer_structure.vnle, equalizer_structure.vnle_pf_mlse] + + % figure('Name',string([char(eqs),''])); + % hold on + + for pre_emph = [0,1] + + dbmode_filtered = dataTable(dataTable.db_mode == ~pre_emph,:); eq_choice = equalizer_structure(eqs); @@ -37,50 +46,238 @@ for pre_emph = [1] end eq_filtered = dbmode_filtered(dbmode_filtered.equalizer_structure == eq_choice,:); - - symbolrate_sorted = sortrows(eq_filtered,{'symbolrate','min_BER_precoded'}, 'ascend'); - [~, ia] = unique(symbolrate_sorted.symbolrate, 'first'); + + % ===== NEW: compute averages + per-row keep masks (robust filtering) ===== + [Tav, keepMask, keepMaskP] = avgBerBySymbolrate(eq_filtered); % <= NEW + + % x-values (bitrate) for raw points (same mapping as your lines) + M = unique(eq_filtered.pam_level); % (assumes single PAM per curve) + if show_bitrate + x_raw = eq_filtered.symbolrate.*1e-9 .* floor(log2(M)*10)/10; + else + x_raw = eq_filtered.symbolrate.*1e-9; + end + + % ===== NEW: scatter kept raw BER points (hidden from legend) ===== + cols = cbrewer2('Paired',12); + thisColor = cols((2*eqs)+1+pre_emph,:); + scatter(x_raw(keepMask), ... % kept points + eq_filtered.BER(keepMask), ... + 14, thisColor, 'filled', ... + 'MarkerFaceAlpha', 0.35, ... + 'MarkerEdgeAlpha', 0.35, ... + 'HandleVisibility','off'); + + if showPrecoded + scatter(x_raw(keepMaskP), ... % kept precoded points + eq_filtered.BER_precoded(keepMaskP), ... + 14, thisColor, 'filled', ... + 'Marker', 'square', ... + 'MarkerFaceAlpha', 0.35, ... + 'MarkerEdgeAlpha', 0.35, ... + 'HandleVisibility','off'); + end + + % ===== NEW: optionally show filtered-out points in red ===== + if showFiltered + bad = ~keepMask; + if any(bad) + scatter(x_raw(bad), eq_filtered.BER(bad), ... + 18, 'r', 'x', 'LineWidth', 1.2, ... + 'HandleVisibility','off'); + end + badp = ~keepMaskP; + if any(badp) + scatter(x_raw(badp), eq_filtered.BER_precoded(badp), ... + 18, 'r', '+', 'LineWidth', 1.2, ... + 'HandleVisibility','off'); + end + end + + % Keep your sorting and one-per-symbolrate behavior (using Tav) + symbolrate_sorted = sortrows(Tav,{'symbolrate','avg_BER_calc'}, 'ascend'); + [~, ia] = unique(symbolrate_sorted.symbolrate, 'first'); symbolrate_sorted = symbolrate_sorted(ia, :); - - % Example data (replace these with your real vectors) - symbolrate = symbolrate_sorted.symbolrate.*1e-9; % in baud - bitrate = symbolrate * floor(log2(M)*10)/10; - ber = symbolrate_sorted.min_BER; % BER - ber_precoded = symbolrate_sorted.min_BER_precoded; % BER - cols = cbrewer2('Paired',12); - cols = [0.4660 0.6740 0.1880 ; 0.9290 0.6940 0.1250 ; 0 0.4470 0.7410; 0.4940 0.1840 0.5560]; %VNLE; PF ; DFE ; DB tgt - - dname = [char(eq_choice)]; - if pre_emph + + if show_bitrate + % Bitrate for the averaged curves (unchanged) + xraw = symbolrate_sorted.symbolrate.*1e-9 .* floor(log2(M)*10)/10; + else + xraw = symbolrate_sorted.symbolrate.*1e-9; + end + % Use the MATLAB-averaged BERs + ber = symbolrate_sorted.avg_BER_calc; + ber_precoded = symbolrate_sorted.avg_BER_precoded_calc; + + + dname = strrep([char(eq_choice)],'_',' '); + if pre_emph dname = [dname,' with pre-emph.']; else dname = [dname,' w/o pre-emph.']; end - plot(bitrate, ber, 'LineWidth', 1.5, 'MarkerSize', 5,'Marker','o','LineStyle','-','Color',cols((2*eqs)+1+pre_emph,:),'MarkerEdgeColor',cols((2*eqs)+1+pre_emph,:),'MarkerFaceColor',[1,1,1],'DisplayName',[dname]); - plot(bitrate, ber_precoded, 'LineWidth', 1.5, 'MarkerSize', 5,'Marker','square','LineStyle',':','Color',cols((2*eqs)+1+pre_emph,:),'MarkerEdgeColor',cols((2*eqs)+1+pre_emph,:),'MarkerFaceColor',[1,1,1],'DisplayName',[dname,'; pre-coded']); - grid on; - - % Axis labels and title - xlabel('Baud Rate GBaud', 'FontSize', 12); - ylabel('BER', 'FontSize', 12); - title('BER vs. Baud Rate', 'FontSize', 14, 'FontWeight', 'bold'); - - % Improve tick formatting - set(gca, 'XScale', 'linear', ... - 'YScale', 'log', ... - 'TickLabelInterpreter', 'latex', ... - 'FontSize', 11); - legend + plot(xraw, ber, ... + 'LineWidth', 1.5, 'MarkerSize', 5, ... + 'Marker','o','LineStyle','-', ... + 'Color',thisColor,'MarkerEdgeColor',thisColor,'MarkerFaceColor',[1,1,1], ... + 'DisplayName', dname); - xticks(bitrate); + if showPrecoded + plot(xraw, ber_precoded, ... + 'LineWidth', 1.5, 'MarkerSize', 5, ... + 'Marker','square','LineStyle',':', ... + 'Color',thisColor,'MarkerEdgeColor',thisColor,'MarkerFaceColor',[1,1,1], ... + 'DisplayName', [dname,'; pre-coded']); + end + + grid on; + + if show_bitrate + xlabel('Net bitrate [GBps]', 'FontSize', 12); + else + xlabel('Symbol rate [GBd]', 'FontSize', 12); + end + ylabel('BER', 'FontSize', 12); + title('BER vs. Baud Rate','FontSize', 14, 'FontWeight', 'bold'); + + set(gca, 'XScale', 'linear', ... + 'YScale', 'log', ... + 'TickLabelInterpreter', 'latex', ... + 'FontSize', 11); - % Optional: tighten axis limits - xlim([min(bitrate), max(bitrate)]); - ylim([5e-5, 0.5]); + xticks(xraw); + if show_bitrate + xticks(200:25:500); + xlim([350 500]); + else + xlim([min(xraw), max(xraw)]); + end + + ylim([5e-4, 0.3]); end - + yline([2.2e-4, 4.85e-3, 2e-2],'LineWidth',1,'LineStyle','--','HandleVisibility','off'); +end + +function [Tav, keepAll, keepAllP] = avgBerBySymbolrate(T, ZT, MIN_G) +% Minimal robust averaging of BER per symbolrate (+ masks for kept points). +% Usage: [Tav, keepAll, keepAllP] = avgBerBySymbolrate(T, ZT, MIN_G) +% Defaults: ZT=3 (MAD z-thresh in log10), MIN_G=2 (min points to filter) + + if nargin < 2, ZT = 5; end + if nargin < 5, MIN_G = 0; end + + hasP = ismember('BER_precoded', T.Properties.VariableNames); + hasNB = ismember('numBits', T.Properties.VariableNames); + + [G,~,idx] = unique(T.symbolrate); + nG = numel(G); + + avgBER = nan(nG,1); + avgBERp = nan(nG,1); + keepAll = false(height(T),1); + keepAllP = false(height(T),1); + + for gi = 1:nG + r = idx==gi; + + x = T.BER(r); + nb = hasNB * T.numBits(r) + ~hasNB; % if missing, nb==1 (scalar expansion ok) + + [avgBER(gi), keepAll(r)] = rmeanBer(x, nb, ZT, MIN_G); + + if hasP + xp = T.BER_precoded(r); + [avgBERp(gi), keepAllP(r)] = rmeanBer(xp, nb, ZT, MIN_G); + end + end + + Tav = table(G, avgBER, avgBERp, ... + 'VariableNames', {'symbolrate','avg_BER_calc','avg_BER_precoded_calc'}); +end + +function [mu, keep] = rmeanBer(x, nb, ZT, MIN_G, onlyHighOutliers, minKeepThreshold) +% Robust arithmetic mean of BER with log-domain MAD filtering (returns keep mask) +% +% Params: +% x : BER values +% nb : numBits (for floor) +% ZT : MAD z-threshold +% MIN_G : min group size before filtering +% onlyHighOutliers : (bool) if true, only discard values above mean +% minKeepThreshold : values below this BER are always kept +% +% Returns: +% mu : robust mean +% keep : logical mask of kept samples + + if nargin < 5, onlyHighOutliers = false; end + if nargin < 6, minKeepThreshold = 0; end + + x(~isfinite(x)) = NaN; + + if ~isscalar(nb), nb(~isfinite(nb)) = NaN; end + if isscalar(nb) && ~isfinite(nb), nb = 1; end + + floorVal = realmin; + if ~isscalar(nb) || (isscalar(nb) && isfinite(nb) && nb~=1) + fv = 0.5 ./ max(nb, eps); % rule-of-three style floor + if isscalar(fv), floorVal = fv; else, floorVal = fv; end + end + + xAdj = x; + bad = ~isfinite(xAdj) | xAdj <= 0; + if isscalar(floorVal) + xAdj(bad) = floorVal; + else + xAdj(bad) = floorVal(bad); + end + + valid = isfinite(xAdj) & xAdj > 0; + keep = false(size(xAdj)); + + if nnz(valid)==0 + mu = NaN; return + end + if nnz(valid) < MIN_G + mu = mean(xAdj(valid),'omitnan'); keep(valid)=true; return + end + + lx = log10(xAdj(valid)); + med = median(lx,'omitnan'); + mad = median(abs(lx-med),'omitnan'); + + if mad<=0 || ~isfinite(mad) + keep(valid) = true; + mu = mean(xAdj(valid),'omitnan'); + return + end + + sigma = 1.4826*mad; + ksel = abs(lx-med) <= ZT*sigma; + + % convert to linear indices + vIdx = find(valid); + + % === Extension A: only drop high outliers === + if onlyHighOutliers + logMean = mean(lx,'omitnan'); + highIdx = lx > logMean; + ksel = ksel | ~highIdx; % always keep values below/equal to mean + end + + % === Extension B: always keep values below minKeepThreshold === + belowThr = xAdj(valid) < minKeepThreshold; + ksel = ksel | belowThr; + + keep(vIdx(ksel)) = true; + + if any(keep) + mu = mean(xAdj(keep),'omitnan'); + else + mu = mean(xAdj(valid),'omitnan'); + keep(valid) = true; + end end -yline([4.85e-3, 2e-2],'LineWidth',1,'LineStyle','--','HandleVisibility','off'); beautifyBERplot(); \ No newline at end of file diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m index ad75339..2a82dc6 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m @@ -1,6 +1,6 @@ % === SETTINGS === -dsp_options.append_to_db = 0; -dsp_options.max_occurences = 1; +dsp_options.append_to_db = 1; +dsp_options.max_occurences = 15; experiment = "highspeed_2024"; dsp_options.mode = "load_run_id"; % 'simulate' & 'load_files' @@ -39,29 +39,23 @@ end % === Get Run ID's === fp = QueryFilter(); -% fp.where('Runs', 'run_id','EQUALS', 987); -fp.where('Runs', 'pam_level','EQUALS', 4); -fp.where('Runs', 'bitrate','EQUALS', 360e9); -fp.where('Runs', 'fiber_length','EQUALS', 2); +% fp.where('Runs', 'run_id','EQUALS', 987); +M = 6; +% fp.where('Runs', 'pam_level','EQUALS', M); +% fp.where('Runs', 'bitrate','EQUALS', 480e9); +% fp.where('Runs', 'symbolrate','EQUALS', 162e9); +fp.where('Runs', 'fiber_length','EQUALS', 1); fp.where('Runs', 'is_mpi','EQUALS', 0); % fp.where('Runs', 'interference_path_length','EQUALS', 1000); % fp.where('Runs', 'loop_id','GREATER_THAN', 11); % fp.where('Runs', 'sir','EQUALS',18); -fp.where('Runs', 'wavelength','EQUALS', 1310); -fp.where('Runs', 'db_mode','EQUALS', 0); -fp.where('Runs', 'rop_attenuation','EQUALS', 0); +fp.where('Runs', 'wavelength','LESS_THAN', 1311); +% fp.where('Runs', 'db_mode','EQUALS', 0); +% fp.where('Runs', 'rop_attenuation','NOT_EQUAL', 0); % fp.where('Runs', 'power_pd_in','GREATER_THAN', 7); - [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); - -% Keep only the rows corresponding to the first occurrence of each 'sir' value -% [~, unique_indices] = unique(dataTable.sir, 'first'); -% dataTable = dataTable(unique_indices, :); - -% dataTable = dataTable(1,:); - % === Set LOOPS & Initialize DataStorage === dsp_options.parameters = struct(); % dsp_options.parameters.pf_ncoeffs = [1,2];%[0,logspace(-4,0,10)]; @@ -75,150 +69,281 @@ wh.addStorage("dbenc_package"); % === RUN IT === -[results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "serial", 'wh', wh, 'waitbar', true); +[results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "parallel", 'wh', wh, 'waitbar', true); + +results_db = results(dataTable.db_mode==1); +results_nodb = results(dataTable.db_mode==0); + +for i = 1:numel(results_db) + + % VNLE (from results_nodb) + gmi_v = cellfun(@(c) c.metrics.GMI, results_nodb{1,i}.vnle_package); + ber_v = cellfun(@(c) c.metrics.BER, results_nodb{1,i}.vnle_package); + air_v = cellfun(@(c) c.metrics.AIR, results_nodb{1,i}.vnle_package); + snr_v = cellfun(@(c) c.metrics.SNR, results_nodb{1,i}.vnle_package); + [BER_VNLE(i), idx_ber] = min(ber_v); + GMI_VNLE(i) = gmi_v(idx_ber); + AIR_VNLE(i) = air_v(idx_ber); + SNR_VNLE(i) = max(snr_v); + idx_gmi_min_vnle(i) = find(gmi_v == min(gmi_v), 1); + idx_air_max_vnle(i) = find(air_v == max(air_v), 1); + + % MLSE (from results_db) + gmi_m = cellfun(@(c) c.metrics.GMI, results_nodb{1,i}.mlse_package); + ber_m = cellfun(@(c) c.metrics.BER, results_nodb{1,i}.mlse_package); + air_m = cellfun(@(c) c.metrics.AIR, results_nodb{1,i}.mlse_package); + [BER_MLSE(i), idx_ber] = min(ber_m); + GMI_MLSE(i) = gmi_m(idx_ber); + AIR_MLSE(i) = air_m(idx_ber); + idx_gmi_min_mlse(i) = find(gmi_m == min(gmi_m), 1); + idx_air_max_mlse(i) = find(air_m == max(air_m), 1); + + % DB (from results_db, BER_precoded) + gmi_db = cellfun(@(c) c.metrics.GMI, results_db{1,i}.dbtgt_package); + ber_db = cellfun(@(c) c.metrics.BER, results_db{1,i}.dbtgt_package); + ber_db_prec = cellfun(@(c) c.metrics.BER_precoded, results_db{1,i}.dbtgt_package); + air_db = cellfun(@(c) c.metrics.AIR, results_db{1,i}.dbtgt_package); + [BER_DB(i), idx_ber] = min(ber_db); + [BER_DB_PREC(i), idx_ber] = min(ber_db_prec); + + GMI_DB(i) = gmi_db(idx_ber); + AIR_DB(i) = air_db(idx_ber); + idx_gmi_min_db(i) = find(gmi_db == min(gmi_db), 1); + idx_air_max_db(i) = find(air_db == max(air_db), 1); + + % metadata + bitrate(i) = dataTable.bitrate(i); + baudrate(i) = dataTable.symbolrate(i); +end + +STYLE_BASE = 2; % adjust this single number to scale markers & lines +MARKER_SIZE = STYLE_BASE; % marker size (MATLAB MarkerSize) +LINE_WIDTH = max(2, STYLE_BASE/3); % line width (keeps lines reasonable when STYLE_BASE large) + +% --- color map / method -> color assignment (keeps colors consistent) --- +cols = cbrewer2('Paired',8); +cols = linspecer(6); +d = 0; +cm.VNLE = cols(1 + d, :); +cm.MLSE = cols(2 + d, :); +cm.DB_precode = cols(3 + d, :); +cm.DB = cols(4 + d, :); % duobinary + +% prepare x values in GBd +xGHz = baudrate .* 1e-9; +xticks_vals = xGHz; +xtick_labels = arrayfun(@(v) sprintf('%d', round(v)), xticks_vals, 'UniformOutput', false); + +% common marker settings (filled, same face+edge color) +mk.VNLE = {'Marker','o','MarkerFaceColor',cm.VNLE,'MarkerEdgeColor',cm.VNLE,'MarkerSize',MARKER_SIZE}; +mk.MLSE = {'Marker','*','MarkerFaceColor',cm.MLSE,'MarkerEdgeColor',cm.MLSE,'MarkerSize',MARKER_SIZE}; +mk.DB_precode = {'Marker','^','MarkerFaceColor',cm.DB_precode,'MarkerEdgeColor',cm.DB_precode,'MarkerSize',MARKER_SIZE}; +mk.DB = {'Marker','d','MarkerFaceColor',cm.DB,'MarkerEdgeColor',cm.DB,'MarkerSize',MARKER_SIZE}; -% wh.getStoValue('ffe_package',0.005); -% wh.getStoValue('mlse_package',0.005); +% ---------------- FIGURE : BER ---------------- +figure(112+M); clf; hold on; +plot(xGHz, BER_VNLE, ... + 'DisplayName','VNLE', ... + mk.VNLE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +plot(xGHz, BER_MLSE, ... + 'DisplayName','MLSE', ... + mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); +plot(xGHz, BER_DB, ... + 'DisplayName','DB tgt.', ... + mk.DB{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.DB); +plot(xGHz, BER_DB_PREC, ... + 'DisplayName','Diff. Precode + DB tgt.', ... + mk.DB{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.DB); +yline(4.85e-3,'LineWidth',1,'HandleVisibility','off'); +yline(2.2e-4,'LineWidth',1,'HandleVisibility','off'); +xlabel('Baudrate in GBd'); +ylabel('BER'); +set(gca, 'yscale', 'log'); +set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); +grid on; +legend('Location','best'); -% [dataTable,~] = db.queryDB(fp, [db.getTableFieldNames('Runs');db.getTableFieldNames('Results');db.getTableFieldNames('Equalizer')]); + +% ---------------- FIGURE 15 : GMI ---------------- +figure(113+M); clf; hold on; +plot(xGHz, GMI_VNLE, ... + 'DisplayName','VNLE', ... + mk.VNLE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +plot(xGHz, GMI_MLSE, ... + 'DisplayName','MLSE', ... + mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); +plot(xGHz, GMI_DB, ... + 'DisplayName','DB tgt.', ... + mk.DB{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.DB); + +ylim([log2(M)-1, log2(M)]); +xlabel('Baudrate in GBd'); +ylabel('GMI'); +set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); +grid on; +legend('Location','best'); + + + +% ---------------- FIGURE 15 : AIR ---------------- +m = floor(log2(M)*10)/10; +figure(114+M); clf; hold on; +plot(xGHz, GMI_VNLE.*xGHz, ... + 'DisplayName','AIR VNLE', ... + mk.VNLE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +% duobinary has only one GMI curve (DB output) +plot(xGHz, GMI_MLSE.*xGHz, ... + 'DisplayName','AIR MLSE', ... + mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); +% MLSE symbol-wise (if present) +plot(xGHz, GMI_DB.*xGHz, ... + 'DisplayName','AIR DB tgt.', ... + mk.DB{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.DB); + +% ylim([log2(M)-1, log2(M)]); +xlabel('Baudrate in GBd'); +ylabel('AIR in Gbps'); +set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); +grid on; +legend('Location','best'); + + + + + + +% ---------------- FIGURE 15 : Information Rates ---------------- +tp = TransmissionPerformance; + + +m = floor(log2(M)*10)/10; +figure(213+M); clf; hold on; + +netrates_vnle = tp.calculateNetRate(baudrate.* m, ... + 'NGMI', GMI_VNLE./m, ... + 'BER', BER_VNLE); % -% dataTable = cleanUpTable(dataTable); -% wh_analyze = wh_adap; -% % wh_analyze = wh_dcremoval_old2; -% % wh_analyze = wh_dcremoval_old2; +% plot(xGHz, GMI_VNLE.*xGHz, ... +% 'DisplayName','GMI*R VNLE', ... +% mk.VNLE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.VNLE); % -% res = cell(wh_analyze.parameter.mu_dc.length,wh_analyze.parameter.dc_buffer_len.length); -% ber_mean = zeros(wh_analyze.parameter.mu_dc.length,wh_analyze.parameter.dc_buffer_len.length); -% for m = 1:wh_analyze.parameter.mu_dc.length -% for b = 1:wh_analyze.parameter.dc_buffer_len.length -% res{m,b}=wh_analyze.getStoValue("ffe_package",wh_analyze.parameter.mu_dc.values(m),wh_analyze.parameter.dc_buffer_len.values(b)); -% try -% cells = res{m,b}{1}; -% idx = cellfun(@(c) ~isempty(c), cells); -% ber = cellfun(@(c) c.metrics.BER, cells(idx)); -% ber_mean(m,b) = mean(ber); -% catch -% ber_mean(m,b) = NaN; -% end -% end -% end -% -% figure() -% hold on -% ber_fix = cellfun(@(c) c.ffe_package{1}.metrics.BER, results_fix); -% ber_adap_m2 = cellfun(@(c) c.ffe_package{1}.metrics.BER, results_adap); -% ber_adap_method1 = cellfun(@(c) c.ffe_package{1}.metrics.BER, results); -% plot(dataTable.sir,ber_fix,'LineWidth',1,'DisplayName',sprintf('DCt; fix mu = 0.5; p=1024'),'Marker','.','MarkerSize',10); -% plot(dataTable.sir,ber_adap_method1,'LineWidth',1,'DisplayName',sprintf('DCt; adap mu 1; p=1024'),'Marker','.','MarkerSize',10); -% plot(dataTable.sir,ber_adap_m2,'LineWidth',1,'DisplayName',sprintf('DCt; adap mu 2; p=1024'),'Marker','.','MarkerSize',10); -% xlabel('BER'); -% xlabel('SIR'); -% yline([4.85e-3,2e-2],'HandleVisibility', 'off','LineWidth',1,'LineStyle','--'); -% % ylim([9e-4, 0.5]); -% set(gca, 'YScale', 'log'); % BER is usually plotted log-scale -% legend('show', 'Location', 'best'); -% grid on; +% plot(xGHz, netrates_vnle.SDHD.NetRate.*1e-9, ... +% 'DisplayName','SD+HD VNLE', ... +% mk.VNLE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +% plot(xGHz, netrates_vnle.HD.NetRate.*1e-9, ... +% 'DisplayName','Staircase VNLE', ... +% mk.VNLE{:}, 'LineStyle','-.','LineWidth',LINE_WIDTH,'Color',cm.VNLE); % % +% % +% % MLSE symbol-wise (if present) +% plot(xGHz, GMI_MLSE.*xGHz, ... +% 'DisplayName','GMI*R MLSE', ... +% mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); % -% -% figure; clf -% -% % Create meshgrid for contourf -% [X, Y] = meshgrid(dsp_options.parameters.mu_dc, dsp_options.parameters.dc_buffer_len); -% -% % Create contour plot -% contourf(X, Y, ber_mean', 20); % 20 contour levels, adjust as needed -% -% % Set axes to logarithmic scale -% set(gca, 'XScale', 'log', 'YScale', 'log'); -% -% colormap("parula"); -% c = colorbar; -% c.Label.String = 'BER'; -% -% xlabel('\mu_{dc}'); -% ylabel('dc\_buffer\_len'); -% title('BER Optimization over \mu_{dc} and dc\_buffer\_len'); -% -% % Make plot prettier -% grid on -% set(gca, 'Layer', 'top'); % Put grid lines on top of contours -% -% -% % === Look at it === -% y_var = 'BER_precoded'; -% x_var = 'bitrate'; -% fixedVars = {'equalizer_structure', x_var}; -% -% [dataTableClean, outliersTable] = removeGroupOutliers(dataTable, fixedVars, y_var); -% -% % --- Group and aggregate --- -% dataTableGrpd_mean = groupIt(fixedVars, dataTableClean, @mean); -% dataTableGrpd_min = groupIt(fixedVars, dataTableClean, @min); -% dataTableGrpd_max = groupIt(fixedVars, dataTableClean, @max); -% -% % Choose a color map -% cols = linspecer(numel(unique(dataTableGrpd_mean.equalizer_structure))); -% -% figure; -% hold on; -% -% % Get unique equalizer structures for grouping -% unique_eq = unique(dataTableGrpd_mean.equalizer_structure); -% -% for i = 1:numel(unique_eq) -% eq_val = unique_eq(i); -% -% % Filter grouped data for this equalizer structure -% filt = dataTableGrpd_mean.equalizer_structure == eq_val; -% -% x = dataTableGrpd_mean.(x_var)(filt); -% y_mean = dataTableGrpd_mean.(y_var)(filt); -% y_min = dataTableGrpd_min.(y_var)(filt); -% y_max = dataTableGrpd_max.(y_var)(filt); -% -% % Bounds for boundedline (distance from mean) -% y_lower = y_mean - y_min; -% y_upper = y_max - y_mean; -% y_bounds = [y_lower, y_upper]; -% -% % --- Bounded line (mean ± min/max) --- -% if exist('boundedline', 'file') -% [hl, hp] = boundedline(x, y_mean, y_bounds, ... -% 'alpha', 'transparency', 0.1, ... -% 'cmap', cols(i,:), ... -% 'nan', 'fill', ... -% 'orientation', 'vert'); -% set(hl, 'LineWidth', 1.2, 'DisplayName', sprintf('Eq %s', eq_val)); -% set(hp, 'HandleVisibility', 'off'); -% else -% % If boundedline is not available, use errorbar -% errorbar(x, y_mean, y_lower, y_upper, ... -% 'o-', 'Color', cols(i,:), 'LineWidth', 1.2, ... -% 'DisplayName', sprintf('Eq %d', eq_val),'HandleVisibility', 'off'); -% end -% -% % --- Normal line (mean only) --- -% plot(x, y_mean, '-', 'Color', cols(i,:), 'LineWidth', 1.5, ... -% 'DisplayName', sprintf('Mean Eq %s', eq_val),'HandleVisibility', 'off'); -% -% % --- Scatter plot for individual points (from original data) --- -% % Filter original data for this group -% orig_filt = dataTableClean.equalizer_structure == eq_val; -% x_scatter = dataTableClean.(x_var)(orig_filt); -% y_scatter = dataTableClean.(y_var)(orig_filt); -% -% scatter(x_scatter, y_scatter, 10,cols(i,:), 'filled', ... -% 'MarkerFaceAlpha', 0.5, 'DisplayName', sprintf('Scatter Eq %s', eq_val),'HandleVisibility', 'off'); -% end -% -% yline([2.2e-4,4.85e-3,2e-2],'HandleVisibility', 'off','LineWidth',1,'LineStyle','--'); -% set(gca, 'YScale', 'log'); % BER is usually plotted log-scale -% xlabel(x_var, 'Interpreter', 'none'); -% ylabel(y_var, 'Interpreter', 'none'); -% legend('show', 'Location', 'best'); -% grid on; -% title(sprintf('%s vs. %s', y_var, x_var), 'Interpreter', 'none'); -% hold off; \ No newline at end of file +% netrates_mlse = tp.calculateNetRate(baudrate.* m, ... +% 'NGMI', GMI_MLSE./m, ... +% 'BER', BER_MLSE); +% plot(xGHz, netrates_mlse.SDHD.NetRate.*1e-9, ... +% 'DisplayName','SD+HD MLSE', ... +% mk.MLSE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.MLSE); +% plot(xGHz, netrates_mlse.HD.NetRate.*1e-9, ... +% 'DisplayName','Staircase MLSE', ... +% mk.MLSE{:}, 'LineStyle','-.','LineWidth',LINE_WIDTH,'Color',cm.MLSE); + + +% duobinary has only one GMI curve (DB output) +figure(1111); clf; hold on; +plot(xGHz, GMI_DB.*xGHz, ... + 'DisplayName','GMI*R DB tgt.', ... + mk.DB{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.DB); + +netrates_db = tp.calculateNetRate(baudrate.* m, ... + 'NGMI', GMI_DB./m, ... + 'BER', BER_DB_PREC); + +plot(xGHz, netrates_db.SDHD.NetRate.*1e-9, ... + 'DisplayName','SD+HD DB', ... + mk.DB{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.DB); +plot(xGHz, netrates_db.STAIR.NetRate.*1e-9, ... + 'DisplayName','Staircase DB', ... + mk.DB_precode{:}, 'LineStyle','-.','LineWidth',LINE_WIDTH,'Color',cm.DB_precode); +plot(xGHz, netrates_db.O_FEC.NetRate.*1e-9, ... + 'DisplayName','O-FEC DB', ... + mk.DB_precode{:}, 'LineStyle','--','LineWidth',LINE_WIDTH,'Color',cm.DB_precode); +plot(xGHz, netrates_db.KP4_hamming.NetRate.*1e-9, ... + 'DisplayName','KP4 Hamming DB', ... + mk.DB_precode{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.DB_precode); + +% ylim([log2(M)-1, log2(M)]); +xlabel('Baudrate in GBd'); +ylabel('AIR in Gbps'); +set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); +grid on; +legend('Location','best'); +% xlim([1, 256]) + + + +figure(2222); clf; hold on; +plot(xGHz, GMI_MLSE.*xGHz, ... + 'DisplayName','GMI*R DB tgt.', ... + mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); + +netrates_mlse = tp.calculateNetRate(baudrate.* m, ... + 'NGMI', GMI_MLSE./m, ... + 'BER', BER_MLSE); + +plot(xGHz, netrates_mlse.SDHD.NetRate.*1e-9, ... + 'DisplayName','SD+HD DB', ... + mk.MLSE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.MLSE); +plot(xGHz, netrates_mlse.STAIR.NetRate.*1e-9, ... + 'DisplayName','Staircase DB', ... + mk.VNLE{:}, 'LineStyle','-.','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +plot(xGHz, netrates_mlse.O_FEC.NetRate.*1e-9, ... + 'DisplayName','O-FEC DB', ... + mk.VNLE{:}, 'LineStyle','--','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +plot(xGHz, netrates_mlse.KP4_hamming.NetRate.*1e-9, ... + 'DisplayName','KP4 Hamming DB', ... + mk.VNLE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.VNLE); + +% ylim([log2(M)-1, log2(M)]); +xlabel('Baudrate in GBd'); +ylabel('AIR in Gbps'); +set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); +grid on; +legend('Location','best'); +% xlim([1, 256]) + + + +figure(3333); clf; hold on; +plot(xGHz, GMI_VNLE.*xGHz, ... + 'DisplayName','GMI*R DB tgt.', ... + mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); + +netrates_vnle = tp.calculateNetRate(baudrate.* m, ... + 'NGMI', GMI_VNLE./m, ... + 'BER', BER_VNLE); + +plot(xGHz, netrates_vnle.SDHD.NetRate.*1e-9, ... + 'DisplayName','SD+HD DB', ... + mk.MLSE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.MLSE); +plot(xGHz, netrates_vnle.STAIR.NetRate.*1e-9, ... + 'DisplayName','Staircase DB', ... + mk.VNLE{:}, 'LineStyle','-.','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +plot(xGHz, netrates_vnle.O_FEC.NetRate.*1e-9, ... + 'DisplayName','O-FEC DB', ... + mk.VNLE{:}, 'LineStyle','--','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +plot(xGHz, netrates_vnle.KP4_hamming.NetRate.*1e-9, ... + 'DisplayName','KP4 Hamming DB', ... + mk.VNLE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.VNLE); + +% ylim([log2(M)-1, log2(M)]); +xlabel('Baudrate in GBd'); +ylabel('AIR in Gbps'); +set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); +grid on; +legend('Location','best'); +% xlim([1, 256]) \ No newline at end of file diff --git a/projects/IMDD_base_system/simulation_bwl.m b/projects/IMDD_base_system/simulation_bwl.m index b040f1e..40ccbb2 100644 --- a/projects/IMDD_base_system/simulation_bwl.m +++ b/projects/IMDD_base_system/simulation_bwl.m @@ -1,7 +1,8 @@ %%% Run parameters % TX M = 6; -fsym = 112e9; +m = floor(log2(M)*10)/10; +fsym = 224e9; apply_pulsef = 1; fdac = 256e9; @@ -9,20 +10,20 @@ fadc = 256e9; random_key = 2; rcalpha = 0.05; -kover = 16; +kover = 8; vbias_rel = 0.5; -u_pi = 2.9; +u_pi = 3.2; vbias = -vbias_rel*u_pi; -laser_wavelength = 1290; -laser_linewidth = 0; +laser_wavelength = 1310; +laser_linewidth = 1e6; % Channel -link_length = 10000; +link_length = 0; vnle_order1 = 50; -vnle_order2 = 3; -vnle_order3 = 3; +vnle_order2 = 0; +vnle_order3 = 0; vnle_order=[vnle_order1,vnle_order2,vnle_order3]; dfe_order = [0 0 0]; @@ -43,20 +44,22 @@ mu_dfe = 0.0004; dfe_ = sum(dfe_order)>0; doub_mode = db_mode.no_db; - -rop = [-5]; +cols = linspecer(6); +rop = [-6]; bwl = [0.5:0.1:1.5]; -fsym = [72:8:170].*1e9; -ber_vnle = []; +fsym = [120:8:256].*1e9; +fsym =150e9; + +ber_vnle = []; ber_mlse = []; ber_viterbi = []; ber_db = []; ber_db_diff_precoded = []; -gmi_vnle = []; +gmi_vnle_bitwise = []; gmi_mlse = []; gmi_mlse_db = []; -parfor r = 1:length(fsym) +for r = 1:length(fsym) Pform = Pulseformer("fsym",fsym(r),"fdac",4*fsym(r),"pulse","rc","pulselength",16,"alpha",rcalpha); @@ -74,29 +77,94 @@ parfor r = 1:length(fsym) "db_precode",db_precode,"db_encode",db_encode,... "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process(); - El_sig = AWG("fdac",fdac,"f_cutoff",fsym(r),"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig); + % El_sig = AWG("fdac",fdac,"f_cutoff",fsym(r),"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig); + El_sig = M8199B("kover",kover).process(Digi_sig); + % AWG("fdac",fdac,"f_cutoff",fsym(r),"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig); %%%%% Low-pass el. components %%%%%% - tx_bwl = 50e9; - El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); + % tx_bwl = 100e9; + % El_sig = Filter('filtdegree',3,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); %%%%% Electrical Driver Amplifier %%%%%% El_sig = El_sig.normalize("mode","oneone"); + % El_sig = El_sig.setPower(1,"dBm"); + % figure;histogram(El_sig.signal); - %%%%% MODULATE E/O CONVERSION %%%%%% + %%%%% MODULATE E/O CONVERSION %%%%% + u_pi = 3.2; + vbias = -u_pi*0.5; [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); - + + figure(15); + hold on + scatter(El_sig.signal(1:100000)+vbias,(abs(Opt_sig.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') + xlabel('Input in V') + ylabel('abs(Eopt)2 in mW','Interpreter','latex') + ylim([0 2]); + xlim([-3.2 0]); + + Opt_sig.eye(fsym(r),M,"fignum",103837); + %%%%%% Fiber %%%%%% Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); %%%%%% ROP %%%%%% Opt_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); - + + % Opt_sig.eye(fsym(r),M,"fignum",103838); + + % % Opt_sig.signal = Opt_sig.signal + 5*abs(mean(Opt_sig.signal)); + % Opt_sig.move_it_spectrum("displayname",'Opt Sig after Amp','fignum',1223323); + % Pc = abs(mean(Opt_sig.signal)).^2; % carrier power + % Ptot = mean(abs(Opt_sig.signal).^2); % total power + % Ps = max(Ptot - Pc, eps); + % Pcdb = 10*log10(Pc); + % Psdb = 10*log10(Ps); + % + % cspr_dB = 10*log10(Pc / Ps); + % + % % Minimal in-place CSPR set (real, nonnegative field constraint) + % E = Opt_sig.signal; % real field samples + % target_cspr_dB = 20; % <-- set your target CSPR (dB) + % + % % Decompose into DC + zero-mean waveform + % m = mean(E); + % x0 = E - m; % zero-mean modulation + % Ps0 = mean(x0.^2); % sideband power (fixed if shape kept) + % + % % Current CSPR (for reference) + % Pc_cur = m^2; + % Ptot_cur = mean(E.^2); + % Ps_cur = max(Ptot_cur - Pc_cur, eps); + % cspr_in = 10*log10(Pc_cur / Ps_cur); + % + % % Bias needed for target CSPR, and minimal bias to keep E>=0 + % R_tgt = 10^(target_cspr_dB/10); % Pc/Ps + % a_req = sqrt(R_tgt * Ps0); % required DC bias + % a_min = -min(x0); % to avoid negatives everywhere + % a = max(a_req, a_min); % if infeasible, lands at CSPR_min + % + % % Apply bias (preserves waveform shape) + % E_new = a + x0; + % + % % Achieved CSPR + % Pc_new = mean(E_new)^2; + % Ptot_new = mean(E_new.^2); + % Ps_new = max(Ptot_new - Pc_new, eps); + % cspr_out = 10*log10(Pc_new / Ps_new); + % + % % (Optional) show feasibility info + % cspr_min = 10*log10((a_min^2)/max(Ps0,eps)); + % disp(table(cspr_in, target_cspr_dB, cspr_min, cspr_out)); + % + % % Use E_new as your adjusted field + % Opt_sig.signal = E_new; + %%%%%% PD Square Law %%%%%% PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",random_key).process(Opt_sig); - + %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% - rx_bwl = 50e9; + rx_bwl = 70e9; PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig); % %%%%%% Low-pass Scope %%%%%% @@ -113,21 +181,29 @@ parfor r = 1:length(fsym) [~, Scpe_cell, ~, found_sync] = Scpe_sig_2sps.tsynch("reference", Symbols, "fs_ref", fsym(r), "debug_plots", 0); Rx_sig = Scpe_cell{1}; + Rx_sig = Rx_sig.normalize("mode","rms"); - if 1 + if 0 %Duobinary Targeting eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"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); - mlse_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); - + db_ref_sequence = Duobinary().encode(Symbols); db_ref_constellation = unique(db_ref_sequence.signal); [eq_signal, eq_noise] = eq_.process(Rx_sig,db_ref_sequence); - - mlse_.DIR = [1,1]; - [mlse_sig_sd,LLR,gmi_mlse_db(r)] = mlse_.process(eq_signal,Symbols); - + + viterbi = 0; + if viterbi + mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + mlse_.DIR = [1,1]; + [mlse_sig_sd] = mlse_.process(eq_signal); + else + mlse_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).get_levels ./ PAMmapper(M,0).get_scaling); + mlse_.DIR = [1,1]; + [mlse_sig_sd,LLR,gmi_mlse_db(r)] = mlse_.process(eq_signal,Symbols); + end + mlse_sig_hd = PAMmapper(M,0,"eth_style",0).quantize(mlse_sig_sd); mlse_sig_hd_precoded = Duobinary().encode(mlse_sig_hd,"M",M); mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded,"M",M); @@ -138,48 +214,84 @@ parfor r = 1:length(fsym) tx_bits_precoded = PAMmapper(M,0,"eth_style",0).demap(tx_symbols_precoded); rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd_precoded); - [~,errors_db_diff_precoded,ber_db_diff_precoded(r),~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + [~,errors_db_diff_precoded,ber_db_diff_precoded(r),a] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + burst_db_pre(r,:) = count_error_bursts(a, 15)./numel(Tx_bits.signal); %B) Just determine BER rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); - [bits_mlse,errors_mlse,ber_db(r),~] = calc_ber(rx_bits_mlse.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + [bits_mlse,errors_db,ber_db(r),a] = calc_ber(rx_bits_mlse.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + burst_db(r,:) = count_error_bursts(a, 15)./numel(Tx_bits.signal); + + fprintf('BER ber_db_diff_precoded: %.2e \n',ber_db_diff_precoded(r)); + fprintf('BER Vber_dbNLE: %.2e \n',ber_db(r)); + % figure();hold on;stem(1:15,burst_db(r,:),'LineWidth',1,'Color',cols(1,:));stem(1:15,burst_db_pre(r,:),'LineWidth',1,'Color',cols(2,:));set(gca, 'yscale', 'log'); + end + % FFE or VNLE - eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"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); + eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.00,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + % eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[50,2,2],"sps",2,"decide",0); + [eq_signal_sd, eq_noise] = eq_.process(Rx_sig, Symbols); - - [gmi_gomez(r)] = calc_air(eq_signal_sd, Symbols, "skip_front", 100, "skip_end", 100); - [gmi_vnle(r)] = calc_ngmi(eq_signal_sd,Symbols); - [gmi_bitwise(r)] = calc_gmi_bitwise(eq_signal_sd,Symbols); - - snr_vnle(r) = calc_snr(Symbols, eq_signal_sd-Symbols); - eq_signal_sd.plot("displayname",'bla','fignum',118); + showEQNoisePSD(eq_noise, "fignum",1273876,"displayname",'noise after EQ'); + [mi_gomez(r)] = calc_air(eq_signal_sd, Symbols, "skip_front", 100, "skip_end", 100); + [gmi_vnle_bitwise(r)] = calc_ngmi(eq_signal_sd,Symbols); + [gmi_bitwise_2(r)] = calc_gmi_bitwise(eq_signal_sd,Symbols); + snr_vnle(r) = calc_snr(Symbols, eq_signal_sd-Symbols); + + eq_signal_sd.plot("displayname",'bla','fignum',199); + eq_signal_sd.eye(fsym(r),M,"fignum",103837); + % Hard decision on VNLE output eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); rx_bits = PAMmapper(M,0,"eth_style",0).demap(eq_signal_hd); - [~,~,ber_vnle(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + [~,tot_err,ber_vnle(r),a] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + burst_vnle(r,:) = count_error_bursts(a, 10)./tot_err; + + showLevelConfusionMatrix(eq_signal_hd,Symbols,"M",M,"fignum",200,"displayname",'bla'); + showLevelScatter(eq_signal_sd,Symbols,"displayname",'VNLE Out','f_sym',fsym(r),'fignum',201); + show2Dconstellation(eq_signal_sd,Symbols,"displayname",'VNLE Out','fignum',2241); + fprintf('BER VNLE: %.2e \n',ber_vnle(r)); - fprintf('NGMI VNLE: %.2f \n',gmi_vnle(r)./log2(M)); + fprintf('NGMI VNLE: %.2f \n',gmi_vnle_bitwise(r)./m); + if 1 % Process through postfilter and MLSE pf_ncoeffs = 1; - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + if fsym(r) < 200e9 + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1,"coefficients",[1,0.1]); + else + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1,"coefficients",[1,0.85]); + end + + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).get_levels ./ PAMmapper(M,0).get_scaling); [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); + mlse_.DIR = pf_.coefficients; alpha(r) = pf_.coefficients(2); [signalclass_hd,LLR,gmi_mlse(r)] = mlse_.process(mlse_sig_sd,Symbols); mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(signalclass_hd); rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); - [~,~,ber_mlse(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('BER: %.2e \n',ber_mlse(r)); - fprintf('GMI MLSE: %.5f \n',gmi_mlse(r)); - + [~,tot_err,ber_mlse(r),a] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + burst_mlse(r,:) = count_error_bursts(a, 10); + + showLevelConfusionMatrix(mlse_sig_hd,Symbols,"M",M,"fignum",300,"displayname",'bla'); + + fprintf('BER MLSE: %.2e \n',ber_mlse(r)); + fprintf('NGMI MLSE: %.5f \n',gmi_mlse(r)./m); + + levels = sort(unique(Symbols.signal(:)).'); % 1×6 + pairs = reshape(mlse_sig_hd.signal,2,[]).'; + isedge = ismember(pairs, [levels(1) levels(end)]); + isforbidden = sum(isedge,2)==2; + fprintf('Found %d forbidden transitions (even→odd edges).\n', nnz(isforbidden)); + + % Process through postfilter and MLSE pf_ncoeffs = 1; pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); @@ -198,36 +310,166 @@ parfor r = 1:length(fsym) end -cols = cbrewer2('Paired',8); -d = 1; -figure(11);hold on -plot(fsym.*1e-9,alpha,'DisplayName','VNLE','Marker','x','LineStyle','-','Color',cols(1+d,:)); + + +% --- style control (one variable controls both marker size and linewidth) --- +STYLE_BASE = 2; % adjust this single number to scale markers & lines +MARKER_SIZE = STYLE_BASE; % marker size (MATLAB MarkerSize) +LINE_WIDTH = max(1.5, STYLE_BASE/3); % line width (keeps lines reasonable when STYLE_BASE large) + +% --- color map / method -> color assignment (keeps colors consistent) --- +cols = cbrewer2('Paired',8); +cols = linspecer(6); +d = 0; +cm.VNLE = cols(1 + d, :); +cm.MLSE = cols(2 + d, :); +cm.DB_precode = cols(3 + d, :); +cm.DB = cols(4 + d, :); % duobinary + +% prepare x values in GBd +xGHz = fsym .* 1e-9; +xticks_vals = xGHz; +xtick_labels = arrayfun(@(v) sprintf('%d', round(v)), xticks_vals, 'UniformOutput', false); + +% common marker settings (filled, same face+edge color) +mk.VNLE = {'Marker','none','MarkerFaceColor',cm.MLSE,'MarkerEdgeColor',cm.VNLE,'MarkerSize',MARKER_SIZE}; +mk.MLSE = {'Marker','none','MarkerFaceColor',cm.MLSE,'MarkerEdgeColor',cm.MLSE,'MarkerSize',MARKER_SIZE}; +mk.DB_precode = {'Marker','none','MarkerFaceColor',cm.DB_precode,'MarkerEdgeColor',cm.DB_precode,'MarkerSize',MARKER_SIZE}; +mk.DB = {'Marker','none','MarkerFaceColor',cm.DB,'MarkerEdgeColor',cm.DB,'MarkerSize',MARKER_SIZE}; + +% ---------------- FIGURE 11 : alpha (VNLE) ---------------- +figure(110+M); clf; hold on; +plot(xGHz, alpha, ... + 'DisplayName','VNLE', ... + mk.VNLE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.VNLE); xlabel('Baudrate in GBd'); ylabel('alpha'); +set(gca, 'XTick', xticks_vals, 'XTickLabel', xtick_labels); +grid on; +legend('Location','best'); -figure(15);hold on -plot(fsym.*1e-9,gmi_gomez,'DisplayName','gomez','Marker','x','LineStyle','-','Color',cols(1+d,:)); -plot(fsym.*1e-9,gmi_vnle,'DisplayName','vnle other','Marker','x','LineStyle','-','Color',cols(3+d,:)); -% plot(fsym.*1e-9,gmi_bitwise,'DisplayName','bitwise','Marker','x','LineStyle','--','Color',cols(5+d,:)); -plot(fsym.*1e-9,gmi_mlse,'DisplayName','MLSE','Marker','*'); -plot(fsym.*1e-9,gmi_mlse_db,'DisplayName','DB Output','Marker','*'); -ylim([log2(M)-1 log2(M)]); +% ---------------- FIGURE 15 : GMI ---------------- +figure(111+M); clf; hold on; +plot(xGHz, mi_gomez, ... + 'DisplayName','MI VNLE', ... + mk.VNLE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +plot(xGHz, gmi_vnle_bitwise, ... + 'DisplayName','GMI VNLE', ... + mk.VNLE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +% duobinary has only one GMI curve (DB output) +plot(xGHz, gmi_mlse_db, ... + 'DisplayName','GMI DB tgt.', ... + mk.DB{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.DB); +% MLSE symbol-wise (if present) +plot(xGHz, gmi_mlse, ... + 'DisplayName','GMI MLSE', ... + mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); + +ylim([log2(M)-1, log2(M)]); xlabel('Baudrate in GBd'); ylabel('GMI'); +set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); +grid on; +legend('Location','best'); +% xlim([184, 256]) + +% ---------------- FIGURE 13 : BER ---------------- +figure(312+M); hold on; +plot(xGHz, ber_vnle, ... + 'DisplayName','VNLE', ... + mk.VNLE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +plot(xGHz, ber_mlse, ... + 'DisplayName','MLSE', ... + mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); +plot(xGHz, ber_viterbi, ... + 'DisplayName','Viterbi', ... + mk.MLSE{:}, 'LineStyle','--','LineWidth',LINE_WIDTH,'Color',cm.MLSE); + +yline(4.85e-3,'LineWidth',1,'HandleVisibility','off'); +yline(2.2e-4,'LineWidth',1,'HandleVisibility','off'); + +plot(xGHz, ber_db, ... + 'DisplayName','DB tgt.', ... + mk.DB{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.DB); + +plot(xGHz, ber_db_diff_precoded, ... + 'DisplayName','Prec. + DB tgt.', ... + mk.DB{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.DB); -figure(13);hold on -plot(fsym.*1e-9,ber_vnle,'DisplayName','VNLE','Marker','x','LineStyle','-','Color',cols(1+d,:)); -plot(fsym.*1e-9,ber_mlse,'DisplayName','MLSE','Marker','x','LineStyle','-','Color',cols(3+d,:)); -plot(fsym.*1e-9,ber_viterbi,'DisplayName','Viterbi','Marker','x','LineStyle','--','Color',cols(5+d,:)); -plot(fsym.*1e-9,ber_db_diff_precoded,'DisplayName','MLSE db diff','Marker','.','MarkerSize',15,'LineStyle','-'); -plot(fsym.*1e-9,ber_db,'DisplayName','MLSE db','Marker','.','MarkerSize',15,'LineStyle','-'); xlabel('Baudrate in GBd'); ylabel('BER'); set(gca, 'yscale', 'log'); -% ylim([1e-6 0.1]); -legend +set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); +grid on; +legend('Location','best'); +% xlim([184, 256]) +% ---------------- FIGURE 15 : Information Rates ---------------- +tp = TransmissionPerformance; + + +m = floor(log2(M)*10)/10; +figure(113+M); clf; hold on; + +netrates_vnle = tp.calculateNetRate(fsym.* m, ... + 'NGMI', gmi_vnle_bitwise./m, ... + 'BER', ber_vnle); +% +plot(xGHz, gmi_vnle_bitwise.*xGHz, ... + 'DisplayName','GMI*R VNLE', ... + mk.VNLE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.VNLE); + +plot(xGHz, netrates_vnle.SDHD.NetRate.*1e-9, ... + 'DisplayName','SD+HD VNLE', ... + mk.VNLE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.VNLE); +plot(xGHz, netrates_vnle.HD.NetRate.*1e-9, ... + 'DisplayName','Staircase VNLE', ... + mk.VNLE{:}, 'LineStyle','-.','LineWidth',LINE_WIDTH,'Color',cm.VNLE); + + +% +% MLSE symbol-wise (if present) +plot(xGHz, gmi_mlse.*xGHz, ... + 'DisplayName','GMI*R MLSE', ... + mk.MLSE{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.MLSE); + +netrates_mlse = tp.calculateNetRate(fsym.* m, ... + 'NGMI', gmi_mlse./m, ... + 'BER', ber_mlse); +plot(xGHz, netrates_mlse.SDHD.NetRate.*1e-9, ... + 'DisplayName','SD+HD MLSE', ... + mk.MLSE{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.MLSE); +plot(xGHz, netrates_mlse.HD.NetRate.*1e-9, ... + 'DisplayName','Staircase MLSE', ... + mk.MLSE{:}, 'LineStyle','-.','LineWidth',LINE_WIDTH,'Color',cm.MLSE); + + +% duobinary has only one GMI curve (DB output) +plot(xGHz, gmi_mlse_db.*xGHz, ... + 'DisplayName','GMI*R DB tgt.', ... + mk.DB{:}, 'LineStyle','-','LineWidth',LINE_WIDTH,'Color',cm.DB); + +netrates_db = tp.calculateNetRate(fsym.* m, ... + 'NGMI', gmi_mlse_db./m, ... + 'BER', ber_db); + +plot(xGHz, netrates_db.SDHD.NetRate.*1e-9, ... + 'DisplayName','SD+HD DB', ... + mk.DB{:}, 'LineStyle',':','LineWidth',LINE_WIDTH,'Color',cm.DB); +plot(xGHz, netrates_db.HD.NetRate.*1e-9, ... + 'DisplayName','Staircase DB', ... + mk.DB{:}, 'LineStyle','-.','LineWidth',LINE_WIDTH,'Color',cm.DB); + + + +% ylim([log2(M)-1, log2(M)]); +xlabel('Baudrate in GBd'); +ylabel('AIR in Gbps'); +set(gca, 'XTick', xticks_vals(1:2:end), 'XTickLabel', xtick_labels(1:2:end)); +grid on; +legend('Location','best'); +xlim([184, 256]) % Auxiliary nested helper for numerically stable log-sum-exp function s = logsumexp(a) diff --git a/projects/WDM/WDM_model.m b/projects/WDM/WDM_model.m new file mode 100644 index 0000000..347897a --- /dev/null +++ b/projects/WDM/WDM_model.m @@ -0,0 +1,268 @@ +%%% Run parameters +% TX + +M = 4; +m = floor(log2(M)*10)/10; +fsym = 224e9; + +apply_pulsef = 1; +fdac = 2*fsym; +fadc = 2*fsym; +random_key = 2; + +rcalpha = 0.05; +kover = 8; +vbias_rel = 0.5; +u_pi = 3.2; +vbias = -vbias_rel*u_pi; + +laser_linewidth = 0e6; + + +% Channel +link_length = 2; +rop = -5; + +vnle_order1 = 50; +vnle_order2 = 0; +vnle_order3 = 0; + +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; +dfe_order = [0 0 0]; + +alpha = 0; + +len_tr = 4096*2; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +% mu_dc = 0; + +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; + +dfe_ = sum(dfe_order)>0; + +doub_mode = db_mode.no_db; +cols = linspecer(6); + +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + +db_precode = 0; +db_encode = 0; +duob_mode = db_mode.no_db; +apply_pulsef = 0; + + +% AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig); + +%%%%% Low-pass el. components %%%%%% +% tx_bwl = 100e9; +% El_sig = Filter('filtdegree',3,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); + +wavelengthplan = calcWavelengthPlan(16,400e9,1310); +wavelengthplan = [1295,1305,1315,1325]; +N = numel(wavelengthplan); +f_plan = physconst('lightspeed')./(wavelengthplan.*1e-9); +margin = 25e12; % some THz left and right +f_span = (max(f_plan)+margin)-(min(f_plan)-margin); +f_nyq = f_span/2; + +upsample_required = f_nyq./(fdac*kover/2); +upsample_pow = 2^nextpow2(upsample_required); +upsample_ceil = ceil(upsample_required); + +f_opt = fdac*kover*upsample_pow; +f_opt_nyq = f_opt/2; + +signal_cell = {}; +Symbols = {}; +Tx_bits = {}; + +rop = linspace(-11,0,6); %12 workers when parallel + +num_realiz = 10; +gmi_vnle_bitwise = NaN(length(wavelengthplan),length(rop),num_realiz); +snr_vnle= NaN(length(wavelengthplan),length(rop),num_realiz); +ber_vnle= NaN(length(wavelengthplan),length(rop),num_realiz); + +for realiz = 1:num_realiz + + for l = 1:N + + [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource(... + "fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key+l+realiz,... + "db_precode",db_precode,"db_encode",db_encode,... + "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process(); + + % Digi_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0); + Lp_awg = Filter('filtdegree',3,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); + El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig); + % El_sig = M8199B("kover",kover).process(Digi_sig); + % El_sig.spectrum("fignum",101,"displayname",'bla','normalizeTo0dB',0,'lambda0_nm',1310,'useWavelengthAxis',0); + + %%%%% Electrical Driver Amplifier %%%%%% + El_sig = El_sig.normalize("mode","oneone"); + % El_sig = El_sig.setPower(1,"dBm"); + % figure;histogram(El_sig.signal); + + %%%%% MODULATE E/O CONVERSION %%%%% + Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",wavelengthplan(l),"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+l+realiz).process(El_sig); + + signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",30).process(Eml_out); + end + + Opt_sig_wdm = Optical_Multiplex("fs_in",signal_cell{l}.fs,"fs_out",upsample_pow*Eml_out.fs,... + "lambda_center",1310,"random_key",0,"filtype",1,"B",200e9).process(signal_cell); + + Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",3+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',1,'max_num_lines',2); + + %%%%%% Fiber %%%%%% + Opt_sig_wdm_fib=Opt_sig_wdm; + + nSegments = 2; + zdw = 1310; + D_local = 0; %if ~=0, simulation uses "segmented fiber with d+,d-) + randomize_D = true; + Dvec = getDispersionVector(nSegments, D_local, zdw, randomize_D, random_key+realiz); + for s = 1:nSegments + + Opt_sig_wdm_fib = DP_Fiber("L",link_length/nSegments,"D",Dvec(s),"Dpmd",0.1,"Ds",0.06,... + "beat_len",10,"corr_len",100,"dz",1,"manakov",0,... + "gamma",0.0023,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01,... + "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib); + + end + + Opt_sig_wdm_fib.spectrum("fignum",realiz,"displayname",'bla','lambda0_nm',1310,'useWavelengthAxis',0); + + % Opt_sig_wdm_fib.move_it_spectrum("fignum",100212,"displayname",'bla'); + + % Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig) + + parfor ri = 1:length(rop) + + %%%%%% ROP %%%%%% + Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop(ri)+10*log10(N)).process(Opt_sig_wdm_fib); + + Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1,"fs_out",Opt_sig_wdm_rx.fs/upsample_pow,"fs_in",Opt_sig_wdm_rx.fs,"lambda_center",1310).process(Opt_sig_wdm_rx); + + PD_cell = {}; + for l = 1:N + + %%%%%% PD Square Law %%%%%% + assert(fdac*kover==Opt_sig_wdm_demux{l}.fs,'Sampling Frequencies do not match! Check previous steps'); + PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11,"randomkey",random_key+l+realiz).process(Opt_sig_wdm_demux{l}); + + PD_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1); + + %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% + rx_bwl = 100e9; + PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(PD_sig); + + % %%%%%% Low-pass Scope %%%%%% + Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + + %%%%%% Scope %%%%%% + Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(PD_sig); + + Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym); + % Scpe_sig.spectrum("fignum",222,"displayname",'bla','normalizeTo0dB',1); + + [~, Scpe_cell, ~, found_sync] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 1); + Rx_sig = Scpe_cell{1}; + Rx_sig = Rx_sig.normalize("mode","rms"); + + % FFE or VNLE + eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.00,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + + [eq_signal_sd, eq_noise] = eq_.process(Rx_sig, Symbols{l}); + showEQNoisePSD(eq_noise, "fignum",1273876,"displayname",'noise after EQ'); + [mi_gomez] = calc_air(eq_signal_sd, Symbols{l}, "skip_front", 100, "skip_end", 100); + [gmi_vnle_bitwise(l,ri,realiz)] = calc_ngmi(eq_signal_sd,Symbols{l}); + % [gmi_bitwise_2] = calc_gmi_bitwise(eq_signal_sd,Symbols{l}); + snr_vnle(l,ri,realiz) = calc_snr(Symbols{l}, eq_signal_sd-Symbols{l}); + + % eq_signal_sd.plot("displayname",'bla','fignum',199); + % eq_signal_sd.eye(fsym,M,"fignum",103837); + + eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); + rx_bits = PAMmapper(M,0,"eth_style",0).demap(eq_signal_hd); + [~,tot_err,ber_vnle(l,ri,realiz),a] = calc_ber(rx_bits.signal,Tx_bits{l}.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + burst_vnle = count_error_bursts(a, 10)./tot_err; + + % showLevelConfusionMatrix(eq_signal_hd,Symbols{l},"M",M,"fignum",200,"displayname",'bla'); + % showLevelScatter(eq_signal_sd,Symbols{l},"displayname",'VNLE Out','f_sym',fsym,'fignum',201); + % show2Dconstellation(eq_signal_sd,Symbols{l},"displayname",'VNLE Out','fignum',2241); + + fprintf('CH %d :BER VNLE: %.2e \n',l,ber_vnle(l,ri,realiz)); + fprintf('CH %d :NGMI VNLE: %.2f \n',l,gmi_vnle_bitwise(l,ri,realiz)./m); + + end + + end + +end + + +figure();hold on; +cols = linspecer(N); +for l = 1:N + % plot(rop,mean(squeeze(ber_vnle(l,:,:)),2,'omitnan'),'Marker','*','DisplayName',sprintf('Ch: %d',wavelengthplan(l))) + plot(rop,squeeze(ber_vnle(l,:,:)),'Marker','*','DisplayName',sprintf('Ch: %d',wavelengthplan(l)),'Color',cols(l,:),'HandleVisibility','on') +end +yline([3.8e-3,2.2e-4],'HandleVisibility','off'); +ylabel('BER'); +xlabel('ROP') +title('BER vs. ROP'); +set(gca, 'XScale', 'linear', ... + 'YScale', 'log', ... + 'TickLabelInterpreter', 'latex', ... + 'FontSize', 11); + + +function dispersion_vector = getDispersionVector(N, D, ref_zdw, randomize_ZDW, randomkey) +% MATLAB version of the Python generator shown above. +% Returns an N×1 vector (ps/(nm·km)). +% +% D is the nominal dispersion magnitude. For D>0 the link is segmented with +% alternating sign (+D, -D, +D, …). For D==0 it is flat (0) except for +% ZDW randomization. The ZDW detuning is ~N(0, 2 nm) around 1310 nm and is +% converted to dispersion via 0.09 ps/(nm·km) per nm. + + % constants (matching the Python code) + meanLambda_nm = 1310; % center wavelength + sigma_nm = 2; % ZDW sigma + Dslope = 0.09; % ps/(nm·km) per nm detuning + + % random ZDW-induced dispersion offset + if randomize_ZDW + rng(randomkey, 'twister'); + rand_zdws_nm = meanLambda_nm + sigma_nm .* randn(N,1); + rand_D = (rand_zdws_nm - ref_zdw) .* Dslope; % ps/(nm·km) + else + rand_D = zeros(N,1); + end + + % nominal segmented pattern (match Python intent; keep length N) + if D > 0 + base = (-1) .^ ((0:N-1).'); % +1,-1,+1,-1,... + else % D == 0 (or anything else) + base = ones(N,1); + end + + dispersion_vector = base .* D + rand_D; % ps/(nm·km) +end diff --git a/test/bcjr_pam.m b/test/bcjr_pam.m new file mode 100644 index 0000000..674b55f --- /dev/null +++ b/test/bcjr_pam.m @@ -0,0 +1,555 @@ +classdef bcjr_pam < handle + %MLSE calculates the most probable sequence for an input signal with given/ known channel impulse response of any length + + properties(Access=public) + M %PAM-M + DIR + trellis_states + duobinary_output + end + + methods (Access=public) + + function obj = bcjr_pam(options) + %NAME Construct an instance of this class + % Detailed explanation goes here + + arguments + options.M double = 4; + options.DIR double = [1]; + options.trellis_states double = [-3 -1 1 3]; + options.duobinary_output logical = false; + + end + + % + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end + end + + function [VITERBI_ESTIMATION_SYMBOLS,LLR_exact,GMI] = process(obj,data_in,data_ref,tx_bits,bit_mapping) + + + debug = 0; + + % States should match the target states of the prev. EQ (EQ's job was to reduce the error between signal and the target) + trellis_state_mode = 2; + % 0 = use provided states (MUST provide the correct states); + % 1 = normalize to = 1 rms; + % 2 = use target symbols; + % 3 = use statistical levels + % 3 analyzes avg of rx signal levels - can help with nonlinear impairments + + trellis_exclusion = 1; % PAM-6 only (only if data is NOT precoded!) + + % Additional scaling between states, expected output (noiseless_received) and the noisy, filtered input signal + scale_mode = 2; % scale_mode: + % 0 = no scaling, + % 1 = use RMS to scale MODEL, + % 2 = use MMSE/time-corr to scale MODEL, -> This best to get the GMI right -> sometimes the LLP's are not centered around zero... + % 3 = use RMS to scale DATA, + % 4 = use MMSE/time-corr to scale DATA + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%% PREPARATIONS %%%%%%%% + + % remove unnecessary zeros at start of impulse response to keep + % number of trellis states minimal + DIR_nonzero = find(obj.DIR ~= 0); + if DIR_nonzero(1) > 1 + obj.DIR(1:DIR_nonzero(1)-1) = []; + end + + if isscalar(obj.DIR) + obj.DIR = [0 obj.DIR]; + end + + % impulse respnse to remove from signal + obj.DIR = flip(obj.DIR); %i.e. -0.2676 -0.0478 1.0000 + + % Trellis States + obj.trellis_states = reshape(obj.trellis_states,1,[]); + if trellis_state_mode == 1 % Normalize the Trellis states to =1 RMS + + obj.trellis_states = obj.trellis_states ./ rms(obj.trellis_states); + + elseif trellis_state_mode == 2 %simply use the states from the ref signal (should be a robust option) + + obj.trellis_states = reshape(unique(data_ref),size(obj.trellis_states)); + + elseif trellis_state_mode == 3 %use_statistical_levels + + %%%% Separate the equalized signal into the respective levels based on the actually transmitted level + constellation = unique(data_ref); + + % find actual levels from rx signal + symbols_for_lvl = NaN(numel(constellation),length(data_ref)); + for l = 1:numel(constellation) + level_amplitude = constellation(l); + symbols_for_lvl(l,data_ref==level_amplitude) = data_in(data_ref==level_amplitude); + end + + %replace the trellis states + avg_levels = mean(symbols_for_lvl,2,'omitnan'); + obj.trellis_states = sort(avg_levels)'; + + %also replace the whole ref signal (PAM-M) levels + [~, idx] = ismember(data_ref, unique(data_ref)); + data_ref = avg_levels(idx); + + end + + + % seems to be the only way to use combvec for a flexible amount + % of vectors. 'combs' contains all trellis states + pre_comb_mat = repmat(obj.trellis_states,length(obj.DIR)-1,1); + pre_comb_cell = mat2cell(pre_comb_mat,ones(1,size(pre_comb_mat,1)),size(pre_comb_mat,2)); + combs = fliplr(combvec(pre_comb_cell{:}).'); + first_sym = combs(:,1); % das ist das älteste/ trailing Symbol aus der sequenz + last_sym = combs(:,end); %hiermit wird entschieden/ das ist das cursor symbol am ende der sequenz + nStates = length(last_sym); + + % % Calculate all possible input symbols for the desired impulse + % % response. Row number is the index of the previous state, + % % column number is the index of the next state + % % noise free received == branch metrics + % assumes: last_sym = combs(:,end); % already defined earlier + levels = sort(unique(obj.trellis_states(:)).'); + edges = [levels(1) levels(end)]; % edge levels (0 and 5 in PAM6) + + noise_free_received = inf(nStates,nStates); % rows: to, cols: from + edge_edge_mask = false(nStates,nStates); % rows: to, cols: from + + for from = 1:nStates + for to = 1:nStates + % valid transition if shift-register overlap holds + if all(combs(to,2:end) == combs(from,1:end-1)) + % noiseless sample for the 'to' state reached from 'from' + noise_free_received(to,from) = ... + dot(combs(to,:), obj.DIR(end:-1:2)) + last_sym(from)*obj.DIR(1); + + % mark edge→edge candidate (to be excluded only on even→odd steps) + edge_edge_mask(to,from) = ... + (last_sym(from)==edges(1) || last_sym(from)==edges(2)) && ... + (last_sym(to) ==edges(1) || last_sym(to) ==edges(2)); + end + end + end + + h = flip(obj.DIR(:)).'; + data_in = data_in(:); + y_ideal = conv(data_ref(:), h, "same"); + + switch scale_mode + case 0 + g = 1; b = 0; + case 1 % RMS: scale model to data + g = rms(data_in)/rms(y_ideal); b = mean(data_in) - g*mean(y_ideal); + case 2 % MMSE/time-corr: scale states to data + [c,lags] = xcorr(data_in(:), y_ideal, 64); + [~,ix] = max(abs(c)); + lag = lags(ix); + y_ideal = circshift(y_ideal, lag); + mu_y = mean(data_in(:)); + mu_i = mean(y_ideal); + y_c = data_in(:)-mu_y; + yi_c = y_ideal-mu_i; + g = (yi_c'*y_c)/(yi_c'*yi_c); + b = mu_y - g*mu_i; + case 3 % RMS flipped: scale data to model + gd = rms(y_ideal)/rms(data_in); bd = mean(y_ideal) - gd*mean(data_in); + data_in = gd*data_in + bd; + g = 1; b = 0; + case 4 % MMSE/time-corr flipped: scale data to states + [c,lags] = xcorr(data_in(:), y_ideal(:), 64); + [~,ix] = max(abs(c)); + lag = lags(ix); + y_ideal = circshift(y_ideal(:), lag); + mu_y = mean(data_in(:)); + mu_i = mean(y_ideal); + y_c = data_in(:) - mu_y; % data_in centered + yi_c = y_ideal - mu_i; % ideal centered + g = (y_c' * yi_c) / (y_c' * y_c); + b = mu_i - g * mu_y; + data_in = g * data_in(:) + b; + g = 1; b = 0; + end + + % apply (g,b) to states/ expected values + noise_free_received = g*noise_free_received + b; + last_sym = g*last_sym + b; + + % calculate noise power + sigma2 = mean(abs(data_in - (g*y_ideal + b)).^2); %noise = mean(abs((RX Signal - IDEAL Signal)))^2 + inv2s2 = 1/(2*sigma2); + + if debug + figure(100); clf; hold on + obj.showLevelScatter_(data_in, data_ref); + yline(noise_free_received(:), 'DisplayName','Transition States','Color','red','HandleVisibility','off'); + yline(obj.trellis_states(:), 'DisplayName','Transition States','Color','green','LineWidth',2,'HandleVisibility','off') + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% FORWARD PASS (VITERBI -Alpha's) %%%%% + + % Initialize the output vector + pm = zeros(nStates,nStates); + bm_fw = zeros(nStates,nStates,length(data_in)); + + % first start is evaluated without ISI/ wihout the full Impulse response + % so simply use the constellation here + bm = -(data_in(1) - last_sym).^2 * inv2s2; + pm = pm + bm; + [alpha(:,1),pm_survivor_fw_idx(:,1)] = max(pm,[],2); + pm = repmat(alpha(:,1).',nStates,1); + bm_fw(:,:,1) = pm; + + % Forward Recursion (FSM Computation) + for n = 2:length(data_in) + + bm = -(data_in(n) - noise_free_received).^2 * inv2s2; + + % exclude edge to edge transitions only for even->odd steps && PAM-6 + if mod(n,2) == 0 && obj.M == 6 && trellis_exclusion + bm(edge_edge_mask) = -Inf; + end + + pm = pm + bm; + [alpha(:,n),pm_survivor_fw_idx(:,n)] = max(pm,[],2); % choose lowest path metric as new state (get min distance for all state transitions towards a new state) + pm = repmat(alpha(:,n).',nStates,1); % update pm (chosen state to 2nd dimension -> FROM state) + + bm_fw(:,:,n) = bm; + + end + + % we can now get the best path as min + viterbi_path = NaN(1,length(data_in)); + + % find ideal trellis path by going through the trellis backwards + [~,viterbi_path(length(data_in))] = max(alpha(:,length(data_in))); + for n = length(data_in):-1:2 + viterbi_path(n-1) = pm_survivor_fw_idx(viterbi_path(n),n); + end + + + if debug + alpha_ = alpha - min(alpha) + eps; + figure();hold on; + n = 10; + scatter(1:n,obj.trellis_states(repmat([1:numel(obj.trellis_states)]',1,n)),abs(alpha_(:,end-n+1:end)),'Marker','o','LineWidth',1); + scatter(1:n,obj.trellis_states(viterbi_path(end-n+1:end)),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','green'); + % scatter(1:n,data_ref(end-n+1:end),500,'Marker','x','LineWidth',1,'MarkerEdgeColor','red'); + yticks(obj.trellis_states); + ylim([min(obj.trellis_states)-1 max(obj.trellis_states)+1]); + end + + VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = first_sym(viterbi_path); + VITERBI_ESTIMATION_SYMBOLS = reshape(VITERBI_ESTIMATION_SYMBOLS,size(data_in)); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% BACKWARD (Beta's) %%%%% + + % Initialize the output vector + pm = zeros(nStates,nStates); + beta = zeros(nStates,length(data_in)); + pm_survivor_bw_idx = zeros(nStates,length(data_in)); + bm_bw = zeros(nStates,nStates,length(data_in)); + + % starting with the state that has the lowest sum path + % metric, follow the stored information about the + % predecessor + for h = length(data_in)-1:-1:1 + + bm = -(data_in(h+1) - noise_free_received).^2 * inv2s2; + + % exclude edge to edge transitions for even->odd steps && PAM-6 + if mod(h+1, 2) == 0 && obj.M == 6 && trellis_exclusion + bm(edge_edge_mask) = -Inf; + end + + pm = pm + bm.'; + [beta(:,h),pm_survivor_bw_idx(:,h)] = max(pm,[],2); % choose lowest path metric as new state + pm = repmat(beta(:,h).',nStates,1); % update pm (chosen state to 2nd dimension -> FROM state) + + bm_bw(:,:,h) = bm; + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% FORWARD (Combine Alpha and Beta to yield LLP's) %%%%% + + %calc the log probabilities (llp's) + + for k = 1:length(data_in) + + if k == 1 + + alpha_ = repmat(alpha(:,k)',[nStates,1])'; + beta_ = beta(:,k); + + LLP(:,k) = max(alpha_ + beta_,[],2); + + else + + alpha_ = repmat(alpha(:,k-1)',[nStates,1])'; + gamma_ = bm_fw(:,:,k)'; + beta_ = beta(:,k); + + LLP(:,k) = max(alpha_ + gamma_,[],1) + beta_'; + end + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% Calc LLR's %%%%% + + % These are interchangeable... + nml_LLP = LLP - max(LLP); %subtract highest value for better numerical stability, LLP's are not always close to zero + expLLP = exp(nml_LLP); + state_prob = expLLP ./ sum(expLLP); % sums to one (or numerically close to one) + + % compute symbol‐posteriors from LLP in the log‐domain: + amax = max(LLP,[],1); + logZ = amax + log(sum(exp(LLP - amax), 1)); + logPstate = LLP - logZ; % still in log‐domain + state_prob = exp(logPstate); % exact, sums to 1 + + if obj.M == 6 + + num_bits = 5; + + % all possible transitions (for now 36, including the "edges" + % of the QAM 32 constellation) + states = [-5 -3 -1 1 3 5]; + pam6transitions = combvec(states,states)'; % pam6transitions = + % [-5 -5; + % -3 -5; + % -1 -5; ... + + [~, idx_sym_1] = ismember(pam6transitions(:,1), states); + [~, idx_sym_2] = ismember(pam6transitions(:,2), states); + pam6ind = [idx_sym_1, idx_sym_2]; + + numPairs = floor(size(LLP,2)/2); + LLR_exact = zeros(numPairs,5); + LLR_maxlogmap = zeros(numPairs,5); + + for k = 1:numPairs + symbol1 = 2*k-1; + symbol2 = 2*k; + + LLP1 = LLP(:,symbol1); + LLP2 = LLP(:,symbol2); + prob1 = state_prob(:,symbol1); + prob2 = state_prob(:,symbol2); + + % All 36 Combinations: M = LLP Symbol 1 + LLP Symbol 2 + Mij = LLP1(pam6ind(:,1)) + LLP2(pam6ind(:,2)); + pij = prob1(pam6ind(:,1)) .* prob2(pam6ind(:,2)); + + % for each of the 5 bits sum exact-probs or max-log + for b = 1:num_bits + idx_sym_1 = bit_mapping(:,b)==1; + idx_bit_1 = bit_mapping(:,b)==0; + + % exact LLR from probabilities + P1 = sum(pij(idx_sym_1)); %prob that bit == 1 + P0 = sum(pij(idx_bit_1)); + LLR_exact(k,b) = log(P1./P0); %ratio by multiplication + + % max-log: + LLR_maxlogmap(k,b) = max( Mij(idx_sym_1) ) - max( Mij(idx_bit_1) ); % ratio by subtraction + end + end + + % GMI calc includes the Tx-bitstream + tx_bits_pam6_reshaped = reshape(tx_bits',5,[])'; % N x 5 + MI = zeros(1, num_bits); + for k = 1:num_bits + + idx_bit_1 = (tx_bits_pam6_reshaped(:,k) == 0); %wo sind die 1en + idx_sym_1 = (tx_bits_pam6_reshaped(:,k) == 1); %wo sind die 0en + + %LLR's for all actually transmitted ones or zeros + llr0 = LLR_exact(idx_bit_1,k); + llr1 = LLR_exact(idx_sym_1,k); + + % Calculate mutual information for bit position k + I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1 + I1 = mean(log2(1 + exp(-llr1))); % exp(-+LLR) = exp(negative) < 1 + MI(k) = 1 - 0.5 * (I0 + I1); + end + + GMI = sum(MI); % Total mutual information per symbol + GMI = GMI/2; % GMI per single symbol not per two symbols + + else + + % Number of symbols and bits per symbol + num_bits = log2(length(obj.trellis_states)); % 2 bits per symbol + + % bit_mapping = PAMmapper(length(obj.trellis_states),0,"eth_style",0).showBitMapping; + + % Initialize LLR storage + LLR_maxlogmap = zeros(length(data_in),num_bits); + LLR_exact = zeros(length(data_in),num_bits); + + % Compute bit-wise LLRs + for bit_idx = 1:num_bits + + % Find indices where bit is 0 and where it is 1 + idx_bit_0 = bit_mapping(:,bit_idx) == 0; + idx_bit_1 = bit_mapping(:,bit_idx) == 1; + + % Sum over log-probabilities + % Max-Log approximation uses the single max LLP value + % instead of sum over all LLP's + LLR_maxlogmap(:,bit_idx) = max(LLP(idx_bit_1,:), [], 1) - max(LLP(idx_bit_0,:), [], 1); + + % Sum probabilities over states for which the bit is 1 and 0, respectively. + P0 = sum(state_prob(idx_bit_0, :),1); + P1 = sum(state_prob(idx_bit_1, :),1); + LLR_exact(:,bit_idx) = log(P1./P0); % N x num_bits + + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% CALC NGMI %%%%% + + MI = zeros(1, num_bits); + for k = 1:num_bits + + idx_bit_0 = (tx_bits(:,k) == 0); %wo sind die 1en + idx_bit_1 = (tx_bits(:,k) == 1); %wo sind die 0en + + %LLR's for all actually transmitted ones or zeros + llr0 = LLR_exact(idx_bit_0,k); + llr1 = LLR_exact(idx_bit_1,k); + + % mutual information for bit position k + I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1 + I1 = mean(log2(1 + exp(-llr1))); % exp(-+LLR) = exp(negative) < 1 + MI(k) = 1 - 0.5 * (I0 + I1); % assumes equally distributed ones and zeros + end + + GMI = sum(MI); % Total bitwise mutual information + + end + + + if debug + %%% DEBUG PLOT LIKELIHOOD RATIOS %%% + figure(115);clf + subplot(2,1,1) + for bit = 1:num_bits + hold on; + histogram(LLR_exact(:,bit),1000,"DisplayName",sprintf('Actual LLR of Bit Pos %d',bit),'LineStyle','none','FaceAlpha',0.4); + end + legend + + subplot(2,1,2) + for bit = 1:num_bits + hold on; + histogram(LLR_maxlogmap(:,bit),1000,"DisplayName",sprintf('Max Log LLR of Bit Pos %d',bit),'LineStyle','none','FaceAlpha',0.4); + end + legend + + if obj.M == 6 + pairs = reshape(VITERBI_ESTIMATION_SYMBOLS,2,[]).'; + levels = sort(unique(VITERBI_ESTIMATION_SYMBOLS)); + isedge = ismember(pairs, [levels(1) levels(end)]); + isforbidden = sum(isedge,2)==2; + fprintf('Found %d forbidden transitions (even -> odd ; edge -> edge).\n', nnz(isforbidden)); + end + + end + + + end + + function [symbols_for_lvl,avg_for_lvl] = showLevelScatter_(~,eq_signal,ref_symbols) + + figure() + + rx_symbols = eq_signal; %./ rms(eq_signal); + correct_symbols = ref_symbols; + + % col = cbrewer2('Paired',numel(unique(correct_symbols))*2); + col = ... + [0.6510 0.8078 0.8902; ... + 0.1216 0.4706 0.7059; ... + 0.6980 0.8745 0.5412; ... + 0.2000 0.6275 0.1725; ... + 0.9843 0.6039 0.6000; ... + 0.8902 0.1020 0.1098; ... + 0.9922 0.7490 0.4353; ... + 1.0000 0.4980 0; ... + 0.7922 0.6980 0.8392; ... + 0.4157 0.2392 0.6039; ... + 1.0000 1.0000 0.6000; ... + 0.6941 0.3490 0.1569; ... + 0.6510 0.8078 0.8902; ... + 0.1216 0.4706 0.7059; ... + 0.6980 0.8745 0.5412; ... + 0.2000 0.6275 0.1725]; + ccnt = -1; + + levels = unique(correct_symbols); + symbols_for_lvl = NaN(numel(levels),length(correct_symbols)); + start = 1; + ende = length(correct_symbols); + + for l = 1:numel(levels) + ccnt = ccnt+2; + + level_amplitude = levels(l); + + symbols_for_lvl(l,correct_symbols==level_amplitude) = rx_symbols(correct_symbols==level_amplitude); + std_lvl(l) = std(symbols_for_lvl(l,:),'omitnan'); + xax = 1:length(correct_symbols); + + scatter(xax(start:ende),symbols_for_lvl(l,start:ende),10,'.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5,'MarkerEdgeColor',col(ccnt,:)); + hold on; + + + end + + std_lvl = round(std_lvl,2); + + ccnt = 0; + avg_for_lvl = NaN(numel(levels),length(correct_symbols)); + % Add the windowed/ smoothed curves + for l = 1:numel(levels) + ccnt = ccnt+2; + level_amplitude = levels(l); + + L = 500; + movmean = 1/L .* movsum(rx_symbols(correct_symbols==level_amplitude),[L/2,L/2], 'Endpoints', 'fill'); + + avg_for_lvl(l,correct_symbols==level_amplitude) = movmean; + + nanx = isnan(avg_for_lvl(l,:)); + t = 1:numel(avg_for_lvl(l,:)); + avg_for_lvl(l,nanx) = interp1(t(~nanx), avg_for_lvl(l,~nanx), t(nanx)); + + plot(xax(start:ende),avg_for_lvl(l,start:ende),'Color',col(ccnt,:)); + + hold on + end + + % yline(levels); + xlabel('Samples'); + ylabel('Amplitude'); + ylim([-3 3]); + + end + + + end +end diff --git a/test/duobinary_emulation_minimal_example.m b/test/duobinary_emulation_minimal_example.m index 3f808b3..4bb75ca 100644 --- a/test/duobinary_emulation_minimal_example.m +++ b/test/duobinary_emulation_minimal_example.m @@ -1,65 +1,49 @@ -useprbs = 1; + + M = 6; -randkey = 1; -datarate = 224e9; -fsym = round(datarate / log2(M)) ; +apply_precode = 1; -db_pre = 1; +bitpattern = []; +s = RandStream('twister','Seed',1); +for i = 1:log2(M) + N = 2^(17-1); %length of prbs + bitpattern(:,i) = randi(s,[0 1], N, 1); +end -Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",0.05); +if M == 6 + bitpattern = reshape(bitpattern',[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); +end -[d,Symbols,Bits] = PAMsource(... - "fsym",fsym,"M",M,"order",17,"useprbs",1,... - "fs_out",fsym,... - "applyclipping",0,"clipfactor",1.5,... - "applypulseform",0,"pulseformer",Pform,... - "randkey",1,... - "db_precode",db_pre,"db_encode",0,... - "mrds_code",0,"mrds_blocklength",512).process(); +bits = Informationsignal(bitpattern); -%%%CHANNEL +symbols = PAMmapper(M,0).map(bits); -% s = RandStream('twister','Seed',2); -% start = 10000; -% burstwidth = 100; -% d_burst = d; -% for pos = start:start+burstwidth -% lvls = 1.5 .* PAMmapper(M,0).levels / rms(PAMmapper(M,0).levels); -% d_burst.signal(pos) = d.signal(pos)+randn(s,1,1); -% end +bits_rx = PAMmapper(M,0).demap(symbols); +[~,~,ber_direct,~] = calc_ber(bits.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); -d_resample = d.resample("fs_out",2.*fsym); +if apply_precode + symbols_tx = Duobinary().precode(symbols); +else + symbols_tx = symbols; +end -eq_ffe = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",1024,"mu_dd",0.0004,"mu_tr",0,"order",25,"sps",2,"decide",1); -d_eq = eq_ffe.process(d_resample,Symbols); +show2Dconstellation(symbols_tx,symbols_tx,"displayname",'VNLE Out','fignum',2241); + -% s = RandStream('twister','Seed',2); -% start = 10000; -% burstwidth = 100; -% d_burst = d_eq; -% for pos = start:start+burstwidth -% lvls = 1.5 .* PAMmapper(M,0).levels / rms(PAMmapper(M,0).levels); -% d_burst.signal(pos) = d_eq.signal(pos)+randn(s,1,1); -% end -% -% d_burst = PAMmapper(M,0).decide_pamlevel(d_burst); - -if db_pre +if apply_precode % Entschiedene Symbole codieren: d_DB(n) = d(n) + d(n-1) (im Fall von PAM4 7 level [0 1 2 3 4 5 6]) - d_db = Duobinary().encode(d); + symbols_db = Duobinary().encode(symbols_tx); % Entschiedene codierte Symbole decodieren: d_dec(n) = d_DB(n) mod4 - d_dec = Duobinary().decode(d_db); + symbols_rx = Duobinary().decode(symbols_db); else - d_dec = d_burst; + symbols_rx = symbols_tx; end % Vergleichen von b(n) und d_dec(n) -Rx_bits = PAMmapper(M,0).demap(d_dec); - -Tx_bits = Bits; - -[~,error_num,ber,error_pos] = calc_ber(Tx_bits.signal,Rx_bits.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); +bits_rx = PAMmapper(M,0).demap(symbols_rx); +[~,~,ber,~] = calc_ber(bits.signal,bits_rx.signal,"skip_front",10,"skip_end",10,"returnErrorLocation",1); disp(['BER: ',sprintf('%.1E',ber),' - - PAM-',num2str(M)]); diff --git a/test/minimal_example_bcjr.m b/test/minimal_example_bcjr.m new file mode 100644 index 0000000..a692d31 --- /dev/null +++ b/test/minimal_example_bcjr.m @@ -0,0 +1,171 @@ + +M_format = [2,4,6,8]; + +for m = 1:length(M_format) + % --- Parameters --- + M = M_format(m); % PAM order (e.g., 2,4,8) + Nsym = 1e5; % number of symbols + h = [1, 0.5]; % Impulse response to remove + + b = log2(M); + if M == 6 b = 5; end + rng(1); + bits_tx = logical(randi([0 1], Nsym, b, 'uint8')); + + tx_symbols = pammap(bits_tx,M); + + if M == 6 + states = unique(tx_symbols); + pam6transitions = combvec(states',states')'; % pam6transitions = + bitmapping = pamdemap(reshape(pam6transitions',1,[])',M); + else + bitmapping = pamdemap(unique(tx_symbols),M); + end + + scaling = sqrt(sum(unique(tx_symbols).^2)/numel(unique(tx_symbols))); + tx_symbols = tx_symbols ./ scaling; + + % apply impulse response to signal + y_filt = filter(h, 1, tx_symbols); + + sir = 10:25; + for s = 1:length(sir) + + % apply noise + y = awgn(y_filt,sir(s),"measured",1); + + % apply bcjr + BCJR = bcjr_pam("DIR",h,"duobinary_output",0,"M",M,"trellis_states",unique(tx_symbols)); + [viterbi_estimate,LLR,GMI(m,s)] = BCJR.process(y,tx_symbols,bits_tx,bitmapping); + + % decode LLR's + bits_LLR = LLR > 0; + + % demap viterbi symbols sequence + rx_symbols = viterbi_estimate .* scaling; + bits_rx = pamdemap(rx_symbols,M); + + % BER calc + BER_vit(m,s) = nnz(bits_tx ~= bits_LLR) / numel(bits_tx); + fprintf('BER LLR = %.2e \n', BER_vit); + + BER_llr(m,s) = nnz(bits_tx ~= bits_rx) / numel(bits_tx); + fprintf('BER = %.2e \n', BER_llr); + end +end + +figure();hold on +for m = 1:length(M_format) + plot(sir,BER_llr(m,:),'DisplayName',sprintf('PAM %d',M_format(m))) + % plot(sir,BER_vit(m,:),'DisplayName',sprintf('PAM %d',M_format(m)),'LineStyle',':','LineWidth',0.1,'HandleVisibility','off'); +end +ylabel('BER'); +xlabel('SNR') +title('BER vs. SNR'); +set(gca, 'XScale', 'linear', ... + 'YScale', 'log', ... + 'TickLabelInterpreter', 'latex', ... + 'FontSize', 11); + + +figure();hold on +for m = 1:length(M_format) + plot(sir,GMI(m,:),'DisplayName',sprintf('GMI PAM %d',M_format(m))) +end +ylabel('GMI'); +xlabel('SNR') +title('GMI vs. SNR'); +set(gca, 'XScale', 'linear', ... + 'YScale', 'linear', ... + 'TickLabelInterpreter', 'latex', ... + 'FontSize', 11); + +function symbols = pammap(bits,M) +bits = logical(bits); +if M == 2 + symbols = bits; +elseif M == 4 + symbols= 2*bits(:,1) + (bits(:,1)==bits(:,2)); + symbols=2*symbols-3; + +elseif M == 6 + + m = 1; + + if size(bits,2)>size(bits,1) + bits = bits'; %vector aufrecht stellen + end + bits = reshape(bits',1,[])'; + thres = [-3 5;-1 5;-3 -5;-1 -5;-5 3;-5 1;-5 -3;-5 -1;-1 3;-1 1;-1 -3;-1 -1;-3 3;-3 1;-3 -3;-3 -1;3 5;1 5;3 -5;1 -5;5 3;5 1;5 -3;5 -1;1 3;1 1;1 -3;1 -1;3 3;3 1;3 -3;3 -1]; + % LUT based mapping + for k = 1:5:fix(length(bits)/5)*5 + symbols(m:m+1,1) = thres(bin2dec(int2str(bits(k:k+4)'))+1,:); + m = m+2; + end + +elseif M == 8 + x1 = bits(:,1); + x2 = (bits(:,1)==bits(:,3)); + x3 = x2~=bits(:,2); + + symbols = 4*x1 + 2*x2 + x3; + symbols=2*symbols-7; +end +end + +function bits = pamdemap(symbols,M) + +if M == 2 + thres=0; +elseif M == 4 + thres=[-2,0,2]; +elseif M == 6 + thres = [-3 5;-1 5;-3 -5;-1 -5;-5 3;-5 1;-5 -3;-5 -1;-1 3;-1 1;-1 -3;-1 -1;-3 3;-3 1;-3 -3;-3 -1;3 5;1 5;3 -5;1 -5;5 3;5 1;5 -3;5 -1;1 3;1 1;1 -3;1 -1;3 3;3 1;3 -3;3 -1]; +elseif M == 8 + thres=-6:2:6; +end + +if M ~= 6 + symbols = symbols'; + a = squeeze(repmat(real(symbols),[1 1 length(thres)])); %Eingangssignal in 3 spalten + b = squeeze(repmat(reshape(thres(:).',[1 1 length(thres)]),[1 length(symbols) 1])); %Threshold in 3 Spalten + comp_real = a > b; %check for each symbol/ sampling if it exeeds the obj.thresholdseshold 1, 2 or 3 + comp_real=repmat(real(symbols),[1 1 length(thres)]) > repmat(reshape(thres(:).',[1 1 length(thres)]),[1 length(symbols) 1]); + s1=size(comp_real,1); + s2=size(comp_real,2); +end + +if M == 2 + data_out=abs(comp_real(:,:,1)); +elseif M == 4 + data_out=[comp_real(:,:,2); ones(s1,s2) - comp_real(:,:,1) + comp_real(:,:,3)]; +elseif M == 6 + + if size(symbols,2) > 1 + symbols = symbols.'; + end + + if length(symbols)/2 ~= round(length(symbols)/2) + symbols = [symbols;0]; + end + + m = 1; + for n = 1:2:length(symbols) + dist = sqrt((symbols(n)-thres(:,1)).^2+(symbols(n+1)-thres(:,2)).^2); + [~,dd_idx] = min(dist); + % dec_out(n:n+1) = LUT(dd_idx,:); + data_out(m:m+4) = bitget(dd_idx-1,5:-1:1); + m = m+5; + end + + data_out = reshape(data_out',5,[]); + +elseif M == 8 + data_out=[comp_real(:,:,4); + comp_real(:,:,1)-comp_real(:,:,3)+comp_real(:,:,5)-comp_real(:,:,7); + 1-comp_real(:,:,2)+comp_real(:,:,6)]; +end + +bits = data_out'; + +end \ No newline at end of file diff --git a/test/pam_6_differential_code_understand.m b/test/pam_6_differential_code_understand.m new file mode 100644 index 0000000..3f13c5a --- /dev/null +++ b/test/pam_6_differential_code_understand.m @@ -0,0 +1,79 @@ +M = 6; +data = [1,2,3,4,5,6]; + +M = 6; + +bitpattern = []; +s = RandStream('twister','Seed',1); +for i = 1:log2(M) + N = 2^(12-1); %length of prbs + bitpattern(:,i) = randi(s,[0 1], N, 1); +end + +if M == 6 + bitpattern = reshape(bitpattern',[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); +end + +bits = Informationsignal(bitpattern); + +symbols = PAMmapper(M,0).map(bits); +symbols_tx_prec = Duobinary().precode(symbols); + +% all possible transitions (for now 36, including the "edges" +% of the QAM 32 constellation) +states = PAMmapper(6,0,"eth_style",0).levels; +pam6transitions = combvec(states,states)'; % pam6transitions = +% [-5 -5; +% -3 -5; +% -1 -5; ... +pam6transitions_serial = reshape(pam6transitions',[],1); + +data = pam6transitions_serial; +data = round(data); +b = min(data); +data = data - b; +data = data ./ 2; +% THIS WAS USED! +bk = zeros(size(data)); +for k = 2:numel(data) + bk(k) = mod(data(k)-bk(k-1),M); +end + + +%% State Analysis +x = bk;%symbols_tx_prec.signal; +levels = sort(unique(x)).'; % or provide known 1x6 level values + +[~,ix] = min(abs(x - levels),[],2); +x = levels(ix); % snapped/quantized + +%% TRANSITION COUNTS & PROBABILITIES +K = numel(levels); +% map to state indices 1..K +[tf, idx] = ismember(x, levels); +idx = idx(:); +from = idx(1:end-1); +to = idx(2:end); +from = idx(1:2:end); +to = idx(2:2:end); + +% counts C(from,to) +C = accumarray([from,to], 1, [K K], @sum, 0); +% row-stochastic transition matrix P(to|from) +rowSums = sum(C,2); +P = C ./ max(rowSums,1); + +%% 1) HEATMAP (which transitions are more probable?) +figure('Name','Transition Probabilities (to | from)'); +h = heatmap(levels, levels, P, 'Colormap', parula, 'ColorbarVisible','on'); +colormap(gca,[[1,1,1];flip(cbrewer2('Spectral',100))]);clim([0,ceil(max(P(:))*10)/10]); +h.XLabel = 'From state (level)'; +h.YLabel = 'To state (level)'; +h.Title = 'P(to | from)'; + +%% 2) WEIGHTED TRANSITION GRAPH +% Use dtmc if you have Econometrics Toolbox: +mc = dtmc(P, 'StateNames', string(levels)); +figure('Name','Markov Graph (dtmc)'); +gp = graphplot(mc, 'ColorEdges',true, 'LabelEdges',true); \ No newline at end of file diff --git a/test/pam_6_states_analysis.m b/test/pam_6_states_analysis.m new file mode 100644 index 0000000..a92d3ca --- /dev/null +++ b/test/pam_6_states_analysis.m @@ -0,0 +1,212 @@ + + +M = 6; + +bitpattern = []; +s = RandStream('twister','Seed',1); +for i = 1:log2(M) + N = 2^(17-1); %length of prbs + bitpattern(:,i) = randi(s,[0 1], N, 1); +end + +if M == 6 + bitpattern = reshape(bitpattern',[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); +end + +bits = Informationsignal(bitpattern); + +symbols = PAMmapper(M,0).map(bits); + +bits_rx = PAMmapper(M,0).demap(symbols); +[~,~,ber_direct,~] = calc_ber(bits.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); +assert(ber_direct==0,'Mapping is wrong'); + +nBursts = 0; +% No Precoding %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%SEND DIRECTLY +symbols_tx = symbols; + +symbols_rx = introduce_symbol_errors(symbols_tx, 1, 10, nBursts, 42); + +%RECEIVE BRANCH (do nothing special) +bits_rx = PAMmapper(M,0).demap(symbols_rx); + +[~,~,ber,errpos] = calc_ber(bits.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); +disp(['BER normal: - ',sprintf('%.1E',ber),' - - PAM-',num2str(M)]); +bursts_normal = count_error_bursts(errpos, 20); + + +% Precode Emulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%SEND DIRECTLY +symbols_tx = symbols; + +symbols_rx = introduce_symbol_errors(symbols_tx, 1, 10, nBursts, 42); + +%REFERENCE BRACH +symbols_db = Duobinary().encode(symbols_tx); +symbols_tx_emu = Duobinary().decode(symbols_db); +bits_tx_emu = PAMmapper(M,0).demap(symbols_tx_emu); + +% symbols_rx = introduce_symbol_errors(symbols_tx, 1, 10, 200, 42); + +%RECEIVE BRANCH +symbols_db = Duobinary().encode(symbols_rx); +symbols_rx_emu = Duobinary().decode(symbols_db); +bits_rx = PAMmapper(M,0).demap(symbols_rx_emu); + +[~,~,ber_precode_emulation,errpos_precode_emulation] = calc_ber(bits_tx_emu.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); +disp(['BER precode emulation: ',sprintf('%.1E',ber_precode_emulation),' - - PAM-',num2str(M)]); +bursts_precode_emulation = count_error_bursts(errpos_precode_emulation, 20); + + + + +% Precode at Tx %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%SEND PRECODED DATA +symbols_tx_prec = Duobinary().precode(symbols); + +symbols_rx_prec = introduce_symbol_errors(symbols_tx_prec, 1, 10, nBursts, 42); + +%RECEIVE BRANCH +symbols_db = Duobinary().encode(symbols_rx_prec); +symbols_rx_prec = Duobinary().decode(symbols_db); +bits_rx = PAMmapper(M,0).demap(symbols_rx_prec); + +[~,~,ber_precoded,errpos_precoded] = calc_ber(bits.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); +disp(['BER precoded: ',sprintf('%.1E',ber_precoded),' - - PAM-',num2str(M)]); +burst_precoded = count_error_bursts(errpos_precoded, 20); + + +% Precode at Tx but omit at Rx %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%SEND PRECODED DATA +symbols_tx_prec = Duobinary().precode(symbols); +bits_tx_prec = PAMmapper(M,0).demap(symbols_tx_prec); + +symbols_rx_omit = introduce_symbol_errors(symbols_tx_prec, 1, 10, nBursts, 42); + +%RECEIVE BRANCH +bits_rx = PAMmapper(M,0).demap(symbols_rx_omit); + +[~,~,ber_omit,errpos_omit] = calc_ber(bits_tx_prec.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); +disp(['BER (omit precode): ',sprintf('%.1E',ber_omit),' - - PAM-',num2str(M)]); +burst_omit = count_error_bursts(errpos_omit, 20); + +if 0 + cols = linspecer(8); + figure();hold on; + stem(1:20,bursts_normal,'LineWidth',2,'Color',cols(4,:),'Marker','_','DisplayName','w/o diff. precoder'); + stem(1:20,bursts_precode_emulation,'LineWidth',2,'Color',cols(3,:),'Marker','.','LineStyle','-','DisplayName','emulated precoder'); + stem(1:20,burst_precoded,'LineWidth',1,'Color',cols(6,:),'Marker','_','DisplayName','w/ diff. precoder'); + stem(1:20,burst_omit,'LineWidth',1,'Color',cols(5,:),'Marker','.','LineStyle',':','DisplayName','omit precoder'); + xlabel('Bit Error Burst Length') + ylabel('Occurence') + set(gca, 'yscale', 'log'); +end + + +%% State Analysis +signal_to_analyze = symbols_tx_emu; +x = signal_to_analyze.signal(:); +levels = sort(unique(x)).'; % or provide known 1x6 level values + +[~,ix] = min(abs(x - levels),[],2); +x = levels(ix); % snapped/quantized + +%% TRANSITION COUNTS & PROBABILITIES +K = numel(levels); +% map to state indices 1..K +[tf, idx] = ismember(x, levels); +idx = idx(:); +from = idx(1:end-1); +to = idx(2:end); +from = idx(1:2:end); +to = idx(2:2:end); + +% counts C(from,to) +C = accumarray([from,to], 1, [K K], @sum, 0); +% row-stochastic transition matrix P(to|from) +rowSums = sum(C,2); +P = C ./ max(rowSums,1); + +%% 1) HEATMAP (which transitions are more probable?) +figure('Name','Transition Probabilities (to | from)'); +h = heatmap(levels, levels, P, 'Colormap', parula, 'ColorbarVisible','on'); +colormap(gca,[[1,1,1];flip(cbrewer2('Spectral',100))]);clim([0,ceil(max(P(:))*10)/10]); +h.XLabel = 'From state (level)'; +h.YLabel = 'To state (level)'; +h.Title = 'P(to | from)'; + +%% 2) WEIGHTED TRANSITION GRAPH +% Use dtmc if you have Econometrics Toolbox: +mc = dtmc(P, 'StateNames', string(levels.*PAMmapper(M,0).get_scaling)); +figure('Name','Markov Graph (dtmc)'); +gp = graphplot(mc, 'ColorEdges',true, 'LabelEdges',true); + + +function symbols = introduce_symbol_errors(symbols, j, maxBurstLen, nBursts, seed) +%INTRODUCE_SYMBOL_ERRORS injects bursty level errors into symbols.signal. +% symbols.signal : column/row vector of quantized levels (exactly one of 6 values) +% j : max level step per sample (default 1) +% maxBurstLen : maximum burst length (default 8) +% nBursts : number of bursts to insert (default ~1% of length) +% seed : RNG seed (optional) + + if nargin < 2 || isempty(j), j = 1; end + if nargin < 3 || isempty(maxBurstLen), maxBurstLen = 8; end + x = symbols.signal(:); + N = numel(x); + if nargin < 4 || isempty(nBursts), nBursts = max(1, round(0.01*N)); end + if nargin >= 5 && ~isempty(seed), rng(seed); end + + % known levels and index mapping + lvls = sort(unique(x)).'; + K = numel(lvls); + + [~, idx] = ismember(x, lvls); % idx in 1..6 + + used = false(N,1); % avoid overlapping bursts + burst_ranges = zeros(nBursts,2); + + for b = 1:nBursts + % pick start not inside an existing burst + s = randi(N); + while used(s), s = randi(N); end + L = randi(maxBurstLen); + e = min(N, s+L-1); + + % mark used range + used(s:e) = true; + burst_ranges(b,:) = [s e]; + + % choose one direction for the whole burst: -1 (down) or +1 (up) + dir = randi([0 1])*2 - 1; + + % apply level errors within the burst + for t = s:e + k = idx(t); % current level index (1..6) + + % force inward movement at edges; prevents "flipping" to opposite edge + if k == 1 && dir == -1, dir = +1; end + if k == K && dir == +1, dir = -1; end + + step = randi([1 j]); % 1..j steps + kNew = k + dir*step; + + % clamp to [1,K], no wrap-around + if kNew < 1, kNew = 1; elseif kNew > K, kNew = K; end + + % if clamped to the same edge repeatedly, flip direction to keep changing + if kNew == k + dir = -dir; + kNew = max(1, min(K, k + dir*step)); + end + + idx(t) = kNew; + end + end + + x_err = lvls(idx); + symbols.signal = reshape(x_err, size(symbols.signal)); % preserve original shape + +end