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

own EQ implemented and returns same performance like Toms

Browse Source
This commit is contained in:
Silas Oettinghaus
2023-07-14 14:42:15 +02:00
parent 3c11e4b2e5
commit 8b3bc688dd
7 changed files with 1287 additions and 164 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
Classes/00_signals/Signal.m
View File

@@ -291,28 +291,10 @@ classdef Signal
delay_n = round(delay_t .* obj.fs); delay_n = round(delay_t .* obj.fs);
% % build "long" hann window to fade the signal in and out
% % -> prevent hard step in the signal!
% hann_wind = hann(200);
% ones_wind = ones(size(obj.signal));
% ones_wind(1:100) = hann_wind(1:100);
% ones_wind(end-100:end) = hann_wind(end-100:end);
%
% % subtract average
% mu = mean(obj.signal,"all");
%
% obj.signal = obj.signal - mu;
%
% %apply hann
% obj.signal = obj.signal .* ones_wind;
%
% %add average again
% obj.signal = obj.signal + mu;
% finally circshift the signal % finally circshift the signal
obj.signal=circshift(obj.signal,delay_n); % obj.signal=circshift(obj.signal,delay_n);
% obj.signal=[obj.signal(delay_n:end); zeros(delay_n-1,1)]; obj.signal=[zeros(delay_n,1); obj.signal(1:end-delay_n) ];
end end
end end

1
Classes/01_transmit/PAMmapper.m
View File

@@ -108,7 +108,6 @@ classdef PAMmapper
end end
thres = thres .* 1/sqrt(5); thres = thres .* 1/sqrt(5);
case 3 case 3
% 8-ASK % 8-ASK
if obj.unipolar==0 if obj.unipolar==0

82
Classes/04_DSP/EQ.m
View File

@@ -102,10 +102,10 @@ classdef EQ
end end
function signalclass_out = process(obj,signalclass_in, reference_signalclass_in) function [signalclass_out,error_log] = process(obj,signalclass_in, reference_signalclass_in)
% actual processing of the signal (steps 1. - 3.) % actual processing of the signal (steps 1. - 3.)
signalclass_in.signal = obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); [signalclass_in.signal,error_log] = obj.process_(signalclass_in.signal', reference_signalclass_in.signal');
signalclass_in.signal = signalclass_in.signal'; signalclass_in.signal = signalclass_in.signal';
% append to logbook % append to logbook
@@ -117,7 +117,7 @@ classdef EQ
end end
function yout = process_(obj,data_in,ref_in) function [yout,error_log] = process_(obj,data_in,ref_in)
%METHOD1 Summary of this method goes here %METHOD1 Summary of this method goes here
% Detailed explanation goes here % Detailed explanation goes here
@@ -211,35 +211,54 @@ classdef EQ
epsilon_ = diag([ones(1,obj.Ne(1))*obj.epsilon(1) ones(1,N2)*obj.epsilon(2) ones(1,N3)*obj.epsilon(3)]); epsilon_ = diag([ones(1,obj.Ne(1))*obj.epsilon(1) ones(1,N2)*obj.epsilon(2) ones(1,N3)*obj.epsilon(3)]);
end end
obj.k0 = obj.delay; % input delay compared to training sequence< obj.k0 = obj.delay; % input delay compared to training sequence
error_log = [];
if 1 % obj.active if 1 % obj.active
%% Calculation of the filter coefficients in training based LMS mode %% Calculation of the filter coefficients in training based LMS mode
e_ = zeros(obj.Ne(1)+N2+N3,1); % initialization of filter coefficients e_ = zeros(obj.Ne(1)+N2+N3,1); % initialization of filter coefficients
% e(ceil(obj.Ne(1)/2)) = 1; % set central tap to 1 (better starting point since it's closer to the expected solution)
b_ = zeros(obj.Nb(1)+Nb2+Nb3,1); b_ = zeros(obj.Nb(1)+Nb2+Nb3,1);
e_dc = mean(data_in); % initilaization of the dc tap with the mean value of the data e_dc = mean(data_in); % initilaization of the dc tap with the mean value of the whole data
% e_save = NaN(361,8.6e5);
% save_ind = 1;
for trainloops = 1:obj.training_loops for trainloops = 1:obj.training_loops
cnt = 1;
m = obj.k0+1; % starting symbol index at the delay compared to the training sequence m = obj.k0+1; % starting symbol index at the delay compared to the training sequence
error_log = [];
for n = obj.K*obj.k0+1:obj.K:obj.K*obj.training_length
m = m+1;
X_1 = data(obj.Ne(1)+n+(obj.K-1):-1:n+obj.K).';
% X_1(X_1~=0) = X_1(X_1~=0) - mean(X_1(X_1~=0));
% n => index in rx data sequence
%Step From: Oversampling(=2) * Startdelay + 1
%Step Width: Oversampling(=2)
%Step To: Oversampling(=2) * Training Length
for n = obj.K*obj.k0+1:obj.K:obj.K*obj.training_length
% m => index in reference sequence
m = m+1;
%
%dc_ = mean(data(obj.Ne(1)+n+(obj.K-1):-1:n+obj.K).');
% cut symbols from rx data sequence
X_1 = data(obj.Ne(1)+n+(obj.K-1):-1:n+obj.K).';
[X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx); [X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx);
input_vec = [X_1;X_2;X_3];
% cut symbols from desired data sequence (correct symbols in training mode)
D_1 = ref(obj.Nb(1)-obj.k0+m-2:-1:m-obj.k0-1).'; D_1 = ref(obj.Nb(1)-obj.k0+m-2:-1:m-obj.k0-1).';
[D_2,D_3] = obj.calc_nl_vecs(D_1,ind_mat_DFE_2nd,ind_mat_DFE_3rd,norm_fac_DFE2,norm_fac_DFE3,delta_DFE2,delta_DFE3,cplx); [D_2,D_3] = obj.calc_nl_vecs(D_1,ind_mat_DFE_2nd,ind_mat_DFE_3rd,norm_fac_DFE2,norm_fac_DFE3,delta_DFE2,delta_DFE3,cplx);
input_vec = [X_1;X_2;X_3];
reference_vec = [D_1;D_2;D_3]; reference_vec = [D_1;D_2;D_3];
error = e_dc + e_.'*input_vec - b_.'*reference_vec - ref_in(m-obj.k0); % error = e_dc + e.'*input_vec - b.'*reference_vec - ref_in(m-obj.k0);
e_ffe = e_.'*input_vec;
e_dfe = b_.'*reference_vec;
error = e_dc + e_ffe - e_dfe - ref_in(m-obj.k0);
%error = e_dc + e_.'*input_vec - b_.'*reference_vec - ref_in(m-obj.k0);
if real(obj.FFEmu) if real(obj.FFEmu)
if obj.l1act if obj.l1act
@@ -247,7 +266,7 @@ classdef EQ
sgn_e(e_~=0) = e_(e_~=0)./abs(e_(e_~=0)); sgn_e(e_~=0) = e_(e_~=0)./abs(e_(e_~=0));
e_ = e_ - obj.rho*sgn_e./(1+epsilon_*abs(e_)) - error*input_vec*obj.FFEmu; e_ = e_ - obj.rho*sgn_e./(1+epsilon_*abs(e_)) - error*input_vec*obj.FFEmu;
else else
e_ = e_ - error*conj(input_vec)*obj.FFEmu; %e = e - error*conj(input_vec)*obj.FFEmu; %e = e - error*input_vec*obj.FFEmu; e_ = e_ - error*conj(input_vec)*obj.FFEmu; %classic LMS gradient decay
end end
else else
if obj.l1act if obj.l1act
@@ -258,21 +277,24 @@ classdef EQ
e_ = e_ - error*input_vec/(input_vec.'*input_vec); e_ = e_ - error*input_vec/(input_vec.'*input_vec);
end end
end end
% e_save(:,save_ind) = e;
% save_ind = save_ind+1;
e_dc = e_dc - obj.DCmu*error;
error_log(end+1) = e_dc;
e_dc = e_dc - obj.DCmu*error;
error_log(cnt,trainloops) = e_dc;
cnt = cnt+1;
if obj.Nb(1) > 0 if obj.Nb(1) > 0
b_ = b_ + obj.DFEmu*error*reference_vec; % Seems like normalized DFE has worse performance b_ = b_ + obj.DFEmu*error*reference_vec; % Seems like normalized DFE has worse performance
end end
% figure(111);stem((e_),'Markersize',2);ylim([-1 1]);title('FFE Filter Taps');
%
% figure(222);
% stem(input_vec);ylim([-3 3]);
% hold on;
% stem(reference_vec);ylim([-3 3]);
% yline(error,'LineWidth',2); title('Input Vector');
% hold off
end end
end end
@@ -355,6 +377,7 @@ classdef EQ
end end
for dd_loop = 1:obj.dd_loops for dd_loop = 1:obj.dd_loops
cnt = obj.training_length+1;
m = 0; m = 0;
output_vec = zeros(1,floor(length(data_in)/obj.K)); % initilaization of the output vector output_vec = zeros(1,floor(length(data_in)/obj.K)); % initilaization of the output vector
dd_DFE = zeros(obj.Nb(1),1); dd_DFE = zeros(obj.Nb(1),1);
@@ -384,10 +407,9 @@ classdef EQ
for k = 1:obj.K:length(data_in) for k = 1:obj.K:length(data_in)
m=m+1; % Symbol index m=m+1; % Symbol index
X_1 = data(obj.Ne(1)+k-1:-1:k).'; X_1 = data(obj.Ne(1)+k-1:-1:k).';
[X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx); [X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx);
if obj.l1act if obj.l1act
@@ -403,6 +425,7 @@ classdef EQ
dd_out(k) = constellation_in_(dd_idx); dd_out(k) = constellation_in_(dd_idx);
end end
if obj.Nb(1) > 0 if obj.Nb(1) > 0
dd_DFE(2:end) = dd_DFE(1:end-1); dd_DFE(2:end) = dd_DFE(1:end-1);
dd_DFE(1) = dd_out(k); dd_DFE(1) = dd_out(k);
@@ -424,13 +447,14 @@ classdef EQ
if mu_mat ~= 0 if mu_mat ~= 0
e_dc = e_dc - obj.DCmu*error; e_dc = e_dc - obj.DCmu*error;
error_log(end+1) = e_dc; error_log(cnt,dd_loop) = e_dc;
cnt = cnt+1;
end end
end end
end end
%figure(2023);plot(error_log(:,1))
% shifting the output sequence by k0 symbols % shifting the output sequence by k0 symbols
yout = (circshift(output_vec.',-(obj.k0))).'; %(circshift(dd_out.',-(obj.k0))).'; yout = (circshift(output_vec.',-(obj.k0))).'; %(circshift(dd_out.',-(obj.k0))).';

727
Classes/04_DSP/EQ_copy.m Normal file
View File

@@ -0,0 +1,727 @@
classdef EQ_copy
%EQ Summary of this class goes here
% Detailed explanation goes here
properties
Ne %Number of feed forward coefficients (1st, 2nd and 3rd order)
Nb %Number of decision feedback coefficients (1st, 2nd and 3rd order)
K %Number of samples per symbol
delay %Delay of incoming signal
training_length %Number of training symbols
training_loops %Number of loops through sequence for training mode
ideal_dfe %Error free DFE decisions
DB_aim %Aim at duobinary output sequence
M %Order of the PAM constellation (only relevant in case of DB aim)
FFEmu %mu parameter for FFE part in training mode (0 means normalized LMS)
DFEmu % mu parameter for DFE part in training mode
dd_loops % Number of loops through sequence for DD mode
DDmu % mu parameters for DD mode (individual value for each order)
DCmu % mu parameter for the dc tap
l1act %Activate/deactive l1 regularization
rho %Parameter for speed of coeff shrinking
epsilon %Reciprocal value of the magnitude of the coeff to converge to zero (1st,2nd,3rd order)
thres %Theshold for neglecting coefficienties (1st,2nd,3rd order)
static_act %Activate/deactive static coefficient reduction
mode2nd %0: no reduction | 1: polynomial | 2: restricted to interval
len_2nd %length of the interval (only for 2nd order mode = 2)
mode3rd %0: no reduction | 1: polynomial | 2: restricted to interval
len_3rd %length of the interval (only for 3rd order mode = 3/4)
plottrain
plotfinal
load_decisions
save_taps
%during simulation
k0
b
b2
b3
e
e2
e3
coeff_number
constellation_in
end
methods
function obj = EQ_copy(options)
%EQ Construct an instance of this class
% Detailed explanation goes here
arguments(Input)
options.Ne = [10 0 0] %Number of feed forward coefficients (1st, 2nd and 3rd order)
options.Nb = [10 0 0]%Number of decision feedback coefficients (1st, 2nd and 3rd order)
options.K = 1 %Number of samples per symbol
options.delay = 0 %Delay of incoming signal
options.training_length = 1024 %Number of training symbols
options.training_loops = 1 %Number of loops through sequence for training mode
options.ideal_dfe = 0 %Error free DFE decisions
options.DB_aim %Aim at duobinary output sequence
options.M = 1 %Order of the PAM constellation (only relevant in case of DB aim)
options.FFEmu = 0 %mu parameter for FFE part in training mode (0 means normalized LMS)
options.DFEmu = 0.005 % mu parameter for DFE part in training mode
options.dd_loops = 1% Number of loops through sequence for DD mode
options.DDmu = [0.0004 0.0004 0.0004 0.0004 ] % mu parameters for DD mode (individual value for each order)
options.DCmu = 0.005 % mu parameter for the dc tap
options.l1act = 0 %Activate/deactive l1 regularization
options.rho = 5e-4%Parameter for speed of coeff shrinking
options.epsilon = [10 100 1000] %Reciprocal value of the magnitude of the coeff to converge to zero (1st,2nd,3rd order)
options.thres = [5e-3 4e-3 5e-4]%Theshold for neglecting coefficienties (1st,2nd,3rd order)
options.static_act = 0 %Activate/deactive static coefficient reduction
options.mode2nd = 1%0: no reduction | 1: polynomial | 2: restricted to interval
options.len_2nd = 1 %length of the interval (only for 2nd order mode = 2)
options.mode3rd = 1%0: no reduction | 1: polynomial | 2: restricted to interval
options.len_3rd = 1%length of the interval (only for 3rd order mode = 3/4)
options.plottrain = 0
options.plotfinal = 0
options.load_decisions = 0
options.save_taps = 0
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
end
function [signalclass_out,error_log] = process(obj,signalclass_in, reference_signalclass_in)
% actual processing of the signal (steps 1. - 3.)
[signalclass_in.signal,error_log] = obj.process_(signalclass_in.signal', reference_signalclass_in.signal');
signalclass_in.signal = signalclass_in.signal';
% append to logbook
lbdesc = ['EQ '];
signalclass_in = signalclass_in.logbookentry(lbdesc);
% write to output
signalclass_out = signalclass_in;
end
function [yout,error_log] = process_(obj,data_in,ref_in)
%METHOD1 Summary of this method goes here
% Detailed explanation goes here
if obj.DB_aim
ref_DB = zeros(size(ref_in));
for k = 1:length(ref_in)
if k == 1
ref_DB(k) = ref_in(k);
else
ref_DB(k) = ref_in(k) + ref_in(k-1);
end
end
ref_in = ref_DB;
end
ref = [zeros(1,obj.Nb(1)-1) ref_in zeros(1,obj.Nb(1))];
if isreal(ref)
cplx = 0;
else
cplx = 1;
end
if obj.static_act
obj.mode2nd = obj.mode2nd + 1;
obj.mode3rd = obj.mode3rd + 1;
else
obj.mode2nd = 1;
obj.mode3rd = 1;
end
if obj.mode2nd == 1
N2 = (obj.Ne(2)*(obj.Ne(2)+1))/2; % Number of coefficients for second order
elseif obj.mode2nd == 2
N2 = obj.Ne(2);
elseif obj.mode2nd == 3
N2 = (obj.len_2nd+1)*(2*obj.Ne(2)-obj.len_2nd)/2;
elseif obj.mode2nd == 4
N2 = (ceil(obj.Ne(2)/2)+1)*(2*obj.Ne(2)-ceil(obj.Ne(2)/2))/2;
end
if obj.mode3rd == 1
if cplx
N3 = obj.Ne(3)^2*(obj.Ne(3)+1)/2;
else
N3 = obj.Ne(3)*(obj.Ne(3)+1)*(obj.Ne(3)+2)/6; % Number of coefficients for third order
end
elseif obj.mode3rd == 2
N3 = obj.Ne(3);
elseif obj.mode3rd == 3
N3 = obj.Ne(3)^2;
elseif obj.mode3rd == 4
N3 = round(1/6*(obj.len_3rd+1)*(obj.len_3rd+2)*(3*obj.Ne(3)-2*obj.len_3rd));
elseif obj.mode3rd == 5
N3 = 2*obj.Ne(3)*obj.len_3rd-obj.len_3rd*(obj.len_3rd+1)+obj.Ne(3);
end
Nb2 = (obj.Nb(2)*(obj.Nb(2)+1))/2;
Nb3 = obj.Nb(3)*(obj.Nb(3)+1)*(obj.Nb(3)+2)/6;
data_in = data_in/sqrt(mean(abs(data_in).^2)); % power normalization of input sequence
if obj.FFEmu == 0
norm_fac2 = sqrt(mean(abs(data_in.^2).^2)); % power normalization for second and third order terms
norm_fac3 = sqrt(mean(abs(data_in.^3).^2)); % (not necessary, but seems to be more stable if applied --> same as different mu values for linear and nl terms)
else
norm_fac2 = 1;
norm_fac3 = 1;
end
norm_fac_DFE2 = sqrt(mean(abs(ref_in.^2).^2)); % same for DFE input (reference)
norm_fac_DFE3 = sqrt(mean(abs(ref_in.^3).^2));
data = [zeros(1,floor(obj.Ne(1)/2)) data_in zeros(1,obj.Ne(1))];
delta_2 = round((obj.Ne(1)-obj.Ne(2))/2);
delta_3 = round((obj.Ne(1)-obj.Ne(3))/2);
delta_DFE2 = 1;%round((obj.Nb-obj.Nb(2))/2);
delta_DFE3 = 1;%round((obj.Nb-obj.Nb(3))/2);
% calculate the indices for the combination of second and third order symbols
% - done in advance because it's the same for each iteration, so time
% can be saved
[ind_mat_2nd,ind_mat_3rd] = obj.calc_ind(obj.Ne(2),N2,obj.Ne(3),N3,obj.mode2nd,obj.mode3rd,obj.len_2nd,obj.len_3rd,cplx);
[ind_mat_DFE_2nd,ind_mat_DFE_3rd] = obj.calc_DFE_ind(obj.Nb(2),Nb2,obj.Nb(3),Nb3);
if obj.l1act
epsilon_ = diag([ones(1,obj.Ne(1))*obj.epsilon(1) ones(1,N2)*obj.epsilon(2) ones(1,N3)*obj.epsilon(3)]);
end
obj.k0 = obj.delay; % input delay compared to training sequence
error_log = [];
if 1 % obj.active
%% Calculation of the filter coefficients in training based LMS mode
e_ = zeros(obj.Ne(1)+N2+N3,1); % initialization of filter coefficients
b_ = zeros(obj.Nb(1)+Nb2+Nb3,1);
e_dc = mean(data_in); % initilaization of the dc tap with the mean value of the whole data
for trainloops = 1:obj.training_loops
cnt = 1;
m = obj.k0+1; % starting symbol index at the delay compared to the training sequence
% n => index in rx data sequence
%Step From: Oversampling(=2) * Startdelay + 1
%Step Width: Oversampling(=2)
%Step To: Oversampling(=2) * Training Length
for n = obj.K*obj.k0+1:obj.K:obj.K*obj.training_length
% m => index in reference sequence
m = m+1;
%
%dc_ = mean(data(obj.Ne(1)+n+(obj.K-1):-1:n+obj.K).');
% cut symbols from rx data sequence
X_1 = data(obj.Ne(1)+n+(obj.K-1):-1:n+obj.K).';
[X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx);
input_vec = [X_1;X_2;X_3];
% cut symbols from desired data sequence (correct symbols in training mode)
D_1 = ref(obj.Nb(1)-obj.k0+m-2:-1:m-obj.k0-1).';
[D_2,D_3] = obj.calc_nl_vecs(D_1,ind_mat_DFE_2nd,ind_mat_DFE_3rd,norm_fac_DFE2,norm_fac_DFE3,delta_DFE2,delta_DFE3,cplx);
reference_vec = [D_1;D_2;D_3];
e_ffe = e_.'*input_vec;
e_dfe = b_.'*reference_vec;
error = e_dc + e_ffe - e_dfe - ref_in(m-obj.k0);
%error = e_dc + e_.'*input_vec - b_.'*reference_vec - ref_in(m-obj.k0);
if real(obj.FFEmu)
if obj.l1act
sgn_e = e_;
sgn_e(e_~=0) = e_(e_~=0)./abs(e_(e_~=0));
e_ = e_ - obj.rho*sgn_e./(1+epsilon_*abs(e_)) - error*input_vec*obj.FFEmu;
else
e_ = e_ - error*conj(input_vec)*obj.FFEmu; %classic LMS gradient decay
end
else
if obj.l1act
sgn_e = e_;
sgn_e(e_~=0) = e_(e_~=0)./abs(e_(e_~=0));
e_ = e_ - obj.rho*sgn_e./(1+epsilon_*abs(e_)) - error*input_vec/(input_vec.'*input_vec);
else
e_ = e_ - error*input_vec/(input_vec.'*input_vec);
end
end
e_dc = e_dc - obj.DCmu*error;
error_log(cnt,trainloops) = error;
cnt = cnt+1;
if obj.Nb(1) > 0
b_ = b_ + obj.DFEmu*error*reference_vec; % Seems like normalized DFE has worse performance
end
% figure(111);stem((e_),'Markersize',2);ylim([-1 1]);title('FFE Filter Taps');
%
% figure(222);
% stem(input_vec);ylim([-3 3]);
% hold on;
% stem(reference_vec);ylim([-3 3]);
% yline(error,'LineWidth',2); title('Input Vector');
% hold off
end
end
%%
% Plot the intermediate coefficients after training mode
obj.b = b_(1:obj.Nb(1));
obj.b2 = b_(obj.Nb(1)+1:obj.Nb(1)+Nb2);
obj.b3 = b_(obj.Nb(1)+Nb2+1:end);
obj.e = e_(1:obj.Ne(1));
obj.e2 = e_(obj.Ne(1)+1:obj.Ne(1)+N2);
obj.e3 = e_(obj.Ne(1)+N2+1:end);
if obj.plottrain
figure(8052)
sgtitle('Training Coeff')
subplot(2,3,1); stem((obj.e),'Markersize',2);
title('FFE coeff linear')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,2); stem(obj.e2,'Markersize',2);
title('FFE coeff nl 2nd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,3); stem(obj.e3,'Markersize',2);
title('FFE coeff nl 3rd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,4);stem(obj.b,'Markersize',2);
title('DFE coeff linear')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,5);stem(obj.b2,'Markersize',2);
title('DFE coeff nl 2nd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,6);stem(obj.b3,'Markersize',2);
title('DFE coeff nl 3rd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
%set(gcf,'Position',[200 500 700 400])
end
if obj.l1act
neg_lin = find(abs(obj.e) < obj.thres(1));
neg_2nd = find(abs(obj.e2) < obj.thres(2));
neg_3rd = find(abs(obj.e3) < obj.thres(3));
neg = [neg_lin;neg_2nd+obj.Ne(1);neg_3rd+obj.Ne(1)+N2]; % indices of the neglected coefficients
rel_lin = find(abs(obj.e) >= obj.thres(1));
rel_2nd = find(abs(obj.e2) >= obj.thres(2));
rel_3rd = find(abs(obj.e3) >= obj.thres(3));
obj.coeff_number = length(rel_lin)+2*length(rel_2nd)+3*length(rel_3rd);
rel = [rel_lin;rel_2nd+obj.Ne(1);rel_3rd+obj.Ne(1)+N2]; % indices of the relevant coefficients
e_(neg) = 0;
ind_mat_2nd(neg_2nd,:) = [];
ind_mat_3rd(neg_3rd,:) = [];
end
%% decision directed mode
if ~obj.DB_aim
constellation_in_ = unique(ref_in); % getting the symbol constellation from reference data
else
if obj.M == 2
constellation_in_ = [-3 -2 -1 0 1 2 3]/sqrt(5)*2;
elseif obj.M == 2.5
constellation_in_ = [-5 -4 -3 -2 -1 0 1 2 3 4 5]/sqrt(10)*2;
elseif obj.M == 3
constellation_in_ = [-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7]/sqrt(21)*2;
else
constellation_in_ = unique(ref_in);
end
end
obj.constellation_in = constellation_in_;
if obj.l1act
coeff = [e_(rel);b_]; % combine FFE and DFE coefficient vectors for DD mode
else
coeff = [e_;b_];
end
for dd_loop = 1:obj.dd_loops
cnt = obj.training_length+1;
m = 0;
output_vec = zeros(1,floor(length(data_in)/obj.K)); % initilaization of the output vector
dd_DFE = zeros(obj.Nb(1),1);
D_2 = zeros(Nb2,1);
D_3 = zeros(Nb3,1);
if all(obj.DDmu == obj.DDmu(1))
mu_mat = obj.DDmu(1);
else
if obj.l1act
mu_mat = diag([ones(1,length(rel_lin))*obj.DDmu(1) ones(1,length(rel_2nd))*obj.DDmu(2) ones(1,length(rel_3rd))*obj.DDmu(3) ones(1,obj.Nb)*obj.DDmu(4)]);
else
mu_mat = diag([ones(1,obj.Ne(1))*obj.DDmu(1) ones(1,N2)*obj.DDmu(2) ones(1,N3)*obj.DDmu(3) ones(1,obj.Nb(1)+Nb2+Nb3)*obj.DDmu(4)]);
end
end
if obj.load_decisions
pathn = evalin('base','modeldir');
temp = load([pathn, 'MLSE_out', '.mat']) ;
%eval(['dd_out_vals = temp.', 'a', ';']) ;
dd_out_vals=temp.a;
dd_out = zeros(size(data_in));
dd_out(1:2:length(data_in)) = dd_out_vals;
else
dd_out = zeros(size(data_in));
end
for k = 1:obj.K:length(data_in)
m=m+1; % Symbol index
X_1 = data(obj.Ne(1)+k-1:-1:k).';
[X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx);
if obj.l1act
input_vec = [X_1(rel_lin);X_2;X_3;-dd_DFE;-D_2;-D_3];
else
input_vec = [X_1;X_2;X_3;-dd_DFE;-D_2;-D_3];
end
output_vec(m) = e_dc + input_vec.'*coeff;
if ~obj.load_decisions
[~,dd_idx] = min(abs(output_vec(m) - constellation_in_)); % decision for closest constellation point
dd_out(k) = constellation_in_(dd_idx);
end
if obj.Nb(1) > 0
dd_DFE(2:end) = dd_DFE(1:end-1);
dd_DFE(1) = dd_out(k);
if obj.ideal_dfe && m > obj.k0
dd_DFE(1) = ref_in(m-obj.k0);
end
[D_2,D_3] = obj.calc_nl_vecs(dd_DFE,ind_mat_DFE_2nd,ind_mat_DFE_3rd,norm_fac_DFE2,norm_fac_DFE3,delta_DFE2,delta_DFE3,cplx);
end
% if dd_loop ~= 21
error = output_vec(m) - dd_out(k);
% else
% error = 0;
% end
coeff = coeff - mu_mat*error*conj(input_vec);
% e_save(:,save_ind) = coeff;
% save_ind = save_ind+1;
if mu_mat ~= 0
e_dc = e_dc - obj.DCmu*error;
error_log(cnt,dd_loop) = error;
cnt = cnt+1;
end
end
end
%figure(2023);plot(error_log(:,1))
% shifting the output sequence by k0 symbols
yout = (circshift(output_vec.',-(obj.k0))).'; %(circshift(dd_out.',-(obj.k0))).';
e_ = coeff(1:end-obj.Nb(1)-Nb2-Nb3);
b_ = coeff(end-obj.Nb(1)-Nb2-Nb3+1:end);
if obj.l1act
obj.e = e_(1:length(rel_lin));
obj.e2 = e_(length(rel_lin)+1:length(rel_lin)+length(rel_2nd));
obj.e3 = e_(length(rel_lin)+length(rel_2nd)+1:end);
else
obj.e = e_(1:obj.Ne(1));
obj.e2 = e_(obj.Ne(1)+1:obj.Ne(1)+N2);
obj.e3 = e_(obj.Ne(1)+N2+1:end);
end
obj.b = b_(1:obj.Nb(1));
obj.b2 = b_(obj.Nb(1)+1:obj.Nb(1)+Nb2);
obj.b3 = b_(obj.Nb(1)+Nb2+1:end);
% plot the final coefficients after DD mode
if obj.plotfinal
figure(8054)
if obj.l1act
sgtitle('Final Coeff')
subplot(2,3,1); stem(rel_lin,obj.e,'Markersize',2);
title('FFE coeff linear')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,2); stem(rel_2nd,obj.e2,'Markersize',2);
title('FFE coeff nl 2nd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,3); stem(rel_3rd,obj.e3,'Markersize',2);
title('FFE coeff nl 3rd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,4);stem(obj.b,'Markersize',2);
title('DFE coeff linear')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,5);stem(obj.b2,'Markersize',2);
title('DFE coeff nl 2nd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,6);stem(obj.b3,'Markersize',2);
title('DFE coeff nl 3rd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
else
sgtitle('Final Coeff')
subplot(2,3,1); stem(obj.e/max(e_),'Markersize',2);
title('FFE coeff linear')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,2); stem(obj.e2,'Markersize',2);
title('FFE coeff nl 2nd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,3); stem(obj.e3,'Markersize',2);
title('FFE coeff nl 3rd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,4);stem(obj.b,'Markersize',2);
title('DFE coeff linear')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,5);stem(obj.b2,'Markersize',2);
title('DFE coeff nl 2nd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
subplot(2,3,6);stem(obj.b3,'Markersize',2);
title('DFE coeff nl 3rd')
xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12)
end
set(gcf,'Position',[1000 500 700 400])
end
% save frequency response to the work space
if obj.save_taps
% save the FFE coefficients to the work space
% pathn = evalin('base','modeldir');
% eval([obj.field_ffe, ' = obj.e ;']) ;
% eval([obj.field_dfe, ' = b ;']) ;
% eval(['save(''', pathn, '\',obj.filen,''', ''', obj.field_ffe,''', ''',obj.field_dfe,''') ;']) ;
save("coefficients",obj.e, obj.b);
end
else
yout = data_in;
end
end
function [X_2,X_3] = calc_nl_vecs(obj,X_1,ind_mat_2,ind_mat_3,norm_fac2,norm_fac3,delta_2,delta_3,cplx)
% calculation of the vectors containing all combinations of input symbols
% of second and third order based on the linear symbols
if ind_mat_2(1) > 0
input_vec_se = X_1(delta_2:end)/norm_fac2;%(K*(k0-1):end)
X_2 = input_vec_se(ind_mat_2(:,1)).*input_vec_se(ind_mat_2(:,2));
else
X_2 = [];
end
if ind_mat_3(1) > 0
if cplx
input_vec_th = X_1(delta_3:end)/norm_fac3;
X_3 = input_vec_th(ind_mat_3(:,1)).*input_vec_th(ind_mat_3(:,2)).*conj(input_vec_th(ind_mat_3(:,3)));
else
input_vec_th = X_1(delta_3:end)/norm_fac3;
X_3 = input_vec_th(ind_mat_3(:,1)).*input_vec_th(ind_mat_3(:,2)).*input_vec_th(ind_mat_3(:,3));
end
else
X_3 = [];
end
end
function [ind_mat_2nd,ind_mat_3rd] = calc_ind(obj,Ne2,N2,Ne3,N3,mode2nd,mode3rd,len_2nd,len_3rd,cplx)
if Ne2 > 0
ind_mat_2nd = NaN(N2,2);
count=1;
if mode2nd == 1
for t = 1:Ne2
for u = t:Ne2
ind_mat_2nd(count,:) = [t u];
count = count + 1 ;
end
end
elseif mode2nd == 2
for t = 1:Ne2
ind_mat_2nd(t,:) = [t t];
end
elseif mode2nd == 3
for t = 1:Ne2
for u = t:Ne2
if u-t<=len_2nd
ind_mat_2nd(count,:) = [t u];
count = count + 1 ;
end
end
end
elseif mode2nd == 4
for t = 1:Ne2
for u = t:Ne2
if u-t<=ceil(Ne2/2)
ind_mat_2nd(count,:) = [t u];
count = count + 1 ;
end
end
end
end
else
ind_mat_2nd = 0;
end
if Ne3 > 0
ind_mat_3rd = NaN(N3,3);
count=1;
if mode3rd == 1
if cplx
for t = 1:Ne3
for u = t:Ne3
for v = 1:Ne3
ind_mat_3rd(count,:) = [t u v];
count = count + 1 ;
end
end
end
else
for t = 1:Ne3
for u = t:Ne3
for v = u:Ne3
ind_mat_3rd(count,:) = [t u v];
count = count + 1 ;
end
end
end
end
elseif mode3rd == 2
for t = 1:Ne3
ind_mat_3rd(t,:) = [t t t];
end
elseif mode3rd == 3
for t = 1:Ne3
for u = t:Ne3
ind_mat_3rd(count,:) = [t t u];
if t ~= u
count = count + 1;
ind_mat_3rd(count,:) = [t u u];
end
count = count + 1;
end
end
elseif mode3rd == 4
for t = 1:Ne3
for u = t:Ne3
for v = u:Ne3
if u-t<=len_3rd && v-t<=len_3rd
ind_mat_3rd(count,:) = [t u v];
count = count + 1 ;
end
end
end
end
elseif mode3rd == 5
for t = 1:Ne3
for u = t:Ne3
if u-t<=len_3rd
ind_mat_3rd(count,:) = [t t u];
if t ~= u
count = count + 1;
ind_mat_3rd(count,:) = [t u u];
end
count = count + 1;
end
end
end
ind_mat_3rd2 = NaN(N3,3);
count = 1;
% for t = 1:Ne3
% ind_mat_3rd2(count,:) = [t t t];
% count = count + 1;
% end
for t = 1:Ne3
% ind_mat_3rd2(count,:) = [t t t];
% count = count + 1;
for u = t:min(Ne3,t+len_3rd)
for v = unique([t u])
ind_mat_3rd2(count,:) = [t v u];
count = count + 1;
% ind_mat_3rd2(count,:) = [t u u];
% count = count + 1;
end
end
end
end
else
ind_mat_3rd = 0;
end
end
function [ind_mat_2nd,ind_mat_3rd] = calc_DFE_ind(obj,Ne2,N2,Ne3,N3)
if Ne2 > 0
ind_mat_2nd = NaN(N2,2);
count=1;
for t = 1:Ne2
for u = t:Ne2
ind_mat_2nd(count,:) = [t u];
count = count + 1 ;
end
end
else
ind_mat_2nd = 0;
end
if Ne3 > 0
ind_mat_3rd = NaN(N3,3);
count=1;
for t = 1:Ne3
for u = t:Ne3
for v = u:Ne3
ind_mat_3rd(count,:) = [t u v];
count = count + 1 ;
end
end
end
else
ind_mat_3rd = 0;
end
end
end
end

369
Classes/04_DSP/EQ_silas.m Normal file
View File

@@ -0,0 +1,369 @@
classdef EQ_silas < handle
%EQ_SILAS FFE and DFE Equalizer Playground
properties
% Important Signals
x_in %Input Sequence to be equalized
x_length
x_norm
d %reference signal
d_norm
d_constellation %constellation points of the reference
y_out %equalizer output signal
% FFE coefficients always named with "e"
Ne
Ce %memory length FFE
Ie1 %Indice Combination of 1nd order FFE
Ie2 %Indice Combination of 2nd order FFE
Ie3 %Indice Combination of 3nd order FFE
e %coefficients for FFE
% DFE coefficients always named with "b"
Nb
Cb %memory length DFE
Ib1 %Indice Combination of 1nd order DFE
Ib2 %Indice Combination of 2nd order DFE
Ib3 %Indice Combination of 3nd order DFE
b %coefficients for DFE
error
e_ffe
e_dfe
e_dc
error_log
mu_dc_train
mu_ffe_train
mu_dfe_train
mu_dc_dd
mu_combined_dd
delay
trainlength
sps
trainloops
ddloops
end
methods
function obj = EQ_silas(options)
%EQ_SILAS Construct an instance of this class
arguments(Input)
options.Ne = [50 5 0] %Number of FFE coefficients (1st, 2nd and 3rd order)
options.Nb = [30 5 3] %Number of DFE coefficients (1st, 2nd and 3rd order)
options.trainloops = 2;
options.trainlength = 4096;
options.ddloops = 2;
options.delay = 0;
options.sps = 2;
options.mu_dc_train = 0.01;
options.mu_ffe_train = 0.005;
options.mu_dfe_train = 0.005;
options.mu_dc_dd = 0.01;
options.mu_combined_dd = [0.0004 0.0005 0.0006 0.0007 ];
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
% Generate helpful vectors and initialize the filters with
% correct length:
obj.Ce = obj.calcVNLEMemoryLength(obj.Ne);
[obj.Ie2,obj.Ie3] = obj.calcIndiceVectors(obj.Ne);
obj.e = zeros(sum(obj.Ce),1);
obj.Cb = obj.calcVNLEMemoryLength(obj.Nb);
[obj.Ib2,obj.Ib3] = obj.calcIndiceVectors(obj.Nb);
obj.b = zeros(sum(obj.Cb),1);
end
function [signalclass_out,error_log] = process(obj,signalclass_in, reference_signalclass_in)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
signalclass_in = signalclass_in.normalize("mode","rms");
% Process the EQ optimization
obj.process_(signalclass_in.signal', reference_signalclass_in.signal');
signalclass_in.signal = obj.y_out';
% append to logbook
lbdesc = ['EQ von Silas ist gelaufen '];
signalclass_in = signalclass_in.logbookentry(lbdesc);
% write to output
signalclass_out = signalclass_in;
end
function process_(obj,x_in,d_in)
% 1) prepare signals
obj.e_dc = mean(x_in);
% 1.1) Input Signal
obj.x_in = [zeros(1,floor(obj.Ne(1)/2)) x_in zeros(1,obj.Ne(1))];
obj.x_length = length(x_in);
obj.x_norm = obj.calcPowerNormalization(x_in);
% 1.2 Reference Signal // Constellation
obj.d = [zeros(1,obj.Nb(1)-1) d_in zeros(1,obj.Nb(1))];
obj.d_constellation = unique(d_in);
obj.d_norm = obj.calcPowerNormalization(d_in);
% 1.3 Training
obj.trainingMode();
% 1.4 Decision Directed Mode
obj.decisionDirectedMode();
end
%% Adaptive Equalization Modes
function trainingMode(obj)
for tloop = 1:obj.trainloops
m = 1+obj.delay;
for n = obj.sps*obj.delay+1:obj.sps:obj.sps*obj.trainlength
m = m+1;
%get Sigal input vectors with correct length for VNLE
x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).';
x_in_vnle_format = obj.calcVNLENonlinVecs(x_in_block,obj.Ie2,obj.Ie3,obj.Ne,[1,1,1]);
%get Reference input vectors with correct length for VNLE
d_block = obj.d(obj.Nb(1)-obj.delay+m-2:-1:m-obj.delay-1).';
d_vnle_format = obj.calcVNLENonlinVecs(d_block,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm);
% Calculate the Error
obj.e_ffe = obj.e.' * x_in_vnle_format;
obj.e_dfe = obj.b.' * d_vnle_format;
obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay);
%update FFE coefficients with LMS
obj.e = obj.e - obj.error*conj(x_in_vnle_format)*obj.mu_ffe_train;
%update DFE coefficients with LMS
obj.b = obj.b + obj.mu_dfe_train*obj.error*d_vnle_format;
%update DC error
obj.e_dc = obj.e_dc - obj.error .* obj.mu_dc_train;
end
end
end
function decisionDirectedMode(obj)
%start the dd mode with coefficients from training
coeff = [obj.e;obj.b];
for ddloop = 1:obj.ddloops
m = 0;
if all(obj.mu_combined_dd == obj.mu_combined_dd(1))
mu_mat = obj.mu_combined_dd(1);
else
mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_combined_dd(1)... %1st order ffe
ones(1,obj.Ce(2))*obj.mu_combined_dd(2)... %2nd order ffe
ones(1,obj.Ce(3))*obj.mu_combined_dd(3)... %3rd order ffe
ones(1,sum(obj.Cb))*obj.mu_combined_dd(4)]); %all order dfe
end
y = zeros(1,floor(obj.x_length/obj.sps));
d_feedback = zeros(obj.Cb(1),1);
d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm);
d_hat = zeros(obj.x_length,1);
for k = 1:obj.sps:obj.x_length
m=m+1;
%get Sigal input vectors with correct length for VNLE
x = obj.x_in(obj.Ne(1)+k-1:-1:k).';
x_vnle = obj.calcVNLENonlinVecs(x,obj.Ie2,obj.Ie3,obj.Ne,[1,1,1]);
%combine FFE with DFE to one vector (cursor between the two sequences)
x_d = [x_vnle;-d_vnle];
%Apply filter
y(m) = obj.e_dc + x_d.'* coeff;
%Decision
[~,symbol_idx] = min(abs(y(m) - obj.d_constellation)); % decision for closest constellation point
d_hat(k) = obj.d_constellation(symbol_idx);
%Error between FFE & DFE filtered signal and Decision
obj.error = y(m) - d_hat(k);
%Update coefficients (both FFE and DFE)
coeff = coeff - mu_mat*obj.error*conj(x_d);
if 1 %mu_mat ~= 0
obj.e_dc = obj.e_dc - obj.mu_dc_dd * obj.error;
obj.error_log(ddloop,m) = obj.e_dc.^2;
end
% Append new decision to decision feedback
if obj.Nb(1) > 0
%shift up one index
d_feedback(2:end) = d_feedback(1:end-1);
%replace 1st index with current estimation
d_feedback(1) = d_hat(k);
%build memorylike VNLE version
d_vnle = obj.calcVNLENonlinVecs(d_feedback,obj.Ib2,obj.Ib3,obj.Nb,obj.d_norm);
end
end
end
obj.y_out = (circshift( y.' ,-(obj.delay))).';
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)
x1 = x_in_block;
x2 = [];
x3 = [];
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
x2 = input_vec_se(I_2(:,1)).*input_vec_se(I_2(:,2));
end
if N_(3) > 0
delta_3 = round((N_(1)-N_(3))/2);
input_vec_th = x_in_block(delta_3:end)/norm_(3);
x3 = input_vec_th(I_3(:,1)).*input_vec_th(I_3(:,2)).*input_vec_th(I_3(:,3));
end
x_in_vnle_format = [x1;x2;x3];
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
for order = 2:numel(N)
n = N(order);
v = 1:n; % Ursprünglicher Vektor
row = 1;
% Schleifen zur Generierung des Indize Vektors
switch order
case 2
indvec2nd = zeros(n*(n+1)/2, order);
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

25
Functions/calc_ber.m
View File

@@ -1,13 +1,30 @@
function [bits,errors,BER] = calc_ber(data_in,data_ref,skip) function [bits,errors,ber,errorLocation] = calc_ber(data_in,data_ref,options)
arguments(Input)
data_in;
data_ref;
options.skip = 0;
options.returnErrorLocation = 0;
end
bits = 0;
errors= 0;
ber= 0;
errorLocation= [];
data_ref = logical(data_ref)'; data_ref = logical(data_ref)';
data_in = logical(data_in)'; data_in = logical(data_in)';
bits = 0; bits = 0;
bits = bits+numel(data_in(:,skip+1:end)); bits = bits+numel(data_in(:,options.skip+1:end));
try try
errors = sum( data_in(:,skip+1:end) ~= data_ref(:,skip+1:end),"all" ); if options.returnErrorLocation == 0
errors = sum( data_in(:,options.skip+1:end) ~= data_ref(:,options.skip+1:end),"all" );
else
errorLocation = sum(data_in(:,options.skip+1:end) ~= data_ref(:,options.skip+1:end),1);
errors = sum(errorLocation ,"all" );
end
catch catch
%warning('BER calculation not optimal: Arrays have incompatible sizes for this operation.') %warning('BER calculation not optimal: Arrays have incompatible sizes for this operation.')
@@ -15,6 +32,6 @@ catch
end end
% Determine BER % Determine BER
BER = sum(errors)/sum(bits); ber = sum(errors)/sum(bits);
end end

215
setup_simulation.m
View File

@@ -1,23 +1,24 @@
sir_loop = -16; sir_loop = -22;
lw = [0.1e6 1e6 10e6];
for lp = 1 for lp = 1
dc_mu = 0.05; dc_mu = 0.05;
sir = sir_loop(lp); sir = sir_loop(1);
rng(10); rng(lp);
%% Set Simulation Variables %% Set Simulation Variables
sir = sir; sir = sir;
delay = 40; %mpi delay in meter delay = 20; %mpi delay in meter
fiblen = 0; %main link in km fiblen = 0; %main link in km
laser_linewidth = 1e6; laser_linewidth =1e6;
O = 17; %order of prbs O = 17; %order of prbs
N = 2^(O-1); %length of prbs N = 2^(O-1); %length of prbs
@@ -62,7 +63,7 @@ for lp = 1
reflectionpoint = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",sir); reflectionpoint = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",sir);
reflectionprop = Fiber("fsimu",fdac*kover,"fiber_length",2*delay/1000,"alpha",0.2,"D",16,"lambda0",1550,"gamma",0); reflectionprop = Fiber("fsimu",fdac*kover,"fiber_length",2*delay/1000,"alpha",0.2,"D",16,"lambda0",1550,"gamma",0);
opticatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",0); opticatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",-8);
phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20); phdiode = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20);
@@ -74,11 +75,13 @@ for lp = 1
"adcresolution",6,"quantbuffer",0.1,'block_dc',1); "adcresolution",6,"quantbuffer",0.1,'block_dc',1);
eq = EQ("K",2,"plottrain",0,"plotfinal",0,... eq = EQ("K",2,"plottrain",0,"plotfinal",0,...
"training_length",2048,"training_loops",5,... "training_length",4096,"training_loops",2,...
"Ne",[50,0,0],"Nb",[0,0,0],... "Ne",[50,5,3],"Nb",[3,2,2],...
"DCmu",dc_mu,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0.005,... "DCmu",dc_mu,"DDmu",[0.0004 0.0005 0.0006 0.0007 ],"DFEmu",0.005,"FFEmu",0.005,...
"dd_loops",2,"epsilon",[10 100 1000 ],"M",2,... "dd_loops",2,"epsilon",[10 100 1000 ],"M",2,...
"thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",1,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0);
eq2 = EQ_silas("Ne",[50,5,0],"Nb",[3,2,0],"trainlength",4096,"mu_dc_dd",0.001,"mu_dc_train",0.05,"mu_ffe_train",0.005,"mu_combined_dd",[0.0004 0.0006 0.0003 0.005],"ddloops",5);
%% PROCESS %% PROCESS
@@ -87,8 +90,6 @@ for lp = 1
for i = 1:log2(M) for i = 1:log2(M)
[bitpattern(:,i),seed] = prbs(O,N,seed); [bitpattern(:,i),seed] = prbs(O,N,seed);
end end
bitpattern = [bitpattern ; flip(bitpattern)];
%bitpattern = prms_out';
% Build Inf. signal class % Build Inf. signal class
bits = Informationsignal(bitpattern); bits = Informationsignal(bitpattern);
@@ -103,9 +104,10 @@ for lp = 1
X = pulseform.process(mod_out); X = pulseform.process(mod_out);
% AWG -> ELECTRICAL DOMAIN % AWG -> ELECTRICAL DOMAIN
awg_out = awg.process(X); X = awg.process(X);
X = lp_laser.process(X);
X = lp_laser.process(awg_out);
X = lp_laser.process(X); X = lp_laser.process(X);
X = X.normalize("mode","oneone"); X = X.normalize("mode","oneone");
@@ -114,22 +116,26 @@ for lp = 1
% Laser; Modulation -> OPTICAL DOMAIN % Laser; Modulation -> OPTICAL DOMAIN
[X,extmodlaser] = extmodlaser.process(X); [X,extmodlaser] = extmodlaser.process(X);
% Fiber Propagation % Fiber Propagation
X = fib.process(X); X = fib.process(X);
%% Reflect with attenuation %% Reflect with attenuation
R = reflectionpoint.process(X); R = reflectionpoint.process(X);
% Propagate % Propagate back and forth
R = reflectionprop.process(R); R = reflectionprop.process(R);
% disp(['SIR ',num2str(10*log10(X.power/R.power))]);
% Delay % Delay
[R,n] = R.delay("delay_meter",delay); [R,n] = R.delay("delay_meter",delay);
% Add together % Add together
X.signal = X.signal(n:end);
R.signal = R.signal(n:end);
disp(['SIR ',num2str(10*log10(X.power/R.power))]);
X = X+R; X = X+R;
%% %%
@@ -140,7 +146,6 @@ for lp = 1
X = phdiode.process(X); X = phdiode.process(X);
X = lp_diode.process(X); X = lp_diode.process(X);
% Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) % Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation)
X = scp.process(X); X = scp.process(X);
@@ -151,7 +156,15 @@ for lp = 1
Eq_in = X.normalize("mode","rms"); Eq_in = X.normalize("mode","rms");
% Equalizer % Equalizer
Eq_out = eq.process(Eq_in,reference); reference.signal = reference.signal(round(n * fsym/fsimu) : end,:);
tracking_speed_bandwidth = eq.FFEmu .* fsym .* 1e-9;
%% MPI reduction DC removal BEFORE EQ
% wl = 1000; % symbols
% Eq_in.signal = Eq_in.signal - 1/wl .* movsum( Eq_in.signal,[wl/2,wl/2]);
[Eq_out] = eq2.process(Eq_in,reference);
%% MPI reduction DC removal %% MPI reduction DC removal
wl = 1000; % symbols wl = 1000; % symbols
@@ -178,21 +191,12 @@ for lp = 1
yk_lvsm.signal(pre_decision_level_uni==level) = yk_lvsm.signal(pre_decision_level_uni==level) - smoothed(pre_decision_level_uni==level); yk_lvsm.signal(pre_decision_level_uni==level) = yk_lvsm.signal(pre_decision_level_uni==level) - smoothed(pre_decision_level_uni==level);
yk_lvlp.signal(pre_decision_level_uni==level) = yk_lvlp.signal(pre_decision_level_uni==level) - filtered(pre_decision_level_uni==level); yk_lvlp.signal(pre_decision_level_uni==level) = yk_lvlp.signal(pre_decision_level_uni==level) - filtered(pre_decision_level_uni==level);
end end
%
% figure(26);plot(1:length(error_log),error_log(:,end),"LineWidth",0.4);ylim([-0.8 0.8]);ylabel('$\epsilon$'),title(['Linewidth: ',num2str(laser_linewidth*1e-6), 'MHz'])
%
% figure(28);plot(1:length(e),e,"LineWidth",0.4);ylim([-0.8 0.8]);
%% MPI Reduction Tunable Notch Filter
% yk = Eq_in.normalize("mode","oneone");
% yk = yk.signal;
% dk = Eq_out.resample("fs_in",fsym,"fs_out",2*fsym);
% dk = ( digimod.decide_pamlevel(dk) *2-3 ) .* 1/sqrt(5);
% dk = dk ./ max(abs(dk));
%
% xk = yk(4:end)-dk(1:end-3);
%
% ya = hilbert(xk);
% figure;periodogram(abs(ya),[],length(ya),2*fsym,'centered')
% plot(abs(ya));
% grid on;
%% Calc EVM %% Calc EVM
@@ -201,31 +205,7 @@ for lp = 1
evm_lsm = calc_evm(yk_lvsm.signal, pre_decision_level_bi); evm_lsm = calc_evm(yk_lvsm.signal, pre_decision_level_bi);
evm_llp = calc_evm(yk_lvlp.signal, pre_decision_level_bi); evm_llp = calc_evm(yk_lvlp.signal, pre_decision_level_bi);
%figure(1);bar([evm_bm' evm_dcsm' evm_llp' evm_lsm']);ylim([0.01 0.1]);set(gca,'yscale','log'); % figure(1);bar([evm_bm' evm_dcsm' evm_llp' evm_lsm']);ylim([0.01 0.1]);set(gca,'yscale','log');
%% Plot stuff
% figure(200)
% hold on
% for i = 1:4
% subplot(4,1,i)
% scatter(1:length(err_),err_(:,i),4,'Marker','.','MarkerEdgeColor',col(i,:));
% end
%
% figure(12);
% hold on;
% scatter(1:length(yk_sm.signal),yk_sm.signal,4,'Marker','.','MarkerEdgeColor',col(6,:));
% for i = 1:4
% x = err_(:,i);
% nanx = isnan(x);
% t = 1:numel(x);
% x(nanx) = interp1(t(~nanx), x(~nanx), t(nanx));
%
% X(i,:) = x;
%
% plot(x(2000:end-2000,1),'LineWidth',2)
% end
%% MPI reduction A3
% Digi Demod % Digi Demod
d_bm = digimod.demap(Eq_out); d_bm = digimod.demap(Eq_out);
@@ -234,21 +214,85 @@ for lp = 1
d_lvlp = digimod.demap(yk_lvlp); d_lvlp = digimod.demap(yk_lvlp);
% BER % BER
bitpattern = bitpattern(round(n * fsym/fsimu):end,:);
dbit = length(d_bm.signal)-length(bitpattern); dbit = length(d_bm.signal)-length(bitpattern);
[~,errors_bm,ber_bm] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); [~,errors_bm,ber_bm,loc] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",1);
[~,errors_dcsm,ber_dcsm] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); [~,errors_dcsm,ber_dcsm] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0);
[~,errors_lvsm,ber_lvsm] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); [~,errors_lvsm,ber_lvsm] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0);
[~,errors_lvlp,ber_lvlp] = calc_ber(d_lvlp.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); [~,errors_lvlp,ber_lvlp] = calc_ber(d_lvlp.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0);
% Display % Display
disp(['BER benchmark: ', sprintf('%2E',ber_bm), ' ERRORS: ' ,num2str(sum(errors_bm))]); disp(['BER benchmark: ', sprintf('%2E',ber_bm), ' ERRORS: ' ,num2str(sum(errors_bm))]);
disp(['BER dc smooth: ', sprintf('%2E',ber_dcsm), ' ERRORS: ' ,num2str(sum(errors_dcsm))]);
disp(['BER lv smooth: ', sprintf('%2E',ber_lvsm), ' ERRORS: ' ,num2str(sum(errors_lvsm))]); % disp(['BER dc smooth: ', sprintf('%2E',ber_dcsm), ' ERRORS: ' ,num2str(sum(errors_dcsm))]);
disp(['BER lv lowpas: ', sprintf('%2E',ber_lvlp), ' ERRORS: ' ,num2str(sum(errors_lvlp))]); % disp(['BER lv smooth: ', sprintf('%2E',ber_lvsm), ' ERRORS: ' ,num2str(sum(errors_lvsm))]);
% disp(['BER lv lowpas: ', sprintf('%2E',ber_lvlp), ' ERRORS: ' ,num2str(sum(errors_lvlp))]);
%% Generate some Plots %% Generate some Plots
if 1
% SCATTER
col = cbrewer2('Paired',8);
xax = 1:Eq_out.length;
figure(3)
sgtitle('')
subplot(1,2,1)
hold on
eq_decision = digimod.decide_pamlevel(Eq_out);
true_symbols = digimod.decide_pamlevel(mod_out);
xindices = 1:Eq_in.length;
errorpos = find(loc~=0)*2;
correct = find(loc==0)*2;
scatter(xindices(correct),Eq_in.signal(correct),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After EQ');
hold on
scatter(xindices(errorpos),Eq_in.signal(errorpos),6,'x','MarkerEdgeColor',col(6,:),'DisplayName','Wrong Decision');
hold off
xlim([1, xindices(end)]);
ylim([-3 3]);
xlabel('Sampling Index')
ylabel('Amplitude')
a = legend;
a.Location = "best";
subplot(1,2,2)
xindices = 1:Eq_out.length;
scatter(xindices(loc==0),Eq_out.signal(loc==0),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After EQ');
hold on
scatter(xindices(loc~=0),Eq_out.signal(loc~=0),8,'x','MarkerEdgeColor',col(6,:),'DisplayName','Wrong Decision');
hold off
xlim([1, xax(end)]);
ylim([-2 2]);
xlabel('Sampling Index')
ylabel('Amplitude')
legend
a = legend;
a.Location = "best";
hold off
% ERROR IN EQ
figure(13)
plot(diff(eq2.error_log(1,:)));
hold on;
for i = 1:size(eq2.error_log,1)
plot(diff(eq2.error_log(i,:)));
end
loc(loc==0) = NaN;
stem(loc.*mean((eq2.error_log(end,:))))
hold off
end
if 0 if 0
col = cbrewer2('Paired',8); col = cbrewer2('Paired',8);
@@ -258,7 +302,7 @@ for lp = 1
subplot(2,1,1) subplot(2,1,1)
hold on hold on
plot(reference.signal(4150:4175),'DisplayName','Tx','Color',col(1,:),'LineWidth',3); plot(reference.signal(4150:4175),'DisplayName','Tx','Color',col(1,:),'LineWidth',3);
plot(eq_out.signal(4150:4175),'DisplayName','Rx after EQ','Color',col(6,:),'LineWidth',1); plot(Eq_out.signal(4150:4175),'DisplayName','Rx after EQ','Color',col(6,:),'LineWidth',1);
title('Modulated Sequence Zoom'); title('Modulated Sequence Zoom');
legend legend
hold off hold off
@@ -271,45 +315,6 @@ for lp = 1
hold off hold off
end end
if 0
col = cbrewer2('Paired',8);
xax = 1:Eq_out.length;
figure(3)
sgtitle('')
% subplot(1,4,1:2)
% scatter(1:4:X.length,X.signal(1:4:end),4,'.','MarkerEdgeColor',col(6,:),'DisplayName','Before EQ');
% xlim([1, xax(end)]);
% %ylim([-2 2]);
% xlabel('Sampling Index')
% ylabel('Amplitude')
% legend
subplot(1,3,1)
scatter(1:Eq_in.length,Eq_in.signal,4,'.','MarkerEdgeColor',col(6,:),'DisplayName','After EQ');
xlim([1, xax(end)]);
ylim([-2 2]);
xlabel('Sampling Index')
ylabel('Amplitude')
legend
subplot(1,3,2)
scatter(xax,Eq_out.signal,4,'.','MarkerEdgeColor',col(6,:),'DisplayName','After EQ');
xlim([1, xax(end)]);
ylim([-2 2]);
xlabel('Sampling Index')
ylabel('Amplitude')
legend
subplot(1,3,3)
scatter(xax,yk_lvsm.signal,4,'.','MarkerEdgeColor',col(6,:),'DisplayName','After EQ');
xlim([1, xax(end)]);
ylim([-2 2]);
xlabel('Sampling Index')
ylabel('Amplitude')
legend
end
end end