diff --git a/Classes/04_DSP/Equalizer/ML_MLSE.m b/Classes/04_DSP/Equalizer/ML_MLSE.m index 3f584c3..168b837 100644 --- a/Classes/04_DSP/Equalizer/ML_MLSE.m +++ b/Classes/04_DSP/Equalizer/ML_MLSE.m @@ -1,3 +1,499 @@ +% classdef ML_MLSE < handle +% % ALGORITHM DESCRIBED IN: +% % W. Lanneer and Y. Lefevre, “Machine Learning-Based Pre-Equalizers for +% % Maximum Likelihood Sequence Estimation in High-Speed PONs,” +% % in 2023 31st European Signal Processing Conference +% +% % Further ML Refs: +% % https://machinelearningmastery.com/cross-entropy-for-machine-learning/ +% % https://docs.pytorch.org/docs/stable/generated/torch.nn.CrossEntropyLoss.html +% +% % The central idea is to overcome the (white-) noise assumption within the previously described +% % Viterbi algorithm, more precisely a closed-loop optimization is proposed that finds a suitable +% % filter-set to directly compute the branch metrics c_k (s,s^' ). These can directly be used to +% % carry out the conventional Viterbi algorithm. The system consists of S^L S=F linear FIR filters, +% % combined with one bias coefficient respectively. These filters take the received input samples to +% % compute the branch metrics estimates (c_k ) ̂(s,s^' ) according toThe central idea is to overcome +% % the (white-) noise assumption within the previously described Viterbi algorithm, more precisely +% % a closed-loop optimization is proposed that finds a suitable filter-set to directly compute the +% % branch metrics c_k (s,s^' ). These can directly be used to carry out the conventional Viterbi +% % algorithm. The system consists of S^L S=F linear FIR filters, combined with one bias coefficient +% % respectively. These filters take the received input samples to compute the branch metrics +% % estimates. Finally, the usual Viterbi is carried out... +% +% % Recommended Settings and some findings: +% +% % Requires many training epochs. According to ML people, 100,200 or +% % even up to 1000 epochs are normal for ML-convergence +% +% % The mu parameter _can_ be adaptive - using the cross entropy and when +% % analyzing the isolated training it looks very promisig. However, is +% % later use I found this is not as stable as a fixed learning rate. +% % mu = 0.1 worked good for me +% +% % Longer orders/ filter length are not always better. For me order=11 +% % was good. +% +% % Delay factor (delta) is good when the order is also increased. With +% % order = 11, a delta of =4 shows good results +% +% properties +% sps % usually 2 +% order +% e +% e_tr +% error +% +% len_tr +% mu_tr +% epochs_tr +% +% dd_mode % 1 or 0 to set DD-mode on or off +% mu_dd %weight update in dd mode +% epochs_dd +% +% adaptive_mu +% +% constellation +% +% L %viterbi memory length +% +% alpha +% DIR +% DIR_flip +% trellis_states +% +% traceback_depth +% +% % --- Added internal class variables used later --- +% S +% Nf +% delta +% nStates +% nFeasible +% combs +% first_sym +% last_sym +% valid +% valid_to_idx +% valid_from_idx +% w +% +% % --- New: fast state lookup --- +% true_to_state_idx +% state_dict % containers.Map: key(sequence)->state index +% key_fmt = '%.8g_'; % key format for sequence strings +% nSym % |constellation| +% +% ber = [] +% ce = ones(1,1); +% end +% +% methods +% function obj = ML_MLSE(options) +% arguments(Input) +% +% options.sps = 2; +% options.order = 15; +% +% options.len_tr = 4096; +% options.mu_tr = 0; +% options.epochs_tr = 5; +% +% options.dd_mode = 1; +% options.mu_dd = 1e-5; +% options.epochs_dd = 5; +% +% options.adaptive_mu = 1; +% +% options.delta = 0; +% options.traceback_depth = 1024; +% +% options.L = 1 +% +% end +% +% fn = fieldnames(options); +% for n = 1:numel(fn) +% obj.(fn{n}) = options.(fn{n}); +% end +% +% obj.e = zeros(obj.order,1); +% obj.error = 0; +% end +% +% function [X,X_viterbi] = process(obj, X, D) +% +% % actual processing of the signal (steps 1. - 3.) +% % 1 normalize RMS +% X = X.normalize("mode","rms"); +% +% % Use sorted constellation for deterministic mapping +% obj.constellation = sort(unique(D.signal),'ascend'); +% obj.nSym = numel(obj.constellation); +% +% if length(X)/length(D) ~= obj.sps +% warning('Signal length does not fit to reference!'); +% end +% +% % ============================================================== +% % INITIALIZATION (only before final epoch and detection mode) +% % ============================================================== +% +% % --- Parameters +% obj.S = numel(obj.constellation); % alphabet size +% obj.Nf = obj.order*obj.sps; % filter length +% % obj.delta = 3;%ceil(obj.Nf/2); % delay parameter +% obj.nStates = obj.S^obj.L; +% obj.nFeasible = obj.nStates*obj.S; +% +% % --- Trellis mapping +% obj.trellis_states = reshape(obj.constellation,1,[]); +% pre_comb_mat = repmat(obj.trellis_states, obj.L, 1); +% pre_comb_cell = mat2cell(pre_comb_mat, ones(1,obj.L), size(pre_comb_mat,2)); +% obj.combs = fliplr(combvec(pre_comb_cell{:}).'); % rows: states, columns: [x_k, x_{k-1}, ...] +% obj.first_sym = obj.combs(:,1); +% obj.last_sym = obj.combs(:,end); +% obj.nStates = size(obj.combs,1); +% +% % --- Valid transitions +% obj.valid = false(obj.nStates); +% for from = 1:obj.nStates +% for to = 1:obj.nStates +% if all(obj.combs(to,2:end) == obj.combs(from,1:end-1)) +% obj.valid(to,from) = true; +% end +% end +% end +% [obj.valid_to_idx, obj.valid_from_idx] = find(obj.valid); +% +% % --- Allocate vectors and weights +% % !! IF SHAPE FIT, then we already have smth there an we want +% % to start with the existing fitler-set +% if isempty(obj.w) || any(size(obj.w) ~= [obj.Nf+1,obj.nFeasible]) +% obj.w = zeros(obj.Nf+1,obj.nFeasible); % filter weights per transition + bias tap +% obj.w = randn(obj.Nf+1,obj.nFeasible); +% end +% +% % --- Precompute dictionary for fast state lookup (sequence -> state) +% keys = cell(obj.nStates,1); +% for i = 1:obj.nStates +% keys{i} = obj.seq_key(obj.combs(i,:)); % combs row is already [x_k, x_{k-1}, ...] +% end +% obj.state_dict = containers.Map(keys, 1:obj.nStates); +% +% % ============================================================== +% % TRAINING +% % ============================================================== +% +% % Training Mode +% n = obj.len_tr; +% training = 1; +% obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,n,training); +% obj.e_tr = obj.e; +% +% % ============================================================== +% % DD-Mode / Fixed Mode +% % ============================================================== +% +% % Decision Directed Mode +% n = X.length; +% training = 0; +% [y,y_vit]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training); +% +% X_viterbi = X; +% +% X.signal = y; +% X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym +% lbdesc = [num2str(obj.order),' tap FFE']; +% X = X.logbookentry(lbdesc); % append to logbook +% +% X_viterbi.signal = y_vit; +% X_viterbi.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym +% lbdesc = [num2str(obj.order),'order FFE + PF + Viterbi']; +% X_viterbi = X_viterbi.logbookentry(lbdesc); % append to logbook +% end +% +% function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training) +% % ============================================================== +% % FFE + Whitening + ML-Based Branch Metric Estimation + Viterbi +% % ============================================================== +% debug = 1; +% showPlots = 1; +% +% % --- Input padding and preallocation +% y = zeros(N,1); +% +% % number of symbol steps in this block +% nSymbols = ceil(N/obj.sps); +% +% for epoch = 1:epochs +% +% % state metrics (log-domain costs): keep as column [nStates×1] +% pm = zeros(obj.nStates,1); % v_{k-1}(s′) +% c_hat = zeros(1,obj.nFeasible); +% v_tilde = zeros(1,obj.nFeasible); +% pred = zeros(nSymbols, obj.nStates, 'uint32'); +% pm_sto = nan(obj.nStates, nSymbols,'like',pm); +% CE_accum = 0; +% +% +% %%% START IDX +% if training +% max_start = length(x) - ( (ceil(N/obj.sps)-1)*obj.sps + 1 ); +% max_start = max(1, max_start); % safety +% start_sample = randi([1, max_start], 1); %rnd training; not really good +% start_sample = 1; +% end_sample = start_sample + (ceil(N/obj.sps)-1)*obj.sps; +% else +% start_sample = 1;%obj.len_tr; +% end_sample = N; +% end +% +% start_symbol = 1 + floor((start_sample - 1)/obj.sps); % ABSOLUTE symbol index +% +% if numel(d) >= obj.L && start_symbol >= obj.L +% init_seq = d(start_symbol-obj.L+1 : start_symbol); % [d_k-L+1 ... d_k] +% true_to_state_idx = obj.state_dict(obj.seq_key(flip(init_seq))); % [d_k ... d_k-L+1] +% else +% % Not enough history – fall back to state 1 +% true_to_state_idx = uint32(1); +% end +% +% symbol = 0; +% for sample = start_sample:obj.sps:end_sample +% symbol = symbol + 1; +% k = symbol; +% sym_idx = start_symbol + (symbol - 1); +% +% % --- Build Δ-delayed observation window y_k +% i1 = sample - obj.Nf + 1 + obj.delta; +% i2 = sample + obj.delta; +% buf = x(max(1,i1):min(length(x),i2)); +% padL = max(0,1 - i1); +% padR = max(0,i2 - length(x)); +% yk = [zeros(padL,1); buf(:); zeros(padR,1)]; % Nf×1 +% yk = [yk;1]; +% +% % --- Predict branch metrics for all feasible transitions: c_hat +% c_hat = (yk.' * obj.w); % [1×nFeasible] +% c_hat = c_hat.'; % [nFeasible×1] +% +% % --- Extended path metrics: v_tilde = pm(from) + c_hat +% % normalize pm to avoid growth (invariant to additive const) +% pm = pm - min(pm); +% v_tilde = pm(obj.valid_from_idx) + c_hat; % [nFeasible×1] +% +% % ===== Gradient update (Algorithm 1) ===== +% +% if 1 %training +% % --- allocate storage once +% if epoch == 1 && symbol == 1 +% obj.true_to_state_idx = ones(ceil(N/obj.sps),1,'uint32'); +% end +% +% % --- previous "to" becomes current "from" +% if symbol > 1 +% true_from_state_idx = obj.true_to_state_idx(symbol-1); +% else +% true_from_state_idx = 1; +% end +% +% % --- compute or reuse "to" state +% if epoch == 1 +% % only compute in first epoch +% if sym_idx >= obj.L +% key_to = obj.seq_key(flip(d(sym_idx-obj.L+1 : sym_idx))); +% if isKey(obj.state_dict, key_to) +% obj.true_to_state_idx(symbol) = obj.state_dict(key_to); +% else +% obj.true_to_state_idx(symbol) = true_from_state_idx; +% end +% else +% obj.true_to_state_idx(symbol) = true_from_state_idx; +% end +% end +% +% % --- reuse cached state from second epoch onward +% true_to_state_idx = obj.true_to_state_idx(symbol); +% +% % --- ensure valid (from,to) +% dirac = zeros(obj.nFeasible,1); +% mask = obj.valid_from_idx==true_from_state_idx & ... +% obj.valid_to_idx ==true_to_state_idx; +% if any(mask) +% dirac(mask) = 1; +% else +% idx = find(obj.valid_from_idx==true_from_state_idx,1,'first'); +% dirac(idx) = 1; +% obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx); +% end +% +% +% +% +% % softmax over -v_tilde (numerically safe shift) +% v_shift = -(v_tilde - min(v_tilde)); % shift to small positive numbers +% v_shift = min(v_shift, 100); % clamp exponent argument (≈ exp(50)=3e21) +% expv = exp(v_shift); +% p = expv ./ (sum(expv) + eps); +% +% % for logging only: +% CE_symbol(symbol) = -log(p(dirac==1) + eps); +% +% if sym_idx > obj.L +% CE_smooth(symbol) = 0.01*CE_symbol(symbol) + 0.99*CE_smooth(symbol-1); +% else +% if epoch > 1 +% CE_smooth(symbol) = obj.ce(end); %use ce from last epoch or =1 for very first round?! +% else +% CE_smooth(symbol) = CE_symbol(symbol); +% end +% end +% +% CE_accum = CE_symbol(symbol) + CE_accum; +% +% +% % gradient term (t - p) +% dmp = (dirac - p)'; % 1×nFeasible +% +% % Per-feature gradient; implicit expansion gives (Nf+1)×nFeasible +% dL_Dw = (yk) .* dmp; +% +% % Start updates only when the ABSOLUTE symbol index has ≥ L history +% if sym_idx >= obj.L +% if obj.adaptive_mu +% mu_eff = CE_smooth(sym_idx); +% mu_eff = max(min(mu_eff, 0.2), 1e-4); +% else +% mu_eff = mu; +% end +% +% obj.w = obj.w - mu_eff .* dL_Dw; % (Nf+1)×nFeasible +% end +% +% % if debug && epoch > 2 +% % figure(100); +% % subplot(4,1,1); +% % heatmap(p'); +% % title('Probs') +% % subplot(4,1,2); +% % heatmap(dmp); +% % title('Update') +% % subplot(4,1,3); +% % heatmap(dL_Dw); +% % title('Update') +% % subplot(4,1,4); +% % heatmap(bj.w); +% % title('Update') +% % +% % end +% +% end +% +% +% +% % --- Compare-Select (matrix form, min of costs) +% v_tilde_mat = inf(obj.nStates, obj.nStates); +% v_tilde_mat(obj.valid) = v_tilde; +% [pm_next, pred(k,:)] = min(v_tilde_mat, [], 2); +% +% % re-center to keep metrics bounded (decision-invariant) +% pm_next = pm_next - min(pm_next); +% +% pm = pm_next; +% pm_sto(:,symbol) = pm; +% end +% +% % --- Traceback (full; you can window with traceback_depth if desired) +% [~, s_end] = min(pm); +% viterbi_path = zeros(symbol,1,'uint32'); +% viterbi_path(symbol) = s_end; +% for n = symbol:-1:2 +% viterbi_path(n-1) = pred(n, viterbi_path(n)); +% end +% +% y_ref = d(start_symbol:end); +% y = obj.first_sym(viterbi_path); +% +% if debug && training +% sym_start = start_symbol; +% sym_end = start_symbol + symbol - 1; +% ref_slice = d(sym_start : sym_end); +% err = sum(y ~= ref_slice(1:numel(y))); +% +% try +% ref_bits = PAMmapper(obj.S,0).demap(ref_slice); +% eq_bits = PAMmapper(obj.S,0).demap(y); +% [~, ~, ber, ~] = calc_ber(ref_bits, eq_bits, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); +% fprintf('Epoch: %d - BER: %.1e \n',epoch, ber); +% obj.ber(epoch) = ber; +% catch +% ser = err./length(y); +% fprintf('Epoch: %d - SER: %.1e \n',epoch, ser); +% end +% +% obj.ce(epoch) = CE_accum./symbol; +% +% if showPlots +% figure(10);clf +% subplot(3,2,1:2); +% heatmap(obj.w); +% title('Filter') +% +% subplot(3,2,3); +% v_tildemat = NaN(obj.nStates, obj.nStates); +% v_tildemat(obj.valid) = v_tilde; % log-domain scores +% heatmap(v_tildemat); +% title('Path Metrics (v_tilde)') +% +% subplot(3,2,4); +% scatter(1:symbol,pm_sto,1,'.') +% title('Path Metric Winners') +% +% subplot(3,2,5);hold on +% scatter(1:symbol,CE_symbol,1,'.'); +% scatter(1:symbol,CE_smooth,1,'.') +% title('Cross Entropy') +% +% subplot(3,2,6); hold on +% +% % Left y-axis: Cross Entropy (linear) +% yyaxis left +% scatter(1:length(obj.ce), obj.ce, 10, 's', 'filled') +% ylabel('Cross Entropy') +% +% % Right y-axis: BER (logarithmic) +% yyaxis right +% scatter(1:length(obj.ber), obj.ber, 10, 'd', 'filled') +% set(gca, 'YScale', 'log') +% ylabel('BER (log scale)') +% +% xlim([1, epochs]) +% xlabel('Epoch') +% title('Cross Entropy // BER') +% grid on +% +% drawnow +% end +% end +% end +% end +% end +% +% methods (Access=private) +% function k = seq_key(obj, seq) +% % Build a stable key string for a sequence row vector in the *same order as combs rows* ([x_k, x_{k-1}, ...]) +% % Use rounding via sprintf to avoid floating-point issues. +% % seq must be a row vector. +% k = sprintf(obj.key_fmt, seq); +% end +% end +% end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + classdef ML_MLSE < handle % --------------------------------------------------------------------- % W. Lanneer and Y. Lefevre, @@ -101,7 +597,7 @@ classdef ML_MLSE < handle obj.S = obj.nSym; obj.Nf = obj.order * obj.sps; obj.nStates = obj.S^obj.L; - obj.nFeasible = obj.nStates * obj.S; + obj.nFeasible = obj.nStates * obj.S; %feasible state transitions % --- Trellis mapping obj.trellis_states = reshape(obj.constellation,1,[]); @@ -164,8 +660,8 @@ classdef ML_MLSE < handle % EQUALIZE % ============================================================== function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training) - debug = 0; - showPlots = 0; + debug = 1; + showPlots = 1; y = zeros(N,1); nSymbols = ceil(N/obj.sps); @@ -241,6 +737,19 @@ classdef ML_MLSE < handle dirac(trans_idx)=1; end + % --- ensure valid (from,to) + if ~any(dirac) + mask = obj.valid_from_idx==true_from_state_idx & ... + obj.valid_to_idx ==true_to_state_idx; + if any(mask) + dirac(mask) = 1; + else + idx = find(obj.valid_from_idx==true_from_state_idx,1,'first'); + dirac(idx) = 1; + obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx); + end + end + % =================================================================== % TRAINING MODE (weight update) % =================================================================== diff --git a/Functions/Theory/Dissertation/mach_zehnder_modulator.m b/Functions/Theory/Dissertation/mach_zehnder_modulator.m new file mode 100644 index 0000000..7b379c2 --- /dev/null +++ b/Functions/Theory/Dissertation/mach_zehnder_modulator.m @@ -0,0 +1,140 @@ +% Minimal MZM transfer-function demo (sinusoidal drive) — aligned with your notation +% +% Implements exactly: +% E_out(t) = E0 * exp(j*w0*t) * exp(-j*w0*L*n_eff/c0) * 1/2 * [ exp(-j*phi1(t)) + rho*exp(-j*phi2(t)) ] +% with phi_{1,2}(t) = pi * v_{1,2}(t)/Vpi +% +% Push-pull: +% v1(t) = +v_drive(t)/2 , v2(t) = -v_drive(t)/2 => phi1 = +pi/2 * v_drive/Vpi, phi2 = -pi/2 * v_drive/Vpi +% +% And the ideal TF (rho=1): +% E_out/E_in = exp(-j*w0*L*n_eff/c0) * cos( (pi/2) * v_drive/Vpi ) +% +% Note: E_in(t) = E0 * exp(j*w0*t) in this script. + +% Parameters +c0 = physconst('lightspeed'); % [m/s] +lambda0 = 1310e-9; % [m] +omega0 = 2*pi*c0/lambda0; + +L = 5e-3; % [m] effective phase section length (set as needed) +n_eff = 2.2; % [-] effective index (set as needed) + +E0 = 1; % field amplitude (arbitrary) +Vpi = 3.2; % [V] half-wave voltage (your V_pi) + +% Drive +f0 = 1e9; % [Hz] +fs = 100e9; % [Hz] +Nper = 1; % number of periods +Vpp = 0.5*Vpi; % [V] peak-to-peak of v_drive(t) + +biasV = 2; % [V] differential bias added to v_drive + +% Analytic +v_ = linspace(-1,2, 2001); +% Field transfer function (amplitude) +Field_mzm_analytic = cos((pi/2)*v_); +% Power transfer function (intensity) +P_mzm_analytic = Field_mzm_analytic.^2; + +% Imbalance factor in YOUR notation: +rho = 1; % rho=1 -> ideal balanced MZM (collapses to ideal TF) + +% Time axis + differential drive voltage v_drive(t) +T = Nper/f0; +t = (0:1/fs:T-1/fs).'; + +v_drive = biasV + (Vpp/2)*sin(2*pi*f0*t); % v_drive(t) (peak = Vpp/2) + +% Push-pull branch voltages (consistent with v_drive = v1 - v2) +v1 = +0.5*v_drive; % arm 1 +v2 = -0.5*v_drive; % arm 2 + +% Phases phi1, phi2 +phi1 = pi * v1 / Vpi; +phi2 = pi * v2 / Vpi; + +% Fields: E_in and E_out (exactly your Eq. (mzm_e_field)) +E_in = E0 .* exp(1i*omega0*t); + +common_phase = exp(-1i * (omega0*L*n_eff/c0)); % exp(-j*omega0*L*n_eff/c0) + +E_out = E0 .* exp(1i*omega0*t) .* common_phase .* 0.5 .* ... + ( exp(-1i*phi1) + rho .* exp(-1i*phi2) ); + +% Transfer function (numerical): E_out/E_in +H_num = E_out ./ E_in; + +% Power (normalized) +Pnorm_num = abs(H_num).^2; % since |E_out/E_in|^2 + + + +% Ideal TF (analytic) for comparison (rho=1, push-pull) +H_ideal = common_phase .* cos( (pi/2) * (v_drive./Vpi) ); + +Pnorm_ideal = abs(H_ideal).^2; +Pnorm_math = cos( (pi/2) * (v_drive./Vpi) ).^2; + + + + + +set(groot, 'defaultLegendInterpreter', 'tex'); +set(groot, 'defaultAxesTickLabelInterpreter', 'tex'); +set(groot, 'defaultTextInterpreter', 'tex'); + +% Normalized voltage axis (multiples of Vpi) +v_norm = v_drive./Vpi; + +colfield = [0,0,0]; %is black +colpow = linspecer(2); +colpow = colpow(1,:); +colvdrive = linspecer(2); +colvdrive = colvdrive(2,:); + +%% SIGNAL IN +figure(1); clf +plot(v_norm,t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on; +ylabel('t [ns]'); xlabel('v_{drive}(t)/V_\pi'); +title('Drive voltage (normalized)'); +xlim([min(v_) max(v_)]); + +%% IN/OUT (static transfer) — normalized x-axis + analytic curve +figure(2); clf +plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF +plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color',colpow); hold on;% analytic power TF +% show input time signal +plot(v_norm,-1+t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on; +% show output time signal +plot(2+t*1e9, Pnorm_num, 'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colvdrive); hold on; +plot(2+t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on; +scatter(v_norm, Pnorm_num, 12, '.', 'LineWidth', 1,'MarkerEdgeColor',colvdrive); +scatter(biasV./Vpi,(cos((pi/2)*biasV./Vpi)^2),10,'Marker','o'); +line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--'); +line([max(v_drive) max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--'); +xline([min(v_norm) max(v_norm)]) + +grid on; +xlabel('v_{drive}(t)/V_\pi'); ylabel('|E_{out}/E_{in}|^2'); +% legend +xlim([min(v_) max(v_)+1]); +ylim([-1 1]); + +% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm.tex'); + + +%% +% % FIELD TF (only field here; do not mix power into this figure) +figure(3); clf +% plot(t*1e9, real(H_num), 'LineWidth', 1.0); hold on; +% plot(t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on; +plot(t*1e9, Pnorm_num, 'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colpow); hold on; +grid on; +xlabel('t [ns]'); ylabel('Re\{E_{out}/E_{in}\}'); +legend +mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm_out.tex'); + + + diff --git a/Functions/mat2tikz_improved.m b/Functions/mat2tikz_improved.m new file mode 100644 index 0000000..a36b1d1 --- /dev/null +++ b/Functions/mat2tikz_improved.m @@ -0,0 +1,23 @@ +function mat2tikz_improved(filename) +arguments + % Default to the path in your example if no argument is provided + filename (1,1) string = 'C:\Users\Silas\Documents\Dissertation\00_Examples\tikz\textfig.tikz'; +end +cleanfigure; +matlab2tikz(char(filename), ... + 'width','\fwidth', ... + 'height','\fheight', ... + 'showInfo',false, ... + 'extraAxisOptions',{ ... + 'legend style={font=\footnotesize}', ... + 'xlabel style={font=\color{white!15!black},font=\small},',... + 'ylabel style={font=\color{white!15!black},font=\small},',... + 'legend columns=1', ... + 'every axis/.append style={font=\scriptsize}',... + 'legend columns=1',... + 'legend style={at={(0.02,0.98)},font=\footnotesize,draw=black!60,rounded corners=2pt,inner sep=1pt,fill=white,column sep=6pt,anchor= north west}',... + 'legend style={at={(0.02,0.98)},draw=white!0!white,font=\scriptsize,inner sep=0.1pt,fill=white,column sep=1pt,anchor= north west}',... + 'every axis/.append style={font=\scriptsize}',... + }); + +end \ No newline at end of file diff --git a/Libs/wesanderson_colors/WesPalette.m b/Libs/wesanderson_colors/WesPalette.m index 7cd1dcb..13e976f 100644 --- a/Libs/wesanderson_colors/WesPalette.m +++ b/Libs/wesanderson_colors/WesPalette.m @@ -1,10 +1,13 @@ classdef WesPalette % WESPALETTE Wes Anderson color palettes with auto-completion % Usage: - % cmap = WesPalette.Zissou1.rgb() - % cmap = WesPalette.Zissou1.rgb(3) - - % https://github.com/karthik/wesanderson?tab=readme-ov-file + % cmap = WesPalette.Zissou1.rgb() % full palette + % cmap = WesPalette.Zissou1.rgb(3) % 3 colors (discrete default) + % cmap = WesPalette.Zissou1.rgb(12,"discrete") % any n, no interpolation + % cmap = WesPalette.Zissou1.rgb(256,"continuous") % smooth colormap (Lab interpolation) + % + % Requires: + % - colorspace.m (Pascal Getreuer) on MATLAB path for "continuous" mode enumeration BottleRocket1 @@ -34,21 +37,33 @@ classdef WesPalette end methods - function cmap = rgb(obj, n) + function cmap = rgb(obj, n, mode) % Return palette as Nx3 RGB colormap [0–1] + % + % n : number of requested colors (optional) + % mode : "discrete" (default) or "continuous" - hex = obj.hex(); + base_hex = obj.hex(); + base_rgb = WesPalette.hex2rgb(base_hex); - rgb = hex2rgb(hex); + if nargin < 2 || isempty(n) + cmap = base_rgb; + return; + end + if nargin < 3 || isempty(mode) + mode = "discrete"; + end + mode = lower(string(mode)); - if nargin == 2 - if n > size(rgb,1) - error('Requested %d colors, but only %d available.', ... - n, size(rgb,1)) - end - cmap = rgb(1:n,:); + validateattributes(n, {'numeric'}, {'scalar','integer','positive'}, mfilename, 'n'); + if mode ~= "discrete" && mode ~= "continuous" + error('mode must be "discrete" or "continuous".'); + end + + if mode == "discrete" + cmap = WesPalette.sample_discrete(base_rgb, n); else - cmap = rgb; + cmap = WesPalette.interpolate_continuous_lab(base_rgb, n); end end end @@ -56,78 +71,116 @@ classdef WesPalette methods (Access = private) function hex = hex(obj) % Internal HEX storage - switch obj case WesPalette.BottleRocket1 hex = {'#A42820','#5F5647','#9B110E','#3F5151','#4E2A1E','#550307','#0C1707'}; - case WesPalette.BottleRocket2 hex = {'#FAD510','#CB2314','#273046','#354823','#1E1E1E'}; - case {WesPalette.Rushmore1, WesPalette.Rushmore} hex = {'#E1BD6D','#EABE94','#0B775E','#35274A','#F2300F'}; - case WesPalette.Royal1 hex = {'#899DA4','#C93312','#FAEFD1','#DC863B'}; - case WesPalette.Royal2 hex = {'#9A8822','#F5CDB4','#F8AFA8','#FDDDA0','#74A089'}; - case WesPalette.Zissou1 hex = {'#3B9AB2','#78B7C5','#EBCC2A','#E1AF00','#F21A00'}; - case WesPalette.Zissou1Continuous hex = {'#3A9AB2','#6FB2C1','#91BAB6','#A5C2A3','#BDC881', ... '#DCCB4E','#E3B710','#E79805','#EC7A05','#EF5703','#F11B00'}; - case WesPalette.Darjeeling1 hex = {'#FF0000','#00A08A','#F2AD00','#F98400','#5BBCD6'}; - case WesPalette.Darjeeling2 hex = {'#ECCBAE','#046C9A','#D69C4E','#ABDDDE','#000000'}; - case WesPalette.Chevalier1 hex = {'#446455','#FDD262','#D3DDDC','#C7B19C'}; - case WesPalette.FantasticFox1 hex = {'#DD8D29','#E2D200','#46ACC8','#E58601','#B40F20'}; - case WesPalette.Moonrise1 hex = {'#F3DF6C','#CEAB07','#D5D5D3','#24281A'}; - case WesPalette.Moonrise2 hex = {'#798E87','#C27D38','#CCC591','#29211F'}; - case WesPalette.Moonrise3 hex = {'#85D4E3','#F4B5BD','#9C964A','#CDC08C','#FAD77B'}; - case WesPalette.Cavalcanti1 hex = {'#D8B70A','#02401B','#A2A475','#81A88D','#972D15'}; - case WesPalette.GrandBudapest1 hex = {'#F1BB7B','#FD6467','#5B1A18','#D67236'}; - case WesPalette.GrandBudapest2 hex = {'#E6A0C4','#C6CDF7','#D8A499','#7294D4'}; - case WesPalette.IsleofDogs1 hex = {'#9986A5','#79402E','#CCBA72','#0F0D0E','#D9D0D3','#8D8680'}; - case WesPalette.IsleofDogs2 hex = {'#EAD3BF','#AA9486','#B6854D','#39312F','#1C1718'}; - case WesPalette.FrenchDispatch hex = {'#90D4CC','#BD3027','#B0AFA2','#7FC0C6','#9D9C85'}; - case WesPalette.AsteroidCity1 hex = {'#0A9F9D','#CEB175','#E54E21','#6C8645','#C18748'}; - case WesPalette.AsteroidCity2 hex = {'#C52E19','#AC9765','#54D8B1','#B67C3B','#175149','#AF4E24'}; - case WesPalette.AsteroidCity3 hex = {'#FBA72A','#D3D4D8','#CB7A5C','#5785C1'}; end end end + + methods (Static, Access = private) + function rgb = hex2rgb(hex) + % hex: cellstr like {'#RRGGBB', ...} + if isstring(hex), hex = cellstr(hex); end + n = numel(hex); + rgb = zeros(n,3); + for i = 1:n + h = char(hex{i}); + if startsWith(h,'#'), h = h(2:end); end + if numel(h) ~= 6 + error('Invalid HEX color: %s', hex{i}); + end + rgb(i,1) = hex2dec(h(1:2))/255; + rgb(i,2) = hex2dec(h(3:4))/255; + rgb(i,3) = hex2dec(h(5:6))/255; + end + end + + function cmap = sample_discrete(base_rgb, n) + % No interpolation; allow any n by sampling/repeating. + k = size(base_rgb,1); + + if n <= k + idx = round(linspace(1, k, n)); % spread across palette + idx = max(1, min(k, idx)); + cmap = base_rgb(idx,:); + else + reps = floor(n / k); + rmd = mod(n, k); + cmap = [repmat(base_rgb, reps, 1); base_rgb(1:rmd,:)]; + end + end + + function cmap = interpolate_continuous_lab(base_rgb, n) + % Smooth interpolation in Lab using colorspace(). + % Requires colorspace.m by Pascal Getreuer on MATLAB path. + + k = size(base_rgb,1); + if k == 1 + cmap = repmat(base_rgb, n, 1); + return; + end + + % Convert to Lab, interpolate each channel, convert back + lab = colorspace('Lab<-RGB', base_rgb); + + t_base = linspace(0, 1, k); + t_new = linspace(0, 1, n); + + lab_new = zeros(n,3); + for c = 1:3 + lab_new(:,c) = interp1(t_base, lab(:,c), t_new, 'linear'); + end + + rgb_new = colorspace('RGB<-Lab', lab_new); + + % Clamp to displayable gamut + cmap = min(max(rgb_new, 0), 1); + end + end end diff --git a/Libs/wesanderson_colors/minimal_example_wespalette.m b/Libs/wesanderson_colors/minimal_example_wespalette.m index 9fcd99b..f308f27 100644 --- a/Libs/wesanderson_colors/minimal_example_wespalette.m +++ b/Libs/wesanderson_colors/minimal_example_wespalette.m @@ -5,8 +5,8 @@ y2 = 1e0 ./ (1 + exp(-0.4*(x-12))); % NLPN y3 = 1e-6 * 10.^(0.45*x); % RP on gamma y4 = 1e-2 * 10.^(0.18*(x-8)); % RP on beta2 -cmap = WesPalette.AsteroidCity1.rgb(4); -cmap = linspecer(4); +cmap = WesPalette.AsteroidCity1; +% cmap = linspecer(4); figure1=figure(202998);clf;hold on lw = 0.8; ms = 4; plot(x,y1,'LineWidth',lw,'Color',cmap(1,:),'Marker','o','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',[1,1,1],'MarkerSize',ms); diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m index 21f881f..986fcdb 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/Copy_of_FIGURE_WAVELENGTH.m @@ -6,7 +6,7 @@ db = DBHandler("dataBase", "labor_highspeed", "type", database_type); pam_levels = [4, 6, 8]; % three tiles bitrate_set = 360e9; -fiberL = 10; +fiberL = 2; fields = [ db.getTableFieldNames('power_state_info'); @@ -96,7 +96,7 @@ end %% ============================================================ % PLOT — 1×3 (PAM-4, PAM-6, PAM-8) % ============================================================ -fig = figure(9110); clf; +fig = figure(9112); clf; tiledlayout(1,3,'TileSpacing','compact','Padding','compact'); lw = 1.8; @@ -223,19 +223,19 @@ ylabel(''); end -pos = 1e3.*[2.7770 1.2017 1.4000 0.3200]; -set(fig, 'Position', pos); +% pos = 1e3.*[2.7770 1.2017 1.4000 0.3200]; +% set(fig, 'Position', pos); %% === EXPORT === -outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz'; -matlab2tikz(outfile, ... - 'width','\fwidth', ... - 'height','\fheight', ... - 'showInfo',false, ... - 'extraAxisOptions',{ ... - 'legend style={font=\footnotesize}', ... - 'legend columns=1' ... - }); +% outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz'; +% matlab2tikz(outfile, ... +% 'width','\fwidth', ... +% 'height','\fheight', ... +% 'showInfo',false, ... +% 'extraAxisOptions',{ ... +% 'legend style={font=\footnotesize}', ... +% 'legend columns=1' ... +% }); diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m index f9af8de..147c7d9 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/final/FIGURE_introduction.m @@ -7,113 +7,133 @@ data = readtable(tablename,"Delimiter",';','DecimalSeparator',','); % PLOT % ============================================================ %% 1. DATA EXTRACTION & SETUP -% Extract columns from the table +%% 1. DATA EXTRACTION & SETUP raw_M = data.M; raw_baud = data.BaudRate; raw_net = data.NetRate; -raw_names = string(data.ZoteroCode); -raw_band = string(data.Band); % Spalte 'Band' als String extrahieren +raw_codes = string(data.ZoteroCode); +raw_names = string(data.Name); +raw_band = string(data.Band); -% Filter: Keep only PAM 2, 4, 6, 8 and rows where Rates are not NaN +% Filter Valid Data target_M = [2, 4, 6, 8]; validIdx = ismember(raw_M, target_M) & ~isnan(raw_baud) & ~isnan(raw_net); -% Apply filter to create the working variables Mvals = raw_M(validIdx); baud = raw_baud(validIdx); netrate = raw_net(validIdx); +codes = raw_codes(validIdx); names = raw_names(validIdx); -bands = raw_band(validIdx); % Gefilterte Bands +bands = raw_band(validIdx); -% Setup Lists for Formatting pam_list = target_M; +colors = flip(cbrewer2('SET1',4)); -% Visual Settings -colors = lines(4); % 4 Farben für PAM-2,4,6,8 - -% 2. PLOT +%% 2. PLOT (For Visual Check only) figure; hold on; -ms = 40; % Marker size (etwas größer für bessere Sichtbarkeit) -lw = 1.5; % Line width für Fit +ms = 20; +lw = 0.5; -for k = 1:length(pam_list) % Loop durch PAM-Formate +for k = 1:length(pam_list) M = pam_list(k); - - % Logical mask for current PAM idxPam = (Mvals == M); - % Extract for this PAM x = baud(idxPam); y = netrate(idxPam); + b = bands(idxPam); n = names(idxPam); - b = bands(idxPam); % Bands für diese PAM-Gruppe - - % Get color for this PAM format col = colors(k,:); - % ----- LEGEND FLAG ----- firstLegend = true; - % ---- PLOT ALL POINTS (Marker basierend auf Band) ---- for i = 1:sum(idxPam) - - % === MARKER LOGIC === - currentBand = b(i); - if strcmpi(currentBand, 'O') - marker = 'o'; % Kreis für O-Band - elseif strcmpi(currentBand, 'C') - marker = 'd'; % Diamond für C-Band + + % Marker Logic + ms = 20; + if strcmpi(b(i), 'O') + marker = 'o'; + elseif strcmpi(b(i), 'C') + marker = 'd'; else - marker = 's'; % Fallback (Quadrat), falls andere Bands auftauchen + marker = 's'; + end + if strcmpi(n(i), 'THIS WORK') + marker = 'pentagram'; + ms = 100; end - % Plotting + % Plot Scatter if firstLegend - h = scatter(x(i), y(i), ms, ... - 'Marker', marker, ... - 'MarkerEdgeColor', col, ... - 'MarkerFaceColor', col, ... + scatter(x(i), y(i), ms, 'Marker', marker, ... + 'MarkerEdgeColor', col, 'MarkerFaceColor', col, ... 'DisplayName', sprintf('PAM-%d', M)); firstLegend = false; else - h = scatter(x(i), y(i), ms, ... - 'Marker', marker, ... - 'MarkerEdgeColor', col, ... - 'MarkerFaceColor', col, ... + scatter(x(i), y(i), ms, 'Marker', marker, ... + 'MarkerEdgeColor', col, 'MarkerFaceColor', col, ... 'HandleVisibility','off'); end - - % ====== CUSTOM DATATIP CONTENT ====== - dt = h.DataTipTemplate; - dt.DataTipRows(1).Label = 'Baud rate'; - dt.DataTipRows(2).Label = 'Net rate'; - dt.DataTipRows(end+1) = dataTipTextRow('Band', b(i)); % Zeigt Band an - dt.DataTipRows(end+1) = dataTipTextRow('Paper', n(i)); end - % ---- Fit (PAM-specific) ---- + % Fit lines if length(x) >= 3 [p, S, mu] = polyfit(x, y, 2); xfit = linspace(min(x), max(x), 200); yfit = polyval(p, xfit, S, mu); - - plot(xfit, yfit, ':', ... - 'LineWidth', lw, ... - 'Color', col, ... - 'HandleVisibility', 'off'); + plot(xfit, yfit, '-', 'LineWidth', lw, 'Color', col, 'HandleVisibility', 'off'); end end grid on; box on; xlabel('Baud rate [GBd]'); ylabel('Net rate [Gb/s]'); -title('Transmission Records (Circle=O, Diamond=C)'); -legend('Location','northwest'); -set(gca,'FontSize',11); +% title('Check Command Window for TikZ Code'); +% legend('Location','northwest'); +%% 3. GENERATE TIKZ ANNOTATION CODE +% This prints the manual \draw commands to the console + +%% GENERATE TIKZ ANNOTATION CODE +% This prints the manual \draw commands to the console + +%% GENERATE TIKZ ANNOTATION CODE (Colored Borders + Tiny Font) +%% GENERATE TIKZ ANNOTATION CODE (No Arrow, Close Text) +fprintf('\n\n%% ===========================================================\n'); +fprintf('%% COPY THE FOLLOWING LINES INTO YOUR .TEX FILE \n'); +fprintf('%% (Paste them just before \\end{axis})\n'); +fprintf('%% ===========================================================\n\n'); + +for i = 1:length(baud) + bx = baud(i); + by = netrate(i); + key = codes(i); + M_val = Mvals(i); + + % --- PLACEMENT LOGIC --- + if M_val == 8 + % PAM-8: Place Top-Left + % 'south east' anchor means the text's bottom-right corner touches the coordinate + % shift moves it slightly up and left to clear the marker + anchorStr = 'south east'; + shiftStr = 'shift={(-3pt, 3pt)}'; + else + % Others: Place Bottom-Right + % 'north west' anchor means the text's top-left corner touches the coordinate + % shift moves it slightly down and right + anchorStr = 'north west'; + shiftStr = 'shift={(3pt, -3pt)}'; + end + + % --- PRINT COMMAND --- + % Uses \node directly at the coordinate (axis cs:...) + fprintf('\\node[anchor=%s, %s, font=\\tiny, fill=white, inner sep=1pt] at (axis cs:%.2f, %.2f) {\\cite{%s}};\n', ... + anchorStr, shiftStr, bx, by, key); +end +fprintf('\n') %% === EXPORT === -outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\highspeedresults.tikz'; +outfile = 'C:\Users\Silas\Documents\latex\JLT_400G_submission\media\matlab2tikz\highspeedresults_test.tikz'; + matlab2tikz(outfile, ... 'width','\fwidth', ... 'height','\fheight', ... diff --git a/projects/ML_based_MLSE/analyze_filter_length.m b/projects/ML_based_MLSE/analyze_filter_length.m index bcf2a91..1b31d6b 100644 --- a/projects/ML_based_MLSE/analyze_filter_length.m +++ b/projects/ML_based_MLSE/analyze_filter_length.m @@ -25,10 +25,10 @@ h = abs([0.3 0.9 0.3]); h = h/norm(h); % h = [1 -1.67085330039878 1.17918163282514 -0.805210559745616 0.571564213123367 -0.296337147529674 0.00649773445209780 0.0854177610195952 -0.0576009020965258 0.0520994427061551 -0.0624586034913656 0.0553280962699552 -0.00705582559925755 -0.0336399056707792 0.0706903719452810 -0.0334124287931977 0.0131699455037966 0.0587431373842994 -0.0515902976066452 0.00647904355473619 0.0137506750904990 -0.0547974515885928 0.00994735499340592 -0.0135513582534086 -0.00463322575007739 0.0277311946101940]; % h = h/norm(h); -symbols_filt = Symbols.filter(1,h); +symbols_filt = Symbols.filter(h,1); symbols_noi = symbols_filt; symbols_noi.signal = awgn(symbols_filt.signal,SNR_dB,'measured'); -symbols_noi.spectrum; +symbols_noi.spectrum(); % --- Generate all parameter pairs [O,D] = ndgrid(order_range, delta_range); @@ -42,7 +42,7 @@ ce_vec = nan(size(pairs,1),1); ce_training = nan(size(pairs,1),training_len); % --- Parallel loop over parameter pairs -parfor k = 1:size(pairs,1) +for k = 1:size(pairs,1) order_k = pairs(k,1); delta_k = pairs(k,2); @@ -55,7 +55,7 @@ parfor k = 1:size(pairs,1) try ml = ML_MLSE("epochs_tr",training_len,"epochs_dd",1,"len_tr",2^15, ... "mu_dd",0.1,"mu_tr",0.1,"order",order_k,"sps",1, ... - "traceback_depth",128,"L",2,"delta",delta_k,"adaptive_mu",0); + "traceback_depth",128,"L",3,"delta",delta_k,"adaptive_mu",0); [y_ml,y_ref] = ml.process(symbols_noi,Symbols); ref_bits = PAMmapper(M,0).demap(y_ref);