diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index 6ba8558..cf4df17 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -45,7 +45,7 @@ classdef PAMmapper signal_in = signal_in.logbookentry(lbdesc,obj); out = signal_in; else - out = signal_in; + out = obj.map_(signal_in); end end @@ -328,8 +328,9 @@ classdef PAMmapper if ~obj.eth_style - data_in = data_in/(sqrt(mean(abs(data_in).^2))); + % data_in = data_in/(sqrt(mean(abs(data_in).^2))); data_in = data_in*sqrt(10); + %data_in = data_in*sqrt((5^2 + 3^2 + 1^2 + 5^2 + 3^2 + 1^2)/6); if size(data_in,2) > 1 data_in = data_in.'; diff --git a/Classes/01_transmit/PAMsource.m b/Classes/01_transmit/PAMsource.m index f1f8e5a..dd77f89 100644 --- a/Classes/01_transmit/PAMsource.m +++ b/Classes/01_transmit/PAMsource.m @@ -164,7 +164,7 @@ classdef PAMsource %%%%% Pulse-forming %%%%%% if obj.applypulseform %%% MY CODE - if 0 + if 1 digi_sig = obj.pulseformer.process(symbols); else diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index 54d0b57..5920e89 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -3,24 +3,34 @@ classdef FFE < handle % 1) Training mode (stable performance when you use NLMS) % 2) Decision directed mode - % Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); + %LMS: mu in order of 0.001 for acceptable convergence speed + %NLMS: mu in order of 0.01 for acceptable convergence speed + %RLS: mu is lambda -> 0.99 -> 1 (has a strong dependency on this! use a loop to find out best values) + + % FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption",adaption_method(adaption),"dd_mode",use_dd_mode); properties sps % usually 2 order e + e_tr error + debug_struct len_tr mu_tr epochs_tr - mu_dd + adaption_technique % nlms, lms, rls + dd_mode % 1 or 0 to set DD-mode on or off + mu_dd %weight update in dd mode epochs_dd - constellation + P % covariance matrix of rls - decide + constellation % symbol constellation + + decide %wether to return the (hard) decisions or the result after FFE (soft) end methods @@ -34,6 +44,8 @@ classdef FFE < handle options.mu_tr = 0; options.epochs_tr = 5; + options.adaption_technique adaption_method = adaption_method.lms; + options.dd_mode = 1; options.mu_dd = 1e-5; options.epochs_dd = 5; @@ -59,10 +71,15 @@ classdef FFE < handle obj.constellation = unique(D.signal); + + delta = 0.05; + obj.P = (1/delta) * eye(obj.order); + % Training Mode training = 1; showviz = 0; obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + obj.e_tr = obj.e; % Decision Directed Mode n = X.length; @@ -71,7 +88,7 @@ classdef FFE < handle [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training,showviz); % Output Signal - if obj.decide + if obj.decide X.signal = decision; else X.signal = signal; @@ -83,24 +100,25 @@ classdef FFE < handle Noi = X; Noi = X - D; - + end - function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz) + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) arguments - obj - x - d - mio - epochs - N - training + obj + x + d + mu + epochs + N + training showviz end x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + lambda = mu; if training mask = ones(obj.order,1); @@ -111,8 +129,17 @@ classdef FFE < handle end mask = ones(obj.order,1); + maincursor_pos=ceil(length(obj.e)/2); + always_ideal_decision = 0; + save_debug = 1; + grad =0; + weight = 0; + update = 0; + + if mu == 0 || (~obj.dd_mode && ~training) + epochs = 1; + end - for epoch = 1 : epochs symbol = 0; for sample = 1 : obj.sps : N @@ -124,30 +151,75 @@ classdef FFE < handle y(symbol,1) = (obj.e.*mask).' * U; % Calculating output of LMS __ * | if training - d_hat(symbol,1) = d(symbol); + d_hat(symbol,1) = d(symbol); else - [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point - d_hat(symbol,1) = obj.constellation(symbol_idx); + if ~always_ideal_decision + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + else + d_hat(symbol,1) = d(symbol); + end end - err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + % err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + err(symbol) = d_hat(symbol) - y(symbol); % Instantaneous error true_err(symbol) = y(symbol) - d(symbol); % Instantaneous error - if mio ~= 0 - obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS - else - normalizationfactor = (U.' * U); - obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS + if training || obj.dd_mode + switch obj.adaption_technique + + case adaption_method.lms + + % mu used as update weight (suggestion: 0.001) + weight = mu; + grad = err(symbol) * U; + update = grad * weight; + obj.e = obj.e + update; + + case adaption_method.nlms + + % mu used as update weight (suggestion: 0.01-0.05; bit higher during tr) + normU = ((U.'*U)) + eps; + weight = mu / normU; + grad = err(symbol) * U; + update = grad * weight; + obj.e = obj.e + update; + + case adaption_method.rls + + + % RLS‐Gain: + denom = lambda + U.' * obj.P * U; + k = (obj.P * U) / denom; + + % Gewichtsupdate: + update = k * err(symbol); + obj.e = obj.e + update; + + % P-Matrix‐Update: + obj.P = (1/lambda) * (obj.P - k * (U.' * obj.P)); + + + end end - obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + + if save_debug + obj.debug_struct.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + if training + obj.debug_struct.error_tr(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error + obj.debug_struct.update_tr(epoch,symbol) = update.'*update ./ rms(obj.e); + end + obj.debug_struct.main_cursor(epoch,symbol) = abs(obj.e(maincursor_pos)); + obj.debug_struct.mu_nlms(epoch,symbol) = weight; + obj.debug_struct.update_gradient(epoch,symbol) = grad.'*grad; + obj.debug_struct.update(epoch,symbol) = update.'*update ./ rms(obj.e); + end end end end - end -end - +end \ No newline at end of file diff --git a/Classes/04_DSP/Equalizer/FFE_DCremoval_adaptive_mu.m b/Classes/04_DSP/Equalizer/FFE_DCremoval_adaptive_mu.m index c6b1a43..9265652 100644 --- a/Classes/04_DSP/Equalizer/FFE_DCremoval_adaptive_mu.m +++ b/Classes/04_DSP/Equalizer/FFE_DCremoval_adaptive_mu.m @@ -9,6 +9,7 @@ classdef FFE_DCremoval_adaptive_mu < handle sps % usually 2 order e + e_tr error len_tr @@ -21,6 +22,8 @@ classdef FFE_DCremoval_adaptive_mu < handle mu_dc dc_buffer_len + adaptive_mu_mode + ffe_buffer_len smoothing_buffer_length @@ -50,6 +53,8 @@ classdef FFE_DCremoval_adaptive_mu < handle options.ffe_buffer_len = 1; + options.adaptive_mu_mode = 1; + options.smoothing_buffer_length = 0; options.smoothing_buffer_update = 0; options.decide = false; @@ -90,6 +95,7 @@ classdef FFE_DCremoval_adaptive_mu < handle % Training Mode training = 1; obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training); + obj.e_tr = obj.e; % Decision Directed Mode N = X.length; @@ -124,6 +130,10 @@ classdef FFE_DCremoval_adaptive_mu < handle training % boolean flag: true->training mode, false->DD mode end + if isempty(obj.e) + obj.e = zeros(obj.order,1); + end + % Zero-padding for filter memory x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; @@ -140,7 +150,7 @@ classdef FFE_DCremoval_adaptive_mu < handle err_prev = 0; % previous error sample for VSS correlation gamma_dc = 1e-6; % meta step-size for DC VSS mu_min = 1e-6; % lower bound for mu_dc - mu_max = 1e-1; % upper bound for mu_dc + mu_max = 3e-1; % upper bound for mu_dc % DC removal buffer L = obj.dc_buffer_len; % buffer length @@ -152,7 +162,6 @@ classdef FFE_DCremoval_adaptive_mu < handle if ~training % each column holds one past gradient of length obj.order grad_buf = NaN(obj.order, L_grad); - end smth_buffer = zeros(1, obj.smoothing_buffer_length); @@ -187,7 +196,11 @@ classdef FFE_DCremoval_adaptive_mu < handle %-- 3) error e_val = y(s) - d_hat(s); - err(s,idx) = e_val; + if epoch == epochs + + err(s,idx) = e_val; + true_err(s,idx) = y(s) - d(s); + end %-- 4) tap-weight update: training immediate, DD buffered if training @@ -199,7 +212,7 @@ classdef FFE_DCremoval_adaptive_mu < handle obj.e = obj.e - e_val * U / normU; end else - if 1 + if 0 % buffer gradient if mu_lms ~= 0 grad = e_val * U; @@ -226,22 +239,40 @@ classdef FFE_DCremoval_adaptive_mu < handle %-- 5) DC adaptation if obj.mu_dc ~= 0 - % VSS for mu_dc - delta_mu = gamma_dc * e_val * err_prev * (U.'*U); - obj.mu_dc = min(max(obj.mu_dc + delta_mu, mu_min), mu_max); - err_prev = e_val; + if obj.adaptive_mu_mode + + % VSS for mu_dc + delta_mu = gamma_dc * e_val * err_prev * (U.'*U); + obj.mu_dc = min(max(obj.mu_dc + delta_mu, mu_min), mu_max); + err_prev = e_val; + + % DC buffer update & periodic estimate + P_err = alpha*P_err + (1-alpha)*e_val^2; + mu_dc_norm = obj.mu_dc / (P_err + eps); + + else + + % DC buffer update & periodic estimate + % P_err = alpha*P_err + (1-alpha)*e_val^2; + % mu_dc_norm = obj.mu_dc / (P_err + eps); + + mu_dc_norm = obj.mu_dc; + + end - % DC buffer update & periodic estimate - P_err = alpha*P_err + (1-alpha)*e_val^2; - mu_dc_norm = obj.mu_dc / (P_err + eps); e_dc_buf = circshift(e_dc_buf, 1); e_dc_buf(1) = e_dc_est - mu_dc_norm * e_val; - + if mod(s, L) == 0 e_dc_est = median(e_dc_buf, 'omitnan'); end + P_err_save(s) = P_err; + % Pcorr_save(s) = e_val * err_prev; + Ucorr_save(s) = (U.'*U); + mu_dc_save(s) = mu_dc_norm; e_dc_save(s) = e_dc_est; + end % store instantaneous squared error @@ -250,23 +281,71 @@ classdef FFE_DCremoval_adaptive_mu < handle end % Optional plotting in DD mode (uncomment if needed) - if ~training - figure(342);clf + if 0%~training + + constellation = unique(d); + lvlcol = cbrewer2('Paired', numel(constellation)*2); + lvlcol = lvlcol(2:2:end, :); + + true_err(true_err==0) = NaN; + true_errmoverr = movsum(true_err, 4096, 'omitnan'); + true_errmoverr = true_errmoverr./rms(true_errmoverr); + + moverr = movsum(err, [100,100], 'omitnan'); + moverr = moverr./rms(moverr); + + figure(500); clf hold on - scatter(1:numSymbols, err + obj.constellation', '.', 'SizeData', 1); - - yline(obj.constellation); - + % 1st subplot: true_errmoverr + % subplot(2,2,1); hold on + % for k = 1:4 + % scatter(1:numSymbols, true_errmoverr(:,k), 1, lvlcol(k,:), '.'); + % end + + % scatter(1:numSymbols, Ucorr_save./rms(Ucorr_save), 1, lvlcol(1,:), '.','DisplayName','Ucorr_save'); + % scatter(1:numSymbols, Pcorr_save./rms(Pcorr_save), 1, lvlcol(1,:), '.','DisplayName','P_corr'); + % scatter(1:numSymbols, P_err_save, 1, lvlcol(1,:), '.','DisplayName','P_err'); + % scatter(1:numSymbols, mu_dc_save, 1, lvlcol(2,:), '.','DisplayName','adapted value of $\mu_{DC}$'); + % scatter(1:numSymbols, sum(moverr,2,'omitnan'), 1, lvlcol(1,:), '.','DisplayName','Mov Error $\hat{d}$ - x over all levels'); + scatter(1:numSymbols, sum(e_dc_save,2,'omitnan'), 1, lvlcol(2,:), '.','DisplayName','Est. Error that is subtracted'); + title('Moving Sum Error'); + hold off + legend + + % 2nd subplot: moverr + subplot(2,2,2); hold on + for k = 1:4 + scatter(1:numSymbols, moverr(:,k), 1, lvlcol(k,:), '.'); + end + title('Moving Sum Error'); + hold off + legend + + % 3rd subplot: err + subplot(2,2,3); hold on + for k = 1:4 + scatter(1:numSymbols, err(:,k), 1, lvlcol(k,:), '.'); + end + title('Error'); + hold off + legend + + % 4th subplot: err + obj.constellation' + subplot(2,2,4); hold on + for k = 1:4 + scatter(1:numSymbols, err(:,k) + obj.constellation(k), 1, lvlcol(k,:), '.'); + end + yline(obj.constellation, '--k'); + title('Error + Constellation'); + hold off + legend + + sgtitle('Error Analysis Subplots'); + end end - function mu = update_mu() - - - - end - end end diff --git a/Classes/04_DSP/Equalizer/FFE_DFE.m b/Classes/04_DSP/Equalizer/FFE_DFE.m index 0841f8d..e896500 100644 --- a/Classes/04_DSP/Equalizer/FFE_DFE.m +++ b/Classes/04_DSP/Equalizer/FFE_DFE.m @@ -60,7 +60,7 @@ classdef FFE_DFE < handle end - function [X] = process(obj, X, D) + function [X,Noi] = process(obj, X, D) % actual processing of the signal (steps 1. - 3.) % 1 normalize RMS @@ -88,6 +88,9 @@ classdef FFE_DFE < handle X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym lbdesc = [num2str(obj.ffe_order),' tap FFE']; X = X.logbookentry(lbdesc); % append to logbook + + Noi = X; + Noi = X - D; end @@ -185,8 +188,6 @@ classdef FFE_DFE < handle obj.e = coeff(1:obj.ffe_order); obj.b = coeff(obj.ffe_order+1:end); - - end end diff --git a/Classes/04_DSP/Equalizer/Postfilter.m b/Classes/04_DSP/Equalizer/Postfilter.m index 8d09ef1..659801e 100644 --- a/Classes/04_DSP/Equalizer/Postfilter.m +++ b/Classes/04_DSP/Equalizer/Postfilter.m @@ -31,7 +31,7 @@ classdef Postfilter < handle end - function signalclass_out = process(obj,signalclass_in,noiseclass_in,options) + function [signalclass_out,noiseclass_out] = process(obj,signalclass_in,noiseclass_in,options) arguments obj @@ -55,7 +55,8 @@ classdef Postfilter < handle else end - + + noiseclass_in = noiseclass_in.filter(obj.coefficients,1); signalclass_in = signalclass_in.filter(obj.coefficients,1); % append to logbook @@ -64,7 +65,7 @@ classdef Postfilter < handle % write to output signalclass_out = signalclass_in; - + noiseclass_out = noiseclass_in; end function showFilter(obj,options) diff --git a/Classes/04_DSP/Sequence Detection/MLSE.m b/Classes/04_DSP/Sequence Detection/MLSE.m index 9576524..c453e18 100644 --- a/Classes/04_DSP/Sequence Detection/MLSE.m +++ b/Classes/04_DSP/Sequence Detection/MLSE.m @@ -34,22 +34,28 @@ classdef MLSE < handle end - function [signalclass_hd,signalclass_sd] = process(obj,signalclass,ref_symbolclass) + function [signalclass_hd,LLR,GMI] = process(obj,signalclass,ref_symbolclass) data_in = signalclass.signal; + shape_in = size(data_in); data_ref = ref_symbolclass.signal; - [data_out_hd,data_out_sd] = obj.process_(data_in,data_ref); + [data_out_hd,LLR,GMI] = obj.process_(data_in,data_ref); + try + data_out_hd = reshape(data_out_hd,shape_in(1),shape_in(2)); + catch + warning('output reshaping failed after MLSE'); + end signalclass_hd = signalclass; signalclass_hd.signal = data_out_hd; - signalclass_sd = signalclass; - signalclass_sd.signal = data_out_sd; + % signalclass_sd = signalclass; + % signalclass_sd.signal = data_out_sd; end - function [VITERBI_ESTIMATION_SYMBOLS,soft_decisions] = process_(obj,data_in,data_ref) + function [VITERBI_ESTIMATION_SYMBOLS,LLR_maxlogmap,GMI] = process_(obj,data_in,data_ref) % remove unnecessary zeros at start of impulse response to keep @@ -70,35 +76,33 @@ classdef MLSE < handle %%%% Separate the equalized signal into the respective levels based on the actually transmitted level constellation = unique(data_ref); decisionLevels = (constellation(1:end-1) + constellation(2:end)) / 2; - tx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(data_ref); + N = length(data_in); + tx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(data_ref); - % impulse respnse i.e. [0.5, 1.0000] - obj.DIR = flip(obj.DIR); - - % RMS normalization of input data - data_in = data_in ./ rms(data_in); + % impulse respnse to remove from signal + obj.DIR = flip(obj.DIR); %i.e. -0.2676 -0.0478 1.0000 % seems to be the only way to use combvec for a flexible amount % of vectors. 'combs' contains all trellis states pre_comb_mat = repmat(obj.trellis_states,length(obj.DIR)-1,1); pre_comb_cell = mat2cell(pre_comb_mat,ones(1,size(pre_comb_mat,1)),size(pre_comb_mat,2)); combs = fliplr(combvec(pre_comb_cell{:}).'); + first_sym = combs(:,1); % das ist das älteste/ trailing Symbol aus der sequenz + last_sym = combs(:,end); %hiermit wird entschieden/ das ist das cursor symbol am ende der sequenz - % Save first and last symbol of each state - first_sym = combs(:,1); - last_sym = combs(:,end); states = sum(combs,2); + nStates = length(last_sym); % Calculate all possible input symbols for the desired impulse % response. Row number is the index of the previous state, % column number is the index of the next state % noise free received == branch metrics - noise_free_received = zeros(length(states),length(states)); + noise_free_received = zeros(nStates,nStates); count_row = 1; count_col = 1; - for l1 = 1:length(states) - for l2 = 1:length(states) + for l1 = 1:nStates + for l2 = 1:nStates if sum(combs(l2,2:end) == combs(l1,1:end-1)) == size(combs,2)-1 noise_free_received(count_row,count_col) = sum(combs(l2,:).*obj.DIR(end:-1:2)) + last_sym(l1)*obj.DIR(1); else @@ -110,8 +114,11 @@ classdef MLSE < handle count_row = 1; end - % match amplitude levels of input signal to those of the calculated ideal symbols - % i.e. match the rms values of data_in to noise_free_received + % first: RMS normalization of input data (rms==1) + data_in = data_in ./ rms(data_in); + + % then, match amplitude levels of input signal to those of the calculated ideal symbols + % i.e. match the rms values of data_in to noise_free_received (rms=1.xx) if isreal(data_in) if obj.M == round(obj.M) data_in = data_in * rms(noise_free_received(noise_free_received ~= inf),'all','omitnan'); @@ -119,338 +126,299 @@ classdef MLSE < handle end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%%%% FORWARD PASS %%%%% + %%%%% FORWARD PASS (VITERBI -Alpha's) %%%%% % Initialize the output vector - pm = zeros(length(states),length(states)); - bm_fw = zeros(length(states),length(states),length(data_in)); + pm = zeros(nStates,nStates); + bm_fw = zeros(nStates,nStates,length(data_in)); + + % first start is evaluated without ISI/ wihout the full Impulse response + % so simply use the constellation here + bm = -(data_in(1) - last_sym).^2 ; + pm = pm + bm; + [alpha(:,1),pm_survivor_fw_idx(:,1)] = max(pm,[],2); + pm = repmat(alpha(:,1).',nStates,1); + bm_fw(:,:,1) = pm; % Forward Recursion (FSM Computation) - for n = 1:length(data_in) + for n = 2:length(data_in) - - bm = abs(data_in(n) - noise_free_received).^2; + bm = -(data_in(n) - noise_free_received).^2 ; pm = pm + bm; - [pm_survivor_fw(:,n),pm_survivor_fw_idx(:,n)] = min(pm,[],2); % choose lowest path metric as new state - pm = repmat(pm_survivor_fw(:,n).',length(states),1); % update pm (chosen state to 2nd dimension -> FROM state) + [alpha(:,n),pm_survivor_fw_idx(:,n)] = max(pm,[],2); % choose lowest path metric as new state + pm = repmat(alpha(:,n).',nStates,1); % update pm (chosen state to 2nd dimension -> FROM state) bm_fw(:,:,n) = bm; end % we can now get the best path as min - best_fw_path = NaN(1,length(data_in)+1); + viterbi_path = NaN(1,length(data_in)); % find ideal trellis path by going through the trellis backwards - [~,best_fw_path(length(data_in)+1)] = min(pm_survivor_fw(:,length(data_in))); + [~,viterbi_path(length(data_in))] = max(alpha(:,length(data_in))); + for n = length(data_in):-1:2 + viterbi_path(n-1) = pm_survivor_fw_idx(viterbi_path(n),n); + end + + VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = first_sym(viterbi_path); - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%%%% BACKWARD PASS %%%%% + %%%%% BACKWARD PASS (Beta's) %%%%% % Initialize the output vector - pm = zeros(length(states),length(states)); - pm_survivor_bw = zeros(length(states),length(data_in)+1); - pm_survivor_bw_idx = zeros(length(states),length(data_in)+1); - bm_bw = zeros(length(states),length(states),length(data_in)+1); + pm = zeros(nStates,nStates); + beta = zeros(nStates,length(data_in)); + pm_survivor_bw_idx = zeros(nStates,length(data_in)); + bm_bw = zeros(nStates,nStates,length(data_in)); % starting with the state that has the lowest sum path % metric, follow the stored information about the % predecessor - for h = length(data_in):-1:1 + for h = length(data_in)-1:-1:1 - best_fw_path(h) = pm_survivor_fw_idx(best_fw_path(h+1),h); - - bm = abs(data_in(h) - noise_free_received).^2; + bm = -(data_in(h+1) - noise_free_received).^2 ; pm = pm + bm.'; - [pm_survivor_bw(:,h),pm_survivor_bw_idx(:,h)] = min(pm,[],2); % choose lowest path metric as new state - pm = repmat(pm_survivor_bw(:,h).',length(states),1); % update pm (chosen state to 2nd dimension -> FROM state) + [beta(:,h),pm_survivor_bw_idx(:,h)] = max(pm,[],2); % choose lowest path metric as new state + pm = repmat(beta(:,h).',nStates,1); % update pm (chosen state to 2nd dimension -> FROM state) bm_bw(:,:,h) = bm; end - VITERBI_ESTIMATION_IDX(1:length(data_in)) = first_sym(best_fw_path(2:end)); - VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = constellation(best_fw_path(2:end)); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% FORWARD PASS PAM 2,4,8 (Combine Alpha and Beta to yield LLP's) %%%%% - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%%%% FORWARD PASS %%%%% - - %calc the log probabilities (llp's) + %calc the log probabilities (llp's) - for k = 2:length(data_in) - - fsm = repmat(pm_survivor_fw(:,k-1)',[length(states),1]); %pm_survivor_fw size: length(states)xlength(sequence) - bm = bm_fw(:,:,k); %bm_fw size: length(states)xlength(states)xlength(sequence) - bsm = repmat(pm_survivor_bw(:,k)',[length(states),1]); %pm_survivor_bw size: length(states)xlength(sequence)+1 + for k = 1:length(data_in) - llp(:,k) = min(fsm+bm+bsm,[],2); - - end + if k == 1 - % Compute soft output PAM4 stream from the metric_sym - soft_output = zeros(length(data_in),1); % Expected symbol value per stage - symbol_prob = zeros(length(data_in), length(states)); % Store full probability distribution (optional) - llp = llp'; - for n = 1:length(data_in) - metrics = llp(n, :); % A posteriori metric for each PAM4 candidate - % For numerical stability, subtract the maximum metric before exponentiating - maxMetric = min(metrics); - expMetrics = exp(metrics - maxMetric); - probs = expMetrics / sum(expMetrics); % Normalize to get probabilities - symbol_prob(n, :) = probs; % (Optional) store distribution for analysis - % Compute the soft output as the expected value of the PAM4 symbols - soft_output(n) = sum(probs .* constellation'); - end + alpha_ = repmat(alpha(:,k)',[nStates,1])'; + beta_ = beta(:,k); - % Number of symbols and bits per symbol - num_symbols = constellation; - num_bits = 2; % 2 bits per symbol + LLP(:,k) = max(alpha_ + beta_,[],2); - bit_mapping = PAMmapper(4,0,"eth_style",1).showBitMapping; % Each row corresponds to the symbol above + else - % Initialize LLR storage - llr = zeros(num_bits, length(data_in)); + alpha_ = repmat(alpha(:,k-1)',[nStates,1])'; + gamma_ = bm_fw(:,:,k)'; + beta_ = beta(:,k); - % Compute bit-wise LLRs - for bit_idx = 1:num_bits - % Find indices where bit is 0 and where it is 1 - idx_bit_0 = find(bit_mapping(:,bit_idx) == 0); - idx_bit_1 = find(bit_mapping(:,bit_idx) == 1); - - % Sum over log-probabilities (Max-Log approximation: using min instead of sum) - llr(:,bit_idx) = min(llp(:,idx_bit_1), [], 1) - min(llp(:,idx_bit_0), [], 1); - - end - - % Convert LLR values to a hard-decision bit stream - bit_stream = llr < 0; - [~,~,ber_llr,~] = calc_ber(bit_stream',tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('LLR BER : %.2e \n',ber_llr); - - % directly decide based on lowest LLP index - [~,llp_based_state_seq]=min(llp); - LLP_EST(1:length(data_in)) = constellation(llp_based_state_seq); - rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(LLP_EST'); - [~,~,ber_llp,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('LLP BER : %.2e \n',ber_llp); - % [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % fprintf('LLP BER +1: %.2e \n',ber_llp); - % [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % fprintf('LLP BER -1: %.2e \n',ber_llp); - - %%%% DECIDE based on Viterbi traceback - rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(VITERBI_ESTIMATION_SYMBOLS'); - [~,~,ber_viterbi,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('Viterbi BER: %.2e \n',ber_viterbi); - % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % fprintf('Viterbi BER: %.2e \n',ber_viterbi); - - % directly decide based on the FW path metrics - [~,fw_direct_state_seq]=min(pm_survivor_fw); - FW_EST(1:length(data_in)) = constellation(fw_direct_state_seq); - rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(FW_EST'); - [~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('FW BER: %.2e \n',ber_fw); - - % directly decide based on the BW path metrics - [~,bw_direct_state_seq]=min(pm_survivor_bw); - BW_EST(1:length(data_in)) = constellation(bw_direct_state_seq(2:end)); - rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(BW_EST'); - [~,~,ber_bw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('BW BER: %.2e \n',ber_bw); - % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % fprintf('BW BER: %.2e \n',ber_viterbi); - % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % fprintf('BW BER: %.2e \n',ber_viterbi); - - PAMmapper(4,0,"eth_style",1).showBitMapping - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - tx_symbolpos = zeros(numel(constellation),length(data_ref)); - est_symbolpos = zeros(numel(constellation),length(data_ref)); - for lvl = 1:numel(constellation) - tx_symbolpos(lvl,data_ref==constellation(lvl)) = 1; - est_symbolpos(lvl,VITERBI_ESTIMATION_IDX==first_sym(lvl)) = 1; - est_symbolpos_llm(lvl,llp_based_state_seq==lvl) = 1; - end - - est_symbolpos= logical(est_symbolpos); - tx_symbolpos = logical(tx_symbolpos); - est_symbolpos_llm = logical(est_symbolpos_llm); - - figure(299);clf;hold on; - llp_filt = NaN(length(states),length(data_in)); - llp_ = llp - min(llp,[],1); - for c = 2%:numel(constellation) - for c2 = 1:3%numel(constellation) - - idx = est_symbolpos_llm(c,:); - llp_filt(c2,idx) = (llp(c2,idx)); - - plotidx = 2:200000; - scatter(plotidx,llp_filt(c2,plotidx),1,'o','LineWidth',1,'DisplayName',sprintf('LLP of Symbol %d | %d was transmitted',c2,c)); + LLP(:,k) = max(alpha_ + gamma_,[],1) + beta_'; + end end - end - % Assume llp is a 4 x N matrix (4 PAM-4 symbols, N time steps) - [numSymbols, numTimeSteps] = size(llp); - P_symbols = zeros(numSymbols, numTimeSteps); % To hold the soft output probabilities + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% FORWARD PASS PAM2,4,8 %%%%% - for k = 1:numTimeSteps - - % Compute the exponentials. (Use -llp_shifted if llp's are costs.) - exponents = -llp_(:, k); - - % Normalize to form a probability vector. - P_symbols(:, k) = exponents / sum(exponents); - end - - % Optionally, if you prefer a single soft output value per symbol (an expectation), - % define your PAM-4 constellation levels, e.g.: - soft_output = sum(P_symbols .* constellation, 1); % 1 x N vector of soft outputs + nml_LLP = LLP - max(LLP); + expLLP = exp(nml_LLP); %subtract highst value for better numerical stability, LLP's are not always close to zero + state_prob = expLLP ./ sum(expLLP); - %%%% BITWISE LLR's ??? %%%%% - figure(100); - clf - hold on - title('LLR between inner and outer bits') - bitmapping = PAMmapper(numel(constellation),0).demap(constellation); - for bp = 1:size(bitmapping,2) - pos_bitone = find(bitmapping(:,bp)==1); - pos_bitzero = find(bitmapping(:,bp)~=1); + if obj.M == 6 - llr_bits(bp,:) = min(llp(pos_bitone,:),[],1) - min(llp(pos_bitzero,:),[],1); + num_bits = 5; - subplot(size(bitmapping,2),1,bp) - scatter(1:length(data_ref),llr_bits(bp,:),1,'.'); - end - %%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % all possible transitions (for now 36, including the "edges" + % of the QAM 32 constellation) + pam6transitions = combvec(obj.trellis_states,obj.trellis_states)'; % pam6transitions = + % [-5 -5; + % -3 -5; + % -1 -5; ... - % From: Log-Likelihood Probabilities LLP --> To: Log-Likelihood Ratios LLR - for s = 1:length(states)-1 - llr(s,:) = llp_(s,:) - llp_(s+1,:); % subtract llp's of successive const. points to get llr's - end - - % Build Sets of LLR's that contain only those values at - % timepoint k, where the symbols: - % a) were actually transmitted: set "S" - % b) were decoded by viterbi: set "SC" - S = NaN(length(states)-1,length(data_in)); - SC = NaN(length(states)-1,length(data_in)); + pam6bits = PAMmapper(6,0,"eth_style",0).demap(reshape(pam6transitions',[],1)./sqrt(10)); - llp_inf = llp_; - llp_inf(llp_==0) = Inf; - [~,scnd_idx] = min(llp_inf,[],1) ; - for s = 1:length(states)-1 + pam6bits = reshape(pam6bits',5,[])'; + %pam6bits = + % 0 0 0 1 0 + % 0 0 0 1 0 + % 0 0 0 1 1 + % 1 0 0 1 1 + % 1 0 0 1 0 + % 1 0 0 1 0 + % 0 0 1 1 0 + % .... - %%% SET "S" - % find all time indices k, where const point s was actually transmitted - indice = tx_symbolpos(s,:)==1; - S(s,indice) = llr(s,indice); - %calc mean of all llr's where symbol s was transmitted - K(s,1) = mean(S(s,indice),'omitnan'); - - % find all time indices k, where const. point s+1 was actually transmitted - indice_plusone = tx_symbolpos(s+1,:)==1; - S(s,indice_plusone) = llr(s,indice_plusone); - K(s,2) = mean(S(s,indice_plusone),'omitnan'); + [ok1, idx_bit_0] = ismember(pam6transitions(:,1), obj.trellis_states); + [ok2, idx2] = ismember(pam6transitions(:,2), obj.trellis_states); + assert(all(ok1)&all(ok2), 'Some transition amplitude not found in trellis_states') + pam6ind = [idx_bit_0, idx2]; - - %%% SET "SC" - % find all time indices k, where symbol s was decoded - % idx: find positions in time where a symbol s was decoded - idx = est_symbolpos_llm(s,:); - % idx 2: find positions where the strongest competitor is - % s+1, i.e. the second best llp is at s+1 - idx2 = scnd_idx == s+1; - SC(s,idx&idx2) = llr(s,idx&idx2); - KC(s,1) = mean(SC(s,idx&idx2),'omitnan'); + tx_bits_pam6_reshaped = reshape(tx_bits,5,[])'; - % find all time indices k, where symbol s+1 was decoded - % idx: find positions in time where a symbol s+1 was decoded - idx = est_symbolpos_llm(s+1,:); - idx2 = scnd_idx == s; - % idx 2: find positions where the strongest competitor is - % s, i.e. the second best llp is at s - SC(s,idx&idx2) = llr(s,idx&idx2); - KC(s,2) = mean(SC(s,idx&idx2),'omitnan'); + numPairs = floor(size(LLP,2)/2); + LLR_exact = zeros(numPairs,5); + LLR_maxlogmap = zeros(numPairs,5); - % scale the sets, using the average (?) of the - llrcn(s,:) = SC(s,:) * ( constellation(s+1)-constellation(s) ) ./ ( K(s,2)-K(s,1)) + decisionLevels(s); + for k = 1:numPairs + symbol1 = 2*k-1; + symbol2 = 2*k; - end - - figure(111) - title("Bla") - clf - for s = 1:length(states)-1 - figure(1111) - hold on - title(sprintf('llrcn = llp %d - llp %d',s, s+1,s, s+1)); - scatter(1:length(llrcn),llrcn(s,:),1,'.'); + LLP1 = LLP (:,symbol1); + LLP2 = LLP (:,symbol2); + prob1 = state_prob(:,symbol1); + prob2 = state_prob(:,symbol2); + % 36 joint‐metrics M = log P(i)*P(j) = L1(i)+L2(j) + Mij = LLP1(pam6ind(:,1)) + LLP2(pam6ind(:,2)); + pij = prob1(pam6ind(:,1)) .* prob2(pam6ind(:,2)); - % STUFENARTIGE LLR'S UND LLP'S %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - figure(111) - - subplot(1,length(states),s) - hold on - title(sprintf('llr = llp %d - llp %d \n SC=llr(tx sym = %d or %d)',s, s+1,s, s+1)); - scatter(1:length(llr),llr(s,:),1,'.'); - scatter(1:length(SC),SC(s,:),1,'.'); - yline([K(s,1),K(s,2)]); - low = min(min(llr)); - hi = max(max(llr)); - ylim([low,hi]); + % now for each of the 5 bits do exact-LLR or max-log + for b = 1:num_bits + idx_bit_0 = pam6bits(:,b)==1; + idx_bit_1 = pam6bits(:,b)==0; - subplot(1,length(states),length(states)) - hold on - scatter(1:length(llr),llr(s,:),1,'.'); - ylim([low,hi]); + %--- exact LLR from probabilities + P1 = sum(pij(idx_bit_0)); + P0 = sum(pij(idx_bit_1)); + LLR_exact(k,b) = log(P1./P0); + %--- max-log: + LLR_maxlogmap(k,b) = max( Mij(idx_bit_0) ) - max( Mij(idx_bit_1) ); % N x num_bits + end + end - % HISTOGRAMME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - figure(112) - subplot(length(states),1,s) - hold on - title(sprintf('histogram of llr; between const point %d - %d',s, s+1)); - histogram(llrcn(s,:),10000,'EdgeAlpha',0) - % histogram(SC(s,:),10000,'EdgeAlpha',0) - xlim([-20, 20]); - subplot(length(states),1,length(states)) - hold on - histogram(llrcn(s,:),10000,'EdgeAlpha',0) - % histogram(SC(s,:),10000,'EdgeAlpha',0) - % xlim([-20, 20]); + MI = zeros(1, num_bits); + for k = 1:num_bits + + idx_bit_1 = (tx_bits_pam6_reshaped(:,k) == 0); %wo sind die 1en + idx_bit_0 = (tx_bits_pam6_reshaped(:,k) == 1); %wo sind die 0en + + %LLR's for all actually transmitted ones or zeros + llr0 = LLR_exact(idx_bit_1,k); + llr1 = LLR_exact(idx_bit_0,k); + + % Calculate mutual information for bit position k + I0 = mean(log2(1 + exp(llr0))); % exp(--LLR) = exp(positive) > 1 + I1 = mean(log2(1 + exp(-llr1))); % exp(-+LLR) = exp(negative) < 1 + MI(k) = 1 - 0.5 * (I0 + I1); + end + + GMI = sum(MI); % Total mutual information per symbol + GMI = GMI/2; + + else + + % Number of symbols and bits per symbol + num_bits = log2(length(obj.trellis_states)); % 2 bits per symbol + bit_mapping = PAMmapper(length(obj.trellis_states),0,"eth_style",0).demap(first_sym./rms(first_sym)); + + % Initialize LLR storage + LLR_maxlogmap = zeros(length(data_in),num_bits); + LLR_exact = zeros(length(data_in),num_bits); + + % Compute bit-wise LLRs + for bit_idx = 1:num_bits + + % Find indices where bit is 0 and where it is 1 + idx_bit_0 = bit_mapping(:,bit_idx) == 0; + idx_bit_1 = bit_mapping(:,bit_idx) == 1; + + % Sum over log-probabilities (Max-Log approximation: using max instead of sum) + LLR_maxlogmap(:,bit_idx) = max(LLP(idx_bit_1,:), [], 1) - max(LLP(idx_bit_0,:), [], 1); + + % Sum probabilities over states for which the bit is 1 and 0, respectively. + P0 = sum(state_prob(idx_bit_0, :),1); + P1 = sum(state_prob(idx_bit_1, :),1); + LLR_exact(:,bit_idx) = log(P1./P0); % N x num_bits + + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% CALC NGMI %%%%% + + MI = zeros(1, num_bits); + LLR_exact = LLR_exact; + for k = 1:num_bits + + idx_bit_1 = (tx_bits(:,k) == 0); %wo sind die 1en + idx_bit_0 = (tx_bits(:,k) == 1); %wo sind die 0en + + %LLR's for all actually transmitted ones or zeros + llr0 = LLR_exact(idx_bit_1,k); + llr1 = LLR_exact(idx_bit_0,k); + + % Calculate mutual information for bit position k + I0 = mean(log2(1 + exp(llr0(100:end-100)))); % exp(--LLR) = exp(positive) > 1 + I1 = mean(log2(1 + exp(-llr1(100:end-100)))); % exp(-+LLR) = exp(negative) < 1 + MI(k) = 1 - 0.5 * (I0 + I1); + end + + GMI = sum(MI); % Total mutual information for 2 symbols end - softdecisions = mean(llrcn,1,'omitnan'); + VITERBI_ESTIMATION_SYMBOLS = VITERBI_ESTIMATION_SYMBOLS./rms(VITERBI_ESTIMATION_SYMBOLS); - distance = abs(VITERBI_ESTIMATION_SYMBOLS-llrcn); - [win_cost,win_idx] = min(distance,[],1); - for i = 1:length(data_in) - soft_decisions(i) = llrcn(win_idx(i),i); + debug = 0; + if debug + %%% DEBUG PLOT LIKELIHOOD RATIOS %%% + figure(115);clf + subplot(2,1,1) + for bit = 1:num_bits + hold on; + histogram(LLR_exact(:,bit),1000,"DisplayName",sprintf('Actual LLR of Bit Pos %d',bit),'LineStyle','none','FaceAlpha',0.4); + end + legend + + subplot(2,1,2) + for bit = 1:num_bits + hold on; + histogram(LLR_maxlogmap(:,bit),1000,"DisplayName",sprintf('Max Log LLR of Bit Pos %d',bit),'LineStyle','none','FaceAlpha',0.4); + end + legend end - soft_decisions = max(llrcn,[],1); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% CHECK BER's %%%%% - showLevelHistogram(soft_decisions,data_ref) + if debug + tx_bits = reshape(tx_bits',[],1); - figure() - % scatter(1:length(data_in),VITERBI_ESTIMATION_SYMBOLS,1,'.'); - scatter(1:length(data_in),soft_decisions,1,'.'); + % DECIDE based on Viterbi traceback + VITERBI_ESTIMATION_SYMBOLS = VITERBI_ESTIMATION_SYMBOLS./rms(VITERBI_ESTIMATION_SYMBOLS); + rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(VITERBI_ESTIMATION_SYMBOLS'); + rx_bits = reshape(rx_bits',[],1); + [~,numErr,ber_viterbi,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('Viterbi BER: %.2e \n',ber_viterbi); + fprintf('Viterbi Errors = %d\n', numErr); + % Convert LLR values to a hard-decision bit stream + bit_stream = LLR_maxlogmap > 0; %ratio separates lower or higher than =0 -> simply decode for the negative values + bit_stream = reshape(bit_stream',[],1); + [~,~,ber_llr,~] = calc_ber(bit_stream,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('LLR MaxLogMAP BER : %.2e \n',ber_llr); - MLM_ESTIMATION(1:length(data_in)) = PAMmapper(numel(constellation),0).quantize(soft_decisions)'; - rx_bits = PAMmapper(numel(constellation),0).demap(MLM_ESTIMATION'); - [~,~,ber_mlm,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - fprintf('MLM BER: %.2e \n',ber_mlm); + % Convert LLR values to a hard-decision bit stream + bit_stream = LLR_exact > 0; %ratio separates lower or higher than =0 -> simply decode for the negative value + bit_stream = reshape(bit_stream',[],1); + [~,~,ber_llr,~] = calc_ber(bit_stream,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('LLR BER : %.2e \n',ber_llr); + + % DECIDE based on lowest LLP index and check BER + [~,llp_based_state_seq]=max(LLP); + LLP_EST(1:length(data_in)) = first_sym(llp_based_state_seq); + LLP_EST = LLP_EST./rms(LLP_EST); + rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(LLP_EST'); + rx_bits = reshape(rx_bits',[],1); + [~,~,ber_llp,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('LLP BER : %.2e \n',ber_llp); + + % directly decide based on the FW path metrics + [~,fw_direct_state_seq]=max(alpha); + FW_EST(1:length(data_in)) = first_sym(fw_direct_state_seq); + FW_EST = FW_EST./rms(FW_EST); + rx_bits = PAMmapper(obj.M,0,"eth_style",0).demap(FW_EST'); + rx_bits = reshape(rx_bits',[],1); + [~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('FW BER: %.2e \n',ber_fw); + + end end end diff --git a/Classes/DataBaseHandler/DBHandler.m b/Classes/DataBaseHandler/DBHandler.m index 926d016..73a4321 100644 --- a/Classes/DataBaseHandler/DBHandler.m +++ b/Classes/DataBaseHandler/DBHandler.m @@ -5,7 +5,7 @@ classdef DBHandler < handle properties conn % Database connection object - pathToDB % Path to the SQLite database + dataBase % Path to the SQLite database tableNames % Cell array containing names of all tables in the database tables = struct(); % Structure containing MATLAB tables for each database table distinctValues @@ -21,7 +21,7 @@ classdef DBHandler < handle % obj = DBHandler('pathToDB', 'path/to/database.db'); arguments - options.pathToDB = ""; % Default value for pathToDB if not provided + options.dataBase = ""; % Default value for pathToDB if not provided options.type = "mysql"; end @@ -37,12 +37,17 @@ classdef DBHandler < handle try if options.type == "sqlite" - obj.conn = sqlite(obj.pathToDB); + + obj.conn = sqlite(obj.dataBase); + elseif options.type == "mysql" - datasource = "jdbc:mysql://134.245.243.254:3306/labor"; + + % datasource = "jdbc:mysql://134.245.243.254:3306/labor"; + + obj.conn = database( ... - "labor", ... % Database name + string(obj.dataBase), ... % Database name "silas", ... % Username "silas", ... % Password (or getSecret) "Vendor", "MySQL", ... @@ -68,8 +73,10 @@ classdef DBHandler < handle function obj = refresh(obj) % Get table names and the first rows of each table to understand the structure + warning off obj.getTableNames(); obj.getTables(); + warning on % obj.getDistinctValues(); end @@ -201,20 +208,24 @@ classdef DBHandler < handle if isstruct(newRow) fields = fieldnames(newRow); - % Handle empty fields - emptyFields = structfun(@isempty, newRow); - if sum(emptyFields) > 0 - emptyFieldNames = fields(emptyFields); - for idx = 1:numel(emptyFieldNames) - newRow.(emptyFieldNames{idx}) = NaN; - end - end + % % Handle empty fields + % emptyFields = structfun(@isempty, newRow); + % if sum(emptyFields) > 0 + % emptyFieldNames = fields(emptyFields); + % for idx = 1:numel(emptyFieldNames) + % newRow.(emptyFieldNames{idx}) = NaN; + % end + % end - % Convert arrays to JSON strings - for i = 1:length(fields) - fieldValue = newRow.(fields{i}); - if isnumeric(fieldValue) && length(fieldValue) > 1 - newRow.(fields{i}) = jsonencode(fieldValue); + % Convert any non-scalar numeric (including [] and vectors) into JSON + for i = 1:numel(fields) + name = fields{i}; + value = newRow.(name); + + if isnumeric(value) && ~isscalar(value) + % jsonencode([]) -> "[]" + % jsonencode([a,b,c]) -> "[a,b,c]" + newRow.(name) = jsonencode(value); end end @@ -238,6 +249,7 @@ classdef DBHandler < handle elseif ischar(value) newRow.(colName{1}) = string(value); end + end % Parameters for retry logic @@ -251,7 +263,7 @@ classdef DBHandler < handle if obj.type == "mysql" newRow_ = obj.convertTableToCellStrings(newRow); end - sqlwrite(obj.conn, tableName, newRow); + sqlwrite(obj.conn, tableName, newRow,"Catalog",obj.dataBase); success = true; catch e if contains(e.message, 'database is locked') || contains(e.message, 'cannot rollback transaction') diff --git a/Classes/DataBaseHandler/Metricstruct.m b/Classes/DataBaseHandler/Metricstruct.m index 9161e47..403e90c 100644 --- a/Classes/DataBaseHandler/Metricstruct.m +++ b/Classes/DataBaseHandler/Metricstruct.m @@ -1,6 +1,6 @@ classdef Metricstruct % ResultData - Class to store and manage metric data from signal processing results - + properties result_id (1,1) double {mustBeNumeric} = NaN run_id (1,1) double {mustBeNumeric} = NaN @@ -24,15 +24,15 @@ classdef Metricstruct GMI (1,1) double {mustBeNumeric} = NaN AIR (1,1) double {mustBeNumeric} = NaN - Alpha (:,1) double {mustBeNumeric} = [] + Alpha (:,1) double {mustBeNumeric} = NaN MLSE_dir (:,1) double {mustBeNumeric} = [] end - + methods function obj = Metricstruct(varargin) % Constructor method for ResultData % Can be called empty or with name-value pairs - + % Process name-value pairs if provided if nargin > 0 for i = 1:2:nargin @@ -44,7 +44,7 @@ classdef Metricstruct end end end - + function s = toStruct(obj) % Convert the object to a struct s = struct(); @@ -53,7 +53,7 @@ classdef Metricstruct s.(props{i}) = obj.(props{i}); end end - + function str = toString(obj) % Convert the object to a formatted string s = obj.toStruct(); @@ -70,9 +70,62 @@ classdef Metricstruct end end end + + function print(obj) + % Print method to display key metrics in a formatted way + + % Define the width for formatting + nameWidth = 15; % Width for parameter names + valueWidth = 12; % Width for values + + % Print header + fprintf('\n%s\n', repmat('=', 1, nameWidth + valueWidth)); + fprintf(' Results \n'); + fprintf('%s\n', repmat('=', 1, nameWidth + valueWidth)); + + % Function to format numbers with appropriate precision + function str = formatNumber(value) + if isnan(value) + str = 'N/A'; + elseif abs(value) < 0.1 && value ~= 0 + str = sprintf('%.2e', value); + else + str = sprintf('%.4f', value); + end + end + + % Print each metric + % BER metrics + fprintf('%-*s: %*s\n', nameWidth, 'BER', valueWidth, formatNumber(obj.BER)); + fprintf('%-*s: %*s\n', nameWidth, 'BER (precoded)', valueWidth, formatNumber(obj.BER_precoded)); + + % SNR + fprintf('%-*s: %*s\n', nameWidth, 'SNR', valueWidth, formatNumber(obj.SNR)); + + % Information rates + fprintf('%-*s: %*s\n', nameWidth, 'GMI', valueWidth, formatNumber(obj.GMI)); + fprintf('%-*s: %*s GBd \n', nameWidth, 'AIR', valueWidth, formatNumber(obj.AIR.*1e-9)); + + % Alpha (if it's not empty, show first few values) + if ~isempty(obj.Alpha) + if length(obj.Alpha) <= 3 + alphaStr = sprintf('%.4f ', obj.Alpha); + else + alphaStr = sprintf('%.4f %.4f %.4f...', obj.Alpha(1:3)); + end + fprintf('%-*s: %s\n', nameWidth, 'Alpha', alphaStr); + else + fprintf('%-*s: %*s\n', nameWidth, 'Alpha', valueWidth, '[]'); + end + + % Print footer + fprintf('%s\n', repmat('=', 1, nameWidth + valueWidth)); + end + end - + methods (Static) + function obj = fromStruct(s) % Create a ResultData object from a struct obj = Metricstruct(); @@ -83,5 +136,6 @@ classdef Metricstruct end end end + end end \ No newline at end of file diff --git a/Classes/Warehouse_class/classes/DataStorage.m b/Classes/Warehouse_class/classes/DataStorage.m index 91530cf..f057c54 100644 --- a/Classes/Warehouse_class/classes/DataStorage.m +++ b/Classes/Warehouse_class/classes/DataStorage.m @@ -150,7 +150,7 @@ classdef DataStorage < handle end value{i} = tmp{1} ; else - value{i} = tmp{1} ; + value{i} = tmp ; end else diff --git a/Datatypes/adaption_method.m b/Datatypes/adaption_method.m new file mode 100644 index 0000000..5b4c12c --- /dev/null +++ b/Datatypes/adaption_method.m @@ -0,0 +1,9 @@ +classdef adaption_method < int32 + + enumeration + lms (1) + nlms (2) + rls (3) + end + +end \ No newline at end of file diff --git a/Functions/EQ_structures/dsp_runid.m b/Functions/EQ_structures/dsp_runid.m index c3d6a63..dbb306b 100644 --- a/Functions/EQ_structures/dsp_runid.m +++ b/Functions/EQ_structures/dsp_runid.m @@ -5,9 +5,11 @@ arguments options.append_to_db = 0; options.max_occurences = 4; options.parameters = struct(); - options.database_path - options.database_name + options.database_type + options.dataBase + options.load_file_path = struct(); options.storage_path + options.mode end % Initialize output structures @@ -17,34 +19,72 @@ output.vnle_package = {}; output.dbtgt_package = {}; output.dbenc_package = {}; -% Initialize database connection -database = DBHandler("pathToDB", [options.database_path, options.database_name], "type", "sqlite"); -dataTable = queryRunid(run_id, database); -fsym = dataTable.symbolrate; -M = double(dataTable.pam_level); -duob_mode = db_mode(dataTable.db_mode); - -if database.checkIfRunExists('Results','run_id',run_id) - disp(['Already got at least one reulst for run id: ',num2str(run_id),' ']) - return +if options.mode == "load_run_id" || options.append_to_db + % Initialize database connection + database = DBHandler("dataBase", [options.dataBase], "type", options.database_type ); end +if options.mode == "load_run_id" + + dataTable = queryRunid(run_id, database); + fsym = dataTable.symbolrate; + M = double(dataTable.pam_level); + duob_mode = db_mode(strrep(dataTable.db_mode,'"','')); + + % if database.checkIfRunExists('Results','run_id',run_id) + % disp(['Already got at least one reulst for run id: ',num2str(run_id),' ']) + % return + % end + + % Load and Sync signal data from DB + [Tx_bits, Symbols, Scpe_cell, ~] = loadAndSyncSignalDataFromDb(dataTable, options); + +elseif options.mode == "load_files" + + Tx_bits = load(options.load_file_path.tx_bits_path); + Symbols = load(options.load_file_path.tx_symbols_path); + Scpe_sig_raw = load(options.load_file_path.rx_raw_path); + + Tx_bits = Tx_bits.Bits; + Symbols = Symbols.Symbols; + Scpe_sig_raw = Scpe_sig_raw.Scpe_sig_raw; + + fsym = Symbols.fs; + M = Symbols.logbook.ModifierCopy{1}.M; + duob_mode = Symbols.logbook.ModifierCopy{1}.duobinary_mode; + + Scpe_sig_resampled = Scpe_sig_raw.resample("fs_in", Scpe_sig_raw.fs, "fs_out", 2*fsym); + [~, Scpe_cell, ~, found_sync] = Scpe_sig_resampled.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + +else + + % Run quick Simulation + tx_simulation; + +end % Handle Settings and argument replacement len_tr = 4096*2; -ffe_order = [50, 5, 5]; +ffe_order = [50, 5, 5]; dfe_order = [0, 0, 0]; pf_ncoeffs = 1; mu_ffe = [0.0001, 0.0008, 0.001]; mu_dfe = 0.0004; -mu_dc = 0.00; +mu_dc = 0.005; +dc_buffer_len = 1; + +mu_tr = 0; +mu_dd = 0.05; +adaption= 1; +use_dd_mode = 1; use_ffe = 1; +use_dfe = 1; use_vnle_mlse = 1; use_dbtgt = 1; -use_dbenc = 0; +use_dbenc = 1; addProcessingResultToDatabase = 0; @@ -60,74 +100,148 @@ if ~isempty(paramStruct) end % Configure equalizers -eq_lin = EQ("Ne",[50,0,0],"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); -eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); -pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); -mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); -mlse_db_ = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); -eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); -% Duobinary signaling (db encoded) -mlse_db_enc = MLSE_viterbi("DIR", [1,1], "duobinary_output", 0, "M", M, "trellis_states", PAMmapper(M,0).levels); -eq_db_enc = EQ("Ne", ffe_order, "Nb", dfe_order, "training_length", len_tr, ... - "training_loops", 5, "dd_loops", 5, "K", 2, "DCmu", mu_dc, ... - "DDmu", [mu_ffe mu_dfe], "DFEmu", 0.005, "FFEmu", 0, "plotfinal", 0, "ideal_dfe", 1); +options.max_occurences = min(options.max_occurences,length(Scpe_cell)); +for r = 1:options.max_occurences + + %FFE + % eq_dfe = FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption",adaption_method(adaption),"dd_mode",use_dd_mode); -% Load signal data -[Tx_bits, Symbols, Scpe_cell, ~] = loadSignalData(dataTable, options); + % + eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + % mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0,"adaption_technique","lms"); + % Duobinary signaling (db encoded) + mlse_db_enc = MLSE("DIR", [1,1], "duobinary_output", 0, "M", M, "trellis_states", PAMmapper(M,0).levels); + eq_db_enc = EQ("Ne", ffe_order, "Nb", dfe_order, "training_length", len_tr, ... + "training_loops", 5, "dd_loops", 5, "K", 2, "DCmu", mu_dc, ... + "DDmu", [mu_ffe mu_dfe], "DFEmu", 0.005, "FFEmu", 0, "plotfinal", 0, "ideal_dfe", 1); -for r = 1:numel(Scpe_cell) - % Preprocess signal Scpe_sig = preprocessSignal(Scpe_cell{r}, Symbols, fsym); - + if duob_mode ~= db_mode.db_encoded + if use_ffe - ffe_results = ffe(eq_lin,M,Scpe_sig,Symbols,Tx_bits,... + + ffe_order = [50, 0, 0]; + eq_dfe = 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); + + dfe_results = ffe(eq_dfe,M,Scpe_sig,Symbols,Tx_bits,... "precode_mode",duob_mode,... - 'showAnalysis',1,... + 'showAnalysis',0,... "postFFE",[],... "eth_style_symbol_mapping",0); - output.ffe_package{r} = ffe_results; + + output.ffe_package{r} = dfe_results; + + dfe_results.metrics.print; + if options.append_to_db - database.addProcessingResult(run_id, ffe_results.metrics, ffe_results.config); + database.addProcessingResult(run_id, dfe_results.metrics, dfe_results.config); end + end - + + if use_dfe + + ffe_order = [50, 5, 5]; + eq_dfe = EQ("Ne",ffe_order,"Nb",[2,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); + + dfe_results = ffe(eq_dfe,M,Scpe_sig,Symbols,Tx_bits,... + "precode_mode",duob_mode,... + 'showAnalysis',0,... + "postFFE",[],... + "eth_style_symbol_mapping",0); + + output.ffe_package{r} = dfe_results; + + dfe_results.metrics.print; + + if options.append_to_db + database.addProcessingResult(run_id, dfe_results.metrics, dfe_results.config); + end + + end + if use_vnle_mlse + + pf_ncoeffs = 1; + eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + % mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + [ffe_results, mlse_results] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... "precode_mode", duob_mode,... 'showAnalysis', 1, ... "postFFE", [],... "eth_style_symbol_mapping", 0); - output.mlse_package{r} = mlse_results; - output.vnle_package{r} = ffe_results; + + ffe_results.metrics.print; + mlse_results.metrics.print; + + output.mlse_package{r} = mlse_results; + output.vnle_package{r} = ffe_results; + + if options.append_to_db + database.addProcessingResult(run_id, mlse_results.metrics, mlse_results.config); + database.addProcessingResult(run_id, ffe_results.metrics, ffe_results.config); + end + + pf_ncoeffs = 2; + eq_ = EQ("Ne",ffe_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + % mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + + [ffe_results, mlse_results] = vnle_postfilter_mlse(eq_, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... + "precode_mode", duob_mode,... + 'showAnalysis', 1, ... + "postFFE", [],... + "eth_style_symbol_mapping", 0); + + ffe_results.metrics.print; + mlse_results.metrics.print; + + output.mlse_package{r} = mlse_results; + output.vnle_package{r} = ffe_results; + if options.append_to_db database.addProcessingResult(run_id, mlse_results.metrics, mlse_results.config); database.addProcessingResult(run_id, ffe_results.metrics, ffe_results.config); end end - + if use_dbtgt + dbt_results = duobinary_target(eq_, mlse_db_, M, Scpe_sig, Symbols, Tx_bits, ... "precode_mode", duob_mode, ... 'showAnalysis', 0,... "postFFE", []); - output.dbtgt_package{r} = dbt_results; + + output.dbtgt_package{r} = dbt_results; + if options.append_to_db database.addProcessingResult(run_id, dbt_results.metrics, dbt_results.config); end + end end - + if duob_mode == db_mode.db_encoded + db_results = duobinary_signaling(eq_db_enc, mlse_db_enc, M, Scpe_sig, Symbols, Tx_bits, "precode_mode",duob_mode, "showAnalysis",0,"postFFE",[]); output.dbenc_package{r} = db_results; if options.append_to_db database.addProcessingResult(run_id, db_results.metrics, db_results.config); end + end diff --git a/Functions/EQ_structures/duobinary_signaling.m b/Functions/EQ_structures/duobinary_signaling.m index 013e184..e47bc19 100644 --- a/Functions/EQ_structures/duobinary_signaling.m +++ b/Functions/EQ_structures/duobinary_signaling.m @@ -36,20 +36,21 @@ if ~isempty(options.postFFE) end % Process through MLSE -eq_signal = mlse_.process(eq_signal); +% [mlse_signal] = mlse_.process(eq_signal); +[mlse_signal,~,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols); % Apply duobinary encoding and decoding -eq_signal = Duobinary().encode(eq_signal); -eq_signal = Duobinary().decode(eq_signal); +mlse_signal = Duobinary().encode(mlse_signal); +mlse_signal = Duobinary().decode(mlse_signal); % Demap symbols to bits -rx_bits = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).demap(eq_signal); +rx_bits = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).demap(mlse_signal); %% Calculate BER and metrics [bits_db, errors_db, ber_db, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1); -% Calculate performance metrics -[snr, snr_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal); +% Calculate performance metrics after duobinary FFE! +[snr, snr_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal); %SNR of duobinary sequence - not directly comparable to [gmi] = calc_air(eq_signal, tx_symbols, "skip_front", 10000, "skip_end", 10000); air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M)); [evm_total, evm_lvl] = calc_evm(eq_signal, tx_symbols); diff --git a/Functions/EQ_structures/duobinary_target.m b/Functions/EQ_structures/duobinary_target.m index b5ff07c..22257d7 100644 --- a/Functions/EQ_structures/duobinary_target.m +++ b/Functions/EQ_structures/duobinary_target.m @@ -23,8 +23,8 @@ if ~isempty(options.postFFE) end mlse_.DIR = [1,1]; -mlse_sig_sd = mlse_.process(eq_signal); - +% mlse_sig_sd = mlse_.process(eq_signal); +[mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols); mlse_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd); % precoding to mitigate error propagation, most prominently used in @@ -73,6 +73,10 @@ end rx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); [bits_db,errors_db,ber_db,errorIndice_db] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); +alpha = arburg(eq_noise.signal,1);%pf_.coefficients(2); +alpha = alpha(2); +gmi_mlse = GMI_MLSE; +air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_mlse ./ log2(double(M)); db_results = struct(); db_results.metrics = Metricstruct; @@ -85,13 +89,10 @@ db_results.metrics.numBits = bits_db; db_results.metrics.numBitErr = errors_db; db_results.metrics.BER_precoded = ber_db_diff_precoded; db_results.metrics.numBitErr_precoded = errors_db_diff_precoded; -db_results.metrics.SNR = NaN; -db_results.metrics.SNR_level = NaN; -db_results.metrics.GMI = NaN; -db_results.metrics.AIR = NaN; -db_results.metrics.EVM = NaN; -db_results.metrics.EVM_level = NaN; +db_results.metrics.GMI = gmi_mlse; +db_results.metrics.AIR = air_mlse; db_results.metrics.MLSE_dir = mlse_.DIR; +db_results.metrics.Alpha = alpha; % Create DB results structure db_results.config = Equalizerstruct(); @@ -99,7 +100,7 @@ eq_.e = []; eq_.e2 = []; eq_.e3 = []; db_results.config.eq = jsonencode(eq_); -mlse_.DIR = []; +% mlse_.DIR = []; db_results.config.mlse = jsonencode(mlse_); db_results.config.equalizer_structure = int32(equalizer_structure.vnle_db_mlse); db_results.config.comment = 'function: Duobinary tgt. (VNLE -> MLSE)'; diff --git a/Functions/EQ_structures/ffe.m b/Functions/EQ_structures/ffe.m index f9ff57f..e3ecc04 100644 --- a/Functions/EQ_structures/ffe.m +++ b/Functions/EQ_structures/ffe.m @@ -36,6 +36,11 @@ if ~isempty(options.postFFE) toc end +try + ch_coefficients = arburg(eq_noise.signal,1); + channel_alpha = ch_coefficients(2); +end + % Hard decision on FFE output eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); @@ -66,10 +71,15 @@ end % Create FFE results structure ffe_results = struct(); +try + eq_.e = []; + eq_.e2 = []; + eq_.e3 = []; + eq_.b = []; + eq_.b2 = []; + eq_.b3 = []; +end -eq_.e = []; -eq_.e2 = []; -eq_.e3 = []; ffe_results.config = Equalizerstruct(); ffe_results.config.eq = jsonencode(eq_); ffe_results.config.equalizer_structure = int32(equalizer_structure.ffe); @@ -95,7 +105,7 @@ ffe_results.metrics.GMI = gmi; ffe_results.metrics.AIR = air; ffe_results.metrics.EVM = evm_total; ffe_results.metrics.EVM_level = evm_lvl; - +ffe_results.metrics.Alpha = channel_alpha; end @@ -140,24 +150,40 @@ end end function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, postFFE) -% Display analysis plots and metrics -figure(336); -showEQNoisePSD(eq_noise, "fignum", 336, "displayname", 'Residual Noise after FFE'); + % Display analysis plots and metrics -if ~isempty(postFFE) - showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 338); -end + % Initialize figure handles + % Corrected line - added tx_symbols as second positional argument + showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 100); -showEQfilter(eq_.e, eq_signal_sd.fs.*2); + warning off + showLevelScatter(eq_signal_sd, tx_symbols, "fignum", 101); + figure(gcf);hold on; plot(((1:length(eq_noise.signal)) / eq_noise.fs) * 1e6,movmean(eq_noise.signal,2000,1), 'LineWidth',3,'Color','black') + warning on -figure(341); clf; -showLevelHistogram(eq_signal_sd, tx_symbols, "fignum", 341); + showLevelHistogram(eq_signal_sd, tx_symbols, "fignum", 102); -warning off -figure(400); clf; -showLevelScatter(eq_signal_sd, tx_symbols, "fignum", 400); + showEQNoisePSD(eq_noise, "fignum", 103, "displayname", 'Residual Noise after FFE'); + + % Figure 2: Post-FFE coefficients (if available) + if ~isempty(postFFE) + showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 104); + end + + try + figure(339); + showEQfilter(eq_.e_tr, eq_signal_sd.fs.*2,"displayname",'training','fignum',339); + showEQfilter(eq_.e, eq_signal_sd.fs.*2,"displayname",'dec. directed','fignum',339); + legend on + end + + try + figure(240); hold on; plot(pow2db(movmean(eq_.debug_struct.error_tr',100)));ylim([-30,3]);title('error training'); + + figure(241); hold on; plot(pow2db(movmean(eq_.debug_struct.update_tr',100)));title('update step training'); + + figure(242); hold on; plot(pow2db(movmean(eq_.debug_struct.update',1000)));title('update step dd'); + end + % eq_signal_sd.eye(eq_signal_sd.fs,M,"displayname",'Eye','fignum',105); -figure(401); clf; -showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 401); -warning on end \ No newline at end of file diff --git a/Functions/EQ_structures/vnle_postfilter_mlse.m b/Functions/EQ_structures/vnle_postfilter_mlse.m index 90691d8..6610221 100644 --- a/Functions/EQ_structures/vnle_postfilter_mlse.m +++ b/Functions/EQ_structures/vnle_postfilter_mlse.m @@ -45,9 +45,10 @@ end eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); % Process through postfilter and MLSE -mlse_sig_sd = pf_.process(eq_signal_sd, eq_noise); +[mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); mlse_.DIR = pf_.coefficients; -mlse_sig_sd = mlse_.process(mlse_sig_sd); + +[mlse_sig_sd,LLR,GMI_MLSE] = mlse_.process(mlse_sig_sd,tx_symbols); mlse_sig_hd = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).quantize(mlse_sig_sd); %% Calculate BER based on precoding mode @@ -63,7 +64,17 @@ air_vnle = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_vnle ./ log2(dou [std_rxraw_total, std_rxraw_lvl] = calc_std(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols); % MLSE metrics -alpha = pf_.coefficients(2); +alpha = arburg(eq_noise.signal,1);%pf_.coefficients(2); +alpha = alpha(2); +gmi_mlse = GMI_MLSE; +air_mlse = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_mlse ./ log2(double(M)); + +%% Display analysis if requested +if options.showAnalysis + displayAnalysis(eq_noise,whitened_noise, eq_signal_sd, rx_signal, eq_, pf_, mlse_, tx_symbols, M, options.postFFE); +end + + %% Prepare output structures % Determine postFFE order @@ -94,19 +105,21 @@ ffe_results.metrics.GMI = gmi_vnle; ffe_results.metrics.AIR = air_vnle; ffe_results.metrics.EVM = evm_vnle_total; ffe_results.metrics.EVM_level = evm_vnle_lvl; -ffe_results.metrics.Alpha = []; +ffe_results.metrics.Alpha = alpha; +try + eq_.e = []; + eq_.e2 = []; + eq_.e3 = []; +end -% Create FFE results structure - -eq_.e = []; -eq_.e2 = []; -eq_.e3 = []; ffe_results.config = Equalizerstruct(); ffe_results.config.eq = jsonencode(eq_); ffe_results.config.equalizer_structure = int32(equalizer_structure.vnle); ffe_results.config.comment = 'function: vnle_postfilter_mlse - FFE part'; + + mlse_results = struct(); mlse_results.metrics = Metricstruct; mlse_results.metrics.result_id = NaN; @@ -118,12 +131,11 @@ mlse_results.metrics.numBits = numbits.mlse; mlse_results.metrics.numBitErr = errors.mlse; mlse_results.metrics.BER_precoded = bers.mlse_precoded; mlse_results.metrics.numBitErr_precoded = errors.mlse_precoded; -mlse_results.metrics.SNR = NaN; -mlse_results.metrics.SNR_level = NaN; -mlse_results.metrics.GMI = NaN; -mlse_results.metrics.AIR = NaN; -mlse_results.metrics.EVM = NaN; -mlse_results.metrics.EVM_level = NaN; +% mlse_results.metrics.SNR = NaN; +mlse_results.metrics.GMI = gmi_mlse; +mlse_results.metrics.AIR = air_mlse; +% mlse_results.metrics.EVM = NaN; +% mlse_results.metrics.EVM_level = NaN; mlse_results.metrics.Alpha = alpha; mlse_results.metrics.MLSE_dir = mlse_.DIR; @@ -131,15 +143,12 @@ mlse_results.metrics.MLSE_dir = mlse_.DIR; mlse_results.config = Equalizerstruct(); mlse_results.config.eq = jsonencode(eq_); -mlse_.DIR = []; +mlse_.DIR = length(mlse_.DIR)-1; mlse_results.config.mlse = jsonencode(mlse_); mlse_results.config.equalizer_structure = int32(equalizer_structure.vnle_pf_mlse); mlse_results.config.comment = 'function: vnle_postfilter_mlse - MLSE part'; -%% Display analysis if requested -if options.showAnalysis - displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, pf_, mlse_, tx_symbols, M, options.postFFE); -end + end @@ -206,7 +215,7 @@ end end -function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, pf_, mlse_, tx_symbols, M, postFFE) +function displayAnalysis(eq_noise, whitened_noise, eq_signal_sd, rx_signal, eq_, pf_, mlse_, tx_symbols, M, postFFE) % Display analysis plots and metrics figure(336); showEQNoisePSD(eq_noise, "fignum", 336, "displayname", 'Residual Noise after VNLE', 'postfilter_taps', pf_.coefficients); @@ -228,4 +237,9 @@ showLevelScatter(eq_signal_sd, tx_symbols, "fignum", 400); showLevelScatter(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols, "fignum", 401); drawnow; warning on + + +whitened_noise.spectrum("displayname",'after postfilter','fignum',342); +eq_noise.spectrum("displayname",'before postfilter','fignum',342); + end \ No newline at end of file diff --git a/Functions/EQ_visuals/showEQNoisePSD.m b/Functions/EQ_visuals/showEQNoisePSD.m index 51b12d7..a10fb15 100644 --- a/Functions/EQ_visuals/showEQNoisePSD.m +++ b/Functions/EQ_visuals/showEQNoisePSD.m @@ -41,4 +41,7 @@ end % Ensure a legend is displayed legend('show'); end + + xlim([-eq_noise.fs/2* 1e-9 eq_noise.fs/2* 1e-9]); + end diff --git a/Functions/EQ_visuals/showEQcoefficients.m b/Functions/EQ_visuals/showEQcoefficients.m index a844155..cb747b8 100644 --- a/Functions/EQ_visuals/showEQcoefficients.m +++ b/Functions/EQ_visuals/showEQcoefficients.m @@ -53,8 +53,8 @@ function showEQcoefficients(options) for i = 1:numSubplots subplot(1, numSubplots, i); - stem(coeffs{i}, 'Color', options.color, 'LineWidth', 0.1, ... - 'Marker', '.', 'MarkerSize', 1); + stem(coeffs{i}, 'Color', options.color, 'LineWidth', 1, ... + 'Marker', '.', 'MarkerSize', 10); title(titles{i}); ylim([-1, 1]); % Set y-axis limits to [-1, 1] grid on; diff --git a/Functions/EQ_visuals/showEQfilter.m b/Functions/EQ_visuals/showEQfilter.m index 1d13022..3e0f4b0 100644 --- a/Functions/EQ_visuals/showEQfilter.m +++ b/Functions/EQ_visuals/showEQfilter.m @@ -1,5 +1,12 @@ -function showEQfilter(coefficients,fs) +function showEQfilter(coefficients,fs,options) + +arguments + coefficients + fs + options.fignum (1,1) double = NaN % Default to NaN if not provided + options.displayname (1,:) char = '' % Default to an empty string if not provided +end % Assuming that obj.e contains the final FFE filter coefficients. % Set the number of frequency points and sampling frequency. @@ -12,24 +19,28 @@ function showEQfilter(coefficients,fs) half_nfft = floor(nfft/2) + 1; f = f(1:half_nfft); H = H(1:half_nfft); - - % Plot the magnitude and phase responses. - figure(339); + + % Determine the figure number to use or create a new figure + if isnan(options.fignum) + fig = figure; % Create a new figure and get its handle + else + fig = figure(options.fignum); % Use the specified figure number + end % Magnitude response (in dB) subplot(2,1,1); hold on - plot(f.*1e-9, 20*log10(abs(1./H))); + plot(f.*1e-9, 20*log10(abs(1./H)),'DisplayName',options.displayname); title('(Inverted) Magnitude Response of FFE Filter'); - xlabel('Frequency (Hz)'); + xlabel('Frequency (GHz)'); ylabel('Magnitude (dB)'); grid on; % Phase response subplot(2,1,2); - plot(f.*1e-9, unwrap(angle(H))); + plot(f.*1e-9, unwrap(angle(H)),'DisplayName',options.displayname); title('Phase Response of FFE Filter'); - xlabel('Frequency (Hz)'); + xlabel('Frequency (GHz)'); ylabel('Phase'); grid on; diff --git a/Functions/EQ_visuals/showLevelHistogram.m b/Functions/EQ_visuals/showLevelHistogram.m index 437b65c..0f31bdc 100644 --- a/Functions/EQ_visuals/showLevelHistogram.m +++ b/Functions/EQ_visuals/showLevelHistogram.m @@ -20,10 +20,14 @@ end fig = figure(options.fignum); % Use the specified figure number end + eq_signal = max(min(eq_signal,3),-3); + %%% Separate Classes constellation = unique(ref_symbols); received_sd = NaN(numel(constellation),length(ref_symbols)); - lvlcol = cbrewer2('Set1',numel(constellation)); + lvlcol = cbrewer2('Paired',numel(constellation)*2); + lvlcol = lvlcol(2:2:end,:); + % lvlcol = cbrewer2('Set1',numel(constellation)); for lvl = 1:numel(constellation) %Separate the equalized signal into the %respective levels based on the actually @@ -43,6 +47,7 @@ end histogram(received_sd(lvl,:),1000,"EdgeAlpha",0,'DisplayName',['Lvl ',num2str(lvl),' ; ',num2str(cnt(lvl)),' '],'FaceColor',lvlcol(lvl,:),'Normalization','pdf'); warning on end + xlim([-3 3]); legend grid on diff --git a/Functions/EQ_visuals/showLevelScatter.m b/Functions/EQ_visuals/showLevelScatter.m index a5d34ef..70ddb39 100644 --- a/Functions/EQ_visuals/showLevelScatter.m +++ b/Functions/EQ_visuals/showLevelScatter.m @@ -7,7 +7,7 @@ arguments options.f_sym = []; end -plot_shit = 0; +plot_shit = 1; if isa(eq_signal,'Signal') options.f_sym = eq_signal.fs; @@ -25,6 +25,7 @@ if plot_shit fig = figure; % Create a new figure and get its handle else fig = figure(options.fignum); % Use the specified figure number + clf; end end @@ -127,5 +128,6 @@ if plot_shit yline(levels); xlabel('Time in $\mu$s'); ylabel('Normalized Amplitude'); +ylim([-3 3]); end end diff --git a/Functions/Job_Processing/loadSignalData.m b/Functions/Job_Processing/loadAndSyncSignalDataFromDb.m similarity index 71% rename from Functions/Job_Processing/loadSignalData.m rename to Functions/Job_Processing/loadAndSyncSignalDataFromDb.m index cff4ff8..db2812c 100644 --- a/Functions/Job_Processing/loadSignalData.m +++ b/Functions/Job_Processing/loadAndSyncSignalDataFromDb.m @@ -1,4 +1,4 @@ -function [Bits, Symbols, Scpe_cell, found_sync] = loadSignalData(dataTable, options) +function [Bits, Symbols, Scpe_cell, found_sync] = loadAndSyncSignalDataFromDb(dataTable, options) % LOADSIGNALDATA Loads and synchronizes signal data from storage % % Inputs: @@ -12,22 +12,30 @@ function [Bits, Symbols, Scpe_cell, found_sync] = loadSignalData(dataTable, opti % found_sync - Boolean indicating if synchronization was successful found_sync = 0; -tempLocalStorage = 1; +tempLocalStorage = 0; + +% Define the fixed storage directory relative to the user's MATLAB preferences directory +storage_dir = fullfile(prefdir, 'temp_sync_data'); % Part A: Check and load from local storage if available if tempLocalStorage == 1 - local_filename = sprintf('sync_data_run_%s.mat', num2str(dataTable.run_id)); + local_filename = fullfile(storage_dir, sprintf('sync_data_run_%s.mat', num2str(dataTable.run_id))); if exist(local_filename, 'file') % Load from local storage and return - load(local_filename, 'Bits', 'Symbols', 'Scpe_cell'); - found_sync = 1; + try + load(local_filename, 'Bits', 'Symbols', 'Scpe_cell'); + found_sync = 1; + return + catch + delete(local_filename); + end end end % If not locally saved, load from storage if ~found_sync % Load transmitted bits - Bits = load([options.storage_path, char(dataTable.tx_bits_path)]); + Bits = load(fullfile([options.storage_path, char(dataTable.tx_bits_path)])); Bits = Bits.Bits; % Map bits to symbols @@ -37,7 +45,7 @@ if ~found_sync Symbols_mapped.fs = fsym; % Load original symbols - Symbols = load([options.storage_path, char(dataTable.tx_symbols_path)]); + Symbols = load(fullfile([options.storage_path, char(dataTable.tx_symbols_path)])); Symbols = Symbols.Symbols; found_sync = 0; @@ -45,7 +53,7 @@ if ~found_sync % Try to load pre-synchronized data try - Scpe_load = load([options.storage_path, char(dataTable.rx_sync_path)]); + Scpe_load = load(fullfile([options.storage_path, char(dataTable.rx_sync_path)])); Scpe_cell = Scpe_load.S; [~,~,~,found_sync] = Scpe_cell{2}.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); catch @@ -74,21 +82,21 @@ end % Part B: Save to local storage if data was loaded and synced if tempLocalStorage == 1 && found_sync - local_filename = sprintf('sync_data_run_%s.mat', num2str(dataTable.run_id)); - % Create directory if it doesn't exist - if ~exist('temp_sync_data', 'dir') - mkdir('temp_sync_data'); + if ~exist(storage_dir, 'dir') + mkdir(storage_dir); end +% local_filename = fullfile(storage_dir, sprintf('sync_data_run_%s.mat', num2str(dataTable.run_id))); + % List existing files and remove oldest if more than N max_local_files = 10; % Store up to N files - files = dir('sync_data_run_*.mat'); + files = dir(fullfile(storage_dir, 'sync_data_run_*.mat')); if length(files) >= max_local_files % Sort by date [~, idx] = sort([files.datenum]); % Delete oldest file - delete(files(idx(1)).name); + delete(fullfile(storage_dir, files(idx(1)).name)); end % Save current data diff --git a/Functions/Job_Processing/submitJobs.m b/Functions/Job_Processing/submitJobs.m index a5a45a4..6cc834f 100644 --- a/Functions/Job_Processing/submitJobs.m +++ b/Functions/Job_Processing/submitJobs.m @@ -89,11 +89,13 @@ if submit_mode == processingMode.parallel % Submit job with proper arguments futures{jobCounter} = parfeval(@dsp_runid, 1, run_ids(r), ... - "database_name", dsp_options.database_name, ... + "database_type", dsp_options.database_type,... + "dataBase", dsp_options.dataBase, ... "append_to_db", dsp_options.append_to_db, ... - "database_path", dsp_options.database_path, ... + "load_file_path", dsp_options.load_file_path, ... "max_occurences", dsp_options.max_occurences, ... "storage_path", dsp_options.storage_path, ... + "mode", dsp_options.mode,... "parameters", optionalVars); fprintf('[RunID %d, Job %d] Submitted to pool.\n', run_ids(r), k); @@ -132,11 +134,13 @@ elseif submit_mode == processingMode.serial fprintf('[RunID %d, Job %d] Running in linear mode...\n', run_ids(r), k); results{k,r} = dsp_runid(run_ids(r),... - "database_name", dsp_options.database_name,... + "database_type", dsp_options.database_type,... + "dataBase", dsp_options.dataBase,... "append_to_db", dsp_options.append_to_db,... - "database_path", dsp_options.database_path,... + "load_file_path", dsp_options.load_file_path,... "max_occurences", dsp_options.max_occurences,... "storage_path", dsp_options.storage_path,... + "mode", dsp_options.mode,... "parameters", optionalVars); fprintf('[RunID %d, Job %d] Completed successfully.\n', run_ids(r), k); @@ -172,7 +176,7 @@ wh = submit_options.wh; function setupParallelPool() curpool = gcp('nocreate'); if isempty(curpool) - parpool; + parpool(10,"IdleTimeout",300); else if ~isempty(curpool.FevalQueue.QueuedFutures) || ~isempty(curpool.FevalQueue.RunningFutures) oldq = length(curpool.FevalQueue.QueuedFutures) + length(curpool.FevalQueue.RunningFutures); diff --git a/Functions/Metrics/calc_snr.m b/Functions/Metrics/calc_snr.m index 2dda8b0..d466022 100644 --- a/Functions/Metrics/calc_snr.m +++ b/Functions/Metrics/calc_snr.m @@ -23,8 +23,6 @@ function [snr_all, snr_per_level] = calc_snr(tx_signal, eq_noise) levels = unique(tx_signal); % Preallocate an array to store the SNR for each unique level - snr_per_level = zeros(size(levels)); - % Loop over each unique level to compute the SNR for that level for i = 1:length(levels) % Find indices where tx_signal equals the current level @@ -32,5 +30,7 @@ function [snr_all, snr_per_level] = calc_snr(tx_signal, eq_noise) % Compute the SNR for these indices snr_per_level(i) = snr(tx_signal(idx), eq_noise(idx)); + + histogram(eq_noise(idx)); end end diff --git a/Functions/Theory/analyze_moving_average_filter.m b/Functions/Theory/analyze_moving_average_filter.m new file mode 100644 index 0000000..43bcc5f --- /dev/null +++ b/Functions/Theory/analyze_moving_average_filter.m @@ -0,0 +1,64 @@ +% Parameters +N_values = [10 100 1000 4096]; % Different filter lengths to analyze +fs = 112e9; + +% Create figure +figure; + +% Plot frequency responses +subplot(211) +hold on +grid on +ylabel('Magnitude (dB)') +title('Frequency Response') +yline(-3,'--r') +ylim([-40 5]) + +subplot(212) +hold on +grid on +xlabel('Frequency (GHz)') +ylabel('Phase (rad)') +title('Phase Response') + +% Color map for different lines +colors = cbrewer2('Set1',length(N_values)); + +% Loop through different filter lengths +for i = 1:length(N_values) + N = N_values(i); + + % Filter coefficients + b = ones(1,N)/N; + a = 1; + + % Frequency response + [h,w] = freqz(b,a,4096*8); + freq = (w/(2*pi))*fs; + h_db = 20*log10(abs(h)); + + % Plot magnitude response + subplot(211) + plot(freq/1e9, h_db, 'Color', colors(i,:), 'DisplayName', sprintf('N=%d', N),'LineWidth',0.1) + + % Plot phase response + subplot(212) + plot(freq/1e9, unwrap(angle(h)), 'Color', colors(i,:), 'DisplayName', sprintf('N=%d', N),'LineWidth',0.1) + + % Find -3dB frequency + cutoff_idx = find(h_db <= -3, 1); + f_cutoff = freq(cutoff_idx)/1e9; + fprintf('N=%d: Cutoff frequency (-3dB point): %.2f GHz\n', N, f_cutoff) +end + +% Add legend and adjust axes +subplot(211) +legend('show') +xlim([0 16]) % Adjust x-axis limit to better see the differences + +subplot(212) +legend('show') +xlim([0 16]) % Adjust x-axis limit to better see the differences + +% Analytical approximation +f_3db_approx = 0.443 * fs./N_values ./ 1e9; \ No newline at end of file diff --git a/projects/ECOC_2025/auswertung_algorithms/mpi_dsp_debug.m b/projects/ECOC_2025/auswertung_algorithms/mpi_dsp_debug.m new file mode 100644 index 0000000..d765e30 --- /dev/null +++ b/projects/ECOC_2025/auswertung_algorithms/mpi_dsp_debug.m @@ -0,0 +1,60 @@ + +load("ffe_debug_snapshot.mat"); + + +dc_buffer_len = logspace(0,3,12); +dc_buffer_len = 1024; + +mu_dc = logspace(-3,0,24); + +parfor d = 1:length(mu_dc) + + eq_lin = FFE_DCremoval_adaptive_mu("epochs_tr",5,"epochs_dd",3,"len_tr",4096*2,"mu_dd",... + 0.0002,"mu_tr",0,"order",25,"sps",2,"decide",0,... + "mu_dc",mu_dc(d),... + "dc_buffer_len",1024, ... + "ffe_buffer_len",1,... + "smoothing_buffer_length",0,... + "smoothing_buffer_update",1,... + "adaptive_mu_mode",0); + % + % eq_lin = FFE("epochs_tr",5,"epochs_dd",3,"len_tr",4096*2,"mu_dd",... + % 0.0002,"mu_tr",0,"order",25,"sps",2,"decide",0); + + ffe_results = ffe(eq_lin,M,Scpe_sig,Symbols,Tx_bits,... + "precode_mode",duob_mode,... + 'showAnalysis',0,... + "postFFE",[],... + "eth_style_symbol_mapping",0); + + ffe_results.metrics.print + + % % eq_lin = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",21,"dfe_order",2,"sps",2,"decide",0); + % + % eq_lin = FFE("epochs_tr",5,"epochs_dd",3,"len_tr",4096*2,"mu_dd",... + % 0.0002,"mu_tr",0,"order",25,"sps",2,"decide",0); + % pf_ = Postfilter("ncoeff",2,"useBurg",1); + % mlse_ = MLSE_viterbi("duobinary_output",0,'M',4,'trellis_states',PAMmapper(4,0).levels); + % + % [ffe_results2, mlse_results] = vnle_postfilter_mlse(eq_lin, pf_, mlse_, M, Scpe_sig, Symbols, Tx_bits, ... + % "precode_mode", duob_mode,... + % 'showAnalysis', 1, ... + % "postFFE", [],... + % "eth_style_symbol_mapping", 0); + + ber(d) = ffe_results.metrics.BER; + + +end + +figure(10) +hold on +plot(mu_dc,ber,'LineWidth',1,'DisplayName',sprintf('DC buffer len = 1024'),'Marker','.','MarkerSize',10); +xlabel('BER'); +xlabel('MU DC'); +title('BER Optimization over dc\_buffer\_len'); +yline([4.85e-3,2e-2],'HandleVisibility', 'off','LineWidth',1,'LineStyle','--'); +ylim([9e-4, 0.5]); +set(gca, 'YScale', 'log'); % BER is usually plotted log-scale +legend('show', 'Location', 'best'); +grid on; \ No newline at end of file diff --git a/projects/ECOC_2025/auswertung_algorithms/run_offline_dsp.m b/projects/ECOC_2025/auswertung_algorithms/run_offline_dsp.m new file mode 100644 index 0000000..af3e9d7 --- /dev/null +++ b/projects/ECOC_2025/auswertung_algorithms/run_offline_dsp.m @@ -0,0 +1,118 @@ +% === SETTINGS === + +dsp_options.append_to_db = 0; +dsp_options.max_occurences = 15; +dsp_options.database_path = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +dsp_options.database_name = 'silas_labor_newdsp_newstructure.db'; +dsp_options.storage_path = 'Z:\2024\sioe_labor\'; + +dsp_options.parameters = struct(); +dsp_options.parameters.mu_dc = [0.005]; + +% === Get Run ID's === +db = DBHandler("pathToDB", [dsp_options.database_path, dsp_options.database_name], "type", "sqlite"); +fp = QueryFilter(); +% fp.where('Runs', 'run_id','EQUALS', 5108); +fp.where('Runs', 'is_mpi','EQUALS', 0); +fp.where('Runs', 'fiber_length','EQUALS', 1); +fp.where('Runs', 'wavelength','EQUALS', 1310); +fp.where('Runs', 'db_mode','EQUALS', 1); +fp.where('Runs', 'rop_attenuation','EQUALS', 0); +fp.where('Runs', 'pam_level','EQUALS', 4); +fp.where('Runs', 'bitrate','EQUALS', 360e9); +% fp.where('Runs', 'power_pd_in','GREATER_THAN', 7); +[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); + +% === Initialize DataStorage === +wh = DataStorage(dsp_options.parameters); +wh.addStorage("ffe_package"); +wh.addStorage("mlse_package"); +wh.addStorage("vnle_package"); +wh.addStorage("dbtgt_package"); +wh.addStorage("dbenc_package"); + +% === RUN IT === + +% [results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "serial", 'wh', wh, 'waitbar', true); +% wh.getStoValue('ffe_package',0.005); +% wh.getStoValue('mlse_package',0.005); + +[dataTable,~] = db.queryDB(fp, [db.getTableFieldNames('Runs');db.getTableFieldNames('Results');db.getTableFieldNames('Equalizer')]); + +dataTable = cleanUpTable(dataTable); + +% === Look at it === +y_var = 'BER_precoded'; +x_var = 'bitrate'; +fixedVars = {'equalizer_structure', x_var}; + +[dataTableClean, outliersTable] = removeGroupOutliers(dataTable, fixedVars, y_var); + +% --- Group and aggregate --- +dataTableGrpd_mean = groupIt(fixedVars, dataTableClean, @mean); +dataTableGrpd_min = groupIt(fixedVars, dataTableClean, @min); +dataTableGrpd_max = groupIt(fixedVars, dataTableClean, @max); + +% Choose a color map +cols = linspecer(numel(unique(dataTableGrpd_mean.equalizer_structure))); + +figure; +hold on; + +% Get unique equalizer structures for grouping +unique_eq = unique(dataTableGrpd_mean.equalizer_structure); + +for i = 1:numel(unique_eq) + eq_val = unique_eq(i); + + % Filter grouped data for this equalizer structure + filt = dataTableGrpd_mean.equalizer_structure == eq_val; + + x = dataTableGrpd_mean.(x_var)(filt); + y_mean = dataTableGrpd_mean.(y_var)(filt); + y_min = dataTableGrpd_min.(y_var)(filt); + y_max = dataTableGrpd_max.(y_var)(filt); + + % Bounds for boundedline (distance from mean) + y_lower = y_mean - y_min; + y_upper = y_max - y_mean; + y_bounds = [y_lower, y_upper]; + + % --- Bounded line (mean ± min/max) --- + if exist('boundedline', 'file') + [hl, hp] = boundedline(x, y_mean, y_bounds, ... + 'alpha', 'transparency', 0.1, ... + 'cmap', cols(i,:), ... + 'nan', 'fill', ... + 'orientation', 'vert'); + set(hl, 'LineWidth', 1.2, 'DisplayName', sprintf('Eq %s', eq_val)); + set(hp, 'HandleVisibility', 'off'); + else + % If boundedline is not available, use errorbar + errorbar(x, y_mean, y_lower, y_upper, ... + 'o-', 'Color', cols(i,:), 'LineWidth', 1.2, ... + 'DisplayName', sprintf('Eq %d', eq_val),'HandleVisibility', 'off'); + end + + % --- Normal line (mean only) --- + plot(x, y_mean, '-', 'Color', cols(i,:), 'LineWidth', 1.5, ... + 'DisplayName', sprintf('Mean Eq %s', eq_val),'HandleVisibility', 'off'); + + % --- Scatter plot for individual points (from original data) --- + % Filter original data for this group + orig_filt = dataTableClean.equalizer_structure == eq_val; + x_scatter = dataTableClean.(x_var)(orig_filt); + y_scatter = dataTableClean.(y_var)(orig_filt); + + scatter(x_scatter, y_scatter, 10,cols(i,:), 'filled', ... + 'MarkerFaceAlpha', 0.5, 'DisplayName', sprintf('Scatter Eq %s', eq_val),'HandleVisibility', 'off'); +end + +yline([2.2e-4,4.85e-3,2e-2],'HandleVisibility', 'off','LineWidth',1,'LineStyle','--'); +set(gca, 'YScale', 'log'); % BER is usually plotted log-scale +xlabel(x_var, 'Interpreter', 'none'); +ylabel(y_var, 'Interpreter', 'none'); +legend('show', 'Location', 'best'); +grid on; +title(sprintf('%s vs. %s', y_var, x_var), 'Interpreter', 'none'); +hold off; \ No newline at end of file diff --git a/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m new file mode 100644 index 0000000..c73a50b --- /dev/null +++ b/projects/HighSpeedExperiment_2024/Auswertung_JLT/run_dsp_from_db.m @@ -0,0 +1,224 @@ +% === SETTINGS === +dsp_options.append_to_db = 1; +dsp_options.max_occurences = 1; + +experiment = "highspeed_2024"; +dsp_options.mode = "load_run_id"; % 'simulate' & 'load_files' +dsp_options.load_file_path = struct(); + +if dsp_options.mode == "load_run_id" + + if experiment == "highspeed_2024" + + 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); + + elseif experiment == "mpi_ecoc_2025" + + dsp_options.database_type = 'mysql'; + dsp_options.dataBase = 'labor'; + dsp_options.storage_path = 'Z:\2025\ECOC Silas\ecoc_2025\'; + db = DBHandler("dataBase", [dsp_options.dataBase], "type", dsp_options.database_type); + + end + +elseif dsp_options.mode == "load_files" + + dsp_options.load_file_path.tx_bits_path = "Z:\2025\ECOC Silas\ecoc_2025\mpi_opti_1000m_pam2 4 6 8\20250417_091513_PAM_4_R_112_bits.mat"'; + dsp_options.load_file_path.tx_symbols_path = "Z:\2025\ECOC Silas\ecoc_2025\mpi_opti_1000m_pam2 4 6 8\20250417_091513_PAM_4_R_112_symbols.mat"'; + dsp_options.load_file_path.rx_raw_path = "Z:\2025\ECOC Silas\ecoc_2025\mpi_opti_1000m_pam2 4 6 8\20250417_091525_PAM_4_R_112_rec01_rx_signal_raw.mat"'; + +elseif dsp_options.mode == "simulate" + + error('Not yet implemented') + +end + +% === Get Run ID's === + +fp = QueryFilter(); +% fp.where('Runs', 'run_id','EQUALS', 987); +% fp.where('Runs', 'pam_level','EQUALS', 4); +fp.where('Runs', 'bitrate','LESS_THAN', 310e9); +fp.where('Runs', 'fiber_length','EQUALS', 1); +fp.where('Runs', 'is_mpi','EQUALS', 0); +% fp.where('Runs', 'interference_path_length','EQUALS', 1000); +% fp.where('Runs', 'loop_id','GREATER_THAN', 11); +% fp.where('Runs', 'sir','EQUALS',18); +fp.where('Runs', 'wavelength','EQUALS', 1310); +% fp.where('Runs', 'db_mode','EQUALS', 1); +fp.where('Runs', 'rop_attenuation','EQUALS', 0); +% fp.where('Runs', 'power_pd_in','GREATER_THAN', 7); + + +[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); + + +% Keep only the rows corresponding to the first occurrence of each 'sir' value +% [~, unique_indices] = unique(dataTable.sir, 'first'); +% dataTable = dataTable(unique_indices, :); + +% dataTable = dataTable(1,:); + +% === Set LOOPS & Initialize DataStorage === +dsp_options.parameters = struct(); +% dsp_options.parameters.pf_ncoeffs = [1,2];%[0,logspace(-4,0,10)]; + +wh = DataStorage(dsp_options.parameters); +wh.addStorage("ffe_package"); +wh.addStorage("mlse_package"); +wh.addStorage("vnle_package"); +wh.addStorage("dbtgt_package"); +wh.addStorage("dbenc_package"); + +% === RUN IT === + +[results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "serial", 'wh', wh, 'waitbar', true); + + +% wh.getStoValue('ffe_package',0.005); +% wh.getStoValue('mlse_package',0.005); + +% [dataTable,~] = db.queryDB(fp, [db.getTableFieldNames('Runs');db.getTableFieldNames('Results');db.getTableFieldNames('Equalizer')]); +% +% dataTable = cleanUpTable(dataTable); +% wh_analyze = wh_adap; +% % wh_analyze = wh_dcremoval_old2; +% % wh_analyze = wh_dcremoval_old2; +% +% res = cell(wh_analyze.parameter.mu_dc.length,wh_analyze.parameter.dc_buffer_len.length); +% ber_mean = zeros(wh_analyze.parameter.mu_dc.length,wh_analyze.parameter.dc_buffer_len.length); +% for m = 1:wh_analyze.parameter.mu_dc.length +% for b = 1:wh_analyze.parameter.dc_buffer_len.length +% res{m,b}=wh_analyze.getStoValue("ffe_package",wh_analyze.parameter.mu_dc.values(m),wh_analyze.parameter.dc_buffer_len.values(b)); +% try +% cells = res{m,b}{1}; +% idx = cellfun(@(c) ~isempty(c), cells); +% ber = cellfun(@(c) c.metrics.BER, cells(idx)); +% ber_mean(m,b) = mean(ber); +% catch +% ber_mean(m,b) = NaN; +% end +% end +% end +% +% figure() +% hold on +% ber_fix = cellfun(@(c) c.ffe_package{1}.metrics.BER, results_fix); +% ber_adap_m2 = cellfun(@(c) c.ffe_package{1}.metrics.BER, results_adap); +% ber_adap_method1 = cellfun(@(c) c.ffe_package{1}.metrics.BER, results); +% plot(dataTable.sir,ber_fix,'LineWidth',1,'DisplayName',sprintf('DCt; fix mu = 0.5; p=1024'),'Marker','.','MarkerSize',10); +% plot(dataTable.sir,ber_adap_method1,'LineWidth',1,'DisplayName',sprintf('DCt; adap mu 1; p=1024'),'Marker','.','MarkerSize',10); +% plot(dataTable.sir,ber_adap_m2,'LineWidth',1,'DisplayName',sprintf('DCt; adap mu 2; p=1024'),'Marker','.','MarkerSize',10); +% xlabel('BER'); +% xlabel('SIR'); +% yline([4.85e-3,2e-2],'HandleVisibility', 'off','LineWidth',1,'LineStyle','--'); +% % ylim([9e-4, 0.5]); +% set(gca, 'YScale', 'log'); % BER is usually plotted log-scale +% legend('show', 'Location', 'best'); +% grid on; +% +% +% +% +% figure; clf +% +% % Create meshgrid for contourf +% [X, Y] = meshgrid(dsp_options.parameters.mu_dc, dsp_options.parameters.dc_buffer_len); +% +% % Create contour plot +% contourf(X, Y, ber_mean', 20); % 20 contour levels, adjust as needed +% +% % Set axes to logarithmic scale +% set(gca, 'XScale', 'log', 'YScale', 'log'); +% +% colormap("parula"); +% c = colorbar; +% c.Label.String = 'BER'; +% +% xlabel('\mu_{dc}'); +% ylabel('dc\_buffer\_len'); +% title('BER Optimization over \mu_{dc} and dc\_buffer\_len'); +% +% % Make plot prettier +% grid on +% set(gca, 'Layer', 'top'); % Put grid lines on top of contours +% +% +% % === Look at it === +% y_var = 'BER_precoded'; +% x_var = 'bitrate'; +% fixedVars = {'equalizer_structure', x_var}; +% +% [dataTableClean, outliersTable] = removeGroupOutliers(dataTable, fixedVars, y_var); +% +% % --- Group and aggregate --- +% dataTableGrpd_mean = groupIt(fixedVars, dataTableClean, @mean); +% dataTableGrpd_min = groupIt(fixedVars, dataTableClean, @min); +% dataTableGrpd_max = groupIt(fixedVars, dataTableClean, @max); +% +% % Choose a color map +% cols = linspecer(numel(unique(dataTableGrpd_mean.equalizer_structure))); +% +% figure; +% hold on; +% +% % Get unique equalizer structures for grouping +% unique_eq = unique(dataTableGrpd_mean.equalizer_structure); +% +% for i = 1:numel(unique_eq) +% eq_val = unique_eq(i); +% +% % Filter grouped data for this equalizer structure +% filt = dataTableGrpd_mean.equalizer_structure == eq_val; +% +% x = dataTableGrpd_mean.(x_var)(filt); +% y_mean = dataTableGrpd_mean.(y_var)(filt); +% y_min = dataTableGrpd_min.(y_var)(filt); +% y_max = dataTableGrpd_max.(y_var)(filt); +% +% % Bounds for boundedline (distance from mean) +% y_lower = y_mean - y_min; +% y_upper = y_max - y_mean; +% y_bounds = [y_lower, y_upper]; +% +% % --- Bounded line (mean ± min/max) --- +% if exist('boundedline', 'file') +% [hl, hp] = boundedline(x, y_mean, y_bounds, ... +% 'alpha', 'transparency', 0.1, ... +% 'cmap', cols(i,:), ... +% 'nan', 'fill', ... +% 'orientation', 'vert'); +% set(hl, 'LineWidth', 1.2, 'DisplayName', sprintf('Eq %s', eq_val)); +% set(hp, 'HandleVisibility', 'off'); +% else +% % If boundedline is not available, use errorbar +% errorbar(x, y_mean, y_lower, y_upper, ... +% 'o-', 'Color', cols(i,:), 'LineWidth', 1.2, ... +% 'DisplayName', sprintf('Eq %d', eq_val),'HandleVisibility', 'off'); +% end +% +% % --- Normal line (mean only) --- +% plot(x, y_mean, '-', 'Color', cols(i,:), 'LineWidth', 1.5, ... +% 'DisplayName', sprintf('Mean Eq %s', eq_val),'HandleVisibility', 'off'); +% +% % --- Scatter plot for individual points (from original data) --- +% % Filter original data for this group +% orig_filt = dataTableClean.equalizer_structure == eq_val; +% x_scatter = dataTableClean.(x_var)(orig_filt); +% y_scatter = dataTableClean.(y_var)(orig_filt); +% +% scatter(x_scatter, y_scatter, 10,cols(i,:), 'filled', ... +% 'MarkerFaceAlpha', 0.5, 'DisplayName', sprintf('Scatter Eq %s', eq_val),'HandleVisibility', 'off'); +% end +% +% yline([2.2e-4,4.85e-3,2e-2],'HandleVisibility', 'off','LineWidth',1,'LineStyle','--'); +% set(gca, 'YScale', 'log'); % BER is usually plotted log-scale +% xlabel(x_var, 'Interpreter', 'none'); +% ylabel(y_var, 'Interpreter', 'none'); +% legend('show', 'Location', 'best'); +% grid on; +% title(sprintf('%s vs. %s', y_var, x_var), 'Interpreter', 'none'); +% hold off; \ No newline at end of file diff --git a/projects/IMDD_base_system/imdd_model.m b/projects/IMDD_base_system/imdd_model.m index c15a919..70b7b29 100644 --- a/projects/IMDD_base_system/imdd_model.m +++ b/projects/IMDD_base_system/imdd_model.m @@ -134,6 +134,7 @@ if fsym_ ~= fsym fsym = fsym_; % fprintf('Adapted symbolrate to %d GBd, to match provided bitrate of %d GBit/s using PAM %d \n',fsym.*1e-9,bitrate.*1e-9, M); end + f_nyquist = fsym/2; %%% run the simulation or measurement or ... diff --git a/projects/IMDD_base_system/simulation_bwl.m b/projects/IMDD_base_system/simulation_bwl.m new file mode 100644 index 0000000..e1477a7 --- /dev/null +++ b/projects/IMDD_base_system/simulation_bwl.m @@ -0,0 +1,227 @@ +%%% Run parameters +% TX +M = 4; +fsym = 180e9; + +apply_pulsef = 1; +fdac = 256e9; +fadc = 256e9; +random_key = 1; + +db_precode = 0; +db_encode = 0; + +rcalpha = 0.05; +kover = 16; +vbias_rel = 0.5; +u_pi = 2.9; +vbias = -vbias_rel*u_pi; +laser_wavelength = 1310; +laser_linewidth = 0; +tx_bw_nyquist = 0.6; +rx_bw_nyquist = 0.6; + + +% Channel +link_length = 1; + +% RX +rop = -8; + + +vnle_order1 = 50; +vnle_order2 = 3; +vnle_order3 = 3; + +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; +dfe_order = [0 0 0]; + +alpha = 0; + +len_tr = 4096*2; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +% mu_dc = 0; + +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; + + +dfe_ = sum(dfe_order)>0; + +doub_mode = db_mode.no_db; + + +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + +db_precode = 0; +db_encode = 0; +duob_mode = db_mode.no_db; +apply_pulsef = 1; +[Digi_sig,Symbols,Tx_bits] = PAMsource(... + "fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key,... + "db_precode",db_precode,"db_encode",db_encode,... + "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process(); + +Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); + + +%%%%% AWG +% El_sig = M8199A("kover",kover).process(Digi_sig); +El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); +% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); +% El_sig = El_sig.setPower(0,"dBm"); + +%%%%% Low-pass el. components %%%%%% +f_nyquist = fsym/2; +tx_bwl = tx_bw_nyquist.*f_nyquist; +% tx_bwl = 80e9; +El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); +% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1); + +%%%%% Electrical Driver Amplifier %%%%%% +El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); +El_sig = El_sig.normalize("mode","oneone"); + +%%%%% MODULATE E/O CONVERSION %%%%%% +[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); + +Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); + + + +rop = [-10]; +ber_vnle = []; +ber_mlse = []; +ber_viterbi = []; +ber_db = []; +ber_db_diff_precoded = []; +gmi_vnle = []; +gmi_mlse = []; +gmi_mlse_db = []; + +for r = 1:length(rop) + + + %%%%%% ROP %%%%%% + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop(r)).process(Opt_sig); + + %%%%%% PD Square Law %%%%%% + Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); + + %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% + rx_bwl = rx_bw_nyquist.*f_nyquist; + % rx_bwl = 80e9; + Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); + + % %%%%%% Low-pass Scope %%%%%% + Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + + % Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); + + %%%%%% Scope %%%%%% + Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); + + Scpe_sig_resampled = Scpe_sig.resample("fs_out",2*fsym); + + [~, Scpe_cell, ~, found_sync] = Scpe_sig_resampled.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + + + eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"Nb",[dfe_order],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + mlse_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + + + %Duobinary Targeting + if 0 + db_ref_sequence = Duobinary().encode(Symbols); + db_ref_constellation = unique(db_ref_sequence.signal); + [eq_signal, eq_noise] = eq_.process(Scpe_cell{1},db_ref_sequence); + + mlse_.DIR = [1,1]; + [mlse_sig_sd,LLR,gmi_mlse_db(r)] = mlse_.process(eq_signal,Symbols); + + mlse_sig_hd = PAMmapper(M,0,"eth_style",0).quantize(mlse_sig_sd); + mlse_sig_hd_precoded = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded,"M",M); + + tx_symbols_precoded = Duobinary().encode(Symbols); + tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); + + tx_bits_precoded = PAMmapper(M,0,"eth_style",0).demap(tx_symbols_precoded); + + rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd_precoded); + [~,errors_db_diff_precoded,ber_db_diff_precoded(r),~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + %B) Just determine BER + rx_bits_mlse = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); + [bits_mlse,errors_mlse,ber_db(r),~] = calc_ber(rx_bits_mlse.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + end + + % FFE or VNLE + eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"Nb",[dfe_order],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + [eq_signal_sd, eq_noise] = eq_.process(Scpe_cell{1}, Symbols); + [gmi_vnle(r)] = calc_air(eq_signal_sd, Symbols, "skip_front", 100, "skip_end", 100); + eq_signal_sd.plot("displayname",'bla','fignum',118); + + % Hard decision on VNLE output + eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); + rx_bits = PAMmapper(M,0,"eth_style",0).demap(eq_signal_hd); + [~,~,ber_vnle(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('BER VNLE: %.2e \n',ber_vnle(r)); + fprintf('NGMI VNLE: %.2f \n',gmi_vnle(r)./log2(M)); + + + % Process through postfilter and MLSE + pf_ncoeffs = 1; + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); + mlse_.DIR = pf_.coefficients; + [signalclass_hd,LLR,gmi_mlse(r)] = mlse_.process(mlse_sig_sd,Symbols); + mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(signalclass_hd); + rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); + [~,~,ber_mlse(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('BER: %.2e \n',ber_mlse(r)); + fprintf('NGMI MLSE: %.2f \n',gmi_mlse(r)./log2(M)); + + + % Process through postfilter and MLSE + pf_ncoeffs = 1; + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); + mlse_.DIR = pf_.coefficients; + mlse_sig_sd = mlse_.process(mlse_sig_sd); + mlse_sig_hd = PAMmapper(M, 0, "eth_style", 0).quantize(mlse_sig_sd); + rx_bits = PAMmapper(M,0,"eth_style",0).demap(mlse_sig_hd); + [~,~,ber_viterbi(r),~] = calc_ber(rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('FW BER: %.2e \n',ber_viterbi(r)); + +end + + +figure();hold on +plot(rop,gmi_mlse,'DisplayName','MLSE','Marker','*'); +plot(rop,gmi_vnle,'DisplayName','VNLE','Marker','*'); +plot(rop,gmi_mlse_db,'DisplayName','MLSE DB Tgt','Marker','*'); +ylim([0 2.6]); + +figure();hold on +plot(rop,ber_vnle,'DisplayName','VNLE'); + +plot(rop,ber_mlse,'DisplayName','MLSE','Marker','*'); +plot(rop,ber_viterbi,'DisplayName','Viterbi','Marker','*'); +plot(rop,ber_db_diff_precoded,'DisplayName','MLSE db diff','Marker','*'); +plot(rop,ber_db,'DisplayName','MLSE db','Marker','*'); +set(gca, 'yscale', 'log'); +legend diff --git a/projects/IMDD_base_system/simulation_bwl_2.m b/projects/IMDD_base_system/simulation_bwl_2.m new file mode 100644 index 0000000..8dff4a2 --- /dev/null +++ b/projects/IMDD_base_system/simulation_bwl_2.m @@ -0,0 +1,144 @@ +%%% Run parameters +% TX +M = 4; +fsym = 180e9; + +apply_pulsef = 1; +fdac = 256e9; +fadc = 256e9; +random_key = 1; + +db_precode = 0; +db_encode = 0; + +rcalpha = 0.05; +kover = 16; +vbias_rel = 0.5; +u_pi = 2.9; +vbias = -vbias_rel*u_pi; +laser_wavelength = 1310; +laser_linewidth = 0; +tx_bw_nyquist = 0.7; + +% Channel +link_length = 1; + +% RX +rop = -8; +rx_bw_nyquist = 0.7; + +vnle_order1 = 50; +vnle_order2 = 3; +vnle_order3 = 3; + +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; +dfe_order = [0 0 0]; + +pf_ncoeffs = 1; + +alpha = 0; + +len_tr = 4096*2; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +% mu_dc = 0; + +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; + + +dfe_ = sum(dfe_order)>0; + +doub_mode = db_mode.no_db; + + +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + +db_precode = 0; +db_encode = 0; +duob_mode = db_mode.no_db; +apply_pulsef = 1; +[Digi_sig,Symbols,Tx_bits] = PAMsource(... + "fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key,... + "db_precode",db_precode,"db_encode",db_encode,... + "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode).process(); + +Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); + + + +%%%%% AWG +% El_sig = M8199A("kover",kover).process(Digi_sig); +El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); +% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); +% El_sig = El_sig.setPower(0,"dBm"); + +%%%%% Low-pass el. components %%%%%% +f_nyquist = fsym/2; +tx_bwl = tx_bw_nyquist.*f_nyquist; +% tx_bwl = 80e9; +El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); +% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1); + +%%%%% Electrical Driver Amplifier %%%%%% +El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); +El_sig = El_sig.normalize("mode","oneone"); + +%%%%% MODULATE E/O CONVERSION %%%%%% +[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); + +Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); + +%%%%%% ROP %%%%%% +Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); + +%%%%%% PD Square Law %%%%%% +Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); + +%%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% +rx_bwl = rx_bw_nyquist.*f_nyquist; +% rx_bwl = 80e9; +Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); + +% %%%%%% Low-pass Scope %%%%%% +Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + +% Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); + +%%%%%% Scope %%%%%% +Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); + +Scpe_sig_resampled = Scpe_sig.resample("fs_out",2*fsym); + +[~, Scpe_cell, ~, found_sync] = Scpe_sig_resampled.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1); + +eq_ = EQ("Ne",[vnle_order1,vnle_order2,vnle_order3],"Nb",[dfe_order],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); +pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); +% mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); +mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + +% FFE or VNLE +[eq_signal_sd, eq_noise] = eq_.process(Scpe_cell{1}, Symbols); + +% Hard decision on VNLE output +eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); + +% Process through postfilter and MLSE +[mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); +mlse_.DIR = pf_.coefficients; + +constellation = [-3, -1, 1, 3]; +chatgpt_answer(mlse_sig_sd.signal, Symbols.signal,mlse_.DIR,constellation); + +% mlse_sig_sd = mlse_.process(mlse_sig_sd,Symbols); +% mlse_sig_hd = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).quantize(mlse_sig_sd); diff --git a/projects/IMDD_base_system/tx_simulation.m b/projects/IMDD_base_system/tx_simulation.m index 77a1679..1cd6216 100644 --- a/projects/IMDD_base_system/tx_simulation.m +++ b/projects/IMDD_base_system/tx_simulation.m @@ -1,6 +1,6 @@ - Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rcalpha); + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rcalpha); [Digi_sig,Symbols,Tx_bits] = PAMsource(... "fsym",fsym,"M",M,"order",19,"useprbs",1,... diff --git a/projects/Job_Processing/run_offline_dsp.m b/projects/Job_Processing/run_offline_dsp.m index 2544680..335b57b 100644 --- a/projects/Job_Processing/run_offline_dsp.m +++ b/projects/Job_Processing/run_offline_dsp.m @@ -32,7 +32,7 @@ wh.addStorage("dbenc_package"); % === RUN IT === -% [results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "parallel", 'wh', wh, 'waitbar', true); +[results,wh] = submitJobs(dataTable.run_id(:), dsp_options, "parallel", 'wh', wh, 'waitbar', true); % wh.getStoValue('ffe_package',0.005); % wh.getStoValue('mlse_package',0.005); diff --git a/projects/Job_Processing/run_offline_dsp_script.m b/projects/Job_Processing/run_offline_dsp_script.m new file mode 100644 index 0000000..a75b3bc --- /dev/null +++ b/projects/Job_Processing/run_offline_dsp_script.m @@ -0,0 +1,116 @@ +dsp_options.append_to_db = false; +dsp_options.max_occurences = 1; +dsp_options.mode = "load_run_id"; +experiment = "highspeed_2024"; + +if dsp_options.mode == "load_run_id" + + dsp_options.load_file_path = []; + if experiment == 'highspeed_2024' + + dsp_options.database_path = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\silas_labor_newdsp_newstructure.db'; + dsp_options.storage_path = 'Z:\2024\sioe_labor\'; + db = DBHandler("pathToDB", [dsp_options.database_path], "type", "sqlite"); + + elseif 'mpi_ecoc_2025' + + dsp_options.database_path = []; %sqlite im netzwerk + dsp_options.storage_path = 'Z:\2025\ECOC Silas\ecoc_2025\'; + db = DBHandler("pathToDB", [dsp_options.database_path], "type", "mysql"); + + end + + +elseif dsp_options.mode == "load_files" + + dsp_options.load_file_path.tx_bits_path = "Z:\2025\ECOC Silas\ecoc_2025\mpi_opti_1000m_pam2 4 6 8\20250417_091513_PAM_4_R_112_bits.mat"'; + dsp_options.load_file_path.tx_symbols_path = "Z:\2025\ECOC Silas\ecoc_2025\mpi_opti_1000m_pam2 4 6 8\20250417_091513_PAM_4_R_112_symbols.mat"'; + dsp_options.load_file_path.rx_raw_path = "Z:\2025\ECOC Silas\ecoc_2025\mpi_opti_1000m_pam2 4 6 8\20250417_091525_PAM_4_R_112_rec01_rx_signal_raw.mat"'; + +elseif dsp_options.mode == "simulate" + + error('Not yet implemented') + +end + + +%% === Define Database Filters === + +if experiment == 'highspeed_2024' + + fp = QueryFilter(); + db = db.getDistinctValues; + %fp.where('Runs', 'run_id','EQUALS', 3); + fp.where('Runs', 'pam_level','EQUALS', 4); + allrates = string(num2str(sort(db.distinctValues.Runs.bitrate.*1e-9))); + fp.where('Runs', 'bitrate','EQUALS',390e9); + alllengths = string(num2str(sort(db.distinctValues.Runs.fiber_length))); + fp.where('Runs', 'fiber_length','EQUALS', 1); + fp.where('Runs', 'is_mpi','EQUALS', 0); + % fp.where('Runs', 'interference_path_length','EQUALS', 1000); + % fp.where('Runs', 'loop_id','GREATER_THAN', 11); + % fp.where('Runs', 'sir','EQUALS',18); + fp.where('Runs', 'wavelength','EQUALS', 1310); + fp.where('Runs', 'db_mode','EQUALS', 0); + % fp.where('Runs', 'rop_attenuation','EQUALS', 0); + fp.where('Runs', 'power_pd_in','GREATER_THAN', 7); + +elseif experiment == 'mpi_ecoc_2025' + + fp = QueryFilter(); + %fp.where('Runs', 'run_id','EQUALS', 3); + fp.where('Runs', 'pam_level','EQUALS', 4); + fp.where('Runs', 'bitrate','EQUALS',224e9); + fp.where('Runs', 'is_mpi','EQUALS', 1); + fp.where('Runs', 'interference_path_length','EQUALS', 1000); + % fp.where('Runs', 'loop_id','GREATER_THAN', 11); + fp.where('Runs', 'sir','EQUALS',18); + fp.where('Runs', 'wavelength','EQUALS', 1310); + fp.where('Runs', 'db_mode','EQUALS', 0); + % fp.where('Runs', 'rop_attenuation','EQUALS', 0); + fp.where('Runs', 'power_pd_in','GREATER_THAN', 7); + % fp.where('Runs', 'loop_id','EQUALS', 14); + +end + + +% === Queue from DB === +[dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); + +% === Adapt or filter Matlab Table === +% [~, unique_indices] = unique(dataTable.sir, 'first'); +% dataTable = dataTable(unique_indices, :); + + +% === Set LOOPS & Initialize DataStorage === +dsp_options.parameters = struct(); +% dsp_options.parameters.adaption = [3]; +% dsp_options.parameters.mu_tr = 0.999;%[0,logspace(-4,0,10)]; +% dsp_options.parameters.mu_dd = [logspace(-0.004,0,20)]; +% dsp_options.parameters.use_dd_mode = [1]; + +wh = DataStorage(dsp_options.parameters); +wh.addStorage("ffe_package"); +wh.addStorage("mlse_package"); +wh.addStorage("vnle_package"); +wh.addStorage("dbtgt_package"); +wh.addStorage("dbenc_package"); + + +% === RUN IT === +% run this function: dsp_runid(run_id, options) +[results,wh] = submitJobs(dataTable.run_id, dsp_options, "serial", 'wh', wh, 'waitbar', false); + + +%% +BERs = cellfun(@(c) c.ffe_package{1,1}.metrics.BER, results); +figure +hold on +plot(dsp_options.parameters.mu_dd,BERs,'Marker','+') +yline([2.2e-4,4.85e-3,2e-2],'HandleVisibility', 'off','LineWidth',1,'LineStyle','--'); +set(gca, 'YScale', 'log'); % BER is usually plotted log-scale +xlabel('$\mu$ DD', 'Interpreter', 'latex'); +ylabel('BER', 'Interpreter', 'latex'); +legend('show', 'Location', 'best'); +grid on; +hold off; \ No newline at end of file diff --git a/test/duobinary_minimal_example.m b/test/duobinary_minimal_example.m index 8866e69..c39183e 100644 --- a/test/duobinary_minimal_example.m +++ b/test/duobinary_minimal_example.m @@ -46,66 +46,55 @@ end Tx_bits = Informationsignal(bitpattern); %%%%% Duobinary %%%%%% - -precode = 1; -db_encode = 0; - close all -Symbols_tx = PAMmapper(M,0).map(Tx_bits); +Symbols_tx = PAMmapper(M,0,"eth_style",0).map(Tx_bits); Symbols_tx.fs = fsym; +precode = db_mode.db_precoded; + %%% precode -if precode - Symbols0 = Duobinary().precode(Symbols_tx); -else - Symbols0 = Symbols_tx; +switch precode + case db_mode.db_precoded + Symbols_tx = Duobinary().precode(Symbols_tx); + case db_mode.db_encoded + Symbols_tx = Duobinary().precode(Symbols_tx); + Symbols_tx = Duobinary().encode(Symbols_tx); + case db_mode.no_db + end -figure;histogram(Symbols0.signal); +for n = 10 -for n = 0:200 + Symbols_rx = Symbols_tx; - if db_encode - Symbols1 = Duobinary().encode(Symbols0); - else - Symbols1 = Symbols0; - end - - Symbols2 = Symbols1; pos = 1; if n~=0 for pos = 1:n po = randi(100); - a = Symbols2.signal(100+pos) == Symbols1.signal(100+po); + a = Symbols_rx.signal(100+pos) == Symbols_tx.signal(100+po); while a == 1 po = po+1; po = randi(100); - a = Symbols2.signal(100+pos) == Symbols1.signal(100+po); + a = Symbols_rx.signal(100+pos) == Symbols_tx.signal(100+po); end - Symbols2.signal(100+pos) = Symbols1.signal(100+po); + Symbols_rx.signal(100+pos) = Symbols_tx.signal(100+po); end end - % disp(Symbols2.signal(100:100+pos)==Symbols1.signal(100:100+pos)) + error_positions = ~(Symbols_rx.signal == Symbols_tx.signal); + error_positions = find(error_positions==1); - %%% encode + switch precode + + case db_mode.db_precoded + Symbols_rx = Duobinary().encode(Symbols_rx); + Symbols_rx = Duobinary().decode(Symbols_rx); + case db_mode.db_encoded + Symbols_rx = Duobinary().decode(Symbols_rx); - if db_encode - - % figure;histogram(Symbols2.signal); - - Symbols3 = Duobinary().decode(Symbols2); - % figure;histogram(Symbols3.signal); - % autoArrangeFigures; - elseif precode - Symbols3 = Duobinary().encode(Symbols2); - Symbols3 = Duobinary().decode(Symbols3); - % figure;histogram(Symbols3.signal); - else - Symbols3 = Symbols2; end - - Rx_bits = PAMmapper(M,0).demap(Symbols3); + + Rx_bits = PAMmapper(M,0).demap(Symbols_rx); %%%%% Check BER of Bit Sequence %%%%%% @@ -113,14 +102,10 @@ for n = 0:200 % disp(['BER: ',sprintf('%.1E',ber),sprintf(' - Num. Err: %.1d',error_num(n+1)-2),' - - PAM-',num2str(M)]); fprintf('n: %d - Num. Err: %.1d \n',n,error_num(n+1)); + end -figure() -hold on -scatter(1:length(Symbols3),Symbols3.signal,1,'.'); -scatter(error_pos,Symbols3.signal(error_pos),14,'o'); - figure(3); clf @@ -128,8 +113,8 @@ clf subplot(2,2,1) hold on title('First Bits') -stairs(Tx_bits.signal(1:100,1),'LineStyle','-','LineWidth',2,'DisplayName','Tx Bits'); -stairs(Rx_bits.signal(1:100,1),'LineWidth',2,'DisplayName','Rx Bits','LineStyle',':') +stairs(Tx_bits.signal(100:150,1),'LineStyle','-','LineWidth',2,'DisplayName','Tx Bits'); +stairs(Rx_bits.signal(100:150,1),'LineWidth',2,'DisplayName','Rx Bits','LineStyle',':') legend subplot(2,2,2) @@ -143,12 +128,12 @@ subplot(2,2,3) hold on title('First Symbols Compare') stairs(Symbols_tx.signal(1:100,1),'LineWidth',2,'DisplayName','Tx Symbols','LineStyle','-') -stairs(Symbols.signal(1:100,1),'LineStyle',':','LineWidth',2,'DisplayName','Rx Symbols'); +stairs(Symbols_rx.signal(1:100,1),'LineStyle',':','LineWidth',2,'DisplayName','Rx Symbols'); legend subplot(2,2,4) hold on title('Last Symbols Compare') stairs(Symbols_tx.signal(end-50:end,1),'LineWidth',2,'DisplayName','Tx Symbols','LineStyle','-') -stairs(Symbols.signal(end-50:end,1),'LineStyle',':','LineWidth',2,'DisplayName','Rx Symbols'); +stairs(Symbols_rx.signal(end-50:end,1),'LineStyle',':','LineWidth',2,'DisplayName','Rx Symbols'); legend \ No newline at end of file