From 5c2e27687d3897abb42295ca7b9d681f2267051b Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Wed, 3 Jun 2026 09:05:33 +0200 Subject: [PATCH] Scattered stuff from Silas during Dissertation --- Classes/00_signals/Signal.m | 381 +++++++------ Classes/01_transmit/ChannelFreqResp.m | 2 +- Classes/04_DSP/Equalizer/Copy_of_VNLE.m | 394 ++++++++++++++ Classes/04_DSP/Equalizer/FFE.m | 6 + Classes/04_DSP/Equalizer/FFE_Kalman.m | 245 --------- .../04_DSP/Equalizer/FFE_Kalman_Feedback.m | 197 ------- Classes/04_DSP/Equalizer/ML_MLSE.m | 500 +----------------- Classes/04_DSP/Equalizer/VNLE.m | 97 +++- Datatypes/channel_model.m | 9 + Functions/EQ_structures/dsp_scope_signal.m | 35 +- Functions/EQ_structures/duobinary_target.m | 6 +- Functions/EQ_structures/ffe.m | 3 +- .../EQ_structures/vnle_postfilter_mlse.m | 4 +- Functions/EQ_visuals/show2Dconstellation.m | 158 +++++- Functions/EQ_visuals/showEQNoisePSD.m | 12 +- Functions/EQ_visuals/showEQNoiseSNR.m | 65 ++- Functions/EQ_visuals/showEQfilter.m | 23 +- Functions/EQ_visuals/showLevelHistogram.m | 115 +++- Functions/EQ_visuals/showMarkovDiagram.m | 35 ++ .../EQ_visuals/showTransitionProbabilities.m | 40 ++ Functions/Lab_helper/loadFreqResp.m | 4 +- Functions/mat2tikz_improved.m | 7 +- Theory/Mapping_Coding/pam_6_states_analysis.m | 18 +- .../partial_response/duobinary_back2back.m | 19 +- .../Signal_Processing/mlse_minimal_example.m | 29 +- .../Auswertung_JLT/final/FIGURE_ROP.m | 15 +- .../Auswertung_JLT/run_dsp_from_db.m | 2 +- projects/Diss/investigate_noise_enhancement.m | 257 +++++++++ projects/Diss/peak_distortion.m | 66 +++ projects/IMDD_base_system/imdd_it.m | 179 +++++-- projects/IMDD_base_system/imdd_model.m | 86 +-- projects/IMDD_base_system/minimal_example.m | 11 +- .../minimal_example_huawei/minimal_example.m | 2 +- .../minimal_example_huawei/ml_mlse_pam.m | 4 +- projects/MLSE_ML_based/rop_evaluation.m | 2 +- .../theoretic_channel_evaluation.m | 34 +- .../simulation_nonlin_dsp.m | 2 +- 37 files changed, 1726 insertions(+), 1338 deletions(-) create mode 100644 Classes/04_DSP/Equalizer/Copy_of_VNLE.m delete mode 100644 Classes/04_DSP/Equalizer/FFE_Kalman.m delete mode 100644 Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m create mode 100644 Datatypes/channel_model.m create mode 100644 Functions/EQ_visuals/showMarkovDiagram.m create mode 100644 Functions/EQ_visuals/showTransitionProbabilities.m create mode 100644 projects/Diss/investigate_noise_enhancement.m create mode 100644 projects/Diss/peak_distortion.m diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 3579fa5..6261b91 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -198,44 +198,65 @@ classdef Signal end %% Add signals from one signal to another, the first object will sustain - function Sum = plus(X,y) - - if isa(y,'Signal') - Sum = X; - Sum.signal = X.signal + y.signal; - elseif isnumeric(y) - Sum = X; - Sum.signal = X.signal + y; - end - - end + function Sum = plus(X,y) + + if isa(X,'Signal') && isa(y,'Signal') + Sum = X; + Sum.signal = X.signal + y.signal; + elseif isa(X,'Signal') && isnumeric(y) + Sum = X; + Sum.signal = X.signal + y; + elseif isnumeric(X) && isa(y,'Signal') + Sum = y; + Sum.signal = X + y.signal; + end + + end %% Add signals from one signal to another, the first object will sustain - function Diff = minus(X,y) - - if isa(y,'Signal') - Diff = X; - Diff.signal = X.signal - y.signal; - elseif isnumeric(y) - Diff = X; - Diff.signal = X.signal - y; - end - - end - - function Product = times(X,y) - - if isa(y,'Signal') - Product = X; - Product.signal = X.signal .* y.signal; - elseif isnumeric(y) - Product = X; - Product.signal = X.signal .* y; - end - - end - - %% Display length + function Diff = minus(X,y) + + if isa(X,'Signal') && isa(y,'Signal') + Diff = X; + Diff.signal = X.signal - y.signal; + elseif isa(X,'Signal') && isnumeric(y) + Diff = X; + Diff.signal = X.signal - y; + elseif isnumeric(X) && isa(y,'Signal') + Diff = y; + Diff.signal = X - y.signal; + end + + end + + function Product = times(X,y) + + if isa(X,'Signal') && isa(y,'Signal') + Product = X; + Product.signal = X.signal .* y.signal; + elseif isa(X,'Signal') && isnumeric(y) + Product = X; + Product.signal = X.signal .* y; + elseif isnumeric(X) && isa(y,'Signal') + Product = y; + Product.signal = X .* y.signal; + end + + end + + function Product = mtimes(X,y) + + if (isa(X,'Signal') && isnumeric(y) && isscalar(y)) || ... + (isnumeric(X) && isscalar(X) && isa(y,'Signal')) + Product = times(X,y); + else + error('Signal:mtimes:UnsupportedOperands', ... + 'Use element-wise .* for Signal multiplication, or scalar * Signal for scaling.'); + end + + end + + %% Display length function return_length = length(obj) %METHOD1 Summary of this method goes here % Detailed explanation goes here @@ -325,7 +346,8 @@ classdef Signal else - obj.signal = resample(obj.signal,options.fs_out,options.fs_in,options.n,options.beta); + [p, q] = rat(options.fs_out / options.fs_in); + obj.signal = resample(obj.signal,p,q,options.n,options.beta); desc = ['resample signal from ', num2str(options.fs_in*1e-9), ' GHz to ', num2str(options.fs_out*1e-9), ' GHz' ]; @@ -344,16 +366,17 @@ classdef Signal obj options.fignum = 2025 options.displayname = ""; - options.color = []; - options.linestyle = '-'; - options.normalizeToNyquist = 0; - options.normalizeToSamplingRate = 0; - options.addDCoffset = 0; - options.normalizeToDC = 0; - options.normalizeTo0dB = 0; - options.show_onesided = false; - options.max_num_lines = []; % Leave empty or omit to disable line rotation - options.fft_length = []; + options.color = []; + options.linestyle = '-'; + options.HandleVisibility (1,1) string {mustBeMember(options.HandleVisibility, ["on","off"])} = "on"; + options.normalizeToNyquist = 0; + options.normalizeToSamplingRate = 0; + options.addDCoffset = 0; + options.normalizeToDC = 0; + options.normalizeTo0dB = 0; + options.show_onesided = false; + options.max_num_lines = []; % Leave empty or omit to disable line rotation + options.fft_length = []; % --- NEW options --- options.useWavelengthAxis (1,1) logical = false % plot x-axis in wavelength options.lambda0_nm (1,1) double = 1310 % center wavelength [nm] @@ -364,23 +387,23 @@ classdef Signal end - useSamplingRateAxis = options.normalizeToSamplingRate ~= 0; - useRadPerSampleAxis = options.normalizeToNyquist ~= 0 && ~useSamplingRateAxis; - - if ~useRadPerSampleAxis && ~useSamplingRateAxis - [p_lin,f_Hz] = pwelch(obj.signal, hanning(options.fft_length), ... - options.fft_length/2, options.fft_length, ... - obj.fs, "centered", "power", "mean"); - f_GHz = f_Hz*1e-9; % keep frequency vector for frequency axis - else - [p_lin,f_rad] = pwelch(obj.signal, hanning(options.fft_length), ... - options.fft_length/2, options.fft_length, ... - "centered", "power", "mean"); - % In normalized modes, pwelch returns rad/sample centered on 0. - % Divide by 2*pi for the f/fs axis where Nyquist is 0.5. - end + useSamplingRateAxis = options.normalizeToSamplingRate ~= 0; + useRadPerSampleAxis = options.normalizeToNyquist ~= 0 && ~useSamplingRateAxis; - % p_lin = movmean(p_lin,4); + if ~useRadPerSampleAxis && ~useSamplingRateAxis + [p_lin,f_Hz] = pwelch(obj.signal, hanning(options.fft_length), ... + options.fft_length/2, options.fft_length, ... + obj.fs, "centered", "power", "mean"); + f_GHz = f_Hz*1e-9; % keep frequency vector for frequency axis + else + [p_lin,f_rad] = pwelch(obj.signal, hanning(options.fft_length), ... + options.fft_length/2, options.fft_length, ... + "centered", "power", "mean"); + % In normalized modes, pwelch returns rad/sample centered on 0. + % Divide by 2*pi for the f/fs axis where Nyquist is 0.5. + end + + p_lin = movmean(p_lin,10); if options.normalizeTo0dB p_lin = p_lin ./ max(p_lin); @@ -392,67 +415,67 @@ classdef Signal end % --- If requested, build wavelength axis from frequency offset --- - if options.useWavelengthAxis && ~useRadPerSampleAxis && ~useSamplingRateAxis - c = physconst('LightSpeed'); % [m/s] - lambda0_m = options.lambda0_nm*1e-9; % center wavelength [m] - f_c = c / lambda0_m; % carrier frequency [Hz] + if options.useWavelengthAxis && ~useRadPerSampleAxis && ~useSamplingRateAxis + c = physconst('LightSpeed'); % [m/s] + lambda0_m = options.lambda0_nm*1e-9; % center wavelength [m] + f_c = c / lambda0_m; % carrier frequency [Hz] % exact mapping f_abs = f_c + f_Hz; % absolute frequency [Hz] lambda_m = c ./ f_abs; % wavelength [m] lambda_nm = lambda_m * 1e9; % wavelength [nm] - % assign axis - x_vec = lambda_nm(:); - x_label = "Wavelength [nm]"; - dc_axis = f_Hz; - - % Sort to ensure axis is ascending - [x_vec, sortIdx] = sort(x_vec, 'ascend'); - p_dbm = p_dbm(sortIdx, :); - dc_axis = dc_axis(sortIdx); - else - % Frequency or normalized axes - if ~useRadPerSampleAxis && ~useSamplingRateAxis - x_vec = f_GHz; - x_label = "Frequency in GHz"; - dc_axis = f_GHz; - elseif useSamplingRateAxis - x_vec = f_rad ./ (2*pi); - x_label = "Normalized Frequency f/fs"; - dc_axis = x_vec; - else - x_vec = f_rad; % normalized frequency in rad/sample - x_label = "Normalized Frequency [rad/sample]"; - dc_axis = x_vec; - end - end - - if options.show_onesided - keep_idx = dc_axis >= 0; - x_vec = x_vec(keep_idx); - dc_axis = dc_axis(keep_idx); - p_dbm = p_dbm(keep_idx, :); - end - - figure(options.fignum); - ax = gca; + % assign axis + x_vec = lambda_nm(:); + x_label = "Wavelength [nm]"; + dc_axis = f_Hz; + + % Sort to ensure axis is ascending + [x_vec, sortIdx] = sort(x_vec, 'ascend'); + p_dbm = p_dbm(sortIdx, :); + dc_axis = dc_axis(sortIdx); + else + % Frequency or normalized axes + if ~useRadPerSampleAxis && ~useSamplingRateAxis + x_vec = f_GHz; + x_label = "Frequency in GHz"; + dc_axis = f_GHz; + elseif useSamplingRateAxis + x_vec = f_rad ./ (2*pi); + x_label = "Normalized Frequency f/fs"; + dc_axis = x_vec; + else + x_vec = f_rad; % normalized frequency in rad/sample + x_label = "Normalized Frequency [rad/sample]"; + dc_axis = x_vec; + end + end + + if options.show_onesided + keep_idx = dc_axis >= 0; + x_vec = x_vec(keep_idx); + dc_axis = dc_axis(keep_idx); + p_dbm = p_dbm(keep_idx, :); + end + + figure(options.fignum); + ax = gca; hold on - p_dbm = p_dbm+options.addDCoffset; - - if options.normalizeToDC - [~,min_idx]=min(abs(dc_axis)); - pow_at_dc = p_dbm(min_idx); - p_dbm = p_dbm-pow_at_dc; - end + p_dbm = p_dbm+options.addDCoffset; + + if options.normalizeToDC + [~,min_idx]=min(abs(dc_axis)); + pow_at_dc = p_dbm(min_idx); + p_dbm = p_dbm-pow_at_dc; + end % p_dbm = movmean(p_dbm,10); for s = 1:min(size(p_dbm)) if isempty(options.color) - plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1); + plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1, 'LineStyle', options.linestyle, 'HandleVisibility', options.HandleVisibility); else - plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1, 'Color', options.color,'LineStyle',options.linestyle); + plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1, 'Color', options.color,'LineStyle',options.linestyle, 'HandleVisibility', options.HandleVisibility); end end @@ -468,27 +491,27 @@ classdef Signal % Axis labels and limits xlabel(x_label); - if options.useWavelengthAxis && ~useRadPerSampleAxis && ~useSamplingRateAxis - xlim([min(x_vec) max(x_vec)]); - else - if ~useRadPerSampleAxis && ~useSamplingRateAxis - % Keep your existing freq handling (you can fine-tune as needed) - % xlim([-128 128]); % example for 256 GSa/s if desired - xlim([min(x_vec) max(x_vec)]); - elseif useSamplingRateAxis - if options.show_onesided - xlim([0, 0.5]); - else - xlim([-0.5, 0.5]); - end - else - if options.show_onesided - xlim([0, pi]); - else - xlim([-pi, pi]); - end - end - end + if options.useWavelengthAxis && ~useRadPerSampleAxis && ~useSamplingRateAxis + xlim([min(x_vec) max(x_vec)]); + else + if ~useRadPerSampleAxis && ~useSamplingRateAxis + % Keep your existing freq handling (you can fine-tune as needed) + % xlim([-128 128]); % example for 256 GSa/s if desired + xlim([min(x_vec) max(x_vec)]); + elseif useSamplingRateAxis + if options.show_onesided + xlim([0, 0.5]); + else + xlim([-0.5, 0.5]); + end + else + if options.show_onesided + xlim([0, pi]); + else + xlim([-pi, pi]); + end + end + end ylabel(ylab); @@ -1007,15 +1030,47 @@ classdef Signal elseif mode == 1 % generate eye diagram using histogram - maxA = max(sig(100:end-100))*1.3; - minA = min(sig(100:end-100))*1.3; + finite_eye = eye_mat(isfinite(eye_mat)); + if isempty(finite_eye) + finite_eye = sig(isfinite(sig)); + end + amp_min = min(finite_eye); + amp_max = max(finite_eye); + amp_center = (amp_max + amp_min) / 2; + amp_span = amp_max - amp_min; + if amp_span == 0 + amp_span = max(abs(amp_center),1); + end + amp_margin = 0.08 * amp_span; + maxA = amp_center + amp_span/2 + amp_margin; + minA = amp_center - amp_span/2 - amp_margin; + if ~isa(obj,'Opticalsignal') && minA < 0 && maxA > 0 + targetStep = max(abs([minA maxA])) / 2; + if targetStep > 0 + stepMagnitude = 10^floor(log10(targetStep)); + normalizedStep = targetStep / stepMagnitude; + if normalizedStep <= 1 + tickStep = stepMagnitude; + elseif normalizedStep <= 2 + tickStep = 2 * stepMagnitude; + elseif normalizedStep <= 5 + tickStep = 5 * stepMagnitude; + else + tickStep = 10 * stepMagnitude; + end + axisLimit = 2 * tickStep; + maxA = axisLimit; + minA = -axisLimit; + end + end % maxA = 0.12; % minA = -0.08; difference= maxA-minA; - data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1; + data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1; + data_ind_y = min(max(data_ind_y,1),histpoints); for n=1:size(data_ind_y,1) nn=histcounts(data_ind_y(n,:),1:histpoints+1); @@ -1041,22 +1096,22 @@ classdef Signal if isa(obj,'Opticalsignal') title(['Optical Eye ',options.displayname]) ylabel("Power in mW"); - y_tickstring = string(linspace(maxA.*1e3,minA.*1e3,6)); - min_ = min(abs(obj.signal(100:end-100)).^2); - max_ = abs(max(obj.signal(100:end-100)).^2); - elseif isa(obj,'Electricalsignal') - title(['Electrical Eye ',options.displayname]) - ylabel("Voltage in V"); - y_tickstring = string(linspace(maxA,minA,6)); - min_ = min(obj.signal(100:end-100)); - max_ = abs(max(obj.signal(100:end-100))); - else - title(['Digital Eye ',options.displayname]) - ylabel("Digital Signal Amplitude"); - y_tickstring = string(linspace(maxA,minA,6)); - min_ = min(obj.signal(100:end-100)); - max_ = abs(max(obj.signal(100:end-100))); - end + yTickValues = linspace(maxA.*1e3,minA.*1e3,5); + min_ = min(abs(obj.signal(100:end-100)).^2); + max_ = abs(max(obj.signal(100:end-100)).^2); + elseif isa(obj,'Electricalsignal') + title(['Electrical Eye ',options.displayname]) + ylabel("Voltage in V"); + yTickValues = linspace(maxA,minA,5); + min_ = min(obj.signal(100:end-100)); + max_ = abs(max(obj.signal(100:end-100))); + else + title(['Digital Eye ',options.displayname]) + ylabel("Digital Signal Amplitude"); + yTickValues = linspace(maxA,minA,5); + min_ = min(obj.signal(100:end-100)); + max_ = abs(max(obj.signal(100:end-100))); + end xlabel('Time in ps') @@ -1094,7 +1149,14 @@ classdef Signal hist_interest_smoth = smooth(hist_interest,20); a = scatter(hist_interest_smoth+posxall,1:length(hist_interest_smoth),4,'.','MarkerEdgeColor','red'); - [pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",10,"NPeaks",M,"MinPeakHeight",30,"MinPeakProminence",10); + minPeakDistance = max(10, floor(histpoints / (2*M))); + minPeakProminence = max(3, 0.05 * max(hist_interest_smoth)); + [pk,loc] = findpeaks(hist_interest_smoth, ... + "MinPeakDistance",minPeakDistance, ... + "NPeaks",M, ... + "MinPeakProminence",minPeakProminence, ... + "SortStr","descend"); + loc = sort(loc); scatter(posxall,loc,'red','Marker','x','LineWidth',2); @@ -1177,13 +1239,14 @@ classdef Signal end - yticks(linspace(0,histpoints,6)); - y_tickstring = sprintfc('%.2f', y_tickstring); - yticklabels(y_tickstring); - - xticks(linspace(0,histpoints_horizontal,6)) - x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12); - xticklabels(x_tickstring); + yTickPositions = linspace(1,histpoints,numel(yTickValues)); + yticks(yTickPositions); + yticklabels(sprintfc('%.2f', yTickValues)); + + xTickValues = linspace(0, 2/fsym, 6) .* 1e12; + xticks(linspace(1,histpoints_horizontal,numel(xTickValues))) + x_tickstring = sprintfc('%.2f', xTickValues); + xticklabels(x_tickstring); % end diff --git a/Classes/01_transmit/ChannelFreqResp.m b/Classes/01_transmit/ChannelFreqResp.m index e6c8a2a..21c068e 100644 --- a/Classes/01_transmit/ChannelFreqResp.m +++ b/Classes/01_transmit/ChannelFreqResp.m @@ -254,7 +254,7 @@ classdef ChannelFreqResp < handle xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on; %%% plot for publication - figure(101); + figure(1996); hold all; box on; title('Magnitude Freq. Response'); diff --git a/Classes/04_DSP/Equalizer/Copy_of_VNLE.m b/Classes/04_DSP/Equalizer/Copy_of_VNLE.m new file mode 100644 index 0000000..acc6c52 --- /dev/null +++ b/Classes/04_DSP/Equalizer/Copy_of_VNLE.m @@ -0,0 +1,394 @@ +classdef Copy_of_VNLE < handle + % Implementation of plain and simple FFE. + % 1) Training mode (stable performance when you use NLMS) + % 2) Decision directed mode + + % Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[25,2,2],"sps",2,"decide",1); + % Somehow it is not possible to use only 1 nonlinear order + + properties + sps % usually 2 + order + e + e_dc + error + + len_tr + mu_tr + epochs_tr + + mu_dd + epochs_dd + mu_dc + + constellation + + decide + save_debug = 0; + debug_struct + + optmize_mus = 0; + mu_optimization + mu_optimization_iter = 0; + + x_norm + ce + ie2 + ie3 + end + + methods + function obj = Copy_of_VNLE(options) + arguments(Input) + + options.sps = 2; + options.order = [15,2,2]; + + options.len_tr = 4096; + options.mu_tr = 0; + options.epochs_tr = 5; + + options.mu_dd = 1e-5; + options.epochs_dd = 5; + options.mu_dc = 0; + + options.decide = false; + options.save_debug = 0; + options.optmize_mus = 0; + + end + + fn = fieldnames(options); + for n = 1:numel(fn) + obj.(fn{n}) = options.(fn{n}); + end + + + obj.error = 0; + obj.e_dc = 0; + + end + + function [X,N] = process(obj, X, D) + + % actual processing of the signal (steps 1. - 3.) + % 1 normalize RMS + X = X.normalize("mode","rms"); + + obj.constellation = unique(D.signal); + obj.x_norm = obj.calcPowerNormalization(X.signal); + obj.ce = obj.calcVNLEMemoryLength(obj.order); + [obj.ie2,obj.ie3] = obj.calcIndiceVectors(obj.order); + + obj.e = zeros( sum(obj.ce) ,1); + obj.e_dc = 0; + + if obj.optmize_mus + obj.optimizeMus(X.signal,D.signal); + obj.e = zeros(sum(obj.ce),1); + obj.e_dc = 0; + end + + % Training Mode + training = 1; + showviz = 0; + obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); + + % Decision Directed Mode + N = X.length; + training = 0; + showviz = 0; + [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz); + + % Output Signal + if obj.decide + X.signal = decision; + else + X.signal = signal; + end + X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym + lbdesc = [num2str(obj.order),' tap FFE']; + X = X.logbookentry(lbdesc); % append to logbook + + N = X; + N = X - D; + + + end + + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) + + arguments + obj + x + d + mu + epochs + N + training + showviz + end + + if all(mu == mu(1)) + % mu = mu(1); + mu = diag(ones(1,sum(obj.ce))*mu(1)); + else + mu = diag([ones(1,obj.ce(1))*mu(1) ... + ones(1,obj.ce(2))*mu(2) ... + ones(1,obj.ce(3))*mu(3) ]); + end + + x = [zeros(floor(obj.order(1)/2),1); x; zeros(obj.order(1),1)]; + + if showviz + f = figure(111); + subplot(2,2,1:2); + hold on + a = scatter(1:numel(x),x,1,'.'); + a2 = scatter(1,1,1,'.'); + a3 = scatter(1,1,2,'.'); + a4 = xline(1); + ylim([-3 3]) + xlim([0 length(x)]); + subplot(2,2,3:4) + c = stem(obj.e); + ylim([-1 1]) + drawnow + end + + for epoch = 1 : epochs + symbol = 0; + for sample = 1 : obj.sps : N + + symbol = symbol+1; + + + % x_in = x(obj.order(1)+sample+(obj.sps-1):-1:sample+obj.sps); + x_in = x(obj.order(1)+sample-1:-1:sample); + x_in = obj.calcVNLENonlinVecs(x_in,obj.ie2,obj.ie3,obj.order,obj.x_norm); + + y(symbol,1) = obj.e_dc + obj.e.' * x_in; % Calculating output of LMS __ * | + + if training + err = y(symbol) - d(symbol); % Instantaneous error + else + [~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); + err = y(symbol) - d_hat(symbol); % Instantaneous error + end + + if ~all(mu==0,'all') %mu has not only zeros + obj.e = obj.e - ( (mu * x_in) * err ) ; % Weight update rule of LMS + else + normalizationfactor = (x_in.' * x_in); + obj.e = obj.e - err * x_in / normalizationfactor; % Weight update rule of NLMS + end + if obj.mu_dc ~= 0 + obj.e_dc = obj.e_dc - obj.mu_dc * err; + end + + if mod(sample,100) == 1 && showviz + a2.XData = 1:2*numel(y); + a2.YData = repelem(y, 2); + a3.XData = 1:2*numel(d_hat); + a3.YData = repelem(d_hat, 2); + a4.Value = sample; + % b.YData = x(symbol:symbol+500); + c.YData = obj.e; + drawnow; + end + + obj.error(epoch,symbol) = err * err'; % Instantaneous square error + if obj.save_debug + obj.debug_struct.error(epoch,symbol) = err * err'; + if training + obj.debug_struct.error_tr(epoch,symbol) = err * err'; + end + end + + end + end + + end + + function optimizeMus(obj,x,d) + mu_range = [1e-5, 1e-2]; + mu_dc_range = [1e-5, 1e-1]; + + vars = [optimizableVariable("mu_tr",mu_range,"Transform","log"), ... + optimizableVariable("mu_dd",mu_range,"Transform","log")]; + optimize_mu_dc = obj.mu_dc ~= 0; + if optimize_mu_dc + vars = [vars, optimizableVariable("mu_dc",mu_dc_range,"Transform","log")]; + end + + obj.mu_optimization_iter = 0; + obj.mu_optimization = bayesopt(@(p)obj.muObjective(p,x,d),vars, ... + "MaxObjectiveEvaluations",10, ... + "AcquisitionFunctionName","expected-improvement-plus", ... + "IsObjectiveDeterministic",false, ... + "Verbose",0, ... + "PlotFcn",[]); + obj.mu_tr = obj.mu_optimization.XAtMinObjective.mu_tr; + obj.mu_dd = obj.mu_optimization.XAtMinObjective.mu_dd; + if optimize_mu_dc + obj.mu_dc = obj.mu_optimization.XAtMinObjective.mu_dc; + end + + objective_db = 10*log10(obj.mu_optimization.MinObjective); + if optimize_mu_dc + fprintf("\nVNLE mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... + obj.mu_tr,obj.mu_dd,obj.mu_dc,obj.mu_optimization.MinObjective,objective_db); + else + fprintf("\nVNLE mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... + obj.mu_tr,obj.mu_dd,obj.mu_optimization.MinObjective,objective_db); + end + end + + function objective = muObjective(obj,params,x,d) + old_debug = obj.save_debug; + old_mu_dc = obj.mu_dc; + obj.save_debug = 1; + optimize_mu_dc = ismember("mu_dc",string(params.Properties.VariableNames)); + if optimize_mu_dc + obj.mu_dc = params.mu_dc; + end + + obj.e = zeros(sum(obj.ce),1); + obj.e_dc = 0; + obj.debug_struct = struct(); + obj.equalize(x,d,params.mu_tr,obj.epochs_tr,obj.len_tr,1,0); + obj.equalize(x,d,params.mu_dd,obj.epochs_dd,numel(x),0,0); + + objective = mean(obj.debug_struct.error(end,:),"omitnan"); + if ~isfinite(objective) + objective = inf; + end + + objective_db = 10*log10(objective); + obj.mu_optimization_iter = obj.mu_optimization_iter + 1; + if optimize_mu_dc + fprintf("\rVNLE mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,params.mu_dc,objective,objective_db); + else + fprintf("\rVNLE mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,objective,objective_db); + end + obj.save_debug = old_debug; + obj.mu_dc = old_mu_dc; + end + + %% Functions needed During Adaption + function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) + % These are the second and third order input signal products of the VNLE EQ + % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3) + l1=length(x_in_block); + l2=length(I_2); + l3=length(I_3); + final_length = l1+l2+l3; + + x_in_vnle_format = zeros(final_length,1); + + idx = l1; + x_in_vnle_format(1:idx) = x_in_block; + + if N_(2) > 0 + delta_2 = round((N_(1)-N_(2)) / 2); + input_vec_se = x_in_block(delta_2:end) / norm_(2); %TODO normalization step + + % Extract columns from I_2 + col1 = input_vec_se(I_2(:,1)); + col2 = input_vec_se(I_2(:,2)); + + x2 = col1 .* col2; + x_in_vnle_format(idx+1:idx+l2) = x2; + end + + if N_(3) > 0 + delta_3 = round((N_(1)-N_(3))/2); + input_vec_th = x_in_block(delta_3:end) / norm_(3); + + % Extract columns from I_3 + col1 = input_vec_th(I_3(:,1)); + col2 = input_vec_th(I_3(:,2)); + col3 = input_vec_th(I_3(:,3)); + + % Perform matrix multiplication + x3 = col1 .* col2 .* col3; + + idx = idx+l2; + x_in_vnle_format(idx+1:idx+l3) = x3; + end + + end + + %% Functions needed for Preparation + function [C] = calcVNLEMemoryLength(~,N) + + %calculates the memory length of VNLE + C = zeros(size(N)); + + for o = 1:numel(N) + switch o + case 1 + C(o) = N(o); + case 2 + C(o) = N(o)*(N(o)+1) / 2; + case 3 + C(o) = N(o)*(N(o)+1)*(N(o)+2) / 6; + end + end + + end + + function [indvec2nd, indvec3rd] = calcIndiceVectors(~,N) + + % Init vectors of 2nd and 3rd order coefficient indices -> + % yield combination with + indvec2nd=[]; + indvec3rd=[]; + for o = 2:numel(N) + n = N(o); + v = 1:n; % Ursprünglicher Vektor + row = 1; + + % Schleifen zur Generierung des Indize Vektors + switch o + + case 2 + + indvec2nd = zeros(n*(n+1)/2, o); + for i = 1:n + for j = i:n + indvec2nd(row, :) = [v(i) v(j)]; + row = row + 1; + end + end + + case 3 + + indvec3rd = zeros(n*(n+1)*(n+2)/6, 3); + for i = 1:n + for j = i:n + for k = j:n + indvec3rd(row, :) = [v(i) v(j) v(k)]; + row = row + 1; + end + end + end + end + end + + end + + function powerNorm = calcPowerNormalization(~,v) + + powerNorm(1) = sqrt(mean(abs(v ).^2)); + powerNorm(2) = sqrt(mean(abs(v.^2).^2)); + powerNorm(3) = sqrt(mean(abs(v.^3).^2)); + + end + + end +end + diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index f3f822d..9ac0ebd 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -9,6 +9,12 @@ classdef FFE < handle % 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); + % eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ... + % "mu_dd",1e-1,"mu_tr",0.4,"order",50, ... + % "sps",2,"decide",0,"optmize_mus",1,"dd_mode",1, ... + % "adaption_technique","nlms","mu_dc",1.021e-05); + + properties sps % usually 2 order diff --git a/Classes/04_DSP/Equalizer/FFE_Kalman.m b/Classes/04_DSP/Equalizer/FFE_Kalman.m deleted file mode 100644 index 62e7b33..0000000 --- a/Classes/04_DSP/Equalizer/FFE_Kalman.m +++ /dev/null @@ -1,245 +0,0 @@ -classdef FFE_Kalman < handle - % Implementation of plain and simple FFE. - % 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); - - properties - sps % usually 2 - order - e - error - - len_tr - mu_tr - epochs_tr - - mu_dd - epochs_dd - - constellation - - decide - end - - methods - function obj = FFE_Kalman(options) - arguments(Input) - - options.sps = 2; - options.order = 15; - - options.len_tr = 4096; - options.mu_tr = 0; - options.epochs_tr = 5; - - options.mu_dd = 1e-5; - options.epochs_dd = 5; - - options.decide = false; - - end - - fn = fieldnames(options); - for n = 1:numel(fn) - obj.(fn{n}) = options.(fn{n}); - end - - obj.e = zeros(obj.order,1); - obj.error = 0; - - end - - function [X,Noi] = process(obj, X, D) - - % actual processing of the signal (steps 1. - 3.) - % 1 normalize RMS - X = X.normalize("mode","rms"); - - obj.constellation = unique(D.signal); - - % Training Mode - training = 1; - showviz = 0; - obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); - - % Decision Directed Mode - n = X.length; - training = 0; - showviz = 0; - [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training,showviz); - - % Output Signal - if obj.decide - X.signal = decision; - else - X.signal = signal; - end - - X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym - lbdesc = [num2str(obj.order),' tap FFE']; - X = X.logbookentry(lbdesc); % append to logbook - - Noi = X; - Noi = X - D; - - - end - - function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz) - - arguments - obj - x - d - mio - epochs - N - training - showviz - end - - x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; - - for epoch = 1 : epochs - - - - symbol = 0; - % Initialization of Kalman filter variables - A = 1; % State transition matrix - H = 1; % Observation matrix - Q = 1e-4; % Process noise covariance - R = 1e-1; % Measurement noise covariance - P = 1; % Initial error covariance - mpi_est = 0; % Initial estimate for MPI noise - K = 0; % Kalman gain - subtract_mpi_est = 1; - - for sample = 1 : obj.sps : N - symbol = symbol + 1; - - % Get the current input sample and the equalizer output - U = x(obj.order + sample - 1 : -1 : sample); - - if subtract_mpi_est || ~training - y(symbol,1) = (obj.e.' * U) - mpi_est .* 1 ; % Subtract MPI estimate - else - y(symbol,1) = obj.e.' * U; - end - - % Decision and error calculation - if training - d_hat(symbol,1) = d(symbol); - else - [~, symbol_idx] = min(abs(y(symbol) - obj.constellation)); % Closest constellation point - d_hat(symbol,1) = obj.constellation(symbol_idx); - end - - err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous residual error - - true_err(symbol) = y(symbol) - d(symbol); - - % Kalman filter update to track the MPI noise - % Prediction step - P = A * P * A' + Q; % Update error covariance - K = P * H' / (H * P * H' + R); % Kalman gain - - % Update step - mpi_est_new = mpi_est + K * (err(symbol) - H * mpi_est); % MPI noise estimation - alpha=0; - mpi_est = alpha * mpi_est + (1 - alpha) * mpi_est_new; - - P = (1 - K * H) * P; % Update error covariance - - % Subtract MPI noise from the signal - if subtract_mpi_est || ~training - y(symbol) = y(symbol);% - mpi_est; - else - - end - - % Equalizer weight update (LMS or NLMS) - if mio ~= 0 - obj.e = obj.e - (mio * err(symbol) * U); % LMS weight update - else - normalizationfactor = (U.' * U); - obj.e = obj.e - err(symbol) * U / normalizationfactor; % NLMS weight update - end - - % Store MPI estimate for visualization - mpi_estimates(symbol) = mpi_est; - - Kgain(symbol) = K; - P_(symbol) = P; - H_(symbol) = H; - - - end %symbols - end %epoch - - if ~training - - if 1 - % figure; - % subplot(2,2,1) - % hold on - % scatter(1:numel(y),y,1,'.'); - % plot(1:numel(mpi_estimates), mpi_estimates); - % subplot(2,2,3) - % plot(1:numel(true_err), true_err); - % subplot(2,2,2) - % scatter(1:numel(y),y'-mpi_estimates,1,'.'); - % xlabel('Sample Index'); - % ylabel('MPI Noise Estimate'); - % title('MPI Noise Estimation Over Time (Kalman Filter)'); - % grid on; - - figure(111); - - subplot(3,1,1) - hold on - cla - scatter(1:numel(y),y,1,'.'); - plot(1:numel(mpi_estimates), mpi_estimates,'DisplayName','mpi est'); - ylim([-2,2]); - - subplot(3,1,2) - cla - plot(1:numel(true_err), true_err,'DisplayName','true err'); - ylim([-1,1]); - - subplot(3,1,3) - plot(1:numel(true_err), true_err-mpi_estimates,'DisplayName',['diff rms: ',num2str(rms(true_err-mpi_estimates))]); - ylim([-1,1]); - legend - - figure(333) - hold on - von = 5000; - bis = 15000; - plot(von:bis,true_err(von:bis),'DisplayName','Error') - plot(von:bis,movmean(true_err(von:bis), [30,30]),'DisplayName','Error') - plot(von:bis,mpi_estimates(von:bis),'DisplayName','Est','LineWidth',2); - legend - - figure(222) - hold on - crr = xcorr(mpi_estimates,true_err,'normalized'); - plot(crr); - - end - - % y = y-mpi_estimates'; - - - end - - - - end - - end -end - diff --git a/Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m b/Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m deleted file mode 100644 index 68875aa..0000000 --- a/Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m +++ /dev/null @@ -1,197 +0,0 @@ -classdef FFE_Kalman_Feedback < handle - % Implementation of plain and simple FFE. - % 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); - - properties - sps % usually 2 - order - e - error - - len_tr - mu_tr - epochs_tr - - mu_dd - epochs_dd - - constellation - - decide - end - - methods - function obj = FFE_Kalman_Feedback(options) - arguments(Input) - - options.sps = 2; - options.order = 15; - - options.len_tr = 4096; - options.mu_tr = 0; - options.epochs_tr = 5; - - options.mu_dd = 1e-5; - options.epochs_dd = 5; - - options.decide = false; - - end - - fn = fieldnames(options); - for n = 1:numel(fn) - obj.(fn{n}) = options.(fn{n}); - end - - obj.e = zeros(obj.order,1); - obj.error = 0; - - end - - function [X,Noi] = process(obj, X, D) - - % actual processing of the signal (steps 1. - 3.) - % 1 normalize RMS - X = X.normalize("mode","rms"); - - obj.constellation = unique(D.signal); - - % Training Mode - training = 1; - showviz = 0; - obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training); - - % Decision Directed Mode - n = X.length; - training = 0; - showviz = 0; - [signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training); - - % Output Signal - if obj.decide - X.signal = decision; - else - X.signal = signal; - end - - X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym - lbdesc = [num2str(obj.order),' tap FFE']; - X = X.logbookentry(lbdesc); % append to logbook - - Noi = X; - Noi = X - D; - - - end - - - - - function [y, d_hat] = equalize(obj, x, d, mio, epochs, N, training) - - arguments - obj - x - d - mio - epochs - N - training - end - - % Initialize Kalman filter variables - A = 0.7; % State transition matrix - H = 1; % Observation matrix - Q = 1e-2; % Process noise covariance - R = 1e-3; % Measurement noise covariance - P = 1; % Initial error covariance - mpi_est = 0; % Initial estimate for MPI noise - K = 0; % Kalman gain - - % Error buffers for parallelized structure - parallel_depth = 100; % Depending on your parallelization depth - err_buffer = zeros(parallel_depth, 1); % Store collected errors - mpi_est_buffer = zeros(parallel_depth, 1); % Store MPI estimates - - x = [zeros(floor(obj.order/2), 1); x; zeros(obj.order, 1)]; - - for epoch = 1 : epochs - symbol = 0; - - for sample = 1 : obj.sps : N - symbol = symbol + 1; - - U = x(obj.order + sample - 1 : -1 : sample); % Input window for FFE - - % Subtract the estimated MPI noise from the input signal before FFE - x_mpi_reduced = U; - y(symbol, 1) = obj.e.' * x_mpi_reduced; % Output of FFE - - y(symbol, 1) = y(symbol, 1) - mpi_est; - - % Decision and error calculation - if training - [~, symbol_idx] = min(abs(d(symbol) - obj.constellation)); % Closest constellation point - d_hat(symbol, 1) = d(symbol); - else - [~, symbol_idx] = min(abs(y(symbol) - obj.constellation)); % Closest constellation point - d_hat(symbol, 1) = obj.constellation(symbol_idx); - end - - % Calculate the residual error - err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous residual error - - % Collect error for feedback after 'parallel_depth' samples - err_buffer(mod(symbol, parallel_depth) + 1) = err(symbol); - - % Update Kalman filter based on the accumulated errors every parallel_depth samples - if mod(symbol, parallel_depth) == 0 - % Compute the mean error over the last 'parallel_depth' samples - avg_error = mean(err_buffer); - - % Prediction step (Kalman filter) - P = A * P * A' + Q; % Update the error covariance - K = P * H' / (H * P * H' + R); % Compute Kalman gain - - % Update MPI noise estimate - mpi_est = mpi_est + K * (avg_error - H * mpi_est); % Update based on averaged error - P = (1 - K * H) * P; % Update error covariance - - end - - % Store MPI estimate for visualization or debugging - k_buffer(symbol) = K; - mpi_est_buffer(symbol) = mpi_est; - - % FFE weight update using NLMS - if mio ~= 0 - obj.e = obj.e - (mio * err(symbol) * U); % LMS weight update - else - normalizationfactor = (U.' * U); - obj.e = obj.e - err(symbol) * U / normalizationfactor; % NLMS weight update - end - - end - end - - if ~training - - figure(111); - - subplot(3,1,1) - hold on - cla - scatter(1:numel(y),y,1,'.'); - plot(1:numel(mpi_est_buffer), mpi_est_buffer,'DisplayName','mpi est'); - - end - end - - end - -end - - diff --git a/Classes/04_DSP/Equalizer/ML_MLSE.m b/Classes/04_DSP/Equalizer/ML_MLSE.m index c5c4860..a27712b 100644 --- a/Classes/04_DSP/Equalizer/ML_MLSE.m +++ b/Classes/04_DSP/Equalizer/ML_MLSE.m @@ -1,499 +1,3 @@ -% classdef ML_MLSE < handle -% % ALGORITHM DESCRIBED IN: -% % W. Lanneer and Y. Lefevre, “Machine Learning-Based Pre-Equalizers for -% % Maximum Likelihood Sequence Estimation in High-Speed PONs,” -% % in 2023 31st European Signal Processing Conference -% -% % Further ML Refs: -% % https://machinelearningmastery.com/cross-entropy-for-machine-learning/ -% % https://docs.pytorch.org/docs/stable/generated/torch.nn.CrossEntropyLoss.html -% -% % The central idea is to overcome the (white-) noise assumption within the previously described -% % Viterbi algorithm, more precisely a closed-loop optimization is proposed that finds a suitable -% % filter-set to directly compute the branch metrics c_k (s,s^' ). These can directly be used to -% % carry out the conventional Viterbi algorithm. The system consists of S^L S=F linear FIR filters, -% % combined with one bias coefficient respectively. These filters take the received input samples to -% % compute the branch metrics estimates (c_k ) ̂(s,s^' ) according toThe central idea is to overcome -% % the (white-) noise assumption within the previously described Viterbi algorithm, more precisely -% % a closed-loop optimization is proposed that finds a suitable filter-set to directly compute the -% % branch metrics c_k (s,s^' ). These can directly be used to carry out the conventional Viterbi -% % algorithm. The system consists of S^L S=F linear FIR filters, combined with one bias coefficient -% % respectively. These filters take the received input samples to compute the branch metrics -% % estimates. Finally, the usual Viterbi is carried out... -% -% % Recommended Settings and some findings: -% -% % Requires many training epochs. According to ML people, 100,200 or -% % even up to 1000 epochs are normal for ML-convergence -% -% % The mu parameter _can_ be adaptive - using the cross entropy and when -% % analyzing the isolated training it looks very promisig. However, is -% % later use I found this is not as stable as a fixed learning rate. -% % mu = 0.1 worked good for me -% -% % Longer orders/ filter length are not always better. For me order=11 -% % was good. -% -% % Delay factor (delta) is good when the order is also increased. With -% % order = 11, a delta of =4 shows good results -% -% properties -% sps % usually 2 -% order -% e -% e_tr -% error -% -% len_tr -% mu_tr -% epochs_tr -% -% dd_mode % 1 or 0 to set DD-mode on or off -% mu_dd %weight update in dd mode -% epochs_dd -% -% adaptive_mu -% -% constellation -% -% L %viterbi memory length -% -% alpha -% DIR -% DIR_flip -% trellis_states -% -% traceback_depth -% -% % --- Added internal class variables used later --- -% S -% Nf -% delta -% nStates -% nFeasible -% combs -% first_sym -% last_sym -% valid -% valid_to_idx -% valid_from_idx -% w -% -% % --- New: fast state lookup --- -% true_to_state_idx -% state_dict % containers.Map: key(sequence)->state index -% key_fmt = '%.8g_'; % key format for sequence strings -% nSym % |constellation| -% -% ber = [] -% ce = ones(1,1); -% end -% -% methods -% function obj = ML_MLSE(options) -% arguments(Input) -% -% options.sps = 2; -% options.order = 15; -% -% options.len_tr = 4096; -% options.mu_tr = 0; -% options.epochs_tr = 5; -% -% options.dd_mode = 1; -% options.mu_dd = 1e-5; -% options.epochs_dd = 5; -% -% options.adaptive_mu = 1; -% -% options.delta = 0; -% options.traceback_depth = 1024; -% -% options.L = 1 -% -% end -% -% fn = fieldnames(options); -% for n = 1:numel(fn) -% obj.(fn{n}) = options.(fn{n}); -% end -% -% obj.e = zeros(obj.order,1); -% obj.error = 0; -% end -% -% function [X,X_viterbi] = process(obj, X, D) -% -% % actual processing of the signal (steps 1. - 3.) -% % 1 normalize RMS -% X = X.normalize("mode","rms"); -% -% % Use sorted constellation for deterministic mapping -% obj.constellation = sort(unique(D.signal),'ascend'); -% obj.nSym = numel(obj.constellation); -% -% if length(X)/length(D) ~= obj.sps -% warning('Signal length does not fit to reference!'); -% end -% -% % ============================================================== -% % INITIALIZATION (only before final epoch and detection mode) -% % ============================================================== -% -% % --- Parameters -% obj.S = numel(obj.constellation); % alphabet size -% obj.Nf = obj.order*obj.sps; % filter length -% % obj.delta = 3;%ceil(obj.Nf/2); % delay parameter -% obj.nStates = obj.S^obj.L; -% obj.nFeasible = obj.nStates*obj.S; -% -% % --- Trellis mapping -% obj.trellis_states = reshape(obj.constellation,1,[]); -% pre_comb_mat = repmat(obj.trellis_states, obj.L, 1); -% pre_comb_cell = mat2cell(pre_comb_mat, ones(1,obj.L), size(pre_comb_mat,2)); -% obj.combs = fliplr(combvec(pre_comb_cell{:}).'); % rows: states, columns: [x_k, x_{k-1}, ...] -% obj.first_sym = obj.combs(:,1); -% obj.last_sym = obj.combs(:,end); -% obj.nStates = size(obj.combs,1); -% -% % --- Valid transitions -% obj.valid = false(obj.nStates); -% for from = 1:obj.nStates -% for to = 1:obj.nStates -% if all(obj.combs(to,2:end) == obj.combs(from,1:end-1)) -% obj.valid(to,from) = true; -% end -% end -% end -% [obj.valid_to_idx, obj.valid_from_idx] = find(obj.valid); -% -% % --- Allocate vectors and weights -% % !! IF SHAPE FIT, then we already have smth there an we want -% % to start with the existing fitler-set -% if isempty(obj.w) || any(size(obj.w) ~= [obj.Nf+1,obj.nFeasible]) -% obj.w = zeros(obj.Nf+1,obj.nFeasible); % filter weights per transition + bias tap -% obj.w = randn(obj.Nf+1,obj.nFeasible); -% end -% -% % --- Precompute dictionary for fast state lookup (sequence -> state) -% keys = cell(obj.nStates,1); -% for i = 1:obj.nStates -% keys{i} = obj.seq_key(obj.combs(i,:)); % combs row is already [x_k, x_{k-1}, ...] -% end -% obj.state_dict = containers.Map(keys, 1:obj.nStates); -% -% % ============================================================== -% % TRAINING -% % ============================================================== -% -% % Training Mode -% n = obj.len_tr; -% training = 1; -% obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,n,training); -% obj.e_tr = obj.e; -% -% % ============================================================== -% % DD-Mode / Fixed Mode -% % ============================================================== -% -% % Decision Directed Mode -% n = X.length; -% training = 0; -% [y,y_vit]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,n,training); -% -% X_viterbi = X; -% -% X.signal = y; -% X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym -% lbdesc = [num2str(obj.order),' tap FFE']; -% X = X.logbookentry(lbdesc); % append to logbook -% -% X_viterbi.signal = y_vit; -% X_viterbi.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym -% lbdesc = [num2str(obj.order),'order FFE + PF + Viterbi']; -% X_viterbi = X_viterbi.logbookentry(lbdesc); % append to logbook -% end -% -% function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training) -% % ============================================================== -% % FFE + Whitening + ML-Based Branch Metric Estimation + Viterbi -% % ============================================================== -% debug = 1; -% showPlots = 1; -% -% % --- Input padding and preallocation -% y = zeros(N,1); -% -% % number of symbol steps in this block -% nSymbols = ceil(N/obj.sps); -% -% for epoch = 1:epochs -% -% % state metrics (log-domain costs): keep as column [nStates×1] -% pm = zeros(obj.nStates,1); % v_{k-1}(s′) -% c_hat = zeros(1,obj.nFeasible); -% v_tilde = zeros(1,obj.nFeasible); -% pred = zeros(nSymbols, obj.nStates, 'uint32'); -% pm_sto = nan(obj.nStates, nSymbols,'like',pm); -% CE_accum = 0; -% -% -% %%% START IDX -% if training -% max_start = length(x) - ( (ceil(N/obj.sps)-1)*obj.sps + 1 ); -% max_start = max(1, max_start); % safety -% start_sample = randi([1, max_start], 1); %rnd training; not really good -% start_sample = 1; -% end_sample = start_sample + (ceil(N/obj.sps)-1)*obj.sps; -% else -% start_sample = 1;%obj.len_tr; -% end_sample = N; -% end -% -% start_symbol = 1 + floor((start_sample - 1)/obj.sps); % ABSOLUTE symbol index -% -% if numel(d) >= obj.L && start_symbol >= obj.L -% init_seq = d(start_symbol-obj.L+1 : start_symbol); % [d_k-L+1 ... d_k] -% true_to_state_idx = obj.state_dict(obj.seq_key(flip(init_seq))); % [d_k ... d_k-L+1] -% else -% % Not enough history – fall back to state 1 -% true_to_state_idx = uint32(1); -% end -% -% symbol = 0; -% for sample = start_sample:obj.sps:end_sample -% symbol = symbol + 1; -% k = symbol; -% sym_idx = start_symbol + (symbol - 1); -% -% % --- Build Δ-delayed observation window y_k -% i1 = sample - obj.Nf + 1 + obj.delta; -% i2 = sample + obj.delta; -% buf = x(max(1,i1):min(length(x),i2)); -% padL = max(0,1 - i1); -% padR = max(0,i2 - length(x)); -% yk = [zeros(padL,1); buf(:); zeros(padR,1)]; % Nf×1 -% yk = [yk;1]; -% -% % --- Predict branch metrics for all feasible transitions: c_hat -% c_hat = (yk.' * obj.w); % [1×nFeasible] -% c_hat = c_hat.'; % [nFeasible×1] -% -% % --- Extended path metrics: v_tilde = pm(from) + c_hat -% % normalize pm to avoid growth (invariant to additive const) -% pm = pm - min(pm); -% v_tilde = pm(obj.valid_from_idx) + c_hat; % [nFeasible×1] -% -% % ===== Gradient update (Algorithm 1) ===== -% -% if 1 %training -% % --- allocate storage once -% if epoch == 1 && symbol == 1 -% obj.true_to_state_idx = ones(ceil(N/obj.sps),1,'uint32'); -% end -% -% % --- previous "to" becomes current "from" -% if symbol > 1 -% true_from_state_idx = obj.true_to_state_idx(symbol-1); -% else -% true_from_state_idx = 1; -% end -% -% % --- compute or reuse "to" state -% if epoch == 1 -% % only compute in first epoch -% if sym_idx >= obj.L -% key_to = obj.seq_key(flip(d(sym_idx-obj.L+1 : sym_idx))); -% if isKey(obj.state_dict, key_to) -% obj.true_to_state_idx(symbol) = obj.state_dict(key_to); -% else -% obj.true_to_state_idx(symbol) = true_from_state_idx; -% end -% else -% obj.true_to_state_idx(symbol) = true_from_state_idx; -% end -% end -% -% % --- reuse cached state from second epoch onward -% true_to_state_idx = obj.true_to_state_idx(symbol); -% -% % --- ensure valid (from,to) -% dirac = zeros(obj.nFeasible,1); -% mask = obj.valid_from_idx==true_from_state_idx & ... -% obj.valid_to_idx ==true_to_state_idx; -% if any(mask) -% dirac(mask) = 1; -% else -% idx = find(obj.valid_from_idx==true_from_state_idx,1,'first'); -% dirac(idx) = 1; -% obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx); -% end -% -% -% -% -% % softmax over -v_tilde (numerically safe shift) -% v_shift = -(v_tilde - min(v_tilde)); % shift to small positive numbers -% v_shift = min(v_shift, 100); % clamp exponent argument (≈ exp(50)=3e21) -% expv = exp(v_shift); -% p = expv ./ (sum(expv) + eps); -% -% % for logging only: -% CE_symbol(symbol) = -log(p(dirac==1) + eps); -% -% if sym_idx > obj.L -% CE_smooth(symbol) = 0.01*CE_symbol(symbol) + 0.99*CE_smooth(symbol-1); -% else -% if epoch > 1 -% CE_smooth(symbol) = obj.ce(end); %use ce from last epoch or =1 for very first round?! -% else -% CE_smooth(symbol) = CE_symbol(symbol); -% end -% end -% -% CE_accum = CE_symbol(symbol) + CE_accum; -% -% -% % gradient term (t - p) -% dmp = (dirac - p)'; % 1×nFeasible -% -% % Per-feature gradient; implicit expansion gives (Nf+1)×nFeasible -% dL_Dw = (yk) .* dmp; -% -% % Start updates only when the ABSOLUTE symbol index has ≥ L history -% if sym_idx >= obj.L -% if obj.adaptive_mu -% mu_eff = CE_smooth(sym_idx); -% mu_eff = max(min(mu_eff, 0.2), 1e-4); -% else -% mu_eff = mu; -% end -% -% obj.w = obj.w - mu_eff .* dL_Dw; % (Nf+1)×nFeasible -% end -% -% % if debug && epoch > 2 -% % figure(100); -% % subplot(4,1,1); -% % heatmap(p'); -% % title('Probs') -% % subplot(4,1,2); -% % heatmap(dmp); -% % title('Update') -% % subplot(4,1,3); -% % heatmap(dL_Dw); -% % title('Update') -% % subplot(4,1,4); -% % heatmap(bj.w); -% % title('Update') -% % -% % end -% -% end -% -% -% -% % --- Compare-Select (matrix form, min of costs) -% v_tilde_mat = inf(obj.nStates, obj.nStates); -% v_tilde_mat(obj.valid) = v_tilde; -% [pm_next, pred(k,:)] = min(v_tilde_mat, [], 2); -% -% % re-center to keep metrics bounded (decision-invariant) -% pm_next = pm_next - min(pm_next); -% -% pm = pm_next; -% pm_sto(:,symbol) = pm; -% end -% -% % --- Traceback (full; you can window with traceback_depth if desired) -% [~, s_end] = min(pm); -% viterbi_path = zeros(symbol,1,'uint32'); -% viterbi_path(symbol) = s_end; -% for n = symbol:-1:2 -% viterbi_path(n-1) = pred(n, viterbi_path(n)); -% end -% -% y_ref = d(start_symbol:end); -% y = obj.first_sym(viterbi_path); -% -% if debug && training -% sym_start = start_symbol; -% sym_end = start_symbol + symbol - 1; -% ref_slice = d(sym_start : sym_end); -% err = sum(y ~= ref_slice(1:numel(y))); -% -% try -% ref_bits = PAMmapper(obj.S,0).demap(ref_slice); -% eq_bits = PAMmapper(obj.S,0).demap(y); -% [~, ~, ber, ~] = calc_ber(ref_bits, eq_bits, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1); -% fprintf('Epoch: %d - BER: %.1e \n',epoch, ber); -% obj.ber(epoch) = ber; -% catch -% ser = err./length(y); -% fprintf('Epoch: %d - SER: %.1e \n',epoch, ser); -% end -% -% obj.ce(epoch) = CE_accum./symbol; -% -% if showPlots -% figure(10);clf -% subplot(3,2,1:2); -% heatmap(obj.w); -% title('Filter') -% -% subplot(3,2,3); -% v_tildemat = NaN(obj.nStates, obj.nStates); -% v_tildemat(obj.valid) = v_tilde; % log-domain scores -% heatmap(v_tildemat); -% title('Path Metrics (v_tilde)') -% -% subplot(3,2,4); -% scatter(1:symbol,pm_sto,1,'.') -% title('Path Metric Winners') -% -% subplot(3,2,5);hold on -% scatter(1:symbol,CE_symbol,1,'.'); -% scatter(1:symbol,CE_smooth,1,'.') -% title('Cross Entropy') -% -% subplot(3,2,6); hold on -% -% % Left y-axis: Cross Entropy (linear) -% yyaxis left -% scatter(1:length(obj.ce), obj.ce, 10, 's', 'filled') -% ylabel('Cross Entropy') -% -% % Right y-axis: BER (logarithmic) -% yyaxis right -% scatter(1:length(obj.ber), obj.ber, 10, 'd', 'filled') -% set(gca, 'YScale', 'log') -% ylabel('BER (log scale)') -% -% xlim([1, epochs]) -% xlabel('Epoch') -% title('Cross Entropy // BER') -% grid on -% -% drawnow -% end -% end -% end -% end -% end -% -% methods (Access=private) -% function k = seq_key(obj, seq) -% % Build a stable key string for a sequence row vector in the *same order as combs rows* ([x_k, x_{k-1}, ...]) -% % Use rounding via sprintf to avoid floating-point issues. -% % seq must be a row vector. -% k = sprintf(obj.key_fmt, seq); -% end -% end -% end - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - classdef ML_MLSE < handle % --------------------------------------------------------------------- % W. Lanneer and Y. Lefevre, @@ -661,8 +165,8 @@ classdef ML_MLSE < handle % EQUALIZE % ============================================================== function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training) - debug = 1; - showPlots = 1; + debug = 0; + showPlots = 0; y = zeros(N,1); nSymbols = ceil(N/obj.sps); diff --git a/Classes/04_DSP/Equalizer/VNLE.m b/Classes/04_DSP/Equalizer/VNLE.m index d95bb98..51d96a5 100644 --- a/Classes/04_DSP/Equalizer/VNLE.m +++ b/Classes/04_DSP/Equalizer/VNLE.m @@ -10,6 +10,7 @@ classdef VNLE < handle sps % usually 2 order e + e_dc error len_tr @@ -18,10 +19,17 @@ classdef VNLE < handle mu_dd epochs_dd + mu_dc constellation decide + save_debug = 0; + debug_struct + + optmize_mus = 0; + mu_optimization + mu_optimization_iter = 0; x_norm ce @@ -42,8 +50,11 @@ classdef VNLE < handle options.mu_dd = 1e-5; options.epochs_dd = 5; + options.mu_dc = 0; options.decide = false; + options.save_debug = 0; + options.optmize_mus = 0; end @@ -54,6 +65,7 @@ classdef VNLE < handle obj.error = 0; + obj.e_dc = 0; end @@ -69,6 +81,13 @@ classdef VNLE < handle [obj.ie2,obj.ie3] = obj.calcIndiceVectors(obj.order); obj.e = zeros( sum(obj.ce) ,1); + obj.e_dc = 0; + + if obj.optmize_mus + obj.optimizeMus(X.signal,D.signal); + obj.e = zeros(sum(obj.ce),1); + obj.e_dc = 0; + end % Training Mode training = 1; @@ -148,7 +167,7 @@ classdef VNLE < handle x_in = x(obj.order(1)+sample-1:-1:sample); x_in = obj.calcVNLENonlinVecs(x_in,obj.ie2,obj.ie3,obj.order,obj.x_norm); - y(symbol,1) = obj.e.' * x_in; % Calculating output of LMS __ * | + y(symbol,1) = obj.e_dc + obj.e.' * x_in; % Calculating output of LMS __ * | if training err = y(symbol) - d(symbol); % Instantaneous error @@ -164,6 +183,9 @@ classdef VNLE < handle normalizationfactor = (x_in.' * x_in); obj.e = obj.e - err * x_in / normalizationfactor; % Weight update rule of NLMS end + if obj.mu_dc ~= 0 + obj.e_dc = obj.e_dc - obj.mu_dc * err; + end if mod(sample,100) == 1 && showviz a2.XData = 1:2*numel(y); @@ -177,12 +199,85 @@ classdef VNLE < handle end obj.error(epoch,symbol) = err * err'; % Instantaneous square error + if obj.save_debug + obj.debug_struct.error(epoch,symbol) = err * err'; + if training + obj.debug_struct.error_tr(epoch,symbol) = err * err'; + end + end end end end + function optimizeMus(obj,x,d) + mu_range = [1e-5, 1e-2]; + mu_dc_range = [1e-5, 1e-1]; + + vars = [optimizableVariable("mu_tr",mu_range,"Transform","log"), ... + optimizableVariable("mu_dd",mu_range,"Transform","log")]; + optimize_mu_dc = obj.mu_dc ~= 0; + if optimize_mu_dc + vars = [vars, optimizableVariable("mu_dc",mu_dc_range,"Transform","log")]; + end + + obj.mu_optimization_iter = 0; + obj.mu_optimization = bayesopt(@(p)obj.muObjective(p,x,d),vars, ... + "MaxObjectiveEvaluations",10, ... + "AcquisitionFunctionName","expected-improvement-plus", ... + "IsObjectiveDeterministic",false, ... + "Verbose",0, ... + "PlotFcn",[]); + obj.mu_tr = obj.mu_optimization.XAtMinObjective.mu_tr; + obj.mu_dd = obj.mu_optimization.XAtMinObjective.mu_dd; + if optimize_mu_dc + obj.mu_dc = obj.mu_optimization.XAtMinObjective.mu_dc; + end + + objective_db = 10*log10(obj.mu_optimization.MinObjective); + if optimize_mu_dc + fprintf("\nVNLE mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... + obj.mu_tr,obj.mu_dd,obj.mu_dc,obj.mu_optimization.MinObjective,objective_db); + else + fprintf("\nVNLE mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ... + obj.mu_tr,obj.mu_dd,obj.mu_optimization.MinObjective,objective_db); + end + end + + function objective = muObjective(obj,params,x,d) + old_debug = obj.save_debug; + old_mu_dc = obj.mu_dc; + obj.save_debug = 1; + optimize_mu_dc = ismember("mu_dc",string(params.Properties.VariableNames)); + if optimize_mu_dc + obj.mu_dc = params.mu_dc; + end + + obj.e = zeros(sum(obj.ce),1); + obj.e_dc = 0; + obj.debug_struct = struct(); + obj.equalize(x,d,params.mu_tr,obj.epochs_tr,obj.len_tr,1,0); + obj.equalize(x,d,params.mu_dd,obj.epochs_dd,numel(x),0,0); + + objective = mean(obj.debug_struct.error(end,:),"omitnan"); + if ~isfinite(objective) + objective = inf; + end + + objective_db = 10*log10(objective); + obj.mu_optimization_iter = obj.mu_optimization_iter + 1; + if optimize_mu_dc + fprintf("\rVNLE mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,params.mu_dc,objective,objective_db); + else + fprintf("\rVNLE mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ... + obj.mu_optimization_iter,params.mu_tr,params.mu_dd,objective,objective_db); + end + obj.save_debug = old_debug; + obj.mu_dc = old_mu_dc; + end + %% Functions needed During Adaption function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_) % These are the second and third order input signal products of the VNLE EQ diff --git a/Datatypes/channel_model.m b/Datatypes/channel_model.m new file mode 100644 index 0000000..0109492 --- /dev/null +++ b/Datatypes/channel_model.m @@ -0,0 +1,9 @@ +classdef channel_model < int32 + + enumeration + awgn (1) + awgn_alphad (2) + physical (3) + end + +end \ No newline at end of file diff --git a/Functions/EQ_structures/dsp_scope_signal.m b/Functions/EQ_structures/dsp_scope_signal.m index 0a6520a..31f0bca 100644 --- a/Functions/EQ_structures/dsp_scope_signal.m +++ b/Functions/EQ_structures/dsp_scope_signal.m @@ -47,12 +47,12 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options) use_dd_mode = 1; use_ffe = 1; - use_dfe = 0; - use_vnle_mlse = 0; - use_dbtgt = 0; + use_dfe = 1; + use_vnle_mlse = 1; + use_dbtgt = 1; use_dbenc = 0; use_ml_mlse = 0; - showAnalysis = 1; + showAnalysis = 0; decoding_mode = []; addProcessingResultToDatabase = 0; %#ok @@ -70,8 +70,8 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options) pf_ = Postfilter("ncoeff", pf_ncoeffs, "useBurg", 1); %#ok mlse_ = MLSE("duobinary_output", 0, 'M', M, 'trellis_states', PAMmapper(M,0).levels); %#ok - 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"); %#ok - eq_post = FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption",adaption_method(adaption),"dd_mode",use_dd_mode); %#ok + 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","mu_dc",mu_dc); %#ok + eq_post = FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",mu_dd,"mu_tr",mu_tr,"order",25,"sps",2,"decide",0, "adaption_technique",adaption_method(adaption),"dd_mode",use_dd_mode,"mu_dc",mu_dc); %#ok mlse_db_enc = MLSE("DIR", [1,1], "duobinary_output", 0, "M", M, "trellis_states", PAMmapper(M,0).levels); %#ok eq_db_enc = EQ("Ne", ffe_order_dbtgt, "Nb", dfe_order_dbtgt, "training_length", len_tr, ... @@ -81,21 +81,22 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options) Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, fsym, ... "mode", options.preprocess_mode, ... "tx_pulseformer", options.tx_pulseformer, ... - "debug_plots", options.debug_plots); - Scpe_sig.spectrum("fignum",200,"normalizeTo0dB",1,"displayname",'Rx','addDCoffset',-6); - Scpe_sig.spectrum("fignum",201,"normalizeTo0dB",0,"displayname",'Rx'); + "debug_plots", 0); - ylim([-30,3]); - xlim([-5,100]); + % Scpe_sig.spectrum("fignum",200,"normalizeTo0dB",1,"displayname",'Rx','addDCoffset',-6); + % ylim([-30,3]); + % xlim([-5,100]); Scpe_sig.signal = Scpe_sig.signal(1:2*Symbols.length); Scpe_sig.signal = real(Scpe_sig.signal); if duob_mode ~= db_mode.db_encoded if use_ffe - eq_dfe = EQ("Ne",ffe_order_ffe,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + eq_ffe = EQ("Ne",ffe_order_ffe,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + eq_ffe = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-1,"mu_tr",0.4,"order",ffe_order_ffe(1),... + "sps",2,"decide",0,"optmize_mus",0,"dd_mode",1,"adaption_technique","nlms","mu_dc",1.021e-05); - ffe_results = ffe(eq_dfe, M, Scpe_sig, Symbols, Tx_bits, ... + ffe_results = ffe(eq_ffe, M, Scpe_sig, Symbols, Tx_bits, ... "precode_mode", duob_mode, ... 'showAnalysis', options.debug_plots, ... "postFFE", [], ... @@ -116,7 +117,7 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options) "eth_style_symbol_mapping", 0); dfe_results.config.equalizer_structure = "dfe"; - dfe_results.metrics.print; + dfe_results.metrics.print("description",'DFE'); output.dfe_package = dfe_results; end @@ -144,8 +145,8 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options) "eth_style_symbol_mapping", 0); ffe_results.config.equalizer_structure = "vnle"; - ffe_results.metrics.print; - mlse_results.metrics.print; + ffe_results.metrics.print("description",'VNLE'); + mlse_results.metrics.print("description",'MLSE'); output.mlse_package = mlse_results; output.vnle_package = ffe_results; @@ -190,7 +191,7 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options) "decoding_mode", decoding_mode); end - dbt_results.metrics.print; + dbt_results.metrics.print("description",'Duob. Target'); output.dbtgt_package = dbt_results; end end diff --git a/Functions/EQ_structures/duobinary_target.m b/Functions/EQ_structures/duobinary_target.m index c24ffca..2e42e75 100644 --- a/Functions/EQ_structures/duobinary_target.m +++ b/Functions/EQ_structures/duobinary_target.m @@ -14,6 +14,7 @@ function [db_results] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx options.decoding_mode db_decoder = db_decoder.sequencedetection; end + %Duobinary Targeting db_ref_sequence = Duobinary().encode(tx_symbols); @@ -172,14 +173,15 @@ function [db_results] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx rx_signal.spectrum("normalizeTo0dB",1,"fignum",250,"displayname","Rx Spectrum"); - Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",10,"displayname","DB encoded reference"); + Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",250,"displayname","DB encoded reference"); showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise after Equalization'); fprintf('DB tgt BER: %.2e \n',ber_db); figure(341); clf; - showLevelHistogram(eq_signal, db_ref_sequence, "fignum", 341); + tx_symbols_uncoded = Duobinary().decode(db_ref_sequence); + showLevelHistogram(eq_signal, db_ref_sequence, "fignum", 341,"ref_symbol_uncoded",tx_symbols_uncoded); end diff --git a/Functions/EQ_structures/ffe.m b/Functions/EQ_structures/ffe.m index a58e215..964905e 100644 --- a/Functions/EQ_structures/ffe.m +++ b/Functions/EQ_structures/ffe.m @@ -171,6 +171,8 @@ function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, showEQNoisePSD(eq_noise, "fignum", 103, "displayname", 'Residual Noise after FFE'); + % showEQcoefficients('n1', eq_.e, "displayname", 'Coefficients', 'fignum', 104); + % Figure 2: Post-FFE coefficients (if available) if ~isempty(postFFE) showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 104); @@ -179,7 +181,6 @@ function displayAnalysis(eq_noise, eq_signal_sd, rx_signal, eq_, tx_symbols, M, try figure(339);hold on showEQfilter(eq_.e, 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 diff --git a/Functions/EQ_structures/vnle_postfilter_mlse.m b/Functions/EQ_structures/vnle_postfilter_mlse.m index 803ec70..7a0eaac 100644 --- a/Functions/EQ_structures/vnle_postfilter_mlse.m +++ b/Functions/EQ_structures/vnle_postfilter_mlse.m @@ -48,6 +48,8 @@ eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal_sd); [mlse_sig_sd,whitened_noise] = pf_.process(eq_signal_sd, eq_noise); +[tx_symbols_pr,~] = pf_.process(tx_symbols, eq_noise); + if 0 %tx_symbols.fs > 190e9 if pf_.ncoeff == 1 if pf_.coefficients(2) < 0 @@ -283,7 +285,7 @@ if ~isempty(postFFE) showEQcoefficients('n1', postFFE.e, "displayname", 'Coefficients', 'fignum', 341); end -% showEQcoefficients('n1', eq_.e1,'n2', eq_.e2,'n3', eq_.e3, "displayname", 'Coefficients', 'fignum', 339); +showEQcoefficients('n1', eq_.e,'n2', eq_.e2,'n3', eq_.e3, "displayname", 'Coefficients', 'fignum', 339); showEQfilter(eq_.e, eq_signal_sd.fs.*2); figure(340); clf; diff --git a/Functions/EQ_visuals/show2Dconstellation.m b/Functions/EQ_visuals/show2Dconstellation.m index 1170d66..14ee6df 100644 --- a/Functions/EQ_visuals/show2Dconstellation.m +++ b/Functions/EQ_visuals/show2Dconstellation.m @@ -4,6 +4,12 @@ arguments ref_symbols options.fignum (1,1) double = NaN % Default to NaN if not provided options.displayname (1,:) char = '' % Default to an empty string if not provided + options.color = []; + options.clear = true; + options.rasterize = false; + options.raster_resolution (1,1) double = 1000; + options.raster_markersize (1,1) double = 1; + options.raster_clim = []; end plot_shit = 1; @@ -22,7 +28,9 @@ if plot_shit fig = figure; % Create a new figure and get its handle else fig = figure(options.fignum); % Use the specified figure number - clf; + if options.clear + clf; + end end end @@ -31,6 +39,8 @@ rx_symbols = eq_signal; %./ rms(eq_signal); correct_symbols = ref_symbols; col = cbrewer2('Paired',numel(unique(correct_symbols))*3); +useDefaultColormap = isempty(options.color); +scatterColor = options.color; ccnt = -1; levels = unique(correct_symbols); @@ -66,20 +76,146 @@ D = sqrt(D_even.^2 + D_odd.^2); [X,Y] = meshgrid(levels, levels); [X_,Y_] = meshgrid(statistical_mean, statistical_mean); +plot_xlim = [floor(min(X_(:)))-1, ceil(max(X_(:)))+1]; +plot_ylim = [floor(min(Y_(:)))-1, ceil(max(Y_(:)))+1]; hold on; -scatter(rx_even,rx_odd,5*ones(1,length(D)),D,'.','DisplayName',['Even/Odd PAM-',num2str(M)],'MarkerEdgeColor',col(2,:)); -% colormap(gca,flip(cbrewer2('Spectral',100))) -colormap(gca,'hsv'); -scatter(X_(:), Y_(:), 2, 'x', 'LineWidth', 10, 'MarkerEdgeColor', col(6,:),'DisplayName','Statistical Rx Levels'); -scatter(X(:), Y(:), 2, 'x', 'LineWidth', 10, 'MarkerEdgeColor', col(4,:),'DisplayName','Tx Levels'); -xlim([floor(min(X_(:)))-1, ceil(max(X_(:)))+1]); -ylim([floor(min(Y_(:)))-1, ceil(max(Y_(:)))+1]); +if isempty(options.displayname) + displayName = ['Even/Odd PAM-',num2str(M)]; +else + displayName = options.displayname; +end +if useDefaultColormap + cols = flip(cbrewer2('RdYlBu',100)); + cols=(cbrewer2("RdYlBu",4096)); + % cols(45:55,:) = []; + colormap(gca,cols); + if isempty(options.raster_clim) + cLim = [0 max(D).*0.45]; + else + cLim = options.raster_clim; + end + if cLim(2) <= cLim(1) + cLim(2) = cLim(1) + eps; + end + clim(cLim); + % colormap(gca,'hsv'); + if options.rasterize + nBins = options.raster_resolution; + xEdges = linspace(plot_xlim(1),plot_xlim(2),nBins+1); + yEdges = linspace(plot_ylim(1),plot_ylim(2),nBins+1); + xBin = discretize(rx_even,xEdges); + yBin = discretize(rx_odd,yEdges); + valid = ~isnan(xBin) & ~isnan(yBin) & isfinite(D); -yticks((levels(1:end-1) + levels(2:end)) / 2); -xticks((levels(1:end-1) + levels(2:end)) / 2); + countImg = accumarray([yBin(valid), xBin(valid)],1,[nBins nBins],@sum,0); + dImg = accumarray([yBin(valid), xBin(valid)],D(valid),[nBins nBins],@mean,NaN); + countImgRaw = countImg; + markerRadius = max(0, round((options.raster_markersize - 1) / 2)); + if markerRadius > 0 + kernel = ones(2*markerRadius+1); + countImg = conv2(countImg,kernel,'same'); + weightedD = dImg; + weightedD(~isfinite(weightedD)) = 0; + weightedD = conv2(weightedD .* countImgRaw,kernel,'same'); + dWeight = conv2(isfinite(dImg) .* countImgRaw,kernel,'same'); + dImg = weightedD ./ dWeight; + end + colorIdx = round((dImg-cLim(1)) ./ diff(cLim) * (size(cols,1)-1)) + 1; + colorIdx = min(max(colorIdx,1),size(cols,1)); + + density = log1p(countImg); + if max(density(:)) > 0 + density = density ./ max(density(:)); + end + + rgbImg = ones(nBins,nBins,3); + for rgbIdx = 1:3 + colorPlane = reshape(cols(colorIdx,rgbIdx),nBins,nBins); + colorPlane(~isfinite(dImg)) = 1; + rgbImg(:,:,rgbIdx) = (1-density) + density .* colorPlane; + end + + hImg = image(linspace(plot_xlim(1),plot_xlim(2),nBins), ... + linspace(plot_ylim(1),plot_ylim(2),nBins), ... + rgbImg); + hImg.HandleVisibility = 'off'; + set(gca,'YDir','normal'); + uistack(hImg,'bottom'); + plot(nan,nan,'.','DisplayName',displayName,'MarkerEdgeColor',col(2,:)); + else + sz = ones(1,length(D)).*1; + scatter(rx_even,rx_odd,sz,D,'.','DisplayName',displayName); + end + rxLevelColor = col(6,:); +else + if options.rasterize + nBins = options.raster_resolution; + xEdges = linspace(plot_xlim(1),plot_xlim(2),nBins+1); + yEdges = linspace(plot_ylim(1),plot_ylim(2),nBins+1); + xBin = discretize(rx_even,xEdges); + yBin = discretize(rx_odd,yEdges); + valid = ~isnan(xBin) & ~isnan(yBin); + + countImg = accumarray([yBin(valid), xBin(valid)],1,[nBins nBins],@sum,0); + markerRadius = max(0, round((options.raster_markersize - 1) / 2)); + if markerRadius > 0 + kernel = ones(2*markerRadius+1); + countImg = conv2(countImg,kernel,'same'); + end + density = log1p(countImg); + if max(density(:)) > 0 + density = density ./ max(density(:)); + end + + rgbImg = ones(nBins,nBins,3); + for rgbIdx = 1:3 + rgbImg(:,:,rgbIdx) = (1-density) + density .* scatterColor(rgbIdx); + end + + hImg = image(linspace(plot_xlim(1),plot_xlim(2),nBins), ... + linspace(plot_ylim(1),plot_ylim(2),nBins), ... + rgbImg); + hImg.HandleVisibility = 'off'; + set(gca,'YDir','normal'); + uistack(hImg,'bottom'); + plot(nan,nan,'.','DisplayName',displayName,'MarkerEdgeColor',scatterColor); + else + scatter(rx_even,rx_odd,5*ones(1,length(D)),'.','DisplayName',displayName,'MarkerEdgeColor',scatterColor); + end + rxLevelColor = scatterColor; +end + +scatter(X_(:), Y_(:), 2, 'x', 'LineWidth', 10, 'MarkerEdgeColor', rxLevelColor,'DisplayName',[displayName, ' Rx Levels']); +scatter(X(:), Y(:), 2, 'x', 'LineWidth', 10, 'MarkerEdgeColor', col(4,:),'DisplayName','Tx Levels'); +xlim(plot_xlim); +ylim(plot_ylim); + +decisionBoundaries = (levels(1:end-1) + levels(2:end)) / 2; +ax = gca; +xticks(levels); +yticks(levels); +xticklabels(compose('%g',levels)); +yticklabels(compose('%g',levels)); +try + ax.XAxis.MinorTickValues = decisionBoundaries; + ax.YAxis.MinorTickValues = decisionBoundaries; +end +ax.XMinorTick = 'on'; +ax.YMinorTick = 'on'; +ax.XMinorGrid = 'on'; +ax.YMinorGrid = 'on'; +ax.GridLineStyle = '--'; +ax.MinorGridLineStyle = ':'; +ax.GridColor = [0.65 0.65 0.65]; +ax.MinorGridColor = [0.35 0.35 0.35]; +ax.GridAlpha = 0.35; +ax.MinorGridAlpha = 0.45; legend xlabel('even symbols'); ylabel('odd symbols'); -axis equal; grid on; +axis equal; +grid on; +grid minor; +ax.Layer = 'top'; end diff --git a/Functions/EQ_visuals/showEQNoisePSD.m b/Functions/EQ_visuals/showEQNoisePSD.m index ed7b440..84f4d1e 100644 --- a/Functions/EQ_visuals/showEQNoisePSD.m +++ b/Functions/EQ_visuals/showEQNoisePSD.m @@ -4,7 +4,7 @@ arguments options.postfilter_taps = NaN options.fignum (1,1) double = NaN % Default to NaN if not provided options.displayname (1,:) char = '' % Default to an empty string if not provided - options.color = [0.2157 0.4941 0.7216]; + options.color = []; end % Determine the figure number to use or create a new figure @@ -16,12 +16,14 @@ end hold on ax = gca; % N = numel(ax.Children); - N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children)); - cmap = linspecer(8); - options.color = cmap(mod(N, size(cmap, 1)) + 1, :); + if isempty(options.color) + N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children)); + cmap = linspecer(8); + options.color = cmap(mod(N, size(cmap, 1)) + 1, :); + end % Ensure the figure is ready before calling spectrum - eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 1,"color",options.color); + eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 0,"color",options.color); if ~isnan(options.postfilter_taps) % Hold on to the figure for further plotting diff --git a/Functions/EQ_visuals/showEQNoiseSNR.m b/Functions/EQ_visuals/showEQNoiseSNR.m index cac55eb..f0715d3 100644 --- a/Functions/EQ_visuals/showEQNoiseSNR.m +++ b/Functions/EQ_visuals/showEQNoiseSNR.m @@ -6,7 +6,8 @@ arguments options.fs_rx options.fignum (1,1) double = NaN % Default to NaN if not provided options.displayname (1,:) char = '' % Default to an empty string if not provided - options.color = [0.2157 0.4941 0.7216]; + options.color = []; + options.linestyle = '-'; end % Determine the figure number to use or create a new figure @@ -30,42 +31,52 @@ end ax = gca; % N = numel(ax.Children); - N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children)); - cmap = linspecer(8); - options.color = cmap(mod(N, size(cmap, 1)) + 1, :); - - % Ensure the figure is ready before calling spectrum - title('SNR of received Signal') + if isempty(options.color) + N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children)); + cmap = linspecer(8); + options.color = cmap(mod(N, size(cmap, 1)) + 1, :); + end fft_length = 2^(nextpow2(length(eq_signal))-7); [s_lin,w] = pwelch(eq_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_tx,"centered","psd","mean"); - [n_lin,w] = pwelch(noise_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_rx,"centered","psd","mean"); + [n_lin,w_noise] = pwelch(noise_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_rx,"centered","psd","mean"); + + if numel(w_noise) ~= numel(w) || any(abs(w_noise - w) > max(options.fs_tx,options.fs_rx)*eps) + n_lin = interp1(w_noise,n_lin,w,"linear",NaN); + end w = w.*1e-9; - snr_dbm = 10*log10(s_lin./n_lin); + snr_lin = s_lin./n_lin; + snr_lin = movmean(snr_lin,10); + snr_dbm = 10*log10(snr_lin); - % figure(231) - hold on - plot(w,snr_dbm,'DisplayName','SNR','LineWidth',0.5,'Color',options.color); - yline(mean(snr_dbm),'HandleVisibility','off','Color',options.color); + if isempty(options.displayname) + options.displayname = 'SNR'; + end + + plot(w,snr_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color,'LineStyle',options.linestyle); + % yline(mean(snr_dbm),'HandleVisibility','off','Color',options.color,'LineStyle','--','LineWidth',1); xlabel("Frequency in GHz"); - - edgetick = 2^(nextpow2(options.fs_tx*1e-9)); - ticks = -edgetick:16:edgetick; - xticks(ticks); - - [~,b]=min(abs((-edgetick:16:edgetick)-max(w))); - xlim([-ticks(b+1) ticks(b+1)]); + ylabel("SNR (dB)"); + xlim([min(w) max(w)]); - max_snr = ceil(max(snr_dbm)/10)*10; - min_snr = floor(min(snr_dbm)/10)*10; - ylim([min_snr,max_snr]); - yticks(-200:10:100); + y_min = min(snr_dbm(:)); + y_max = max(snr_dbm(:)); + y_range = y_max - y_min; + if y_range == 0 + y_range = 10; + end + y_margin = 0.05 * y_range; + ylim([y_min - y_margin, y_max + y_margin]); + try + yticks(round(linspace(y_min, y_max, min(10, max(4, ceil(y_range/10)))))); + end - grid on; grid minor; - legend('Interpreter','none'); - title('Noise of soft decision signal (not MLSE)'); + grid on; + if isempty(get(gca, 'Legend')) + legend('Interpreter','none'); + end end diff --git a/Functions/EQ_visuals/showEQfilter.m b/Functions/EQ_visuals/showEQfilter.m index 4b0f403..d0a9f8f 100644 --- a/Functions/EQ_visuals/showEQfilter.m +++ b/Functions/EQ_visuals/showEQfilter.m @@ -6,6 +6,7 @@ arguments 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 + options.color = []; end % Assuming that obj.e contains the final FFE filter coefficients. @@ -34,20 +35,24 @@ end end % Magnitude response (in dB) - subplot(2,1,1); + % subplot(2,1,1); hold on - plot(f.*1e-9, H_db,'DisplayName',options.displayname); + if isempty(options.color) + plot(f.*1e-9, H_db,'DisplayName',options.displayname,'LineWidth',1); + else + plot(f.*1e-9, H_db,'DisplayName',options.displayname,'Color',options.color,'LineWidth',1); + end title('(Inverted) Magnitude Response of FFE Filter'); xlabel('Frequency (GHz)'); ylabel('Magnitude (dB)'); grid on; - % Phase response - subplot(2,1,2); - plot(f.*1e-9, unwrap(angle(H)),'DisplayName',options.displayname); - title('Phase Response of FFE Filter'); - xlabel('Frequency (GHz)'); - ylabel('Phase'); - grid on; + % % Phase response + % subplot(2,1,2); + % plot(f.*1e-9, unwrap(angle(H)),'DisplayName',options.displayname); + % title('Phase Response of FFE Filter'); + % xlabel('Frequency (GHz)'); + % ylabel('Phase'); + % grid on; end \ No newline at end of file diff --git a/Functions/EQ_visuals/showLevelHistogram.m b/Functions/EQ_visuals/showLevelHistogram.m index 5d96767..6ef093d 100644 --- a/Functions/EQ_visuals/showLevelHistogram.m +++ b/Functions/EQ_visuals/showLevelHistogram.m @@ -4,6 +4,8 @@ arguments ref_symbols options.fignum (1,1) double = NaN % Default to NaN if not provided options.displayname (1,:) char = '' % Default to an empty string if not provided + options.ref_symbol_uncoded = [] + options.nbins (1,1) double = 1000 end if isa(eq_signal,'Signal') @@ -12,6 +14,20 @@ end if isa(ref_symbols,'Signal') ref_symbols = ref_symbols.signal; end + if isa(options.ref_symbol_uncoded,'Signal') + options.ref_symbol_uncoded = options.ref_symbol_uncoded.signal; + end + + eq_signal = eq_signal(:); + ref_symbols = ref_symbols(:); + ref_symbol_uncoded = options.ref_symbol_uncoded(:); + + assert(numel(eq_signal) == numel(ref_symbols), ... + 'showLevelHistogram:LengthMismatch', ... + 'eq_signal and ref_symbols must have the same number of samples.'); + assert(isempty(ref_symbol_uncoded) || numel(ref_symbol_uncoded) == numel(eq_signal), ... + 'showLevelHistogram:LengthMismatch', ... + 'options.ref_symbol_uncoded must have the same number of samples as eq_signal.'); % Determine the figure number to use or create a new figure if isnan(options.fignum) @@ -22,31 +38,70 @@ 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('Paired',numel(constellation)*2); - lvlcol = lvlcol(2:2:end,:); - lvlcol = linspecer(numel(constellation)); - % lvlcol = cbrewer2('Set1',numel(constellation)); - for lvl = 1:numel(constellation) - %Separate the equalized signal into the - %respective levels based on the actually - %transmitted level! - received_sd(lvl,ref_symbols==constellation(lvl)) = eq_signal(ref_symbols==constellation(lvl)); - end - - % scaling = PAMmapper(numel(constellation),0).get_scaling; - - %%% FFE histogram + %%% histogram clf - for lvl = 1:numel(constellation) - intermediate = received_sd(lvl,:); - cnt(lvl) = round(numel(intermediate(~isnan(intermediate)))./length(eq_signal),3).*100; - hold on - warning off - histogram(received_sd(lvl,:),1000,"EdgeAlpha",0,'DisplayName',['Lvl ',num2str(lvl),' ; ',num2str(cnt(lvl)),' '],'FaceColor',lvlcol(lvl,:),'Normalization','pdf'); - warning on + if isempty(ref_symbol_uncoded) + % Normal mode: split the received signal by the actually + % transmitted reference level. + constellation = unique(ref_symbols); + received_sd = NaN(numel(constellation),numel(ref_symbols)); + lvlcol = linspecer(numel(constellation)); + + for lvl = 1:numel(constellation) + class_mask = ref_symbols == constellation(lvl); + received_sd(lvl,class_mask) = eq_signal(class_mask); + end + + for lvl = 1:numel(constellation) + intermediate = received_sd(lvl,:); + cnt(lvl) = round(numel(intermediate(~isnan(intermediate)))./numel(eq_signal),3).*100; + hold on + warning off + histogram(received_sd(lvl,:),options.nbins, ... + "EdgeAlpha",0, ... + "DisplayName",['Lvl ',num2str(lvl),' ; ',num2str(cnt(lvl)),' '], ... + "FaceColor",lvlcol(lvl,:), ... + "Normalization","pdf"); + warning on + end + else + % Plot p(y | x_n): y is the noisy observation and x_n is the + % uncoded PAM class. For duobinary this naturally gives multi-modal + % PDFs because one x_n class can map to multiple DB amplitudes. Each + % DB lobe is drawn separately so the DB-level structure stays visible. + db_constellation = unique(ref_symbols); + classes = unique(ref_symbol_uncoded); + lvlcol = linspecer(numel(classes)); + + for db_lvl = 1:numel(db_constellation) + db_mask = ref_symbols == db_constellation(db_lvl); + mapped_class = mode(ref_symbol_uncoded(db_mask)); + [~,class_idx] = min(abs(classes - mapped_class)); + cnt = round(nnz(ref_symbol_uncoded == mapped_class)./numel(eq_signal),3).*100; + + if db_lvl == findFirstMappedDbLevel(ref_symbols,ref_symbol_uncoded,db_constellation,mapped_class) + display_name = ['p(y|x_',num2str(class_idx),') ; ',num2str(cnt),' ']; + display_name = ['p(y|x_',num2str(class_idx),')']; + handle_visibility = "on"; + else + display_name = ['p(y|x_',num2str(class_idx),') ; ',num2str(cnt),' ']; + display_name = ['p(y|x_',num2str(class_idx),')']; + handle_visibility = "off"; + end + + weighted_lobe = NaN(size(eq_signal)); + weighted_lobe(db_mask) = eq_signal(db_mask); + + hold on + warning off + histogram(weighted_lobe,options.nbins, ... + "EdgeAlpha",0, ... + "DisplayName",display_name, ... + "FaceColor",lvlcol(class_idx,:), ... + "HandleVisibility",handle_visibility, ... + "Normalization","pdf"); + warning on + end end xlim([-3 3]); legend @@ -58,3 +113,15 @@ end end +function first_db_lvl = findFirstMappedDbLevel(ref_symbols,ref_symbol_uncoded,db_constellation,mapped_class) + first_db_lvl = NaN; + + for db_lvl = 1:numel(db_constellation) + db_mask = ref_symbols == db_constellation(db_lvl); + if mode(ref_symbol_uncoded(db_mask)) == mapped_class + first_db_lvl = db_lvl; + return + end + end +end + diff --git a/Functions/EQ_visuals/showMarkovDiagram.m b/Functions/EQ_visuals/showMarkovDiagram.m new file mode 100644 index 0000000..49068e3 --- /dev/null +++ b/Functions/EQ_visuals/showMarkovDiagram.m @@ -0,0 +1,35 @@ +function showMarkovDiagram(sequence,M) + + +if isa(sequence,'Signal') + sequence = sequence.signal; +end +sequence = sequence(2+mod(1,length(sequence)):end); %filtered sequences often have one "old" sample at idx=1 + +x = sequence; + +levels = sort(unique(x)).'; % or provide known 1x6 level values + +[~,ix] = min(abs(x - levels),[],2); +x = levels(ix); + +K = numel(levels); +% map to state indices 1..K +[tf, idx] = ismember(x, levels); +idx = idx(:); +from = idx(1:end-1); +to = idx(2:end); +from = idx(1:2:end); +to = idx(2:2:end); + +% counts C(from,to) +C = accumarray([from,to], 1, [K K], @sum, 0); +% row-stochastic transition matrix P(to|from) +rowSums = sum(C,2); +P = C ./ max(rowSums,1); + +mc = dtmc(P, 'StateNames', string(levels.*PAMmapper(M,0).get_scaling)); +figure('Name','Markov Graph (dtmc)'); +gp = graphplot(mc, 'ColorEdges',true, 'LabelEdges',true); + +end \ No newline at end of file diff --git a/Functions/EQ_visuals/showTransitionProbabilities.m b/Functions/EQ_visuals/showTransitionProbabilities.m new file mode 100644 index 0000000..db07af7 --- /dev/null +++ b/Functions/EQ_visuals/showTransitionProbabilities.m @@ -0,0 +1,40 @@ +function showTransitionProbabilities(sequence) + + +if isa(sequence,'Signal') + sequence = sequence.signal; +end + +sequence = sequence(2+mod(1,length(sequence)):end); %filtered sequences often have one "old" sample at idx=1 + +x = sequence; +levels = sort(unique(x)).'; % or provide known 1x6 level values + +[~,ix] = min(abs(x - levels),[],2); +x = levels(ix); + + +K = numel(levels); +% map to state indices 1..K +[tf, idx] = ismember(x, levels); +idx = idx(:); +% from = idx(1:end-1); +% to = idx(2:end); +from = idx(1:2:end); +to = idx(2:2:end); + +% counts C(from,to) +C = accumarray([from,to], 1, [K K], @sum, 0); +% row-stochastic transition matrix P(to|from) +rowSums = sum(C,2); +P = C ./ max(rowSums,1); + +%% 1) HEATMAP (which transitions are more probable?) +figure('Name','Transition Probabilities (to | from)'); +h = heatmap(levels, levels, P, 'Colormap', parula, 'ColorbarVisible','on'); +colormap(gca,[[1,1,1];flip(cbrewer2('Spectral',100))]);clim([0,ceil(max(P(:))*10)/10]); +h.XLabel = 'From state (level)'; +h.YLabel = 'To state (level)'; +h.Title = 'P(to | from)'; + +end \ No newline at end of file diff --git a/Functions/Lab_helper/loadFreqResp.m b/Functions/Lab_helper/loadFreqResp.m index e2c99e8..c6e4da8 100644 --- a/Functions/Lab_helper/loadFreqResp.m +++ b/Functions/Lab_helper/loadFreqResp.m @@ -1,6 +1,6 @@ % Define the precomp path -precomp_path = "C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\precomp"; - +precomp_path = "W:\labdata\sioe_labor\precomp"; +precomp_path = "C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\ECOC_2025_MPI"; % Step 1: Find all valid files (assume .mat files for ChannelFreqResp) fileList = dir(fullfile(precomp_path, '*.mat')); fileNames = {fileList.name}; diff --git a/Functions/mat2tikz_improved.m b/Functions/mat2tikz_improved.m index 552e4af..345f057 100644 --- a/Functions/mat2tikz_improved.m +++ b/Functions/mat2tikz_improved.m @@ -1,9 +1,12 @@ -function mat2tikz_improved(filename) +function mat2tikz_improved(filename,options) arguments % Default to the path in your example if no argument is provided filename (1,1) string = 'C:\Users\Silas\Documents\Dissertation\00_Examples\tikz\textfig.tikz'; + options.cleanfigure = false; +end +if options.cleanfigure + cleanfigure; end -cleanfigure; matlab2tikz(char(filename), ... 'width', '\fwidth', ... 'height', '\fheight', ... diff --git a/Theory/Mapping_Coding/pam_6_states_analysis.m b/Theory/Mapping_Coding/pam_6_states_analysis.m index a92d3ca..7cd09c0 100644 --- a/Theory/Mapping_Coding/pam_6_states_analysis.m +++ b/Theory/Mapping_Coding/pam_6_states_analysis.m @@ -1,20 +1,8 @@ -M = 6; +M = 4; -bitpattern = []; -s = RandStream('twister','Seed',1); -for i = 1:log2(M) - N = 2^(17-1); %length of prbs - bitpattern(:,i) = randi(s,[0 1], N, 1); -end - -if M == 6 - bitpattern = reshape(bitpattern',[],1); - bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); -end - -bits = Informationsignal(bitpattern); +bits = Signalgenerator("form", signalform.prms,"M", M,"order", 16).process(); symbols = PAMmapper(M,0).map(bits); @@ -106,7 +94,7 @@ end %% State Analysis -signal_to_analyze = symbols_tx_emu; +signal_to_analyze = symbols_db; x = signal_to_analyze.signal(:); levels = sort(unique(x)).'; % or provide known 1x6 level values diff --git a/Theory/Mapping_Coding/partial_response/duobinary_back2back.m b/Theory/Mapping_Coding/partial_response/duobinary_back2back.m index 9e5badd..3fccde2 100644 --- a/Theory/Mapping_Coding/partial_response/duobinary_back2back.m +++ b/Theory/Mapping_Coding/partial_response/duobinary_back2back.m @@ -1,5 +1,5 @@ -for M = [2 4 8] - bits = Signalgenerator("form", signalform.prms,"M", M,"order", 16).process(); +for M = [4] + bits = Signalgenerator("form", signalform.prms,"M", M,"order", 18).process(); mapper = PAMmapper(M,0); symbols = mapper.map(bits); @@ -8,6 +8,13 @@ for M = [2 4 8] symbols_pre = db.precode(symbols); symbols_db = db.encode(symbols_pre); + + %%% TEST 1 + symbols_db_ = awgn_channel(symbols_db,"snr_dB",15); + figure(1003042); + showLevelHistogram(symbols_db_,symbols_db,"ref_symbol_uncoded",symbols); + showLevelHistogram(symbols_db_,symbols_db); + symbols_rx = db.decode(symbols_db); bits_rx = PAMmapper(M,0).demap(symbols_rx); [~,~,ber,~] = calc_ber(bits.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); @@ -15,6 +22,12 @@ for M = [2 4 8] symbols_pre_pr = pr1.precode(symbols,"M",M); symbols_db_pr = pr1.encode(symbols_pre_pr,"M",M); + + %%% TEST 2 + symbols_db_pr_ = awgn_channel(symbols_db_pr,"snr_dB",15); + figure(1003043); + % showLevelHistogram(symbols_db_pr_,symbols_db_pr,"ref_symbol_uncoded",symbols); + symbols_rx_pr = pr1.decode(symbols_db_pr,"M",M); bits_rx_pr = mapper.demap(symbols_rx_pr); [~,~,ber_pr,~] = calc_ber(bits.signal,bits_rx_pr.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); @@ -36,6 +49,8 @@ for M = [2 4 8] ' | max abs diff: ',sprintf('%.3g',dec_maxdiff)]); disp('---'); + + f = figure("Name",sprintf("PAM-%d Partial-Response Histograms",M)); t = tiledlayout(f,1,4,"TileSpacing","compact","Padding","compact"); title(t,sprintf("PAM-%d Symbol Histograms",M)); diff --git a/Theory/Signal_Processing/mlse_minimal_example.m b/Theory/Signal_Processing/mlse_minimal_example.m index 04e74db..8266f0f 100644 --- a/Theory/Signal_Processing/mlse_minimal_example.m +++ b/Theory/Signal_Processing/mlse_minimal_example.m @@ -3,29 +3,10 @@ M = 4; randkey = 2; fsym = 112e9; -%%%%% PRBS Generation in correct shape for Modulation Format %%%%%% -O = 18; %order of prbs -N = 2^(O-1); %length of prbs -[~,seed] = prbs(O,1); %initialize first seed of prbs -bitpattern=[]; - -if useprbs - for i = 1:log2(M) - [bitpattern(:,i),seed] = prbs(O,N,seed); - end -else - s = RandStream('twister','Seed',randkey); - for i = 1:log2(M) - bitpattern(:,i) = randi(s,[0 1], N, 1); - end -end - -if M == 6 - bitpattern = reshape(bitpattern,[],1); - bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); -end - -Tx_bits = Informationsignal(bitpattern); +Tx_bits = Signalgenerator( ... + "form", signalform.prms, ... + "M", M, ... + "order", 15).process(); Digi_Mod = PAMmapper(M,0); Symbols_tx = Digi_Mod.map(Tx_bits); @@ -52,7 +33,7 @@ else Symbols.spectrum("fignum",129,"displayname",['coeff:',num2str(coeff)]); - Symbols = MLSE("DIR",coeff,"duobinary_output",0,"trellis_states",Digi_Mod.levels,"M",M).process(Symbols); + Symbols = MLSE("DIR",coeff,"duobinary_output",0,"trellis_states",Digi_Mod.levels,"M",M).process(Symbols,Symbols_tx); Rx_bits = PAMmapper(M,0).demap(Symbols); diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_ROP.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_ROP.m index 5a3c530..9f276d4 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_ROP.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_ROP.m @@ -2,8 +2,17 @@ % LOAD DATA (PAM-4, sweep over ROP) % ============================================================ database_type = 'mysql'; -db = DBHandler("dataBase", "labor_highspeed", "type", database_type); - +% db = DBHandler("dataBase", "labor_highspeed", "type", database_type); + 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.server = "192.168.178.192";% "134.245.243.254"; + dsp_options.user = "silas"; + dsp_options.password = "silas"; + dsp_options.storage_path = 'W:\labdata\sioe_labor\'; + db = DBHandler("dataBase", [dsp_options.dataBase],... + "type", dsp_options.database_type,... + "server", dsp_options.server,... + "user", dsp_options.user, "password", dsp_options.password); pam_level = 4; fiberL = 1; % 1 km wlen = 1310; @@ -59,7 +68,7 @@ for k = 1:numel(curves) precoded = decide_precoded(pam_level,curves(k).eq); cfg = struct; - cfg.x_axis = 'power_mzm'; % ROP axis + cfg.x_axis = 'power_rop'; % ROP axis cfg.y_axis = 'BER'; cfg.agg = 'min'; cfg.outlier = 'none'; diff --git a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/run_dsp_from_db.m b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/run_dsp_from_db.m index e11693f..27b8aca 100644 --- a/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/run_dsp_from_db.m +++ b/projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/run_dsp_from_db.m @@ -48,7 +48,7 @@ fp = QueryFilter(); % fp.where('Runs', 'run_id','EQUALS', 2776); M = 4; fp.where('Runs', 'pam_level','EQUALS', M); -fp.where('Runs', 'bitrate','EQUALS', 360e9);%360,390 +fp.where('Runs', 'bitrate','EQUALS', 420e9);%360,390 % fp.where('Runs', 'symbolrate','EQUALS', 195e9); fp.where('Runs', 'fiber_length','EQUALS', 2); fp.where('Runs', 'is_mpi','EQUALS', 0); diff --git a/projects/Diss/investigate_noise_enhancement.m b/projects/Diss/investigate_noise_enhancement.m new file mode 100644 index 0000000..68dd625 --- /dev/null +++ b/projects/Diss/investigate_noise_enhancement.m @@ -0,0 +1,257 @@ + +% Minimal IMDD model for investigating FFE noise enhancement. +% Parameters mirror the current IMDD_base_system/imdd_it.m, +% IMDD_base_system/imdd_model.m, and the FFE branch in dsp_scope_signal.m. + +clearvars; +close all; + +%% Sweep parameters from imdd_it.m +fsym_values = (152:16:200).*1e9; + +precomp = 0; +laser_wavelength = 1310; +M = 4; +link_length = 2; +channel_alpha = 0.5; +duob_mode = db_mode.db_precoded; +decoding_mode = db_decoder.sequencedetection; +channel_mode = channel_model.physical; +channel_snr_dB = 20; +rop_values = -8:2:4; + +%% Inner model parameters from imdd_model.m +bitrate = []; +apply_pulsef = 0; +fdac = 256e9; +fadc = 256e9; +random_key = 1; + +rcalpha = 0.05; +kover = 16; + +vbias_rel = 0.5; +u_pi = 3; +vbias = -vbias_rel*u_pi; + +laser_linewidth = 0; +eml_alpha = 0; + +len_tr = 4096*2; +pf_ncoeffs = 1; +mu_dc = 0.005; + +nRates = numel(fsym_values); +nRops = numel(rop_values); +eq_names = ["FFE","VNLE"]; +nEqs = numel(eq_names); +emptyResult = struct("fsym",[],"ROP",[],"equalizer",[],"BER",[],"SNR",[],"numBitErr",[], ... + "eq_noise",[],"eq_signal_sd",[],"eq",[],"Rx_sig",[],"Tx_sig",[],"Tx_symbols",[]); +results = repmat(emptyResult,nRates,nRops,nEqs); +cols = flip(cbrewer2('RdYlBu',nRates+4)); +midIdx = floor(size(cols,1)/2) + (-1:2); +cols(midIdx,:) = []; + +for i = 1:numel(fsym_values) + fsym = fsym_values(i); + fprintf("Running %.0f GBd (%d/%d)\n", fsym*1e-9, i, numel(fsym_values)); + + if isempty(bitrate) + bitrateCurrent = fsym * log2(M); + end + + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + [Digi_sig, Symbols, ~] = 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, ... + "duobinary_mode",duob_mode, ... + "mrds_code",0,"mrds_blocklength",512).process(); + + rateResults = repmat(emptyResult,nRops,nEqs); + + + Tx_sig = Digi_sig; + + El_sig = M8199B("kover",kover).process(Digi_sig); + + tx_bwl = 85e9; + El_sig = Filter("filtdegree",4,"f_cutoff",tx_bwl,"fs",fdac*kover, ... + "filterType",filtertypes.bessel_inp,"active",true).process(El_sig); + + El_sig = El_sig.normalize("mode","oneone"); + scaling = 0.7*(u_pi/2-abs(vbias-u_pi/2)); + El_sig = El_sig .* scaling; + + 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, ... + "alpha",eml_alpha).process(El_sig); + + Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length, ... + "alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); + + parfor j = 1:nRops + rop = rop_values(j); + fprintf(" ROP %.1f dBm (%d/%d)\n", rop, j, nRops); + + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ... + "amplification_db",rop).process(Opt_sig); + + Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08, ... + "responsivity",0.6,"temperature",20,"nep",1.8e-11).process(Rx_sig); + + rx_bwl = 90e9; + Rx_sig = Filter("filtdegree",4,"f_cutoff",rx_bwl,"fs",fdac*kover, ... + "filterType",filtertypes.butterworth,"active",true).process(Rx_sig); + + Lp_scpe = Filter("filtdegree",4,"f_cutoff",110e9,"fs",fadc, ... + "filterType",filtertypes.butterworth,"active",true); + + Rx_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); + + txPulseformer = []; + if apply_pulsef + txPulseformer = Pform; + end + + Scpe_sig = preprocessSignal(Rx_sig, Symbols, fsym, ... + "mode","auto", ... + "tx_pulseformer",txPulseformer, ... + "debug_plots",0); + Scpe_sig.signal = Scpe_sig.signal(1:2*Symbols.length); + Scpe_sig.signal = real(Scpe_sig.signal); + + for k = 1:nEqs + switch eq_names(k) + case "FFE" + eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ... + "mu_dd",1e-1,"mu_tr",0.4,"order",50, ... + "sps",2,"decide",0,"optmize_mus",1,"dd_mode",1, ... + "adaption_technique","nlms","mu_dc",1.021e-05); + case "VNLE" + eq_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr, ... + "mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0, ... + "order",[25,2,2],"sps",2,"decide",0, ... + "optmize_mus",1,"mu_dc",0.1); + end + + [eq_signal_sd, eq_noise] = eq_.process(Scpe_sig, Symbols); + + eq_signal_hd = PAMmapper(M,0).quantize(eq_signal_sd); + rx_bits = PAMmapper(M,0).demap(eq_signal_hd); + tx_bits_demapped = PAMmapper(M,0).demap(Symbols); + [~, errors, ber] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, ... + "skip_front",10,"skip_end",10); + + localResult = emptyResult; + localResult.fsym = fsym; + localResult.ROP = rop; + localResult.equalizer = eq_names(k); + localResult.BER = ber; + localResult.SNR = snr(eq_signal_sd.signal, eq_noise.signal); + localResult.numBitErr = errors; + localResult.eq_noise = eq_noise; + localResult.eq_signal_sd = eq_signal_sd; + localResult.eq = eq_; + localResult.Rx_sig = Rx_sig; + localResult.Tx_sig = Tx_sig; + localResult.Tx_symbols = Symbols; + rateResults(j,k) = localResult; + end + end + results(i,:,:) = reshape(rateResults,1,nRops,nEqs); + +end + + + + +%% Plot section +close all +berGrid = reshape([results.BER],nRates,nRops,nEqs); + +figure(400); +clf; +for k = 1:nEqs + subplot(1,nEqs,k); + plot(fsym_values.*1e-9, berGrid(:,:,k), "LineWidth",0.8, ... + "LineStyle","-","Marker",".","MarkerSize",15); + set(gca,"YScale","log"); + grid on; + grid minor; + xlabel("Symbol rate (GBd)"); + ylabel("BER"); + title(eq_names(k) + " BER vs. Symbol Rate"); + legend(compose("ROP %.1f dBm",rop_values),"Location","best"); +end + +%% +figure(500); +clf; +for k = 1:nEqs + subplot(1,nEqs,k); + plot(rop_values, berGrid(:,:,k).', "LineWidth",0.8, ... + "LineStyle","-","Marker",".","MarkerSize",15); + set(gca,"YScale","log"); + grid on; + grid minor; + xlabel("ROP"); + ylabel("BER"); + title(eq_names(k) + " BER vs. ROP"); + legend(compose("%.0f GBd",fsym_values.*1e-9),"Location","best"); +end + +%% +plotResults = reshape(results(:,end,:),nRates,nEqs); +for k = 1:nEqs +for i = 1:nRates + fsym = plotResults(i,k).fsym; + eq_noise = plotResults(i,k).eq_noise; + eq_signal_sd = plotResults(i,k).eq_signal_sd; + Rx_sig = plotResults(i,k).Rx_sig; + Tx_sig = plotResults(i,k).Tx_sig; + Tx_symbols = plotResults(i,k).Tx_symbols; + displayName = sprintf('%s, %d GBd, ROP %.1f dBm',plotResults(i,k).equalizer,fsym.*1e-9,plotResults(i,k).ROP); + + showEQNoisePSD(eq_noise, "fignum", 1, "displayname", displayName,"color",cols(i,:)); + xlim([0, 120]) + % mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/03_DSP_Techniques/tikz/ffe_noise_enhancement/noise_psd.tex") + + Tx_sig.normalize("mode","rms").spectrum("displayname",displayName,'fignum',2,'normalizeTo0dB',0,"color",cols(i,:),"linestyle",'-','HandleVisibility','on'); + xlim([0, 120]) + % mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/03_DSP_Techniques/tikz/ffe_noise_enhancement/tx_signals.tex") + + Rx_sig.normalize("mode","rms").spectrum("displayname",displayName,'fignum',3,'normalizeTo0dB',0,"color",cols(i,:)); + xlim([0, 120]) + % mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/03_DSP_Techniques/tikz/ffe_noise_enhancement/rx_signals.tex") + + showEQNoiseSNR(eq_signal_sd, eq_noise,"displayname",displayName,"color",cols(i,:),"fignum",4); + xlim([0, 120]) + % mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/03_DSP_Techniques/tikz/ffe_noise_enhancement/snr_after_ffe.tex") + + % showEQfilter(plotResults(i,k).eq.e, eq_signal_sd.fs.*2,"displayname",displayName,'fignum',5,"color",cols(i,:)); + % xlim([0, 120]); + + fprintf('%s BER %.2e \n',plotResults(i,k).equalizer,plotResults(i,k).BER) + + if i == 1 || i == nRates-1 || i == nRates + show2Dconstellation(PAMmapper(M,0).get_scaling.*eq_signal_sd.signal(1:end), PAMmapper(M,0).get_scaling.*Tx_symbols.signal(1:end), ... + "displayname",displayName, ... + "fignum",5+(k-1)*nRates+i, ... + "clear",i == 1 && k == 1,"rasterize", false); + mat2tikz_improved(sprintf("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/03_DSP_Techniques/tikz/ffe_noise_enhancement/noise_correlation_%s.tex",displayName),"cleanfigure",true); + eq_signal_sd.eye(fsym,M,"fignum",8+(k-1)*nRates+i,"mode",1); + end +end +end + +%% + + diff --git a/projects/Diss/peak_distortion.m b/projects/Diss/peak_distortion.m new file mode 100644 index 0000000..79a1d69 --- /dev/null +++ b/projects/Diss/peak_distortion.m @@ -0,0 +1,66 @@ +% Build convolution matrix and solve c such that conv(x,c) approx delta +x = [0.05; -0.12; 0.3; 1; 0.25; -0.1; 0.04]; % channel pulse response +Ntaps = 7; +mu = 0.002; +Nit = 500; + +[c, hhist] = peakDistEqualizer(x, Ntaps, mu, Nit); + +h_eff = conv(x,c); + +stem(h_eff) +grid on +xlabel('Sample index') +ylabel('Combined impulse response') + + +function [c, h_hist] = peakDistEqualizer(x, Ntaps, mu, Nit) +% Peak-distortion / Lucky sign-gradient equalizer +% +% x : sampled channel pulse response, main cursor should be near center +% Ntaps : odd number of equalizer taps +% mu : step size +% Nit : number of training iterations +% +% c : equalizer coefficients +% h : combined impulse response h = conv(x,c) + + x = x(:); + assert(mod(Ntaps,2)==1, 'Ntaps must be odd'); + + midTap = (Ntaps+1)/2; + c = zeros(Ntaps,1); + c(midTap) = 1; % start with through connection + + h_hist = cell(Nit,1); + + for it = 1:Nit + h = conv(x,c); + h_hist{it} = h; + + % locate main cursor + [~, mainIdx] = max(abs(h)); + + % normalize main cursor conceptually + h = h / h(mainIdx); + + % update non-center equalizer taps + for j = 1:Ntaps + if j == midTap + continue + end + + % map equalizer tap j to corresponding ISI sample + isiIdx = mainIdx + (j - midTap); + + if isiIdx >= 1 && isiIdx <= length(h) + c(j) = c(j) - mu * sign(h(isiIdx)); + end + end + + % optional: normalize center/main gain + h = conv(x,c); + [~, mainIdx] = max(abs(h)); + c = c / h(mainIdx); + end +end \ No newline at end of file diff --git a/projects/IMDD_base_system/imdd_it.m b/projects/IMDD_base_system/imdd_it.m index 7281bcb..1ae0b54 100644 --- a/projects/IMDD_base_system/imdd_it.m +++ b/projects/IMDD_base_system/imdd_it.m @@ -11,25 +11,25 @@ % -> FFE / DFE / VNLE+MLSE / DB target DSP branches % -> result packages stored in DataStorage -close all; +% close all; if 1 uloops = struct; uloops.precomp = [0]; - uloops.fsym = 180e9; - % uloops.bitrate = 100e9; + uloops.fsym = [152:8:212].*1e9; + % uloops.bitrate = 380e9; % uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4]; uloops.laser_wavelength = [1310]; - uloops.M = [2,4,6,8]; - uloops.link_length = 1; + uloops.M = [4]; + uloops.link_length = 2; % uloops.link_length = [0:2:10]; % 1,2,3,5,6,8,10 - uloops.channel_alpha = [0.8]; - uloops.duob_mode = db_mode.no_db; - uloops.decoding_mode = "memoryless"; - uloops.channel_mode = "awgn"; %awgn_alphad + uloops.channel_alpha = [0.5]; + uloops.duob_mode = db_mode.db_precoded; + uloops.decoding_mode = db_decoder.sequencedetection; + uloops.channel_mode = channel_model.physical; uloops.channel_snr_dB = [20]; - uloops.rop = -6; + uloops.rop = 0; wh = DataStorage(uloops); wh.addStorage("ber"); @@ -39,35 +39,119 @@ if 1 end %% -figure + +wh = load("C:\Users\Silas\Documents\MATLAB\Datensätze\Diss_400G\highspeed_loop.mat"); +wh = wh.wh; +M = 6; +rate = 212e9; +precomp = wh.parameter.precomp.values(1); +wavelength = wh.parameter.laser_wavelength.values(1); +link_length = wh.parameter.link_length.values(1); +alpha = wh.parameter.channel_alpha.values(1); +duob_mode = wh.parameter.duob_mode.values(1); +decoding_mode = wh.parameter.decoding_mode.values(1); +channel_mode = wh.parameter.channel_mode.values(1); +snr_dB = wh.parameter.channel_snr_dB.values(1); +rop = wh.parameter.rop.values; +fsym = wh.parameter.fsym.values; + +%% get all ROP values for M = 4 and rate = 92e9: +a = wh.getStoValue('ber',precomp, rate, wavelength, M, link_length,alpha,duob_mode,decoding_mode,channel_mode,snr_dB,rop); + +ber_ffe = getResultBER(a, "ffe_results"); +ber_dfe = getResultBER(a, "dfe_results"); +ber_vnle = getResultBER(a, "vnle_results"); +ber_mlse = getResultBER(a, "mlse_results"); +ber_dbt = getResultBER(a, "dbt_results"); + +figure(44) +clf hold on -if isfield(uloops, 'fsym') && isfield(uloops, 'bitrate') - error('Define only one of uloops.fsym or uloops.bitrate.'); -elseif isfield(uloops, 'fsym') - rateLookup = uloops.fsym; -elseif isfield(uloops, 'bitrate') - rateLookup = uloops.bitrate; -else - error('Define either uloops.fsym or uloops.bitrate.'); -end - -for alpha = uloops.channel_alpha - a=wh.getStoValue('ber',0, rateLookup, uloops.laser_wavelength, uloops.M, uloops.link_length,alpha,uloops.duob_mode,uloops.decoding_mode,uloops.channel_mode,uloops.channel_snr_dB); - gmi_ffe = cellfun(@(x) x.ffe_results.metrics.GMI, a); - alpha_ffe = cellfun(@(x) x.ffe_results.metrics.Alpha, a); - - gmi_mlse = cellfun(@(x) x.mlse_results.metrics.GMI, a); - - plot(uloops.channel_snr_dB,gmi_ffe,'DisplayName',sprintf(''),'LineStyle','-','HandleVisibility','on'); - plot(uloops.channel_snr_dB,gmi_mlse,'DisplayName',sprintf(''),'LineStyle','-','HandleVisibility','on'); -end - -set(gca, 'YScale', 'log'); -% ylim([5e-5 0.5]); -% yline([3.8e-3, 2e-2],'HandleVisibility','off'); -legend +plot(rop,ber_ffe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE"); +plot(rop,ber_dfe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","DFE"); +plot(rop,ber_vnle,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","VNLE"); +plot(rop,ber_mlse,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","MLSE"); +plot(rop,ber_dbt,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","DB target"); +yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); +xlabel('Received Optical Power (dBm)'); +ylabel('Bit Error Rate (BER)'); +title(sprintf('BER vs. ROP | PAM %d | %.0f GBd',M,rate.*1e-9)); +set(gca,'yscale','log'); +set(gca,'Box','on'); +grid on; +grid minor +legend('Interpreter','none'); beautifyBERplot() -ylabel('BER'); + +%% get all baudrate values for M = 4 and fixed ROP: +rop_ = -6; +a = wh.getStoValue('ber',precomp, fsym, wavelength, M, link_length,alpha,duob_mode,decoding_mode,channel_mode,snr_dB,rop_); + +ber_ffe = getResultBER(a, "ffe_results"); +ber_dfe = getResultBER(a, "dfe_results"); +ber_vnle = getResultBER(a, "vnle_results"); +ber_mlse = getResultBER(a, "mlse_results"); +ber_dbt = getResultBER(a, "dbt_results"); + +figure(45) +clf +hold on +plot(fsym.*1e-9,ber_ffe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","FFE"); +plot(fsym.*1e-9,ber_dfe,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","DFE"); +plot(fsym.*1e-9,ber_vnle,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","VNLE"); +plot(fsym.*1e-9,ber_mlse,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","MLSE"); +plot(fsym.*1e-9,ber_dbt,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName","DB target"); +yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); +xlabel('Baudrate (GBd)'); +ylabel('Bit Error Rate (BER)'); +title(sprintf('BER vs. Baudrate | PAM %d | ROP %g dBm',M,rop_)); +set(gca,'yscale','log'); +set(gca,'Box','on'); +grid on; +grid minor +legend('Interpreter','none'); +beautifyBERplot() + + + + + + + + + + + + +% figure +% hold on +% if isfield(uloops, 'fsym') && isfield(uloops, 'bitrate') +% error('Define only one of uloops.fsym or uloops.bitrate.'); +% elseif isfield(uloops, 'fsym') +% rateLookup = uloops.fsym; +% elseif isfield(uloops, 'bitrate') +% rateLookup = uloops.bitrate; +% else +% error('Define either uloops.fsym or uloops.bitrate.'); +% end + +% for alpha = uloops.channel_alpha +% a=wh.getStoValue('ber',0, rateLookup, uloops.laser_wavelength, uloops.M, uloops.link_length,alpha,uloops.duob_mode,uloops.decoding_mode,uloops.channel_mode,uloops.channel_snr_dB); +% gmi_ffe = cellfun(@(x) x.ffe_results.metrics.GMI, a); +% alpha_ffe = cellfun(@(x) x.ffe_results.metrics.Alpha, a); + +% gmi_mlse = cellfun(@(x) x.mlse_results.metrics.GMI, a); + +% plot(uloops.channel_snr_dB,gmi_ffe,'DisplayName',sprintf(''),'LineStyle','-','HandleVisibility','on'); +% plot(uloops.channel_snr_dB,gmi_mlse,'DisplayName',sprintf(''),'LineStyle','-','HandleVisibility','on'); +% end + +% set(gca, 'YScale', 'log'); +% % ylim([5e-5 0.5]); +% % yline([3.8e-3, 2e-2],'HandleVisibility','off'); +% legend +% beautifyBERplot() +% ylabel('BER'); @@ -585,3 +669,24 @@ ylabel('BER'); % % beautifyBERplot(); % legend('BER contour lines'); + +function ber = getResultBER(results, resultField) +ber = nan(size(results)); +resultField = char(resultField); + +for k = 1:numel(results) + if isempty(results{k}) || ~isstruct(results{k}) || ~isfield(results{k}, resultField) + continue + end + + branch = results{k}.(resultField); + if isempty(branch) || ~isstruct(branch) || ~isfield(branch, 'metrics') + continue + end + + val = branch.metrics.BER; + if ~isempty(val) + ber(k) = val(1); + end +end +end diff --git a/projects/IMDD_base_system/imdd_model.m b/projects/IMDD_base_system/imdd_model.m index 7c3f4ea..1e391fa 100644 --- a/projects/IMDD_base_system/imdd_model.m +++ b/projects/IMDD_base_system/imdd_model.m @@ -3,11 +3,11 @@ function [output] = imdd_model(varargin) simulation_mode = 1; %%% Change folder -curFolder = pwd; -funcFolder=fileparts(mfilename('fullpath')); -if ~isempty(funcFolder) - cd(funcFolder); -end +% curFolder = pwd; +% funcFolder=fileparts(mfilename('fullpath')); +% if ~isempty(funcFolder) +% cd(funcFolder); +% end %%% Run parameters % TX @@ -16,8 +16,8 @@ fsym = 150e9; bitrate = []; apply_pulsef = 0; -fdac = 4*fsym;%256e9; -fadc = 4*fsym;%256e9; +fdac = 256e9; +fadc = 256e9; random_key = 1; rcalpha = 0.05; @@ -29,7 +29,7 @@ vbias = -vbias_rel*u_pi; laser_wavelength = 1293; laser_linewidth = 0; -tx_bw_nyquist = 0.65; +tx_bw_nyquist = 0.65; % Channel link_length = 1; @@ -112,18 +112,26 @@ end f_nyquist = fsym/2; -Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); -[Digi_sig,Symbols,Tx_bits] = PAMsource(... - "fsym",fsym,"M",M,"order",16,"useprbs",0,... - "fs_out",fdac,... - "applyclipping",0,"clipfactor",1.5,... - "applypulseform",apply_pulsef,"pulseformer",Pform,... - "randkey",random_key,... - 'duobinary_mode',duob_mode,... - "mrds_code",0,"mrds_blocklength",512).process(); +%%%%% START SIMULATION -Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); +measure_tf = 0; +if measure_tf + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",70,"f_ref",fadc); + Digi_sig = freqresp.buildOFDM(); +else + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + [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,... + 'duobinary_mode',duob_mode,... + "mrds_code",0,"mrds_blocklength",512).process(); +end + +% Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); channel_mode = lower(string(channel_mode)); @@ -131,26 +139,35 @@ switch channel_mode case "physical" %%%%% AWG - % El_sig = M8199A("kover",kover).process(Digi_sig); - El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0).process(Digi_sig); + El_sig = M8199B("kover",kover).process(Digi_sig); + % El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",0,"H_lpf",H_awg).process(Digi_sig); % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); % El_sig = El_sig.setPower(0,"dBm"); %%%%% Low-pass el. components %%%%%% - 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); + % tx_bwl = tx_bw_nyquist.*f_nyquist; + tx_bwl = 85e9; + El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.bessel_inp,"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"); - scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2)); + scaling = 0.7*(u_pi/2-abs(vbias-u_pi/2)); El_sig = El_sig .* scaling; %%%%% 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,"alpha",eml_alpha).process(El_sig); + % figure(10) + % hold on + % scatter(El_sig.signal(1:100000)+vbias,(abs(Opt_sig.signal(1:100000)).^2)*1e3,0.1,'.','DisplayName','Modulator TF') + % ylim([0 u_pi]); + % xlim([-u_pi/2, u_pi/2]+vbias); + % xlabel('Input in V') + % ylabel('abs(Output) in mW') + + % Opt_sig.spectrum("displayname",'Opt Spectrum','fignum',10,'normalizeTo0dB',1); Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); @@ -161,14 +178,14 @@ switch channel_mode 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_bwl = rx_bw_nyquist.*f_nyquist; + rx_bwl = 90e9; 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); + % Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',10,'normalizeTo0dB',1); %%%%%% Scope %%%%%% Rx_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... @@ -193,6 +210,15 @@ switch channel_mode error('Unknown channel_mode "%s". Supported: physical, awgn, awgn_alphaD.', channel_mode); end +if measure_tf + Rx_sig = Rx_sig.resample("fs_out",freqresp.f_ref); + freqresp.estimate(Rx_sig,"fileName",'','save',false); + freqresp.plot() + return; +end + +Rx_sig.spectrum("displayname",'Sampled Rx Spectrum','fignum',1996,'normalizeTo0dB',1); + txPulseformer = []; if apply_pulsef txPulseformer = Pform; @@ -205,7 +231,7 @@ dspParameters.mu_dfe = mu_dfe; dspParameters.mu_dc = mu_dc; dspParameters.use_ffe = 1; dspParameters.use_dfe = 0; -dspParameters.use_vnle_mlse = 1; +dspParameters.use_vnle_mlse = 0; dspParameters.use_dbtgt = 0; dspParameters.use_dbenc = 0; dspParameters.use_ml_mlse = 0; @@ -213,9 +239,9 @@ dspParameters.pf_ncoeffs = pf_ncoeffs; dspParameters.ffe_order_ffe = [50, 0, 0]; dspParameters.ffe_order_dfe = [50, 0, 0]; dspParameters.dfe_feedback_order = [2, 0, 0]; -dspParameters.ffe_order_vnle = [200, 0, 0]; +dspParameters.ffe_order_vnle = [50, 4, 4]; dspParameters.dfe_order_vnle = dfe_order; -dspParameters.ffe_order_dbtgt = [50, 0, 0]; +dspParameters.ffe_order_dbtgt = [50, 4, 4]; dspParameters.dfe_order_dbtgt = dfe_order; dspParameters.decoding_mode = decoding_mode; diff --git a/projects/IMDD_base_system/minimal_example.m b/projects/IMDD_base_system/minimal_example.m index 4a2d867..2b97947 100644 --- a/projects/IMDD_base_system/minimal_example.m +++ b/projects/IMDD_base_system/minimal_example.m @@ -130,17 +130,20 @@ Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols)); %% -if 0 +if 1 % -------------------- FFE -------------------- - ffe_order = [50, 0, 0]; + ffe_order = [150, 0, 0]; eq_ = 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); - % eq_ = FFE("epochs_tr",4,"epochs_dd",5,"len_tr",4096,"mu_dd",0.01,"mu_tr",0.01,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1); - eq_ = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",512,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",99,"dfe_order",99,"sps",2,"decide",0); + % mu_tr_rls = 9.805e-01; mu_dd_rls = 0.999989348903919; + eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",1.933e-04,"order",ffe_order(1),"sps",2,"decide",0,"optmize_mus",1,"dd_mode",1,"adaption_technique","nlms"); + + % eq_ = FFE("epochs_tr",4,"epochs_dd",5,"len_tr",4096,"mu_dd",0.01,"mu_tr",0.01,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1); + % eq_ = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",512,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",99,"dfe_order",99,"sps",2,"decide",0); output.ffe_results = ffe(eq_,M,Scpe_sig,Symbols,Tx_bits, ... "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ... diff --git a/projects/MLSE_ML_based/minimal_example_huawei/minimal_example.m b/projects/MLSE_ML_based/minimal_example_huawei/minimal_example.m index 504a1b2..67f2b87 100644 --- a/projects/MLSE_ML_based/minimal_example_huawei/minimal_example.m +++ b/projects/MLSE_ML_based/minimal_example_huawei/minimal_example.m @@ -1,5 +1,5 @@ -if 0 +if 1 % A) RUN FULL LOOP M_format = [2,4,6,8]; snr = 10:25; diff --git a/projects/MLSE_ML_based/minimal_example_huawei/ml_mlse_pam.m b/projects/MLSE_ML_based/minimal_example_huawei/ml_mlse_pam.m index 84ebfb4..f33d40d 100644 --- a/projects/MLSE_ML_based/minimal_example_huawei/ml_mlse_pam.m +++ b/projects/MLSE_ML_based/minimal_example_huawei/ml_mlse_pam.m @@ -210,8 +210,8 @@ classdef ml_mlse_pam < handle % ============================================================== % ML-Based Branch Metric Estimation + Viterbi % ============================================================== - debug = 1; - showPlots = 1; + % debug = 1; + % showPlots = 1; nSymbols = ceil(N/obj.sps); diff --git a/projects/MLSE_ML_based/rop_evaluation.m b/projects/MLSE_ML_based/rop_evaluation.m index eb77a3f..fd8677d 100644 --- a/projects/MLSE_ML_based/rop_evaluation.m +++ b/projects/MLSE_ML_based/rop_evaluation.m @@ -35,7 +35,7 @@ for b = 1:numel(baudrates) [Rx_sig_2sps_v1, Symbols_v1, Tx_bits_v1] = standard_link_model( ... "M",M,"fsym",baudrate,"rop",rop,"laser_linewidth",1310, ... - "link_length_m",0,"random_key",1); + "link_length_km",0,"random_key",1); diff --git a/projects/MLSE_ML_based/theoretic_channel_evaluation.m b/projects/MLSE_ML_based/theoretic_channel_evaluation.m index 3eb6f09..cdcdc81 100644 --- a/projects/MLSE_ML_based/theoretic_channel_evaluation.m +++ b/projects/MLSE_ML_based/theoretic_channel_evaluation.m @@ -3,19 +3,24 @@ M = 4; order = 18; randkey = 1; -bitpattern = []; -s = RandStream('twister','Seed',randkey); -for i = 1:log2(M) - N = 2^(order-1); %length of prbs - bitpattern(:,i) = randi(s,[0 1], N, 1); -end +% bitpattern = []; +% s = RandStream('twister','Seed',randkey); +% for i = 1:log2(M) +% N = 2^(order-1); %length of prbs +% bitpattern(:,i) = randi(s,[0 1], N, 1); +% end +% +% if M == 6 +% bitpattern = reshape(bitpattern',[],1); +% bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); +% end +% +% Bits = Informationsignal(bitpattern); -if M == 6 - bitpattern = reshape(bitpattern',[],1); - bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); -end - -Bits = Informationsignal(bitpattern); +Bits = Signalgenerator( ... + "form", signalform.prms, ... + "M", M, ... + "order", order).process(); Symbols = PAMmapper(M,0).map(Bits); Symbols.fs = 200e9; @@ -31,8 +36,7 @@ symbols_filt = Symbols.filter(h,1); %% SHOW FIG 3 in Paper: "ML Base Pre-Eq" -SNR_dB = [20:1:25]; -SNR_db = linspace(12,25,12); +SNR_dB = [15:2:25]; ber_ffe = zeros(size(SNR_dB)); ber_mlse_l5 = zeros(size(SNR_dB)); @@ -43,7 +47,7 @@ ber_ml_mlse_l4 = zeros(size(SNR_dB)); epochs_training = 100; -for i = 1:numel(SNR_dB) +parfor i = 1:numel(SNR_dB) symbols_noi = symbols_filt; symbols_noi.signal = awgn(symbols_filt.signal, SNR_dB(i), 'measured'); % AWGN with given SNR diff --git a/projects/MLSE_nonuniform_states/simulation_nonlin_dsp.m b/projects/MLSE_nonuniform_states/simulation_nonlin_dsp.m index 1896dd4..acfa6e0 100644 --- a/projects/MLSE_nonuniform_states/simulation_nonlin_dsp.m +++ b/projects/MLSE_nonuniform_states/simulation_nonlin_dsp.m @@ -1,7 +1,7 @@ %%% Run parameters % TX M = 4; -m = floor(log2(M)*10)/10; +m = floor(log2(M)*10)/10; fsym = 224e9; apply_pulsef = 1;