From cc4a66adae3a4afc46d4ef74f37c2dea2c987930 Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Fri, 6 Oct 2023 15:19:09 +0200 Subject: [PATCH] EQ update vom Arbeits PC --- Classes/04_DSP/EQ_silas.m | 26 +- Classes/04_DSP/EQ_silas_plain.m | 369 +++++++++++++++++++++++ MPI_offline_analysis/mpi_analysis_pam4.m | 133 +++++--- 3 files changed, 470 insertions(+), 58 deletions(-) create mode 100644 Classes/04_DSP/EQ_silas_plain.m diff --git a/Classes/04_DSP/EQ_silas.m b/Classes/04_DSP/EQ_silas.m index cb33343..d32d875 100644 --- a/Classes/04_DSP/EQ_silas.m +++ b/Classes/04_DSP/EQ_silas.m @@ -76,7 +76,7 @@ classdef EQ_silas < handle options.mu_dc_dd = 0.01; options.mu_combined_dd = [0.0004 0.0005 0.0006 0.0007 ]; - options.dcmode + options.dcmode = 1; end fn = fieldnames(options); @@ -207,18 +207,18 @@ classdef EQ_silas < handle 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 - mu_ffe = [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)]; - mu_dfe = [ones(1,sum(obj.Cb))*obj.mu_combined_dd(4)]; - end + + 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 + + mu_ffe = [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)]; + + mu_dfe = ones(1,sum(obj.Cb))*obj.mu_combined_dd(4); + y = zeros(1,floor(obj.x_length/obj.sps)); d_feedback = zeros(obj.Cb(1),1); diff --git a/Classes/04_DSP/EQ_silas_plain.m b/Classes/04_DSP/EQ_silas_plain.m new file mode 100644 index 0000000..13ccb44 --- /dev/null +++ b/Classes/04_DSP/EQ_silas_plain.m @@ -0,0 +1,369 @@ +classdef EQ_silas_plain < 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_plain(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/MPI_offline_analysis/mpi_analysis_pam4.m b/MPI_offline_analysis/mpi_analysis_pam4.m index bce1b23..bac262e 100644 --- a/MPI_offline_analysis/mpi_analysis_pam4.m +++ b/MPI_offline_analysis/mpi_analysis_pam4.m @@ -1,64 +1,107 @@ clear; col = linspecer(6); -% load dataset -simdata = load('/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km/pam4__loop_10_92Gbd_19092023_1508.mat'); -% get common variables of simulation loop -fdac = simdata.saveStructTemp.common.f_DAC; -fadc = simdata.saveStructTemp.common.f_ADC; -fsym = simdata.saveStructTemp.common.f_sym; +foldername = 'C:\Users\Silas\Documents\MATLAB\Labor_Datensatz_PAM4_MPI\pam4_10km_1km'; -% 0) build RX Signal +allfiles = dir(foldername); -y_rx = Electricalsignal(simdata.saveStructTemp.dp_tsynch_out'); -y_rx.fs = 2.*fsym; - -% 0) build tx reference Signal for eq training -y_digimod = Electricalsignal(simdata.saveStructTemp.digi_mod_out'); -y_digimod.fs = fsym; +parfsym = zeros(1,length(allfiles)); +parsir = zeros(1,length(allfiles)); +ber_wo_reduction = zeros(1,length(allfiles)); +ber_w_dcremoval = zeros(1,length(allfiles)); +mudc = zeros(1,length(allfiles)); -% 1) normlaize -y_rx = y_rx.normalize("mode","rms"); +parfor i = 1:length(allfiles) + + if allfiles(i).bytes ~= 0 + current_filename = allfiles(i).name; + else + continue + end + + % load dataset + simdata = load([foldername,filesep, current_filename]); + % get common variables of simulation loop + + %even if redundant, save stuff to results struct + parsir(i) = simdata.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2; + parfsym(i) = simdata.saveStructTemp.common.f_sym; + + fdac = simdata.saveStructTemp.common.f_DAC; + fadc = simdata.saveStructTemp.common.f_ADC; + fsym = simdata.saveStructTemp.common.f_sym; + + % 0) build RX Signal + + y_rx = Electricalsignal(simdata.saveStructTemp.dp_tsynch_out'); + + y_rx.fs = 2.*fsym; + + % 0) build tx reference Signal for eq training + y_digimod = Electricalsignal(simdata.saveStructTemp.digi_mod_out'); + y_digimod.fs = fsym; -% 2.1) equalize -eq = EQ("K",2,"plottrain",0,"plotfinal",0,... - "training_length",4096,"training_loops",3,... - "Ne",[50,9,3],"Nb",[2,0,0],... - "DCmu",0,"DDmu",[0.0004 0.0005 0.0005 0.0005 ],"DFEmu",0.005,"FFEmu",0.005,... - "dd_loops",3,"epsilon",[10 100 1000 ],"M",2,... - "thres",[0.005 0.0004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); - -eq2 = EQ_silas("Ne",[50,9,9],"Nb",[0,0,0],"trainlength",4096,"mu_dc_dd",0.05,"mu_dc_train",0.05,"mu_ffe_train",0.005,"mu_combined_dd",[0.0004 0.0005 0.0005 0.0005 ],"ddloops",3,"trainloops",3); + % 1) normlaize + y_rx = y_rx.normalize("mode","rms"); -[Eq_out] = eq.process(y_rx,y_digimod); -[Eq_out2] = eq2.process(y_rx,y_digimod); - -figure() -scatter(1:length(y_rx.signal),y_rx.signal,4,'.','MarkerEdgeColor',col(2,:),'DisplayName','After EQ'); -hold on -%scatter(1:length(Eq_out.signal),Eq_out.signal,4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After EQ'); -hold on -scatter(1:length(Eq_out2.signal),Eq_out2.signal,3,'.','MarkerEdgeColor',col(3,:),'DisplayName','After EQ'); + % 2.1) equalize + eq = EQ("K",2,"plottrain",0,"plotfinal",0,... + "training_length",4096,"training_loops",4,... + "Ne",[50,9,9],"Nb",[2,0,0],... + "DCmu",0.05,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,... + "dd_loops",4,"epsilon",[10 100 1000 ],"M",2,... + "thres",[0.005 0.004 0.005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); -% 2.2) MPI reduction blocks + [Eq_out] = eq.process(y_rx,y_digimod); -% 3) digital demodulation -digimod = PAMmapper(2^simdata.saveStructTemp.common.M,0); -d_hat = digimod.demap(Eq_out); -d_hat2 = digimod.demap(Eq_out2); -d_ = digimod.demap(Electricalsignal(simdata.saveStructTemp.digi_mod_out')); -d = simdata.saveStructTemp.prms_out; + % 2.2) MPI reduction blocks + + % 3) digital demodulation + digimod = PAMmapper(2^simdata.saveStructTemp.common.M,0); + + d_hat = digimod.demap(Eq_out); + + d_ = digimod.demap(Electricalsignal(simdata.saveStructTemp.digi_mod_out')); + d = simdata.saveStructTemp.prms_out; + + % 4) BER calculation + [totalbits,errors_bm,ber_dhat,loc] = calc_ber(d_hat.signal(1:end,:) ,d(1:end,:)',"skip",0,"returnErrorLocation",1); + + ber_nodctap = mean(simdata.saveStructTemp.prms_compare_state.BER ); + + disp(['BER_new_eq1: ',num2str(ber_dhat),' BER_old: ',num2str(ber_nodctap)]); + + % 5) save to struct + ber_wo_reduction(i) = ber_nodctap; + ber_w_dcremoval(i) = ber_dhat; + mudc(i) = eq.DCmu; + +end + + +rate = [92e9, 56e9]; +figure(99) + +for r = 1:length(rate) + a = [parsir(parfsym==rate(r));ber_wo_reduction(parfsym==rate(r))]; + a=sort(a,2); + plot(a(1,:),a(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; w/o reduction'],'LineWidth',2,'Marker','o','MarkerSize',3); +%scatter(parsir(parfsym==rate(r)),ber_wo_reduction(parfsym==rate(r)),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; w/o reduction'],'LineWidth',2,'Marker','o') + hold on + b = [parsir(parfsym==rate(r));ber_w_dcremoval(parfsym==rate(r))]; + b=sort(b,2); + plot(b(1,:),b(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; adaptive DC removal'],'LineWidth',2,'Marker','+','MarkerSize',3); +%scatter(parsir(parfsym==rate(r)),ber_w_dcremoval(parfsym==rate(r)),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; adaptive DC removal'],'LineWidth',2,'Marker','+'); +end +set(gca, 'YScale', 'log') +yline(3.8e-3,'LineWidth',1.5); + -% 4) BER calculation -[totalbits,errors_bm,ber_dhat,loc] = calc_ber(d_hat.signal(1:end,:) ,d(1:end,:)',"skip",0,"returnErrorLocation",1); -[totalbits2,errors_bm2,ber_dhat2,loc] = calc_ber(d_hat2.signal(1:end,:) ,d(1:end,:)',"skip",0,"returnErrorLocation",1); -ber_nodctap = mean(simdata.saveStructTemp.prms_compare_state.BER ); -disp(['BER_new_eq1: ',num2str(ber_dhat),' BER_new_eq2: ',num2str(ber_dhat2),' BER_old: ',num2str(ber_nodctap)]);