silas/imdd_silas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Scattered stuff from Silas during Dissertation

Browse Source
This commit is contained in:
Silas Oettinghaus
2026-06-03 09:05:33 +02:00
parent a91da3b97c
commit 5c2e27687d
37 changed files with 1726 additions and 1338 deletions
Download Patch File Download Diff File Expand all files Collapse all files

105
Classes/00_signals/Signal.m
View File

@@ -200,12 +200,15 @@ classdef Signal
%% Add signals from one signal to another, the first object will sustain
function Sum = plus(X,y)
if isa(y,'Signal')
if isa(X,'Signal') && isa(y,'Signal')
Sum = X;
Sum.signal = X.signal + y.signal;
elseif isnumeric(y)
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
@@ -213,24 +216,42 @@ classdef Signal
%% Add signals from one signal to another, the first object will sustain
function Diff = minus(X,y)
if isa(y,'Signal')
if isa(X,'Signal') && isa(y,'Signal')
Diff = X;
Diff.signal = X.signal - y.signal;
elseif isnumeric(y)
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(y,'Signal')
if isa(X,'Signal') && isa(y,'Signal')
Product = X;
Product.signal = X.signal .* y.signal;
elseif isnumeric(y)
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
@@ -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' ];
@@ -346,6 +368,7 @@ classdef Signal
options.displayname = "";
options.color = [];
options.linestyle = '-';
options.HandleVisibility (1,1) string {mustBeMember(options.HandleVisibility, ["on","off"])} = "on";
options.normalizeToNyquist = 0;
options.normalizeToSamplingRate = 0;
options.addDCoffset = 0;
@@ -380,7 +403,7 @@ classdef Signal
% Divide by 2*pi for the f/fs axis where Nyquist is 0.5.
end
% p_lin = movmean(p_lin,4);
p_lin = movmean(p_lin,10);
if options.normalizeTo0dB
p_lin = p_lin ./ max(p_lin);
@@ -450,9 +473,9 @@ classdef Signal
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
@@ -1007,8 +1030,39 @@ 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;
@@ -1016,6 +1070,7 @@ classdef Signal
difference= maxA-minA;
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,19 +1096,19 @@ classdef Signal
if isa(obj,'Opticalsignal')
title(['Optical Eye ',options.displayname])
ylabel("Power in mW");
y_tickstring = string(linspace(maxA.*1e3,minA.*1e3,6));
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");
y_tickstring = string(linspace(maxA,minA,6));
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");
y_tickstring = string(linspace(maxA,minA,6));
yTickValues = linspace(maxA,minA,5);
min_ = min(obj.signal(100:end-100));
max_ = abs(max(obj.signal(100:end-100)));
end
@@ -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,12 +1239,13 @@ classdef Signal
end
yticks(linspace(0,histpoints,6));
y_tickstring = sprintfc('%.2f', y_tickstring);
yticklabels(y_tickstring);
yTickPositions = linspace(1,histpoints,numel(yTickValues));
yticks(yTickPositions);
yticklabels(sprintfc('%.2f', yTickValues));
xticks(linspace(0,histpoints_horizontal,6))
x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12);
xTickValues = linspace(0, 2/fsym, 6) .* 1e12;
xticks(linspace(1,histpoints_horizontal,numel(xTickValues)))
x_tickstring = sprintfc('%.2f', xTickValues);
xticklabels(x_tickstring);
%

2
Classes/01_transmit/ChannelFreqResp.m
View File

@@ -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');

394
Classes/04_DSP/Equalizer/Copy_of_VNLE.m Normal file
View File

@@ -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

6
Classes/04_DSP/Equalizer/FFE.m
View File

@@ -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

245
Classes/04_DSP/Equalizer/FFE_Kalman.m
View File

@@ -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

197
Classes/04_DSP/Equalizer/FFE_Kalman_Feedback.m
View File

@@ -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

500
Classes/04_DSP/Equalizer/ML_MLSE.m
View File

@@ -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);

97
Classes/04_DSP/Equalizer/VNLE.m
View File

@@ -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

9
Datatypes/channel_model.m Normal file
View File

@@ -0,0 +1,9 @@
classdef channel_model < int32
enumeration
awgn (1)
awgn_alphad (2)
physical (3)
end
end

35
Functions/EQ_structures/dsp_scope_signal.m
View File

@@ -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<NASGU>
@@ -70,8 +70,8 @@ function output = dsp_scope_signal(Scpe_sig_raw, Symbols, Tx_bits, options)
pf_ = Postfilter("ncoeff", pf_ncoeffs, "useBurg", 1); %#ok<NASGU>
mlse_ = MLSE("duobinary_output", 0, 'M', M, 'trellis_states', PAMmapper(M,0).levels); %#ok<NASGU>
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<NASGU>
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<NASGU>
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<NASGU>
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<NASGU>
mlse_db_enc = MLSE("DIR", [1,1], "duobinary_output", 0, "M", M, "trellis_states", PAMmapper(M,0).levels); %#ok<NASGU>
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

6
Functions/EQ_structures/duobinary_target.m
View File

@@ -15,6 +15,7 @@ function [db_results] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx
end
%Duobinary Targeting
db_ref_sequence = Duobinary().encode(tx_symbols);
db_ref_constellation = unique(db_ref_sequence.signal);
@@ -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

3
Functions/EQ_structures/ffe.m
View File

@@ -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

4
Functions/EQ_structures/vnle_postfilter_mlse.m
View File

@@ -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;

158
Functions/EQ_visuals/show2Dconstellation.m
View File

@@ -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

12
Functions/EQ_visuals/showEQNoisePSD.m
View File

@@ -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

65
Functions/EQ_visuals/showEQNoiseSNR.m
View File

@@ -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");
ylabel("SNR (dB)");
xlim([min(w) max(w)]);
edgetick = 2^(nextpow2(options.fs_tx*1e-9));
ticks = -edgetick:16:edgetick;
xticks(ticks);
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
[~,b]=min(abs((-edgetick:16:edgetick)-max(w)));
xlim([-ticks(b+1) ticks(b+1)]);
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);
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

23
Functions/EQ_visuals/showEQfilter.m
View File

@@ -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

115
Functions/EQ_visuals/showLevelHistogram.m
View File

@@ -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

35
Functions/EQ_visuals/showMarkovDiagram.m Normal file
View File

@@ -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

40
Functions/EQ_visuals/showTransitionProbabilities.m Normal file
View File

@@ -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

4
Functions/Lab_helper/loadFreqResp.m
View File

@@ -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};

7
Functions/mat2tikz_improved.m
View File

@@ -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', ...

18
Theory/Mapping_Coding/pam_6_states_analysis.m
View File

@@ -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

19
Theory/Mapping_Coding/partial_response/duobinary_back2back.m
View File

@@ -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));

29
Theory/Signal_Processing/mlse_minimal_example.m
View File

@@ -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);

15
projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/final/FIGURE_ROP.m
View File

@@ -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';

2
projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/run_dsp_from_db.m
View File

@@ -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);

257
projects/Diss/investigate_noise_enhancement.m Normal file
View File

@@ -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
%%

66
projects/Diss/peak_distortion.m Normal file
View File

@@ -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

179
projects/IMDD_base_system/imdd_it.m
View File

@@ -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

86
projects/IMDD_base_system/imdd_model.m
View File

@@ -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;

11
projects/IMDD_base_system/minimal_example.m
View File

@@ -130,18 +130,21 @@ 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",[], ...
"eth_style_symbol_mapping",0);

2
projects/MLSE_ML_based/minimal_example_huawei/minimal_example.m
View File

@@ -1,5 +1,5 @@
if 0
if 1
% A) RUN FULL LOOP
M_format = [2,4,6,8];
snr = 10:25;

4
projects/MLSE_ML_based/minimal_example_huawei/ml_mlse_pam.m
View File

@@ -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);

2
projects/MLSE_ML_based/rop_evaluation.m
View File

@@ -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);

34
projects/MLSE_ML_based/theoretic_channel_evaluation.m
View File

@@ -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

2
projects/MLSE_nonuniform_states/simulation_nonlin_dsp.m
View File

@@ -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;