From 1bf17ec6bee092ca4f64960b44beabefc1d77cf2 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Tue, 20 Jan 2026 13:48:15 +0100 Subject: [PATCH 1/4] FSO: - try to analyze spectrum of FSO and use the coeffs in the ML analysis. simple channel coeff + awgn example --- projects/FSO_transmission/first_analysis.m | 16 ++++++++++++++-- projects/ML_based_MLSE/analyze_filter_length.m | 18 +++++++++++------- 2 files changed, 25 insertions(+), 9 deletions(-) diff --git a/projects/FSO_transmission/first_analysis.m b/projects/FSO_transmission/first_analysis.m index d8361bc..9a327ab 100644 --- a/projects/FSO_transmission/first_analysis.m +++ b/projects/FSO_transmission/first_analysis.m @@ -83,8 +83,8 @@ Scope_sig = Electricalsignal(traceData.YData,"fs",fs); % pulse is symmetric, hence we can use pulsef directly as matched filter. % It feels off (bit I think correct) that the fsym is now the output freq.!! % -> output 2 sps to omit timing recovery!? -apply_matched_filter = 0; -k = 4; +apply_matched_filter = 1; +k = 1; if apply_matched_filter Pform = Pulseformer("fsym",fsym,"fdac",k*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1); @@ -96,6 +96,18 @@ else end Rx_matched.spectrum(); + +%% + +coefficients = arburg(Rx_matched.signal,25); + +figure() +[h,w] = freqz(1,coefficients,Rx_matched.length,"whole",Rx_matched.fs); +h = h/max(abs(h)); +hold on +w_ = (w - Rx_matched.fs/2); +plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(coefficients,2)), ' '],'LineWidth',2); + %% Timing Rec apply_timing_rec = 1; if apply_timing_rec diff --git a/projects/ML_based_MLSE/analyze_filter_length.m b/projects/ML_based_MLSE/analyze_filter_length.m index 616b715..bcf2a91 100644 --- a/projects/ML_based_MLSE/analyze_filter_length.m +++ b/projects/ML_based_MLSE/analyze_filter_length.m @@ -5,9 +5,9 @@ M = 4; randkey = 1; % --- Parameter sweep -order_range = 2:3:11; % FFE order -delta_range = 0:2:4; % delta -SNR_dB = 20; +order_range = 5:5:50; % FFE order +delta_range = 0:5:20; % delta +SNR_dB = 30; % --- Prepare bit sequence order_bits = 19; @@ -21,10 +21,14 @@ Symbols = PAMmapper(M,0).map(Bits); Symbols.fs = 200e9; % --- Channel (minimal ISI + AWGN) -h = [0.3 0.9 0.3]; h = h/norm(h); -symbols_filt = Symbols.filter(h,1); +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_noi = symbols_filt; symbols_noi.signal = awgn(symbols_filt.signal,SNR_dB,'measured'); +symbols_noi.spectrum; % --- Generate all parameter pairs [O,D] = ndgrid(order_range, delta_range); @@ -51,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",3,"delta",delta_k,"adaptive_mu",0); + "traceback_depth",128,"L",2,"delta",delta_k,"adaptive_mu",0); [y_ml,y_ref] = ml.process(symbols_noi,Symbols); ref_bits = PAMmapper(M,0).demap(y_ref); @@ -89,7 +93,7 @@ end beautifyBERplot ylabel('BER'); xlabel('Filter Order [N]'); title('BER vs. Filter order'); -ylim([1e-4, 0.1]); +% ylim([1e-4, 0.1]); yline(3.8e-3,'HandleVisibility','off'); yline(2.2e-4,'HandleVisibility','off'); From 7eb33648146b49212f6a8e0b67972ae7c96f2a37 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Thu, 29 Jan 2026 16:49:50 +0100 Subject: [PATCH 2/4] Minimal changes Add Theory plots for Silas Diss --- Classes/04_DSP/Equalizer/ML_MLSE.m | 515 +++++++++++++++++- .../Dissertation/mach_zehnder_modulator.m | 140 +++++ Functions/mat2tikz_improved.m | 23 + Libs/wesanderson_colors/WesPalette.m | 127 +++-- .../minimal_example_wespalette.m | 4 +- .../final/Copy_of_FIGURE_WAVELENGTH.m | 26 +- .../final/FIGURE_introduction.m | 132 +++-- .../ML_based_MLSE/analyze_filter_length.m | 8 +- 8 files changed, 860 insertions(+), 115 deletions(-) create mode 100644 Functions/Theory/Dissertation/mach_zehnder_modulator.m create mode 100644 Functions/mat2tikz_improved.m 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); From 5f3fd61519f9a866d11925c97dc69eee3910ed9d Mon Sep 17 00:00:00 2001 From: "silas (home)" Date: Fri, 30 Jan 2026 12:16:27 +0100 Subject: [PATCH 3/4] add few files from home pc --- projects/Messung_Zürich/minimal_matched.m | 60 +++++++++++++++++++++++ test/bitwise_demapping_pam6.m | 43 ++++++++++++++++ 2 files changed, 103 insertions(+) create mode 100644 projects/Messung_Zürich/minimal_matched.m create mode 100644 test/bitwise_demapping_pam6.m diff --git a/projects/Messung_Zürich/minimal_matched.m b/projects/Messung_Zürich/minimal_matched.m new file mode 100644 index 0000000..a5dfa33 --- /dev/null +++ b/projects/Messung_Zürich/minimal_matched.m @@ -0,0 +1,60 @@ +L = 4; % Oversampling factor +fsym = 1e9/4; +rollOff = 0.5; % Pulse shaping roll-off factor +htx = rcosdesign(rollOff,16, L, 'sqrt'); +hrx = conj(fliplr(htx)); + +signal = zeros(1000,1); +signal(500) = 1; +Digi_sig = Informationsignal(signal,"fs",fsym); + +% Digi_sig = Digi_sig.resample("fs_out",L*fsym); +% +% Digi_sig_tx = Digi_sig.filter(htx,1); +% Digi_sig_rx = Digi_sig_tx.filter(hrx,1); +% +% figure;plot(Digi_sig_tx.signal);hold on;plot(Digi_sig_rx.signal),plot(Digi_sig.signal); +% +% +% +% +pulsef = Pulseformer("alpha",1,"matched",0,"fdac",Digi_sig.fs*L,"fsym",fsym,"pulse","rrc","pulselength",16); +Digi_sig = pulsef.process(Digi_sig); + +pulsef = Pulseformer("alpha",1,"matched",1,"fdac",Digi_sig.fs,"fsym",fsym,"pulse","rrc","pulselength",16); +Digi_sig_matched = pulsef.process(Digi_sig); + +% Filter: +htx = rcosdesign(rollOff,16, L, 'sqrt'); +% Note half of the target delay is used, because when combined +% to the matched filter, the total delay will be achieved. +hrx = conj(fliplr(htx)); + +figure +plot(htx) +title('Transmit Filter') +xlabel('Index') +ylabel('Amplitude') + +figure +plot(hrx) +title('Rx Filter (Matched Filter)') +xlabel('Index') +ylabel('Amplitude') + +p = conv(htx,hrx); + +figure +plot(p) +title('Combined Tx-Rx = Raised Cosine') +xlabel('Index') +ylabel('Amplitude') + +% And let's highlight the zero-crossings +zeroCrossings = NaN*ones(size(p)); +zeroCrossings(1:L:end) = 0; +zeroCrossings((rcDelay)*L + 1) = NaN; % Except for the central index +hold on +plot(zeroCrossings, 'o') +legend('RC Pulse', 'Zero Crossings') +hold off \ No newline at end of file diff --git a/test/bitwise_demapping_pam6.m b/test/bitwise_demapping_pam6.m new file mode 100644 index 0000000..6ea0b3d --- /dev/null +++ b/test/bitwise_demapping_pam6.m @@ -0,0 +1,43 @@ + + +%%%%% SETTINGS %%%%%% +useprbs = 1; +M = 8; +randkey = 1; +fsym = 112e9; +viewresults = 0; + +%%%%% Mapping %%%%% +M = 6; +data = 0:M-1; +bitpersymbol = log2(M); + +s = RandStream('twister','Seed',1); +bitpattern = randi(s,[0 1], 2^18, 1); +bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); + +bits_tx = Informationsignal(bitpattern); + +symbols = PAMmapper(M,0).map(bits); + +pam6transitions = combvec(PAMmapper(M,0).levels,PAMmapper(M,0).levels)'; + +pam6bits = PAMmapper(6,0,"eth_style",0).demap(reshape(pam6transitions',[],1)./sqrt(10)); +symbols_rx = PAMmapper(M,0).map(pam6bits).*sqrt(10); + +pam6bits = reshape(pam6bits',5,[])'; + +figure; hold on +scatter(pam6transitions(:,1), pam6transitions(:,2), 'x', 'LineWidth', 1); +n = size(pam6transitions,1); +labels = cellstr(char(pam6bits + '0')); % -> N x 1 cell array of char rows +text(pam6transitions(:,1), pam6transitions(:,2), labels, ... + 'HorizontalAlignment','left', 'VerticalAlignment','bottom'); + + + +bits_rx = PAMmapper(M,0).demap(symbols); + +[~,error_num,ber,error_pos] = calc_ber(bits_tx.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); + +PAMmapper(8,0).showBitMapping \ No newline at end of file From fba7cbdba2438ed73749cb09b5f28eb5f279220c Mon Sep 17 00:00:00 2001 From: "silas (home)" Date: Sun, 1 Feb 2026 13:53:20 +0100 Subject: [PATCH 4/4] Faster DP_fiber with GPU processing... Run test\gpu_cpu_comparison.m to see difference on your setup --- Classes/02_optical/DP_Fiber.m | 19 +- Classes/02_optical/Optical_Demultiplex.m | 45 +-- Classes/02_optical/Optical_Multiplex.m | 45 ++- Classes/02_optical/dp_fiber_lib/CNLSE_plain.m | 189 ++++++---- Classes/02_optical/dp_fiber_lib/canUseGPU.m | 20 + .../02_optical/dp_fiber_lib/getNLstepsize.m | 40 +- Classes/02_optical/dp_fiber_lib/lin_step.m | 64 ++-- Classes/DataBaseHandler/DBHandler.m | 13 +- .../Dissertation/mach_zehnder_modulator.m | 153 ++++++-- .../Auswertung_JLT/run_dsp_from_db.m | 5 +- projects/WDM/WDM_model_10km_queue.m | 2 +- test/bayesopt_ffe_tuning.m | 219 +++++++++++ test/gpu_cpu_comparison.m | 349 ++++++++++++++++++ test/gpu_processing_dpfiber.m | 265 +++++++++++++ test/gpu_processing_test.m | 70 ++++ 15 files changed, 1281 insertions(+), 217 deletions(-) create mode 100644 Classes/02_optical/dp_fiber_lib/canUseGPU.m create mode 100644 test/bayesopt_ffe_tuning.m create mode 100644 test/gpu_cpu_comparison.m create mode 100644 test/gpu_processing_dpfiber.m create mode 100644 test/gpu_processing_test.m diff --git a/Classes/02_optical/DP_Fiber.m b/Classes/02_optical/DP_Fiber.m index a5f3dce..c1f08ea 100644 --- a/Classes/02_optical/DP_Fiber.m +++ b/Classes/02_optical/DP_Fiber.m @@ -24,6 +24,8 @@ classdef DP_Fiber SS_dzmax % [m] max dz (adaptive SSFM) SS_dzmin % [m] min dz (adaptive SSFM) n_waveplates % number of PMD waveplates + useGPU % GPU acceleration: true, false, or 'auto' (default) + useSingle % Use single precision on GPU (default: false) % ---- Internal state (persistent between calls) ---- state % struct mirroring legacy 'state' @@ -56,6 +58,8 @@ classdef DP_Fiber options.SS_dzmax = 2e4 % m options.SS_dzmin = 100 % m options.n_waveplates = 100 + options.useGPU = 'auto' % 'auto', true, or false + options.useSingle = false % single precision GPU end % Copy provided options into properties @@ -208,7 +212,7 @@ classdef DP_Fiber % 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; + st.brf.simdgd = 0; % cumsum used in legacy only for debug; keep compatibility variable: ~cumsum(st.brf.db0); % no-op to mirror legacy path @@ -228,7 +232,18 @@ classdef DP_Fiber x_in = signal_in(:,1).'; y_in = signal_in(:,2).'; - [x_out, y_out, obj.state] = CNLSE_plain(x_in, y_in, obj.state); + % Determine GPU usage + if ischar(obj.useGPU) || isstring(obj.useGPU) + if strcmpi(obj.useGPU, 'auto') + gpuFlag = []; % Let CNLSE_plain auto-detect + else + error('DP_Fiber:InvalidGPU', 'useGPU must be true, false, or ''auto'''); + end + else + gpuFlag = logical(obj.useGPU); + end + + [x_out, y_out, obj.state] = CNLSE_plain(x_in, y_in, obj.state, gpuFlag, obj.useSingle); obj.state.propagated_length = obj.state.propagated_length + obj.state.L; diff --git a/Classes/02_optical/Optical_Demultiplex.m b/Classes/02_optical/Optical_Demultiplex.m index 7d7c39e..8e272d2 100644 --- a/Classes/02_optical/Optical_Demultiplex.m +++ b/Classes/02_optical/Optical_Demultiplex.m @@ -42,7 +42,7 @@ classdef Optical_Demultiplex < handle function signalclasses_out = process(obj, signalclass_in) - % ---- Infer wavelength: either given or from input total signal + % ---- Infer wavelength: either given or from input total signal if isempty(obj.wavelengthplan) obj.wavelengthplan = signalclass_in.lambda; %meter else @@ -81,7 +81,7 @@ classdef Optical_Demultiplex < handle obj signal_in end - + w = obj.fs_out ./ obj.fs_in ; blocklen_in = length(signal_in); blocklen_out = w*blocklen_in; @@ -119,30 +119,31 @@ classdef Optical_Demultiplex < handle N = size(lo,1); C = size(lo,2); - x_envelopes = zeros(N, C, 'like', signal_in); - y_envelopes = zeros(N, C, 'like', signal_in); + % ---- VECTORIZED: Process all channels in parallel ---- + % Batched FFT operates on each column simultaneously on GPU - s1 = signal_in(:,1); - s2 = signal_in(:,2); + % Extract polarization signals + s1 = signal_in(:,1); % X polarization [N×1] + s2 = signal_in(:,2); % Y polarization [N×1] - % Reusable work buffers (avoid reallocations) - wrk_time = zeros(N,1, 'like', signal_in); - wrk_freq = zeros(N,1, 'like', signal_in); + % Broadcast signal to all channels and multiply with LO + % s1, s2 are [N×1], lo is [N×C] → result is [N×C] + x_mixed = att .* s1 .* lo; % [N×C] + y_mixed = att .* s2 .* lo; % [N×C] + + % Batched FFT: each column computed in parallel + x_freq = fft(x_mixed); % [N×C] + y_freq = fft(y_mixed); % [N×C] + + % Apply filter (H is [N×1], broadcasts across columns) + x_filtered = x_freq .* H; % [N×C] + y_filtered = y_freq .* H; % [N×C] + + % Batched IFFT + x_envelopes = ifft(x_filtered); % [N×C] + y_envelopes = ifft(y_filtered); % [N×C] - for c = 1:C - % ---- X branch ---- - wrk_time(:) = att .* s1 .* lo(:,c); % N×1 - wrk_freq(:) = fft(wrk_time); % N×1 - wrk_freq(:) = wrk_freq .* H; % N×1 - x_envelopes(:,c) = ifft(wrk_freq); % N×1 - % ---- Y branch ---- - wrk_time(:) = att .* s2 .* lo(:,c); - wrk_freq(:) = fft(wrk_time); - wrk_freq(:) = wrk_freq .* H; - y_envelopes(:,c) = ifft(wrk_freq); - end - end end end diff --git a/Classes/02_optical/Optical_Multiplex.m b/Classes/02_optical/Optical_Multiplex.m index 1d722b0..7fc8a33 100644 --- a/Classes/02_optical/Optical_Multiplex.m +++ b/Classes/02_optical/Optical_Multiplex.m @@ -1,10 +1,10 @@ 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 - + % 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... + % 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,... @@ -117,7 +117,7 @@ classdef Optical_Multiplex < handle % 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 @@ -139,25 +139,36 @@ classdef Optical_Multiplex < handle x_envelopes = NaN([blocklen_out N]); y_envelopes = x_envelopes; - + + % ---- OPTIMIZED: Pre-compute all LO phases as [blocklen_out × N] matrix ---- + time_idx = (0:blocklen_out-1).'; % [blocklen_out × 1] + lo_phases = mod(2*pi * time_idx * obj.df_T / obj.fs_out, 2*pi); % [blocklen_out × N] + lo_all = cos(lo_phases) + 1i*sin(lo_phases); % [blocklen_out × N] + + % Collect all resampled signals first (still requires loop due to cell array) + x_signals = zeros(blocklen_out, N); + y_signals = zeros(blocklen_out, N); + 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; - + data_in_resampled = data_in{o}.resample("fs_out", obj.fs_out); + x_signals(:, o) = data_in_resampled.signal(:, 1); + y_signals(:, o) = data_in_resampled.signal(:, 2); end + % ---- VECTORIZED: Batched FFT/IFFT for all channels ---- + % Apply filter to all channels at once + x_filtered = ifft(fft(x_signals) .* H); % [blocklen_out × N] + y_filtered = ifft(fft(y_signals) .* H); % [blocklen_out × N] + + % Apply attenuation and LO shift to all channels + x_envelopes = att .* x_filtered .* lo_all; % [blocklen_out × N] + y_envelopes = att .* y_filtered .* lo_all; % [blocklen_out × N] + 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 diff --git a/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m b/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m index cfa7294..463adca 100644 --- a/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m +++ b/Classes/02_optical/dp_fiber_lib/CNLSE_plain.m @@ -1,41 +1,72 @@ -function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state) - +function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state,useGPU,useSingle) +% GPU auto-detection if not specified +if nargin < 4 || isempty(useGPU) + useGPU = canUseGPU(); +end +if nargin < 5 || isempty(useSingle) + useSingle = false; +end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % pre calculations - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% pre calculations +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - state.common_beta=struct('X',0,'Y',0); +state.common_beta=struct('X',0,'Y',0); - for n=1:2 - % get current polarization name and contrary one - curPol = state.polNames{n}; +for n=1:2 + % get current polarization name and contrary one + curPol = state.polNames{n}; - % extend linear transfer function depending on beta values for the current polarization - % Was ist der Sinn dieser komischen beta notation? zB. state.beta.X = [0.3142 0 -9.1105e-28 5.1068e-41] - for n_beta = 1:length(state.beta.(curPol)) - state.common_beta.(curPol) = state.common_beta.(curPol) + state.beta.(curPol)(n_beta) * (state.omega).^(n_beta-1) / factorial(n_beta-1); - end - - %opt_out_struct.(curPol)=opt_in_struct.(curPol).envelope; + % extend linear transfer function depending on beta values for the current polarization + % Was ist der Sinn dieser komischen beta notation? zB. state.beta.X = [0.3142 0 -9.1105e-28 5.1068e-41] + for n_beta = 1:length(state.beta.(curPol)) + state.common_beta.(curPol) = state.common_beta.(curPol) + state.beta.(curPol)(n_beta) * (state.omega).^(n_beta-1) / factorial(n_beta-1); end - beta_const = state.beta.('X')(1); - beta_1 = state.beta.('X')(2); - beta_2 = state.beta.('X')(3); - beta_3 = state.beta.('X')(4); - deltaomega = state.omega; - beta_x = beta_const + beta_1 * deltaomega + 1/2 * beta_2 * deltaomega.^2 + 1/6 *beta_3 * deltaomega.^3; - -% opt_in_x = gpuArray(opt_in_x); -% opt_in_y = gpuArray(opt_in_y); + %opt_out_struct.(curPol)=opt_in_struct.(curPol).envelope; +end - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Split Step Method - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +beta_const = state.beta.('X')(1); +beta_1 = state.beta.('X')(2); +beta_2 = state.beta.('X')(3); +beta_3 = state.beta.('X')(4); +deltaomega = state.omega; +beta_x = beta_const + beta_1 * deltaomega + 1/2 * beta_2 * deltaomega.^2 + 1/6 *beta_3 * deltaomega.^3; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% GPU Transfer (if enabled) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if useGPU + % Convert to single precision if requested (faster on most GPUs) + if useSingle + opt_in_x = gpuArray(single(opt_in_x)); + opt_in_y = gpuArray(single(opt_in_y)); + state.common_beta.X = gpuArray(single(state.common_beta.X)); + state.common_beta.Y = gpuArray(single(state.common_beta.Y)); + state.brf.db1 = gpuArray(single(state.brf.db1)); + state.brf.db0 = gpuArray(single(state.brf.db0)); + for k = 1:numel(state.brf.matR) + state.brf.matR{k} = gpuArray(single(state.brf.matR{k})); + end + else + opt_in_x = gpuArray(opt_in_x); + opt_in_y = gpuArray(opt_in_y); + state.common_beta.X = gpuArray(state.common_beta.X); + state.common_beta.Y = gpuArray(state.common_beta.Y); + state.brf.db1 = gpuArray(state.brf.db1); + state.brf.db0 = gpuArray(state.brf.db0); + for k = 1:numel(state.brf.matR) + state.brf.matR{k} = gpuArray(state.brf.matR{k}); + end + end +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Split Step Method +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [opt_out_x,opt_out_y] = split_step_loop(state.L,opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin,... % state.lin_z_test,state.corr_length,state.n_plates_done,state.missing_dz,state.brf,state.common_beta,.... @@ -44,35 +75,33 @@ function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state) % [opt_out_x,opt_out_y] = split_step_loop_mex(state.L,opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin,... % state.lin_z_test,state.corr_length,state.n_plates_done,state.missing_dz,state.brf,state.common_beta,.... % state.chi,state.manakov,state.beat_len); - - % get nonlinear step size - [state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); - %[state.dz] = getNLstepsize_original(state,opt_out_struct); - - state.n_step = 0; - state.z_prop = 0; - state.test_dz = []; - state.powers = []; - tic +% get nonlinear step size +[state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); +%[state.dz] = getNLstepsize_original(state,opt_out_struct); - while state.z_prop < state.L +state.n_step = 0; +state.z_prop = 0; +state.test_dz = []; +state.powers = []; - % reduce step length (dz) if we are to overshoot the fiber length - % (L) in the next step - if state.z_prop + state.dz > state.L - state.dz = state.L - state.z_prop; - end - - % update step number (n) - state.n_step=state.n_step+1; +while state.z_prop < state.L - % append current step length to logbook (dzs) - state.dzs(state.n_step)=state.dz; + % reduce step length (dz) if we are to overshoot the fiber length + % (L) in the next step + if state.z_prop + state.dz > state.L + state.dz = state.L - state.z_prop; + end - - % half linear step - [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... + % update step number (n) + state.n_step=state.n_step+1; + + % append current step length to logbook (dzs) + state.dzs(state.n_step)=state.dz; + + + % half linear step + [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,... state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X]... @@ -82,12 +111,12 @@ function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state) state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X); - % complete nonlinear step - - [opt_in_x,opt_in_y] = nl_step(opt_in_x,opt_in_y, state.dz, state.gamma, state.chi, state.manakov, state.beat_len ,state.alpha_lin.X, state.alpha_lin.Y); - - % half linear step - [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... + % complete nonlinear step + + [opt_in_x,opt_in_y] = nl_step(opt_in_x,opt_in_y, state.dz, state.gamma, state.chi, state.manakov, state.beat_len ,state.alpha_lin.X, state.alpha_lin.Y); + + % half linear step + [opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,... state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,... state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X]... @@ -97,24 +126,34 @@ function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state) state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,... state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X); - % get nonlinear step size - [state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); - %[state.dz] = getNLstepsize_original(state,opt_out_struct); + % get nonlinear step size + [state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin); + %[state.dz] = getNLstepsize_original(state,opt_out_struct); - + +end + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% GPU Gather (if enabled) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if useGPU + opt_out_x = gather(opt_in_x); + opt_out_y = gather(opt_in_y); + + % Gather state arrays back to CPU + state.common_beta.X = gather(state.common_beta.X); + state.common_beta.Y = gather(state.common_beta.Y); + state.brf.db1 = gather(state.brf.db1); + state.brf.db0 = gather(state.brf.db0); + for k = 1:numel(state.brf.matR) + state.brf.matR{k} = gather(state.brf.matR{k}); end - - toc +else + opt_out_x = opt_in_x; + opt_out_y = opt_in_y; +end - - - - - opt_out_x = (opt_in_x); - opt_out_y = (opt_in_y); - -% opt_out_x = gather(opt_in_x); -% opt_out_y = gather(opt_in_y); - end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/canUseGPU.m b/Classes/02_optical/dp_fiber_lib/canUseGPU.m new file mode 100644 index 0000000..84666bb --- /dev/null +++ b/Classes/02_optical/dp_fiber_lib/canUseGPU.m @@ -0,0 +1,20 @@ +function canUse = canUseGPU() +%CANUSEGPU Check if a compatible GPU is available for parallel computing +% Returns true if MATLAB Parallel Computing Toolbox is available and +% a CUDA-capable GPU is detected. + +canUse = false; + +% Check if Parallel Computing Toolbox is installed +if ~license('test', 'Distrib_Computing_Toolbox') + return; +end + +% Check for GPU device +try + gpu = gpuDevice(); + canUse = gpu.DeviceSupported; +catch + canUse = false; +end +end diff --git a/Classes/02_optical/dp_fiber_lib/getNLstepsize.m b/Classes/02_optical/dp_fiber_lib/getNLstepsize.m index a70edb9..fd42eba 100644 --- a/Classes/02_optical/dp_fiber_lib/getNLstepsize.m +++ b/Classes/02_optical/dp_fiber_lib/getNLstepsize.m @@ -1,25 +1,25 @@ function [rDZ] = getNLstepsize(ux,uy,gamma,dzmin,dzmax,dphimax,alpha_lin) - - maxPow = max(gamma.*max(real(ux).^2+imag(ux).^2+real(uy).^2+imag(uy).^2)); - - Leff = dphimax/maxPow; - alpha_lin = max([alpha_lin.X alpha_lin.Y]); - nl_att_len_ratio = alpha_lin*Leff; - - if nl_att_len_ratio >= 1 - rDZ = dzmax; +% gather() handles gpuArray inputs - ensures scalar is on CPU for comparisons +maxPow = gather(max(gamma.*max(real(ux).^2+imag(ux).^2+real(uy).^2+imag(uy).^2))); + +Leff = dphimax/maxPow; +alpha_lin = max([alpha_lin.X alpha_lin.Y]); +nl_att_len_ratio = alpha_lin*Leff; + +if nl_att_len_ratio >= 1 + rDZ = dzmax; +else + if alpha_lin == 0 + step = Leff; else - if alpha_lin == 0 - step = Leff; - else - %effective length? - step = -1/alpha_lin*log(1-nl_att_len_ratio); - end - - rDZ = min([step dzmax]); - rDZ = max([rDZ dzmin]); - end - + %effective length? + step = -1/alpha_lin*log(1-nl_att_len_ratio); + end + + rDZ = min([step dzmax]); + rDZ = max([rDZ dzmin]); +end + end \ No newline at end of file diff --git a/Classes/02_optical/dp_fiber_lib/lin_step.m b/Classes/02_optical/dp_fiber_lib/lin_step.m index eea8a1e..cdb3562 100644 --- a/Classes/02_optical/dp_fiber_lib/lin_step.m +++ b/Classes/02_optical/dp_fiber_lib/lin_step.m @@ -22,37 +22,37 @@ n_plates_left = n_plates - n_plates_done; % compute last plate size ( if it fits, it should be 0) if missing_dz > aStepSize - + last_plate = aStepSize; missing_dz = missing_dz-aStepSize; plate_sizes = last_plate; plate_numbers = n_plates; - + else - + last_plate = aStepSize - missing_dz - (n_plates_left-1)*corr_length; - + if missing_dz == 0 missing_dz = []; end - + %build vector of plate lengths with missing plate part from prev. %iterartion , then some normal plates and finally a fraction of a plate %to fit into the step length plate_sizes = [missing_dz corr_length*ones(1,n_plates_left-1) last_plate]; - + if n_plates_done == 0 plate_numbers =[(n_plates_done+1):(n_plates-1) n_plates]; else plate_numbers = [n_plates_done (n_plates_done+1):(n_plates-1) n_plates]; % not wrking yet end - - %remember for next step + + %remember for next step missing_dz = corr_length - last_plate; - + end -plate_steps = repmat(n_step,1,length(plate_sizes)); +plate_steps = repmat(n_step,1,length(plate_sizes)); %#ok % %figure;stem(plate_sizes); @@ -66,50 +66,39 @@ test_plate_numbers = [test_plate_numbers, plate_numbers]; opt_x=fft(opt_x); opt_y=fft(opt_y); -% db1 = gpuArray(brf.db1); -% db0 = gpuArray(brf.db0); -% common_beta_x = gpuArray(common_beta_x); -% common_beta_y = gpuArray(common_beta_y); - -db1 = (brf.db1); -db0 = (brf.db0); -common_beta_x = (common_beta_x); -common_beta_y = (common_beta_y); +% Note: db1, db0, common_beta_x, common_beta_y are already on GPU when +% useGPU=true (transferred in CNLSE_plain) +db1 = brf.db1; +db0 = brf.db0; % process every waveplate with given sizes in plate_sizes for n=1:length(plate_sizes) dz = plate_sizes(n); - - % figure(87);subplot(2,1,1);plot(real(x(900:1150)));subplot(2,1,2);plot(real(y(900:1150))); - % MOV1=[MOV1 getframe(87)]; - + % extract rotation matrix from pre calculated matrices matR = brf.matR{plate_numbers(n)}; - + % transform to eigenvalue of of fiber segment tOpt.X = conj(matR(1,1))*opt_x + conj(matR(2,1))*opt_y; tOpt.Y = conj(matR(1,2))*opt_x + conj(matR(2,2))*opt_y; - + % calculate statistical delta beta for pmd delta_beta = 0.5*(db1+db0(n))/corr_length; - % build transfer function with delta beta - %common.beta = beta1+beta2*omega^2 - + %accumulate delta beta for log... brf.simdgd = brf.simdgd + (db1(length(db1)/2+1)+db0(n))/corr_length; - + h.X = exp(-1j*(common_beta_x-delta_beta)*dz); h.Y = exp(-1j*(common_beta_y+delta_beta)*dz); - % delta_beta has to be added to the transfer function - + % process with transfer function tOpt.X = h.X.*tOpt.X ; tOpt.Y = h.Y.*tOpt.Y ; - + % rotate back opt_x = matR(1,1)*tOpt.X + matR(1,2)*tOpt.Y; opt_y = matR(2,1)*tOpt.X + matR(2,2)*tOpt.Y; - + end lin_z_test = lin_z_test + sum(plate_sizes,2); @@ -117,9 +106,8 @@ lin_z_test = lin_z_test + sum(plate_sizes,2); %update the number of processed plates so far n_plates_done = n_plates_done + n_plates_left; -% attanuate the signal each linear state with alpha -% ( 0.2dB = 4.6052e-05 ) -rOpt_x=ifft(exp(-alpha_lin_x*aStepSize/2).*opt_x); % /2 not sure why (have to find it in formulas) -rOpt_y=ifft(exp(-alpha_lin_y*aStepSize/2).*opt_y); % but not relevant for now +% attenuate the signal each linear step with alpha +rOpt_x=ifft(exp(-alpha_lin_x*aStepSize/2).*opt_x); +rOpt_y=ifft(exp(-alpha_lin_y*aStepSize/2).*opt_y); -end \ No newline at end of file +end diff --git a/Classes/DataBaseHandler/DBHandler.m b/Classes/DataBaseHandler/DBHandler.m index 5569746..b7ba026 100644 --- a/Classes/DataBaseHandler/DBHandler.m +++ b/Classes/DataBaseHandler/DBHandler.m @@ -23,6 +23,12 @@ classdef DBHandler < handle arguments options.dataBase = ""; % Default value for pathToDB if not provided options.type = "mysql"; + options.server = ""; + options.port = 3306; + options.user = ""; + options.password = ""; + + end % Assign values to class properties based on input arguments @@ -46,12 +52,13 @@ classdef DBHandler < handle obj.conn = database( ... string(obj.dataBase), ... % Database name - "silas", ... % Username - "silas", ... % Password (or getSecret) + options.user, ... % Username + options.password, ... % Password (or getSecret) "Vendor", "MySQL", ... - "Server", "134.245.243.254", ... + "Server", options.server, ... "PortNumber", 3306, ... "JDBCDriverLocation", "C:\Users\Silas\Documents\mysql-connector-j-9.3.0\mysql-connector-j-9.3.0.jar"); + end catch e diff --git a/Functions/Theory/Dissertation/mach_zehnder_modulator.m b/Functions/Theory/Dissertation/mach_zehnder_modulator.m index 7b379c2..f9eca58 100644 --- a/Functions/Theory/Dissertation/mach_zehnder_modulator.m +++ b/Functions/Theory/Dissertation/mach_zehnder_modulator.m @@ -1,16 +1,4 @@ -% 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] @@ -25,27 +13,60 @@ 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) +fs = 200e9; % [Hz] +Nper = 2; % number of periods +Vpp = 0.6*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) +biasV = 1.1; % [V] differential bias added to v_drive % 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) + +if 1 + % SINE + v_drive = biasV + (Vpp/2)*sin(2*pi*f0*t); % v_drive(t) (peak = Vpp/2) + +else + + % --- Generate PAM-4 Sequence --- + symbols = linspace(-0.5, 0.5, 4); + num_symbols = 12; % Increased slightly for better visual + rng(44); + random_data = symbols(randi(4, 1, num_symbols)); + + % Create time axis (Note: T is your period from the sine code) + sps = round(T * fs); + t = (0:1/fs:(num_symbols*T)-1/fs).'; + + % Upsample to rectangular waveform + v_pam = repelem(random_data, sps).'; + + % Apply swing and bias: Resulting range is [biasV-Vpp/2, biasV+Vpp/2] + v_drive_rect = biasV + (v_pam * Vpp); + + % --- Round the edges --- + filter_span = round(sps/1.5); % Increased span for smoother "rounding" + window = gausswin(filter_span); + window = window / sum(window); + + % Apply filter (using 'same' to keep vector length, but be aware of edge transients) + v_drive = conv(v_drive_rect, window, 'same'); + +end + + +% 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; % Push-pull branch voltages (consistent with v_drive = v1 - v2) v1 = +0.5*v_drive; % arm 1 @@ -69,8 +90,6 @@ 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) ); @@ -78,9 +97,6 @@ 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'); @@ -100,8 +116,11 @@ 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_)]); +% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\linear_casee\mzm_input_signal.tex'); + %% IN/OUT (static transfer) — normalized x-axis + analytic curve +if 0 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 @@ -123,18 +142,78 @@ xlim([min(v_) max(v_)+1]); ylim([-1 1]); % mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm.tex'); - - +end %% -% % FIELD TF (only field here; do not mix power into this figure) + figure(3); 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 + +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_)]); +ylim([-1 1]); + +% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm_tramsfer_function_matlab.tex'); + + +%% + +figure(4); 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','black'); hold on;% analytic power TF +input_dots = linspace(min(v_drive),max(v_drive),4)./Vpi; +% input_dots = unique(v_drive_rect)./Vpi; +output_dots = (cos((pi/2)*input_dots).^2); +scatter(input_dots,output_dots,'Marker','x','LineWidth',1,'MarkerEdgeColor','black'); +scatter(input_dots,zeros(size(input_dots)),'Marker','^','LineWidth',2,'MarkerEdgeColor','black'); +% scatter(ones(size(input_dots)),output_dots,'Marker','<','LineWidth',2,'MarkerEdgeColor','black'); + +for i = 1:numel(input_dots) + % Draw the dashed projection lines + line([input_dots(i), input_dots(i)], [output_dots(i), 0], 'linewidth', 0.5, 'color', 'black', 'linestyle', '--', 'handlevisibility', 'off'); + line([input_dots(i), 1], [output_dots(i), output_dots(i)], 'linewidth', 0.5, 'color', 'black', 'linestyle', '--', 'handlevisibility', 'off'); + + % Add the level annotation boxes near the output (y-axis) + % Adjust the '1.05' to move the box further right or 'output_dots(i)' for height + j = 3-(i-1)*2; + text(1, output_dots(i), sprintf('Level %d', j), ... + 'FontSize', 8, ... + 'EdgeColor', 'black', ... + 'BackgroundColor', 'white', ... + 'Margin', 2); +end + +xlim([0,1.5]); +ylim([0,1]) + +% line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), 0],'linewidth',0.5,'color','black','linestyle','--'); +% line([max(v_drive), max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), 0],'linewidth',0.5,'color','black','linestyle','--'); + +% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\linear_casee\mzm_tf.tex'); +% xticks(sort(input_dots)); +% yticks(sort(output_dots)); +grid off + + +%% +% % FIELD TF (only field here; do not mix power into this figure) +figure(5); 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'); +yticks(sort(output_dots)); +% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\linear_casee\mzm_output_signal.tex'); 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 ab31451..67a2ff8 100644 --- a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m @@ -10,10 +10,11 @@ if dsp_options.mode == "load_run_id" if experiment == "highspeed_2024" - dsp_options.database_type = 'mysql'; + dsp_options.database_type = "mysql"; dsp_options.dataBase = 'labor_highspeed';%'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\silas_labor_newdsp_newstructure.db'; dsp_options.storage_path = 'Z:\2024\sioe_labor\'; - db = DBHandler("dataBase", [dsp_options.dataBase], "type", dsp_options.database_type); + db = DBHandler("dataBase", [dsp_options.dataBase],... + "type", dsp_options.database_type,"server","192.168.178.192","user","silas","password","silas"); elseif experiment == "mpi_ecoc_2025" diff --git a/projects/WDM/WDM_model_10km_queue.m b/projects/WDM/WDM_model_10km_queue.m index a41cb02..8aae2ea 100644 --- a/projects/WDM/WDM_model_10km_queue.m +++ b/projects/WDM/WDM_model_10km_queue.m @@ -312,7 +312,7 @@ for realiz = 1:s.num_realiz Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... - "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib); + "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",1,"useSingle",1).process(Opt_sig_wdm_fib); % -------- Evaluate at intermediate distance (enqueue jobs) -------- if eval_ptr <= nEval && seg == eval_seg(eval_ptr) diff --git a/test/bayesopt_ffe_tuning.m b/test/bayesopt_ffe_tuning.m new file mode 100644 index 0000000..adedb3d --- /dev/null +++ b/test/bayesopt_ffe_tuning.m @@ -0,0 +1,219 @@ +%% Bayesian Optimization for FFE Parameter Tuning +% This script uses bayesopt to find optimal mu_dd and mu_tr values +% that minimize BER for the FFE equalizer. + +clear; clc; + +%% Setup - Same as gpu_processing_dpfiber.m +s.wavelengthplan = calcWavelengthPlan(4, 400e9, 1310); +link_length = 10; +s.pmd = 0.1; +s.gamma = 0.0023; + +s.M = 4; +fsym = 112e9; +fdac = 2*fsym; +fadc = 120000000000; +s.random_key = 1; + +% Laser / Modulator +vbias_rel = 0.5; +u_pi = 4.6; +vbias = -vbias_rel*u_pi; +laser_linewidth = 0e6; + +duob_mode = db_mode.no_db; +rcalpha = 0.05; +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + +s.chirpalpha = 0; +s.p_launch = 3; +s.p = "co"; + +N = numel(s.wavelengthplan); + +switch s.p + case "co" + pol_rot = 100.*ones(1,N); + d_local = 0; +end + +f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9); +margin = 25e12; +f_span = (max(f_plan)+margin)-(min(f_plan)-margin); +f_nyq = f_span/2; + +kover = 4; +upsample_required = f_nyq./(fdac*kover/2); +upsample_pow = 2^nextpow2(upsample_required); + +s.f_opt = fdac*kover*upsample_pow; +s.f_opt_nyq = s.f_opt/2; + +s.rop = -8; % Fixed ROP for optimization + +%% Generate TX signals (run once) +fprintf('Generating TX signals...\n'); +for l = 1:N + [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ... + "fsym",fsym,"M",s.M,"order",15,"useprbs",0, ... + "fs_out",fdac, ... + "applyclipping",0,"clipfactor",1.5, ... + "applypulseform",1,"pulseformer",Pform, ... + "randkey",s.random_key+l, ... + "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ... + ).process(); + + Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ... + "filterType",filtertypes.gaussian,"active",true); + + El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ... + "bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ... + "H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig); + + El_sig = El_sig.normalize("mode","oneone"); + scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); + El_sig = El_sig .* scaling; + + Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ... + "lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ... + "linewidth",laser_linewidth,"randomkey",s.random_key+l,"alpha",s.chirpalpha).process(El_sig); + + signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out); +end + +%% WDM mux + launch +Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ... + "lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell); + +Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... + "amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm); + +%% Fiber propagation +segment_length = 1; +nSegments = link_length/segment_length; +nSegments = round(nSegments); + +zdw = 1310; +randomize_D = true; +Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key); + +Opt_sig_wdm_fib = Opt_sig_wdm; +fprintf('Running fiber propagation...\n'); +for seg = 1:nSegments + fprintf('Segment %d/%d\n', seg, nSegments); + Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... + "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... + "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... + "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",true).process(Opt_sig_wdm_fib); +end + +%% Pre-process to get Rx_sig (do demux once) +fprintf('Pre-processing receiver chain...\n'); +l = 1; % Use channel 1 for optimization + +Opt_sig_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ... + "fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib); + +Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... + "amplification_db",s.rop).process(Opt_sig_demux{l}); + +PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ... + "nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx); + +rx_bwl = 100e9; +PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ... + "filterType",filtertypes.butterworth,"active",true).process(PD_sig); + +Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); +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',0,'H_lpf',Lp_scpe).process(PD_sig); + +Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym); + +[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0); +Rx_sig = Scpe_cell{1}; +Rx_sig = Rx_sig.normalize("mode","rms"); + +fprintf('Receiver pre-processing complete. Ready for optimization.\n\n'); + +%% Define the objective function for bayesopt +function ber = ffe_objective(params, Rx_sig, Symbols_l, Tx_bits_l, M, duob_mode) +mu_dd = params.mu_dd; +mu_tr = params.mu_tr; + +try + eq_ffe = FFE("epochs_tr", 5, "epochs_dd", 2, "len_tr", 2^13, ... + "mu_dd", mu_dd, "mu_tr", mu_tr, ... + "order", 50, "sps", 2, "decide", 0, ... + "adaption", adaption_method.nlms, "dd_mode", 1); + + ffe_results = ffe(eq_ffe, M, Rx_sig, Symbols_l, Tx_bits_l, ... + "precode_mode", duob_mode, ... + 'showAnalysis', 0, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + + ber = ffe_results.metrics.BER; + + if ber == 0 + ber = 1e-10; + end + + if ~isfinite(ber) + ber = 0.5; + end + + fprintf(' mu_dd=%.4e, mu_tr=%.4e -> BER=%.4e\n', mu_dd, mu_tr, ber); + +catch ME + fprintf(' mu_dd=%.4e, mu_tr=%.4e -> FAILED (%s)\n', mu_dd, mu_tr, ME.message); + ber = 0.5; +end +end + +%% Define optimizable variables +mu_dd_var = optimizableVariable('mu_dd', [1e-5, 0.1], 'Transform', 'log'); +mu_tr_var = optimizableVariable('mu_tr', [1e-5, 0.1], 'Transform', 'log'); + +%% Run Bayesian Optimization +fprintf('========== Starting Bayesian Optimization ==========\n'); +fprintf('Optimizing mu_dd and mu_tr to minimize BER\n'); +fprintf('Search range: mu_dd=[1e-5, 0.1], mu_tr=[1e-5, 0.1]\n\n'); + +objective_fn = @(params) ffe_objective(params, Rx_sig, Symbols{l}, Tx_bits{l}, s.M, duob_mode); + +results = bayesopt(objective_fn, [mu_dd_var, mu_tr_var], ... + 'MaxObjectiveEvaluations', 30, ... + 'AcquisitionFunctionName', 'expected-improvement-plus', ... + 'IsObjectiveDeterministic', false, ... + 'ExplorationRatio', 0.5, ... + 'Verbose', 1, ... + 'PlotFcn', []); + +%% Display Results +fprintf('\n========== FFE Optimization Complete ==========\n'); +fprintf('Best FFE parameters found:\n'); +fprintf(' mu_dd = %.6e\n', results.XAtMinObjective.mu_dd); +fprintf(' mu_tr = %.6e\n', results.XAtMinObjective.mu_tr); +fprintf(' BER = %.6e\n', results.MinObjective); + +%% Verify with optimal parameters +fprintf('\nVerifying optimal FFE parameters...\n'); +best_mu_dd = results.XAtMinObjective.mu_dd; +best_mu_tr = results.XAtMinObjective.mu_tr; + +eq_ffe_best = FFE("epochs_tr", 5, "epochs_dd", 2, "len_tr", 2^13, ... + "mu_dd", best_mu_dd, "mu_tr", best_mu_tr, ... + "order", 50, "sps", 2, "decide", 0, ... + "adaption", adaption_method.nlms, "dd_mode", 1); + +ffe_results_best = ffe(eq_ffe_best, s.M, Rx_sig, Symbols{l}, Tx_bits{l}, ... + "precode_mode", duob_mode, ... + 'showAnalysis', 1, ... + "postFFE", [], ... + "eth_style_symbol_mapping", 0); + +fprintf('\nFinal FFE BER with optimal parameters: %.6e\n', ffe_results_best.metrics.BER); diff --git a/test/gpu_cpu_comparison.m b/test/gpu_cpu_comparison.m new file mode 100644 index 0000000..02d40ba --- /dev/null +++ b/test/gpu_cpu_comparison.m @@ -0,0 +1,349 @@ +%% GPU vs CPU Comparison Test for DP_Fiber +% This script runs the fiber simulation with and without GPU acceleration +% and compares the numerical results. + +% clear; clc; + +%% Setup (same as gpu_processing_dpfiber.m but simplified) +s.wavelengthplan = calcWavelengthPlan(4, 400e9, 1310); +link_length = 10; +s.pmd = 0.1; +s.gamma = 0.0023; + +s.M = 4; +fsym = 112e9; +fdac = 2*fsym; +fadc = 120000000000; +s.random_key = 1; + +% Laser / Modulator +vbias_rel = 0.5; +u_pi = 4.6; +vbias = -vbias_rel*u_pi; +laser_linewidth = 0e6; + +% DB Stuff +duob_mode = db_mode.no_db; + +rcalpha = 0.05; +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + +s.chirpalpha = 0; +s.p_launch = 3; +s.p = "co"; + +N = numel(s.wavelengthplan); + +switch s.p + case "co" + pol_rot = 100.*ones(1,N); + d_local = 0; + case "pair" + pol_rot = repmat([100,100,0,0],1,N/4); + d_local = 0; + case "alt" + pol_rot = repmat([100,0,100,0],1,N/4); + d_local = 0; + case "seg" + pol_rot = 100.*ones(1,N); + d_local = 3; + otherwise + error('Unknown fwm_mitigation_technique: %s', string(s.p)); +end + +f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9); +margin = 5e12; +f_span = (max(f_plan)+margin)-(min(f_plan)-margin); +f_nyq = f_span/2; + +kover = 4; +upsample_required = f_nyq./(fdac*kover/2); +upsample_pow = 2^nextpow2(upsample_required); + +s.f_opt = fdac*kover*upsample_pow; +s.f_opt_nyq = s.f_opt/2; + +%% TX per channel +for l = 1:N + [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ... + "fsym",fsym,"M",s.M,"order",15,"useprbs",0, ... + "fs_out",fdac, ... + "applyclipping",0,"clipfactor",1.5, ... + "applypulseform",1,"pulseformer",Pform, ... + "randkey",s.random_key+l, ... + "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ... + ).process(); + + Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ... + "filterType",filtertypes.gaussian,"active",true); + + El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ... + "bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ... + "H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig); + + clear Digi_sig + + El_sig = El_sig.normalize("mode","oneone"); + scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); + El_sig = El_sig .* scaling; + + Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ... + "lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ... + "linewidth",laser_linewidth,"randomkey",s.random_key+l,"alpha",s.chirpalpha).process(El_sig); + + clear El_sig + + signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out); + + clear Eml_out Lp_awg +end + +disp('Signal generated for all channels.'); + +%% WDM mux + launch +Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ... + "lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell); + +Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... + "amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm); + +% Save input for both runs +Opt_sig_input = Opt_sig_wdm; + +segment_length = 1; +nSegments = link_length/segment_length; +if abs(nSegments - round(nSegments)) > 1e-12 + error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length); +end +nSegments = round(nSegments); + +zdw = 1310; +randomize_D = true; + +if nSegments > 0 + Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key); +else + Dvec = []; +end + +%% Run WITHOUT GPU +fprintf('\n========== Running WITHOUT GPU (CPU) ==========\n'); +Opt_sig_cpu = Opt_sig_input; + +tic; +for seg = 1:nSegments + fprintf('CPU Segment %d/%d \n',seg, nSegments); + + Opt_sig_cpu = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... + "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... + "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... + "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",false).process(Opt_sig_cpu); +end +time_cpu = toc; +fprintf('CPU Time: %.3f seconds\n', time_cpu); + +%% Run WITH GPU (Double Precision) +fprintf('\n========== Running WITH GPU (Double Precision) ==========\n'); +Opt_sig_gpu_double = Opt_sig_input; + +tic; +for seg = 1:nSegments + fprintf('GPU-Double Segment %d/%d \n',seg, nSegments); + + Opt_sig_gpu_double = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... + "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... + "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... + "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",false).process(Opt_sig_gpu_double); +end +time_gpu_double = toc; +fprintf('GPU Double Time: %.3f seconds\n', time_gpu_double); + +%% Run WITH GPU (Single Precision) +fprintf('\n========== Running WITH GPU (Single Precision) ==========\n'); +Opt_sig_gpu_single = Opt_sig_input; + +tic; +for seg = 1:nSegments + fprintf('GPU-Single Segment %d/%d \n',seg, nSegments); + + Opt_sig_gpu_single = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... + "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... + "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... + "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",true).process(Opt_sig_gpu_single); +end +time_gpu_single = toc; +fprintf('GPU Single Time: %.3f seconds\n', time_gpu_single); + +%% Compare Results +fprintf('\n========== Numerical Comparison ==========\n'); + +sig_cpu = Opt_sig_cpu.signal; +sig_gpu_double = Opt_sig_gpu_double.signal; +sig_gpu_single = Opt_sig_gpu_single.signal; + +% Check dimensions +fprintf('CPU signal size: [%d x %d]\n', size(sig_cpu,1), size(sig_cpu,2)); +fprintf('GPU Double signal size: [%d x %d]\n', size(sig_gpu_double,1), size(sig_gpu_double,2)); +fprintf('GPU Single signal size: [%d x %d]\n', size(sig_gpu_single,1), size(sig_gpu_single,2)); + +% CPU vs GPU Double +fprintf('\n--- CPU vs GPU Double ---\n'); +diff_cpu_double = abs(sig_cpu - sig_gpu_double); +max_diff_cpu_double = max(diff_cpu_double(:)); +mean_diff_cpu_double = mean(diff_cpu_double(:)); +rel_diff_cpu_double = max_diff_cpu_double / max(abs(sig_cpu(:))); +fprintf('Max absolute difference: %.6e\n', max_diff_cpu_double); +fprintf('Mean absolute difference: %.6e\n', mean_diff_cpu_double); +fprintf('Max relative difference: %.6e\n', rel_diff_cpu_double); + +% CPU vs GPU Single +fprintf('\n--- CPU vs GPU Single ---\n'); +diff_cpu_single = abs(sig_cpu - sig_gpu_single); +max_diff_cpu_single = max(diff_cpu_single(:)); +mean_diff_cpu_single = mean(diff_cpu_single(:)); +rel_diff_cpu_single = max_diff_cpu_single / max(abs(sig_cpu(:))); +fprintf('Max absolute difference: %.6e\n', max_diff_cpu_single); +fprintf('Mean absolute difference: %.6e\n', mean_diff_cpu_single); +fprintf('Max relative difference: %.6e\n', rel_diff_cpu_single); + +% GPU Double vs GPU Single +fprintf('\n--- GPU Double vs GPU Single ---\n'); +diff_double_single = abs(sig_gpu_double - sig_gpu_single); +max_diff_double_single = max(diff_double_single(:)); +mean_diff_double_single = mean(diff_double_single(:)); +rel_diff_double_single = max_diff_double_single / max(abs(sig_gpu_double(:))); +fprintf('Max absolute difference: %.6e\n', max_diff_double_single); +fprintf('Mean absolute difference: %.6e\n', mean_diff_double_single); +fprintf('Max relative difference: %.6e\n', rel_diff_double_single); + +% Check tolerances +fprintf('\n--- Tolerance Check ---\n'); +tol_double = 1e-10; +tol_single = 1e-5; % Single precision has ~7 significant digits + +if max_diff_cpu_double < tol_double + fprintf('✓ CPU vs GPU Double: EQUIVALENT (diff < %.0e)\n', tol_double); +else + fprintf('✗ CPU vs GPU Double: DIFFER beyond tolerance (%.0e)\n', tol_double); +end + +if max_diff_cpu_single < tol_single + fprintf('✓ CPU vs GPU Single: ACCEPTABLE (diff < %.0e)\n', tol_single); +else + fprintf('⚠ CPU vs GPU Single: Precision loss detected (diff = %.2e, tol = %.0e)\n', max_diff_cpu_single, tol_single); +end + +% Performance comparison +fprintf('\n========== Performance Summary ==========\n'); +fprintf('CPU Time: %.3f s\n', time_cpu); +fprintf('GPU Double Time: %.3f s\n', time_gpu_double); +fprintf('GPU Single Time: %.3f s\n', time_gpu_single); +fprintf('\n'); +fprintf('Speedup (GPU Double vs CPU): %.2fx\n', time_cpu/time_gpu_double); +fprintf('Speedup (GPU Single vs CPU): %.2fx\n', time_cpu/time_gpu_single); +fprintf('Speedup (GPU Single vs GPU Double): %.2fx\n', time_gpu_double/time_gpu_single); + +%% ========== BER Comparison ========== +% Process each fiber output through simplified receiver to check if +% single-precision affects actual BER performance + +fprintf('\n========== BER Comparison ==========\n'); +fprintf('Processing signals through receiver chain...\n'); + +% Receiver parameters +rop = -7; % Received optical power [dBm] +len_tr = 4096; % Training length +mu_dc = 0.005; +mu_ffe = [0.0001 0.0008 0.001]; +mu_dfe = 0.0004; + +% Helper function to process through receiver and get BER +function ber = process_receiver(Opt_sig_fib, l, Symbols, Tx_bits, ... + fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode) + +% Demux single channel +Opt_sig_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ... + "fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_fib); + +% ROP amplifier +Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... + "amplification_db",rop).process(Opt_sig_demux{l}); + +% Photodiode +PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ... + "nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx); + +% Low-pass filter +rx_bwl = 100e9; +PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ... + "filterType",filtertypes.butterworth,"active",true).process(PD_sig); + +% Scope +Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); +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',0,'H_lpf',Lp_scpe).process(PD_sig); + +% Resample to 2 sps +Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym); + +% Time sync +[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0); +Rx_sig = Scpe_cell{1}; +Rx_sig = Rx_sig.normalize("mode","rms"); + +% FFE Equalizer +ffe_order = [50, 0, 0]; +eq_ffe = EQ("Ne",ffe_order,"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.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + +ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols{l},Tx_bits{l}, ... + "precode_mode",duob_mode, ... + 'showAnalysis',0, ... + "postFFE",[], ... + "eth_style_symbol_mapping",0); + +ber = ffe_results.metrics.BER; +end + +% Process each mode for channel 1 +l = 4; % Use first channel for comparison + +fprintf('Processing CPU result...\n'); +ber_cpu = process_receiver(Opt_sig_cpu, l, Symbols, Tx_bits, ... + fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode); + +fprintf('Processing GPU Double result...\n'); +ber_gpu_double = process_receiver(Opt_sig_gpu_double, l, Symbols, Tx_bits, ... + fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode); + +fprintf('Processing GPU Single result...\n'); +ber_gpu_single = process_receiver(Opt_sig_gpu_single, l, Symbols, Tx_bits, ... + fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode); + +% Display BER results +fprintf('\n========== BER Results (Channel %d, ROP = %d dBm) ==========\n', l, rop); +fprintf('CPU: BER = %.4e\n', ber_cpu); +fprintf('GPU Double: BER = %.4e\n', ber_gpu_double); +fprintf('GPU Single: BER = %.4e\n', ber_gpu_single); + +fprintf('\n--- BER Comparison ---\n'); +if ber_cpu == 0 && ber_gpu_double == 0 && ber_gpu_single == 0 + fprintf('✓ All BERs are zero (no errors detected)\n'); +else + ber_diff_double = abs(ber_cpu - ber_gpu_double); + ber_diff_single = abs(ber_cpu - ber_gpu_single); + fprintf('|BER_cpu - BER_gpu_double| = %.4e\n', ber_diff_double); + fprintf('|BER_cpu - BER_gpu_single| = %.4e\n', ber_diff_single); + + if ber_diff_double < 1e-6 && ber_diff_single < 1e-6 + fprintf('✓ BER differences are negligible\n'); + elseif ber_diff_single > ber_diff_double * 10 + fprintf('⚠ Single precision shows measurable BER impact\n'); + else + fprintf('✓ BER differences within acceptable range\n'); + end +end + +fprintf('\n========== Test Complete ==========\n'); diff --git a/test/gpu_processing_dpfiber.m b/test/gpu_processing_dpfiber.m new file mode 100644 index 0000000..317833e --- /dev/null +++ b/test/gpu_processing_dpfiber.m @@ -0,0 +1,265 @@ + + +s.wavelengthplan = calcWavelengthPlan(16, 400e9, 1310); +N = numel(s.wavelengthplan); +link_length = 10; +s.pmd = 0.1;%0.1; +s.gamma = 0.0023; + +s.M = 4; +fsym = 112e9; +fdac = 2*fsym; +fadc = 120000000000; +s.random_key = 1; + +% Laser / s.Modulator +vbias_rel = 0.5; +u_pi = 4.6; +vbias = -vbias_rel*u_pi; +laser_linewidth = 0e6; + +% DB Stuff +duob_mode = db_mode.no_db; + +rcalpha = 0.05; +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + +s.chirpalpha = 0; + +s.p_launch = 3; +s.p = "co"; + +switch s.p + case "co" + pol_rot = 100.*ones(1,N); + d_local = 0; + case "pair" + pol_rot = repmat([100,100,0,0],1,N/4); + d_local = 0; + case "alt" + pol_rot = repmat([100,0,100,0],1,N/4); + d_local = 0; + case "seg" + pol_rot = 100.*ones(1,N); + d_local = 3; + otherwise + error('Unknown fwm_mitigation_technique: %s', string(s.p)); +end + +f_plan = physconst('lightspeed')./(s.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; + +kover = 4; +upsample_required = f_nyq./(fdac*kover/2); +upsample_pow = 2^nextpow2(upsample_required); + +s.f_opt = fdac*kover*upsample_pow; +s.f_opt_nyq = s.f_opt/2; + +s.rop = -10:1:0; + +profile on + +%% ---------- TX per channel ---------- +for l = 1:N + + [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ... + "fsym",fsym,"M",s.M,"order",17,"useprbs",0, ... + "fs_out",fdac, ... + "applyclipping",0,"clipfactor",1.5, ... + "applypulseform",1,"pulseformer",Pform, ... + "randkey",s.random_key+l, ... + "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ... + ).process(); + + Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ... + "filterType",filtertypes.gaussian,"active",true); + + El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ... + "bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ... + "H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig); + + % Digi_sig not needed after AWG + clear Digi_sig + + % Electrical Driver Amplifier + El_sig = El_sig.normalize("mode","oneone"); + scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); + El_sig = El_sig .* scaling; + + % E/O Conversion + Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ... + "lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ... + "linewidth",laser_linewidth,"randomkey",s.random_key+l,"alpha",s.chirpalpha).process(El_sig); + + % El_sig not needed after EML + clear El_sig + + signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out); + + % Eml_out not needed after pol controller + clear Eml_out Lp_awg +end + +disp('Signal generated for all channels.'); + +%% ---------- WDM mux + launch ---------- +Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ... + "lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell); + +Opt_sig_wdm.spectrum(); + +Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... + "amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm); + +Opt_sig_wdm_fib = Opt_sig_wdm; + +segment_length = 1; +nSegments = link_length/segment_length; +if abs(nSegments - round(nSegments)) > 1e-12 + error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length); +end +nSegments = round(nSegments); + +zdw = 1310; +randomize_D = true; + +% Guard for 0 km: avoid calling getDispersionVector(0,...) if it doesn't support it +if nSegments > 0 + Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key); +else + Dvec = []; +end + +for seg = 1:nSegments + + fprintf('Segment %d/%d \n',seg, nSegments); + + Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ... + "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ... + "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ... + "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",1).process(Opt_sig_wdm_fib); + +end + +profile off +profile viewer + +Opt_sig_wdm_fib.spectrum(); + +%% ========== BER Evaluation ========== +fprintf('\n========== BER Evaluation ==========\n'); +fprintf('Processing signals through receiver chain...\n'); + +% Receiver parameters +len_tr = 4096; % Training length +mu_dc = 0.005; +mu_ffe = [0.0001 0.0008 0.001]; +mu_dfe = 0.0004; + +% Helper function to process through receiver and get BER +function ber = process_receiver(Opt_sig_fib, l, Symbols, Tx_bits, ... + fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode) + +% Demux single channel +Opt_sig_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ... + "fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_fib); + +% ROP amplifier +Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... + "amplification_db",rop).process(Opt_sig_demux{l}); + +% Photodiode +PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ... + "nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx); + +% Low-pass filter +rx_bwl = 100e9; +PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ... + "filterType",filtertypes.butterworth,"active",true).process(PD_sig); + +% Scope +Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); +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',0,'H_lpf',Lp_scpe).process(PD_sig); + +% Resample to 2 sps +Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym); + +% Time sync +[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0); +Rx_sig = Scpe_cell{1}; +Rx_sig = Rx_sig.normalize("mode","rms"); + +% FFE Equalizer +ffe_order = [50, 0, 0]; +eq_ffe = EQ("Ne",ffe_order,"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.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); +eq_ffe = FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",6.624e-05,"mu_tr",0.058136,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1); + +ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols{l},Tx_bits{l}, ... + "precode_mode",duob_mode, ... + 'showAnalysis',0, ... + "postFFE",[], ... + "eth_style_symbol_mapping",0); + +ber = ffe_results.metrics.BER; +end + +% Process each ROP value for selected channels using parfor +ber_results = zeros(length(s.rop), N); + +% Flatten loop for parfor: iterate over all (ROP, channel) combinations +num_rop = length(s.rop); +rop_vals = s.rop; +ber_flat = zeros(num_rop * N, 1); + +parfor idx = 1:(num_rop * N) + % Convert linear index to (ri, l) subscripts + ri = ceil(idx / N); + l = mod(idx - 1, N) + 1; + + fprintf('ROP %d dBm, Channel %d/%d\n', rop_vals(ri), l, N); + ber_flat(idx) = process_receiver(Opt_sig_wdm_fib, l, Symbols, Tx_bits, ... + fdac, kover, upsample_pow, fsym, fadc, rop_vals(ri), s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode); +end + +% Reshape back to [num_rop × N] matrix +ber_results = reshape(ber_flat, [N, num_rop]).'; + + +%% Display BER Results +fprintf('\n========== BER Results ==========\n'); +fprintf('ROP [dBm] | '); +for l = 1:N + fprintf('Ch%d | ', l); +end +fprintf('\n'); + +for ri = 1:length(s.rop) + fprintf('%8d | ', s.rop(ri)); + for l = 1:N + fprintf('%.2e | ', ber_results(ri, l)); + end + fprintf('\n'); +end + +% Plot BER vs ROP +figure; +semilogy(s.rop, mean(ber_results, 2), '-o', 'LineWidth', 2); +hold on; +for l = 1:N + semilogy(s.rop, ber_results(:, l), '--', 'LineWidth', 1); +end +hold off; +xlabel('ROP [dBm]'); +ylabel('BER'); +title('BER vs Received Optical Power (GPU Single Precision)'); +legend(['Mean', arrayfun(@(x) sprintf('Ch%d', x), 1:N, 'UniformOutput', false)]); +grid on; + +fprintf('\n========== Test Complete ==========\n'); diff --git a/test/gpu_processing_test.m b/test/gpu_processing_test.m new file mode 100644 index 0000000..b08c03f --- /dev/null +++ b/test/gpu_processing_test.m @@ -0,0 +1,70 @@ +% 1. Setup Data OUTSIDE the timer +d = gpuDevice; +N = 10000; + +fprintf('Preparing data...\n'); +A_cpu_double = rand(N, N); % Create double on CPU +A_cpu_single = single(A_cpu_double); % Create single on CPU + +% Warmup run (wakes up the GPU from idle state) +A_warm = gpuArray.rand(1000, 1000, 'single'); +B_warm = A_warm * A_warm; +wait(d); + +fprintf('------------------------------------------------\n'); + +% TEST 1: Double Precision (The "Slow" way) +% We move data to GPU first so we only measure calculation time +A_gpu_double = gpuArray(A_cpu_double); +wait(d); % Ensure transfer is done before starting timer + +fprintf('Running DOUBLE precision test... '); +tic; +B_gpu = A_gpu_double * A_gpu_double; +wait(d); % FORCE MATLAB TO WAIT FOR GPU +time_double = toc; +fprintf('Done.\n'); +fprintf('Double Precision Time: %.4f seconds\n', time_double); + +% TEST 2: Single Precision (The "Fast" way) +A_gpu_single = gpuArray(A_cpu_single); +wait(d); % Ensure transfer is done + +fprintf('Running SINGLE precision test... '); +tic; +B_gpu = A_gpu_single * A_gpu_single; +wait(d); % FORCE MATLAB TO WAIT FOR GPU +time_single = toc; +fprintf('Done.\n'); +fprintf('Single Precision Time: %.4f seconds\n', time_single); + +% Calculate Speedup +fprintf('------------------------------------------------\n'); +fprintf('Speedup Factor using single precision: %.2fx\n', time_double / time_single); + +% +% Create 10,000 small matrices (10x10) stacked in a 3D array +A_stack = gpuArray.rand(10, 10, 10000, 'single'); +B_stack = gpuArray.rand(10, 10, 10000, 'single'); + +% BAD: Looping (GPU overhead kills you) +tic; +for i=1:10000 + C(:,:,i) = A_stack(:,:,i) * B_stack(:,:,i); +end +wait(d); +loop_time = toc; + +% GOOD: Pagefun (Executes all 10,000 mults simultaneously) +tic; +C_stack = pagefun(@mtimes, A_stack, B_stack); +wait(d); +pagefun_time = toc; + +fprintf('------------------------------------------------\n'); +fprintf('Speedup Factor using pagefun: %.2fx\n', loop_time / pagefun_time); + + +fprintf('Total VRAM: %.2f GB\n', d.TotalMemory / 1e9); +fprintf('Available VRAM: %.2f GB\n', d.AvailableMemory / 1e9); +fprintf('Usage: %.1f%%\n', 100 * (1 - d.AvailableMemory / d.TotalMemory)); \ No newline at end of file