From 025498d120ea8c3f42dc32ac7dd9a6b2fc6cc47e Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Thu, 12 Oct 2023 16:43:54 +0200 Subject: [PATCH] DC removal algorithm now with latency --- Classes/04_DSP/EQ_silas.m | 41 +++-- MPI_offline_analysis/mpi_analysis_pam4.m | 194 ++++++++++++----------- setup_simulation.m | 74 +++++---- 3 files changed, 169 insertions(+), 140 deletions(-) diff --git a/Classes/04_DSP/EQ_silas.m b/Classes/04_DSP/EQ_silas.m index 47b0591..695684e 100644 --- a/Classes/04_DSP/EQ_silas.m +++ b/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 diff --git a/MPI_offline_analysis/mpi_analysis_pam4.m b/MPI_offline_analysis/mpi_analysis_pam4.m index e4683a4..dabfd70 100644 --- a/MPI_offline_analysis/mpi_analysis_pam4.m +++ b/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)); -parsir = zeros(1,length(allfiles)); -ber_wo_reduction = zeros(1,length(allfiles)); -ber_w_dcremoval = zeros(1,length(allfiles)); -mudc = zeros(1,length(allfiles)); +mudc_loop = [0.05 0.01 0.005]; + +parfsym = zeros(length(mudc_loop),length(allfiles)); +parsir = zeros(length(mudc_loop),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) - 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','+'); +for j = 1:length(mudc_loop) + for r = 1:length(rate) + a = [parsir(j,parfsym(j,:)==rate(r));ber_wo_reduction(j,parfsym(j,:)==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(j,parfsym(j,:)==rate(r));ber_w_dcremoval(j,parfsym(j,:)==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 end -set(gca, 'YScale', 'log') +set(gca, 'YScale', 'log'); yline(3.8e-3,'LineWidth',1.5); diff --git a/setup_simulation.m b/setup_simulation.m index 4f89514..b9242c6 100644 --- a/setup_simulation.m +++ b/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,... - "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); + eq = EQ_silas("Ne",[50,5,3],"Nb",[2,0,0],"trainlength",4096,... + "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'); -% xline(delay/(physconst("LightSpeed")/1.4677),'LineWidth',4,'LineStyle','--','HandleVisibility','off') - + scatter(xax_sec(pos_high),abs(X.signal(pos_high).^2)*1000,4,'.','MarkerEdgeColor',col(2,:),'DisplayName','Destructive Interference'); + %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