EQ update vom Arbeits PC

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);
+                mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_combined_dd(1)... %1st order ffe
-                else
+                               ones(1,obj.Ce(2))*obj.mu_combined_dd(2)... %2nd order ffe
-                    mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_combined_dd(1)... %1st order ffe
+                               ones(1,obj.Ce(3))*obj.mu_combined_dd(3)... %3rd order ffe
-                                   ones(1,obj.Ce(2))*obj.mu_combined_dd(2)... %2nd order ffe
+                             ones(1,sum(obj.Cb))*obj.mu_combined_dd(4)]); %all order dfe
-                                   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
-                    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(2))*obj.mu_combined_dd(2)... %2nd order ffe
+                               ones(1,obj.Ce(3))*obj.mu_combined_dd(3)];
-                                   ones(1,obj.Ce(3))*obj.mu_combined_dd(3)];
+                
-                    mu_dfe = [ones(1,sum(obj.Cb))*obj.mu_combined_dd(4)];
+                mu_dfe = ones(1,sum(obj.Cb))*obj.mu_combined_dd(4);
-                end
+
 
                 y               = zeros(1,floor(obj.x_length/obj.sps));
                 d_feedback      = zeros(obj.Cb(1),1);

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

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;
+foldername = 'C:\Users\Silas\Documents\MATLAB\Labor_Datensatz_PAM4_MPI\pam4_10km_1km';
-fadc = simdata.saveStructTemp.common.f_ADC;
-fsym = simdata.saveStructTemp.common.f_sym;
 
-% 0) build RX Signal
+allfiles = dir(foldername);
 
-y_rx = Electricalsignal(simdata.saveStructTemp.dp_tsynch_out');
+parfsym = zeros(1,length(allfiles));
-y_rx.fs = 2.*fsym;
+parsir = zeros(1,length(allfiles));
-
+ber_wo_reduction = zeros(1,length(allfiles));
-% 0) build tx reference Signal for eq training
+ber_w_dcremoval = zeros(1,length(allfiles));
-y_digimod = Electricalsignal(simdata.saveStructTemp.digi_mod_out');
+mudc = zeros(1,length(allfiles));
-y_digimod.fs = fsym;
 
 
-% 1) normlaize
+parfor i = 1:length(allfiles)
-y_rx = y_rx.normalize("mode","rms");
+
+    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
+    % 1) normlaize
-eq          = EQ("K",2,"plottrain",0,"plotfinal",0,...
+    y_rx = y_rx.normalize("mode","rms");
-    "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);
 
 
-[Eq_out] = eq.process(y_rx,y_digimod);
+    % 2.1) equalize
-[Eq_out2] = eq2.process(y_rx,y_digimod);
+    eq          = EQ("K",2,"plottrain",0,"plotfinal",0,...
-
+        "training_length",4096,"training_loops",4,...
-figure()
+        "Ne",[50,9,9],"Nb",[2,0,0],...
-scatter(1:length(y_rx.signal),y_rx.signal,4,'.','MarkerEdgeColor',col(2,:),'DisplayName','After EQ');
+        "DCmu",0.05,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,...
-hold on
+        "dd_loops",4,"epsilon",[10 100 1000 ],"M",2,...
-%scatter(1:length(Eq_out.signal),Eq_out.signal,4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After EQ');
+        "thres",[0.005 0.004 0.005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0);
-hold on
-scatter(1:length(Eq_out2.signal),Eq_out2.signal,3,'.','MarkerEdgeColor',col(3,:),'DisplayName','After EQ');
 
 
-% 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);
+    % 2.2) MPI reduction blocks
-d_hat2 = digimod.demap(Eq_out2);
+
-d_    = digimod.demap(Electricalsignal(simdata.saveStructTemp.digi_mod_out'));
+    % 3) digital demodulation
-d = simdata.saveStructTemp.prms_out;
+    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)]);