From 8b3bc688ddbd664d0980def8138e47b37b30d31d Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Fri, 14 Jul 2023 14:42:15 +0200 Subject: [PATCH] own EQ implemented and returns same performance like Toms --- Classes/00_signals/Signal.m | 22 +- Classes/01_transmit/PAMmapper.m | 1 - Classes/04_DSP/EQ.m | 86 ++-- Classes/04_DSP/EQ_copy.m | 727 ++++++++++++++++++++++++++++++++ Classes/04_DSP/EQ_silas.m | 369 ++++++++++++++++ Functions/calc_ber.m | 25 +- setup_simulation.m | 221 +++++----- 7 files changed, 1287 insertions(+), 164 deletions(-) create mode 100644 Classes/04_DSP/EQ_copy.m create mode 100644 Classes/04_DSP/EQ_silas.m diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 051c6ca..39be34e 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -291,28 +291,10 @@ classdef Signal 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 - 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 diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index e7356d0..50592c7 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -108,7 +108,6 @@ classdef PAMmapper end thres = thres .* 1/sqrt(5); - case 3 % 8-ASK if obj.unipolar==0 diff --git a/Classes/04_DSP/EQ.m b/Classes/04_DSP/EQ.m index e34d62a..3e2837d 100644 --- a/Classes/04_DSP/EQ.m +++ b/Classes/04_DSP/EQ.m @@ -102,10 +102,10 @@ classdef EQ 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.) - 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'; % append to logbook @@ -117,7 +117,7 @@ classdef EQ 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 % Detailed explanation goes here @@ -211,43 +211,62 @@ classdef EQ 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< - + 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 - % 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); - e_dc = mean(data_in); % initilaization of the dc tap with the mean value of the data - % e_save = NaN(361,8.6e5); - % save_ind = 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 - 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); + 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); - input_vec = [X_1;X_2;X_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 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; %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 else if obj.l1act @@ -258,21 +277,24 @@ classdef EQ e_ = e_ - error*input_vec/(input_vec.'*input_vec); 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 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 @@ -355,6 +377,7 @@ classdef EQ 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); @@ -384,9 +407,8 @@ classdef EQ for k = 1:obj.K:length(data_in) 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); @@ -403,6 +425,7 @@ classdef EQ 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); @@ -423,14 +446,15 @@ classdef EQ % save_ind = save_ind+1; if mu_mat ~= 0 - e_dc = e_dc - obj.DCmu*error; - error_log(end+1) = e_dc; + e_dc = e_dc - obj.DCmu*error; + error_log(cnt,dd_loop) = e_dc; + 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))).'; diff --git a/Classes/04_DSP/EQ_copy.m b/Classes/04_DSP/EQ_copy.m new file mode 100644 index 0000000..b2d68e7 --- /dev/null +++ b/Classes/04_DSP/EQ_copy.m @@ -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 + diff --git a/Classes/04_DSP/EQ_silas.m b/Classes/04_DSP/EQ_silas.m new file mode 100644 index 0000000..4ec87c0 --- /dev/null +++ b/Classes/04_DSP/EQ_silas.m @@ -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 + + + + + + + + + + + + + + + + + + + + + + diff --git a/Functions/calc_ber.m b/Functions/calc_ber.m index 2d79f38..6ed7bc3 100644 --- a/Functions/calc_ber.m +++ b/Functions/calc_ber.m @@ -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_in = logical(data_in)'; bits = 0; -bits = bits+numel(data_in(:,skip+1:end)); +bits = bits+numel(data_in(:,options.skip+1:end)); 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 %warning('BER calculation not optimal: Arrays have incompatible sizes for this operation.') @@ -15,6 +32,6 @@ catch end % Determine BER -BER = sum(errors)/sum(bits); +ber = sum(errors)/sum(bits); end \ No newline at end of file diff --git a/setup_simulation.m b/setup_simulation.m index a732aa9..2218f6f 100644 --- a/setup_simulation.m +++ b/setup_simulation.m @@ -1,23 +1,24 @@ -sir_loop = -16; +sir_loop = -22; +lw = [0.1e6 1e6 10e6]; for lp = 1 dc_mu = 0.05; - sir = sir_loop(lp); + sir = sir_loop(1); - rng(10); + rng(lp); %% Set Simulation Variables sir = sir; - delay = 40; %mpi delay in meter + delay = 20; %mpi delay in meter fiblen = 0; %main link in km - laser_linewidth = 1e6; + laser_linewidth =1e6; O = 17; %order 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); 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); @@ -74,11 +75,13 @@ for lp = 1 "adcresolution",6,"quantbuffer",0.1,'block_dc',1); eq = EQ("K",2,"plottrain",0,"plotfinal",0,... - "training_length",2048,"training_loops",5,... - "Ne",[50,0,0],"Nb",[0,0,0],... - "DCmu",dc_mu,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0.005,... + "training_length",4096,"training_loops",2,... + "Ne",[50,5,3],"Nb",[3,2,2],... + "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,... - "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 @@ -87,8 +90,6 @@ for lp = 1 for i = 1:log2(M) [bitpattern(:,i),seed] = prbs(O,N,seed); end - bitpattern = [bitpattern ; flip(bitpattern)]; - %bitpattern = prms_out'; % Build Inf. signal class bits = Informationsignal(bitpattern); @@ -103,9 +104,10 @@ for lp = 1 X = pulseform.process(mod_out); % 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 = X.normalize("mode","oneone"); @@ -114,22 +116,26 @@ for lp = 1 % Laser; Modulation -> OPTICAL DOMAIN [X,extmodlaser] = extmodlaser.process(X); - % Fiber Propagation X = fib.process(X); - + %% Reflect with attenuation R = reflectionpoint.process(X); - % Propagate + % Propagate back and forth R = reflectionprop.process(R); - % disp(['SIR ',num2str(10*log10(X.power/R.power))]); % Delay [R,n] = R.delay("delay_meter",delay); % 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; %% @@ -140,7 +146,6 @@ for lp = 1 X = phdiode.process(X); X = lp_diode.process(X); - % Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) X = scp.process(X); @@ -151,7 +156,15 @@ for lp = 1 Eq_in = X.normalize("mode","rms"); % 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 wl = 1000; % symbols @@ -178,22 +191,13 @@ 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_lvlp.signal(pre_decision_level_uni==level) = yk_lvlp.signal(pre_decision_level_uni==level) - filtered(pre_decision_level_uni==level); 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 evm_bm = calc_evm(Eq_out.signal, pre_decision_level_bi); @@ -201,31 +205,7 @@ for lp = 1 evm_lsm = calc_evm(yk_lvsm.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'); - - %% 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 - + % figure(1);bar([evm_bm' evm_dcsm' evm_llp' evm_lsm']);ylim([0.01 0.1]);set(gca,'yscale','log'); % Digi Demod d_bm = digimod.demap(Eq_out); @@ -234,21 +214,85 @@ for lp = 1 d_lvlp = digimod.demap(yk_lvlp); % BER - + bitpattern = bitpattern(round(n * fsym/fsimu):end,:); dbit = length(d_bm.signal)-length(bitpattern); - [~,errors_bm,ber_bm] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); - [~,errors_dcsm,ber_dcsm] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); - [~,errors_lvsm,ber_lvsm] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); - [~,errors_lvlp,ber_lvlp] = calc_ber(d_lvlp.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,:),"skip",0,"returnErrorLocation",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,:),"skip",0,"returnErrorLocation",0); % Display 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 lv lowpas: ', sprintf('%2E',ber_lvlp), ' ERRORS: ' ,num2str(sum(errors_lvlp))]); + + % 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 lv lowpas: ', sprintf('%2E',ber_lvlp), ' ERRORS: ' ,num2str(sum(errors_lvlp))]); %% 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 col = cbrewer2('Paired',8); @@ -258,7 +302,7 @@ for lp = 1 subplot(2,1,1) hold on 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'); legend hold off @@ -271,45 +315,6 @@ for lp = 1 hold off 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