DC removal algorithm now with latency

Classes/04_DSP/EQ_silas.m

@@ -33,7 +33,6 @@ classdef EQ_silas < handle
         e_ffe
         e_dfe
         e_dc
-        error_log
 
         % coefficients
         mu_dc_train
@@ -52,8 +51,8 @@ classdef EQ_silas < handle
         trainloops
         ddloops
 
-        dcmode
+        eq_blocklength % block lengt of EQ (until now, only the dc subtraction is affected by this)
-
+        eq_updatelatency % time in symbols until the calculated updates reach the signal again (until now, only the dc subtraction is affected by this)
 
 
 
@@ -80,7 +79,8 @@ classdef EQ_silas < handle
                 options.mu_ffe_dd =  [0.0004 0.0005 0.0006];
                 options.mu_dfe_dd =  0.0005;
 
-                options.dcmode = 1;
+                options.eq_blocklength = 1;
+                options.eq_updatelatency = 1;
             end
 
             fn = fieldnames(options);
@@ -106,7 +106,7 @@ classdef EQ_silas < handle
 
         end
 
-        function [signalclass_out,error_log] = process(obj,signalclass_in, reference_signalclass_in)
+        function [signalclass_out] = process(obj,signalclass_in, reference_signalclass_in)
 
             % actual processing of the signal (steps 1. - 3.)
             % 1 normalize RMS
@@ -151,12 +151,15 @@ classdef EQ_silas < handle
         %% Adaptive Equalization Modes
 
         function trainingMode(obj)
+            
+            dc_block = ones(obj.eq_blocklength,1);
 
             for tloop = 1:obj.trainloops
                 m = 1+obj.delay;
-
+                dc_cnt =  0;
                 for n = obj.sps*obj.delay+1:obj.sps:obj.sps*obj.trainlength
                     m = m+1;
+                    dc_cnt = dc_cnt+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).';
@@ -173,9 +176,7 @@ classdef EQ_silas < handle
 
                     % Calculate the Error
                     obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay);
-
+   
-                    err_track(n) = obj.error;
-
 
                     if obj.mu_ffe_train ~= 0
                         %update FFE coefficients with LMS
@@ -189,7 +190,12 @@ classdef EQ_silas < handle
                     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;
+                    dc_block(dc_cnt) = obj.error .* obj.mu_dc_train;
+                    
+                    if dc_cnt == obj.eq_blocklength
+                        obj.e_dc = obj.e_dc - mean(dc_block(dc_cnt));
+                        dc_cnt = 0;
+                    end
 
                 end
 
@@ -203,10 +209,12 @@ classdef EQ_silas < handle
             %start the dd mode with coefficients from training
             coeff = [obj.e;obj.b];
 
+            dc_block = ones(obj.eq_blocklength,1);
+
             for ddloop = 1:obj.ddloops
 
                 m = 0;
-
+                dc_cnt = 0;
 
                 mu_mat = diag([ones(1,obj.Ce(1))*obj.mu_ffe_dd(1)... %1st order ffe
                     ones(1,obj.Ce(2))*obj.mu_ffe_dd(2)... %2nd order ffe
@@ -226,7 +234,7 @@ classdef EQ_silas < handle
                 d_hat           = zeros(obj.x_length,1);
 
                 for k = 1:obj.sps:obj.x_length
-
+                    dc_cnt = dc_cnt+1;
                     m=m+1;
 
                     %get Sigal input vectors with correct length for VNLE
@@ -254,9 +262,14 @@ classdef EQ_silas < handle
 
                     coeff = coeff - mu_mat*obj.error*conj(x_d);
                    
-                    obj.e_dc = obj.e_dc - obj.mu_dc_dd * obj.error;
+                    
+                    %update DC error
+                    dc_block(dc_cnt) = obj.error .* obj.mu_dc_dd;
 
-                    obj.error_log(ddloop,m) = obj.e_dc.^2;
+                    if dc_cnt == obj.eq_blocklength
+                        obj.e_dc = obj.e_dc -  mean(dc_block(dc_cnt));
+                        dc_cnt = 0;
+                    end
 
                     % Append new decision to decision feedback
                     if obj.Nb(1) > 0

MPI_offline_analysis/mpi_analysis_pam4.m

@@ -2,115 +2,117 @@
 clear;
 col = linspecer(6);
 
-foldername = '/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km';
+if ismac
+    foldername = '/Users/silasoettinghaus/Documents/MATLAB/Labor_Datensatz_PAM4_MPI/pam4_10km_1km';
+else
+    foldername = 'C:\Users\Silas\Documents\MATLAB\Labor_Datensatz_PAM4_MPI\pam4_10km_1km';
+end
 
 allfiles = dir(foldername);
 
-parfsym = zeros(1,length(allfiles));
+mudc_loop = [0.05 0.01 0.005];
-parsir = zeros(1,length(allfiles));
+
-ber_wo_reduction = zeros(1,length(allfiles));
+parfsym             = zeros(length(mudc_loop),length(allfiles));
-ber_w_dcremoval = zeros(1,length(allfiles));
+parsir              = zeros(length(mudc_loop),length(allfiles));
-mudc = zeros(1,length(allfiles));
+ber_wo_reduction    = zeros(length(mudc_loop),length(allfiles));
+ber_w_dcremoval     = zeros(length(mudc_loop),length(allfiles));
+mudc                = zeros(length(mudc_loop),length(allfiles));
 
 
-for i = 1:length(allfiles)
+for j = 1:length(mudc_loop)
+    mudc_current = mudc_loop(j);
+    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(j,i)  = simdata.saveStructTemp.awg2scope_keysight_state.eigenlight_mpi - 3.3 + 7.2;
+        parfsym(j,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;
+
+
+        % 1) normlaize
+        y_rx = y_rx.normalize("mode","rms");
+
+        eq = EQ_silas("Ne",[50,9,9],"Nb",[2,0,0],"trainlength",4096,...
+            "sps",2,...
+            "mu_dc_dd",mudc_current,...
+            "mu_dc_train",mudc_current,...
+            "mu_ffe_train",0,...
+            "mu_dfe_train",0.005,...
+            "mu_ffe_dd",[0.0004 0.0004 0.0004],...
+            "mu_dfe_dd",0.005,...
+            "ddloops",4,...
+            "trainloops",4,...
+            "eq_blocklength",10);
+
+        [Eq_out] = eq.process(y_rx,y_digimod);
+
+
+        % 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(['fsym: ',num2str(parfsym(j,i).*1e-9),'GBd, SIR: ',num2str(parsir(j,i)),' ->> BER_new: ',sprintf('%2E',ber_dhat),' BER_old: ',sprintf('%2E',ber_nodctap)]);
+
+        % 5) save to struct
+        
+        ber_wo_reduction(j,i) = ber_nodctap;
+        ber_w_dcremoval(j,i) = ber_dhat;
+        mudc(j,i) = eq.mu_dc_dd;
 
-    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;
-
-
-    % 1) normlaize
-    y_rx = y_rx.normalize("mode","rms");
-
-
-    % 2.1) equalize
-    eq = EQ("K",2,"plottrain",1,"plotfinal",1,...
-        "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);
-
-    eq = EQ_silas("Ne",[50,9,9],"Nb",[2,0,0],"trainlength",4096,...
-        "sps",2,...
-        "mu_dc_dd",0,...
-        "mu_dc_train",0,...
-        "mu_ffe_train",0,...
-        "mu_dfe_train",0.005,...
-        "mu_ffe_dd",[0.0004 0.0004 0.0004],...
-        "mu_dfe_dd",0.005,...
-        "ddloops",3,...
-        "trainloops",4,...
-        "dcmode",1);
-
-    [Eq_out] = eq.process(y_rx,y_digimod);
-
-
-    % 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)
+for j = 1:length(mudc_loop)
-    a = [parsir(parfsym==rate(r));ber_wo_reduction(parfsym==rate(r))];
+    for r = 1:length(rate)
-    a=sort(a,2);
+        a = [parsir(j,parfsym(j,:)==rate(r));ber_wo_reduction(j,parfsym(j,:)==rate(r))];
-    plot(a(1,:),a(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; w/o reduction'],'LineWidth',2,'Marker','o','MarkerSize',3);
+        a=sort(a,2);
-    %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')
+        plot(a(1,:),a(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; w/o reduction'],'LineWidth',2,'Marker','o','MarkerSize',3);
-    hold on
+        %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')
-    b = [parsir(parfsym==rate(r));ber_w_dcremoval(parfsym==rate(r))];
+        hold on
-    b=sort(b,2);
+        b = [parsir(j,parfsym(j,:)==rate(r));ber_w_dcremoval(j,parfsym(j,:)==rate(r))];
-    plot(b(1,:),b(2,:),'DisplayName',[num2str(rate(r).*1e-9),' GBd; PAM4; adaptive DC removal'],'LineWidth',2,'Marker','+','MarkerSize',3);
+        b=sort(b,2);
-    %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','+');
+        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
 end
-set(gca, 'YScale', 'log')
+set(gca, 'YScale', 'log');
 yline(3.8e-3,'LineWidth',1.5);
 
 

setup_simulation.m

@@ -13,9 +13,9 @@ for lp = 1
 
     fiblen = 0; %main link in km
 
-    laser_linewidth =1e6;
+    laser_linewidth =5e6;
 
-    O = 18; %order of prbs
+    O = 19; %order of prbs
     N = 2^(O-1); %length of prbs
     [~,seed] = prbs(O,1); %initialize first seed of prbs
 
@@ -52,7 +52,7 @@ for lp = 1
 
     u_pi = 3.5;
     vbias = (0.5*u_pi)-u_pi;
-    extmodlaser = EML("mode",eml_mode.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth);
+    extmodlaser = EML("mode",eml_mode.im_cosinus,"power",5,"fsimu",fsimu,"lambda",1310,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth);
 
     fib         = Fiber("fsimu",fdac*kover,"fiber_length",fiblen,"alpha",0.2,"D",16,"lambda0",thz2nm(193.1),"gamma",0);
 
@@ -77,8 +77,17 @@ for lp = 1
         "dd_loops",2,"epsilon",[10 100 1000 ],"M",2,...
         "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0);
 
-    eq2 = EQ_silas("Ne",[30,0,0],"Nb",[3,0,0],"trainlength",4096,"mu_dc_dd",0.005,...
+    eq = EQ_silas("Ne",[50,5,3],"Nb",[2,0,0],"trainlength",4096,...
-        "mu_dc_train",0.05,"mu_ffe_train",0.005,"mu_combined_dd",[0.0004 0.0006 0.0003 0.005],"ddloops",3,'dcmode',3);
+        "sps",2,...
+        "mu_dc_dd",0.05,...
+        "mu_dc_train",0.05,...
+        "mu_ffe_train",0,...
+        "mu_dfe_train",0.005,...
+        "mu_ffe_dd",[0.0004 0.0004 0.0004],...
+        "mu_dfe_dd",0.005,...
+        "ddloops",3,...
+        "trainloops",4,...
+        "dcmode",1);
 
 
     %% C) PROCESS TX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -102,7 +111,7 @@ for lp = 1
 
     % 6) Lowpass behavior of laser and hf-cable? why twice?
     X = lp_laser.process(X);
-    X = lp_laser.process(X);
+    % X = lp_laser.process(X);
 
     % 7) Normalize signal
     X = X.normalize("mode","oneone");
@@ -141,38 +150,55 @@ for lp = 1
         [pos_high] = find(abs(phasediff)>thresh);
         [pos_low] = find(abs(phasediff)<=thresh);
 
+
+        figure(6);hold on;histogram(phaseX-phaseR,500,'EdgeAlpha',0)
+
         figure()
         subplot(3,1,1)
-        scatter(xax_sec.',wrapToPi(phaseX),4,'.','MarkerEdgeColor',col(2,:),'DisplayName','Phase of original Signal')
+        scatter(xax_sec.',wrapToPi(phaseX),4,'.','MarkerEdgeColor',col(2,:),'DisplayName','Desired Signal')
         hold on
-        scatter(xax_sec,wrapToPi(phaseR),4,'.','MarkerEdgeColor',col(1,:),'DisplayName','Phase of delayed Signal')
+        scatter(xax_sec,wrapToPi(phaseR),4,'.','MarkerEdgeColor',col(1,:),'DisplayName','Reflected Signal')
         hold on
-%         xline(delay/(physconst("LightSpeed")/1.4677),'LineWidth',4,'LineStyle','--','HandleVisibility','off')
         lg = legend;
         lg.Location = "southwest";
+        xlim([xax_sec(1) xax_sec(end)]);
+        xlabel('time in $\mu$s')
+        ylim([-pi pi]);
+        yticks([-pi 0 pi])
+        yticklabels({'$-\pi$',0, '$\pi$'})
+        ylabel('$\phi$')
 
         subplot(3,1,2)
         scatter(xax_sec(pos_low),phasediff(pos_low),4,'.','MarkerEdgeColor',col(5,:),'DisplayName','Phase Difference')
         hold on
         scatter(xax_sec(pos_high),phasediff(pos_high),4,'.','MarkerEdgeColor',col(2,:),'HandleVisibility','off')
-%         xline(delay/(physconst("LightSpeed")/1.4677),'LineWidth',4,'LineStyle','--','HandleVisibility','off')
+%         yline(thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off')
-        yline(thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off')
+%         yline(-thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off')
-        yline(-thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off')
         lg = legend;
         lg.Location = "southwest";
+        xlabel('time in $\mu$s')
+        xlim([xax_sec(1) xax_sec(end)]);
+        ylim([-pi pi]);
+        yticks([-pi 0 pi])
+        yticklabels({'$-\pi$',0, '$\pi$'})
+        ylabel('$\Delta \phi$')
+
 
         X.signal = X.signal(n:end);
         R.signal = R.signal(n:end);
+        X1 = X;
         X = X+R;
 
         subplot(3,1,3)
-        scatter(xax_sec(pos_low),abs(X.signal(pos_low).^2),4,'.','MarkerEdgeColor',col(5,:),'DisplayName','Constructive Interference');
+        scatter(xax_sec(pos_low),abs(X.signal(pos_low).^2)*1000,4,'.','MarkerEdgeColor',col(5,:),'DisplayName','Constructive Interference');
         hold on
-        scatter(xax_sec(pos_high),abs(X.signal(pos_high).^2),4,'.','MarkerEdgeColor',col(2,:),'DisplayName','Destructive Interference');
+        scatter(xax_sec(pos_high),abs(X.signal(pos_high).^2)*1000,4,'.','MarkerEdgeColor',col(2,:),'DisplayName','Destructive Interference');
-%         xline(delay/(physconst("LightSpeed")/1.4677),'LineWidth',4,'LineStyle','--','HandleVisibility','off')
+        %scatter(xax_sec(pos_high),abs(X1.signal(pos_high).^2)*1000,4,'.','MarkerEdgeColor',col(1,:),'MarkerFaceAlpha',0.3,'DisplayName','Destructive Interference');
-
+        xlim([xax_sec(1) xax_sec(end)]);
+        xlabel('time in $\mu$s')
         lg = legend;
         lg.Location = "southwest";
+        ylabel('Optical power in mW');
 
 
 
@@ -255,8 +281,6 @@ for lp = 1
 
     %     figure(1);bar([evm_bm' evm_dcsm' evm_llp' evm_lsm']);ylim([0.01 0.1]);set(gca,'yscale','log');
 
-
-
     %% PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
     % 1) Digi Demod
@@ -316,7 +340,8 @@ for lp = 1
 
         subplot(1,3,2)
         xindices = 1:Eq_out.length;
-        scatter(xindices(loc==0),Eq_out.signal(loc==0),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After EQ');
+        xax_sec = (1:Eq_out.length) / fsym .* 1e6;
+        scatter(xax_sec(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
@@ -345,17 +370,6 @@ for lp = 1
         hold off
 
 
-                % ERROR IN EQ
-        figure(13)
-        plot((eq2.error_log(1,:)),'DisplayName','First iteration DD mode','Linewidth',0.5);
-        hold on; 
-        for i = 1:size(eq2.error_log,1)
-            plot((eq2.error_log(i,:)),'Linewidth',0.5);
-        end
-        loc(loc==0) = NaN;
-        stem(loc.*mean((eq2.error_log(end,:))),'DisplayName','Error Positions')
-        hold off
-
     end
 
     if 0