From fce274c4cf30ba0ad9236b9eec321e0081f3379f Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Fri, 29 Sep 2023 08:18:22 +0200 Subject: [PATCH] Equalizer with different DC removal positions --- Classes/01_transmit/AWG.m | 2 +- Classes/02_optical/EML.m | 1 - Classes/04_DSP/EQ_silas.m | 63 ++++++-- imddmodel.m | 312 +++++++++++++++++++------------------- run_loop.m | 116 +++++++------- setup_simulation.m | 275 ++++++++++++++++++++------------- 6 files changed, 443 insertions(+), 326 deletions(-) diff --git a/Classes/01_transmit/AWG.m b/Classes/01_transmit/AWG.m index cb8f558..f7a8cd4 100644 --- a/Classes/01_transmit/AWG.m +++ b/Classes/01_transmit/AWG.m @@ -135,7 +135,7 @@ classdef AWG elec_out = repmat(elec_out,obj.repetitions,obj.kover); elec_out = reshape(elec_out',[],1); - % 3. Add skew (not working so far) + % 3. Add skew (not implemented so far) if obj.skew_active elec_out = obj.skew(elec_out); end diff --git a/Classes/02_optical/EML.m b/Classes/02_optical/EML.m index d054822..ebeb60b 100644 --- a/Classes/02_optical/EML.m +++ b/Classes/02_optical/EML.m @@ -111,7 +111,6 @@ classdef EML %scale with noisefactor noi = noi * obj.noisefactor; - %cumsum to accumulate noise over time vector noi = cumsum(noi); diff --git a/Classes/04_DSP/EQ_silas.m b/Classes/04_DSP/EQ_silas.m index 4ec87c0..cb33343 100644 --- a/Classes/04_DSP/EQ_silas.m +++ b/Classes/04_DSP/EQ_silas.m @@ -49,6 +49,8 @@ classdef EQ_silas < handle trainloops ddloops + dcmode + @@ -73,6 +75,8 @@ classdef EQ_silas < handle options.mu_dc_dd = 0.01; options.mu_combined_dd = [0.0004 0.0005 0.0006 0.0007 ]; + + options.dcmode end fn = fieldnames(options); @@ -151,18 +155,32 @@ classdef EQ_silas < handle 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).'; + + if obj.dcmode ~= 3 + x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).'; + elseif obj.dcmode == 3 + x_in_block = obj.x_in(obj.Ne(1)+n+(obj.sps-1):-1:n+obj.sps).' + obj.e_dc; + end + 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); + obj.e_ffe = obj.e.' * x_in_vnle_format; + + % Calculate the Error + if obj.dcmode == 1 + obj.error = obj.e_dc + obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay); + elseif obj.dcmode == 2 + obj.e_ffe = obj.e_ffe + obj.e_dc; + obj.error = obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay); + elseif obj.dcmode == 3 + obj.error = obj.e_ffe - obj.e_dfe - obj.d(obj.Nb(1)-1+m-obj.delay); + end %update FFE coefficients with LMS obj.e = obj.e - obj.error*conj(x_in_vnle_format)*obj.mu_ffe_train; @@ -179,6 +197,7 @@ classdef EQ_silas < handle end + function decisionDirectedMode(obj) %start the dd mode with coefficients from training @@ -195,6 +214,10 @@ classdef EQ_silas < handle 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 y = zeros(1,floor(obj.x_length/obj.sps)); @@ -207,29 +230,47 @@ classdef EQ_silas < handle m=m+1; %get Sigal input vectors with correct length for VNLE - x = obj.x_in(obj.Ne(1)+k-1:-1:k).'; + if obj.dcmode ~= 3 + x = obj.x_in(obj.Ne(1)+k-1:-1:k).'; + elseif obj.dcmode == 3 + x = obj.x_in(obj.Ne(1)+k-1:-1:k).' + obj.e_dc; + end + 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; + if obj.dcmode == 1 + y(m) = obj.e_dc + x_d.'* coeff; + elseif obj.dcmode == 2 || obj.dcmode == 3 +% x_ffe = obj.e.' * x_vnle; +% x_dfe = obj.b.' * d_vnle; +% y(m) = x_ffe - x_dfe; + y(m) = x_d.'* coeff; + end %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); + if obj.dcmode == 1 || obj.dcmode == 3 + obj.error = y(m) - d_hat(k); + elseif obj.dcmode == 2 + obj.error = y(m) - d_hat(k) + obj.e_dc; + end %Update coefficients (both FFE and DFE) + obj.e = obj.e - obj.error * mu_ffe * conj(x_vnle); + obj.b = obj.b + obj.error * mu_dfe * conj(d_vnle); + 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 + obj.e_dc = obj.e_dc - obj.mu_dc_dd * obj.error; + + obj.error_log(ddloop,m) = obj.e_dc.^2; % Append new decision to decision feedback if obj.Nb(1) > 0 diff --git a/imddmodel.m b/imddmodel.m index 3743fc8..7510896 100644 --- a/imddmodel.m +++ b/imddmodel.m @@ -1,206 +1,218 @@ -function output = imddmodel(sir,winlen,linewidth) +function output = imddmodel(sir,dcmode) for realiz = 1:3 rng(realiz); - + %% Set Simulation Variables - - sir = sir; - - delay = 50; %mpi delay in meter - + delay = 0;10+(10*realiz); %mpi delay in meter + fiblen = 0; %main link in km - - laser_linewidth = linewidth; - - O = 17; %order of prbs + + laser_linewidth =0e6; + + O = 18; %order of prbs N = 2^(O-1); %length of prbs [~,seed] = prbs(O,1); %initialize first seed of prbs - + % Modulation M = 4; %PAM-M bitpattern = zeros(N,log2(M)); - + % Symbol Rate fsym = 112e9; - + % DAC Rate fdac = 120e9; - + % Simulation oversampling rate "k"; kover = 16; - + % ADC Rate fadc = 256e9; - + % Simulation frequency in "analog domain" fsimu = kover * fdac ; - - - %% CONSTRUCT ALL CLASSES - + + + + %% B) CONSTRUCT ALL CLASSES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + digimod = PAMmapper(M,0); - - pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.1); - + + pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.027); + awg = AWG('fdac',fdac,'kover',kover,'lpf_active',1,'f_cutoff',56e9,'lpf_type',filtertypes.gaussian,'bit_resolution',5.5); - + lp_laser = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); - + u_pi = 3.5; vbias = (0.5*u_pi)-u_pi; - extmodlaser = EML("mode",eml_mode.im_cosinus,"power",15,"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",1550,"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); - + 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); - - edfaamp = Amplifier("amp_mode","edfa_increase_nase","gain_mode","output_power","amplification_db",0,"nase_mode","generate_ase","noifig",10.5); - + reflectionprop = Fiber("fsimu",fdac*kover,"fiber_length",delay/1000,"alpha",0.2,"D",16,"lambda0",1550,"gamma",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); - + lp_diode = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); - + scp = Scope("fsimu",fdac*kover,"fadc",fadc,... "delay",0,"fixed_delay",0,"lpf_bw",113e9,"filtertype",filtertypes.butterworth,... "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... "adcresolution",6,"quantbuffer",0.1,'block_dc',1); - + eq = EQ("K",2,"plottrain",0,"plotfinal",0,... - "training_length",1024,"training_loops",5,... - "Ne",[50,0,0],"Nb",[0,0,0],... - "DCmu",0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0.00,... - "dd_loops",2,"epsilon",[10 100 1000 ],"M",4,... - "thres",[0.005 0.004 0.0005 ],"l1act",0,"delay",1,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); - - - %% PROCESS - - % PRBS Generation + "training_length",4096,"training_loops",2,... + "Ne",[50,5,3],"Nb",[3,2,2],... + "DCmu",0.05,"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",0,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); + + eq2 = EQ_silas("Ne",[20,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",dcmode); + + + %% C) PROCESS TX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) PRBS Generation 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); - - % Digi Mod - mod_out = digimod.map(bits); - - % merken für EQ training - reference = mod_out; - - % shape shape - X = pulseform.process(mod_out); - - % AWG -> ELECTRICAL DOMAIN - awg_out = awg.process(X); - - X = lp_laser.process(awg_out); - X = lp_laser.process(X); - - X = X.normalize("mode","oneone"); - X.signal = X.signal .* 1.3800; - - % Laser; Modulation -> OPTICAL DOMAIN - X = extmodlaser.process(X); - - % Fiber Propagation - X = fib.process(X); - - %% Reflect with attenuation + + % 2 ) Build Inf. signal class + bits = Informationsignal(bitpattern); + + % 3) Digi modulation -> PAM-M signal + digimod_out = digimod.map(bits); + + % 4) Pulse shaping -> racos + X = pulseform.process(digimod_out); + + % 5) AWG (lowpass, quantization, sample and hold) + X = awg.process(X); + + % 6) Lowpass behavior of laser and hf-cable? why twice? + X = lp_laser.process(X); + X = lp_laser.process(X); + + % 7) Normalize signal + X = X.normalize("mode","oneone"); + X.signal = X.signal .* 1.3800; + + + + %% D) PROCESS OPTICAL CHANNEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Laser; Modulation -> OPTICAL DOMAIN + [X,extmodlaser] = extmodlaser.process(X); + + % 2) Main fiber propagation + X = fib.process(X); + + if delay ~= 0 + + % 3) Reflection + % Reflection is just an attenuation R = reflectionpoint.process(X); - - % Propagate + + % Propagate back and forth (actual fiber propagation) R = reflectionprop.process(R); - - % disp(['SIR ',num2str(10*log10(X.power/R.power))]); - - % Delay - R = R.delay("delay_meter",delay); - + + % Delay the reflected signal + [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; + + % 4) Equalize + % cut reference signal to correct length (nessecary due to MPI delay) + digimod_out.signal = digimod_out.signal(round(n * fsym/fsimu) : end,:); + bitpattern = bitpattern(round(n * fsym/fsimu):end,:); + end - %% - X = opticatten.process(X); - % X = edfaamp.process(X); - - % Photo Diode -> ELECTRICAL DOMAIN - X = phdiode.process(X); - X = lp_diode.process(X); - - - % Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) - X = scp.process(X); - - % Resample to Symbol Rate - X = X.resample("fs_out",2*fsym,"fs_in",fadc); - - % Normalize - Eq_in = X.normalize("mode","rms"); - - % Equalizer - Eq_out = eq.process(Eq_in,reference); - - %% MPI reduction DC removal - wl = winlen; % symbols + % 4) Attenuation + X = opticatten.process(X); + % X = edfaamp.process(X); + + % 5) Photo Diode -> ELECTRICAL DOMAIN + X = phdiode.process(X); + X = lp_diode.process(X); + + + %% E) PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) + X = scp.process(X); + + % 2) Resample to 2x symbol rate + X = X.resample("fs_out",2*fsym,"fs_in",fadc); + + % 3) Normalize + Eq_in = X.normalize("mode","rms"); + + + + % MPI reduction DC removal BEFORE EQ + wl = 3000; % symbols + Eq_in.signal = Eq_in.signal - 1/wl .* movsum( Eq_in.signal,[wl/2,wl/2]); + + % Equalize Signal + [Eq_out] = eq2.process(Eq_in,digimod_out); + + + + %% A1: MPI reduction DC removal + wl = 1000; % symbols yk_dcsm = Eq_out; yk_dcsm.signal = Eq_out.signal - 1/wl .* movsum( Eq_out.signal,[wl/2,wl/2]); - - %% MPI reduction Level wise error removal + + %% A2: MPI reduction Level wise error removal yk_lvsm = Eq_out; yk_lvlp = Eq_out; pre_decision_level_uni = digimod.decide_pamlevel(Eq_out); pre_decision_level_bi = ( pre_decision_level_uni*2-3 ) .* 1/sqrt(5); - + e = Eq_out.signal - pre_decision_level_bi; - - lp_mpi = Filter('filtdegree',1,"f_cutoff",max(10*linewidth,1e6),"fsamp",fsym,"filterType",filtertypes.bessel_inp); - + + lp_mpi = Filter('filtdegree',1,"f_cutoff",2e6,"fsamp",fsym,"filterType",filtertypes.bessel_inp); + filtered = lp_mpi.process(e); - - wl = winlen; % symbols + + wl = 30; % symbols smoothed = ( 1/wl .* movsum(e,[wl/2,wl/2]) ); - + % remove interference for level = 0:3 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 - - %% EVM - [evm_bm(realiz,:),std_bm(realiz,:)] = calc_evm(Eq_out.signal, pre_decision_level_bi); - [evm_dcsm(realiz,:),std_dcsm(realiz,:)] = calc_evm(yk_dcsm.signal, pre_decision_level_bi); - [evm_lsm(realiz,:),std_lsm(realiz,:)] = calc_evm(yk_lvsm.signal, pre_decision_level_bi); - [evm_llp(realiz,:),std_llp(realiz,:)] = calc_evm(yk_lvlp.signal, pre_decision_level_bi); - - %% Digi Demod - d_bm = digimod.demap(Eq_out); - d_dcsm = digimod.demap(yk_dcsm); - d_lvsm = digimod.demap(yk_lvsm); - d_lvlp = digimod.demap(yk_lvlp); - - %% BER - dbit = length(d_bm.signal)-length(bitpattern); - - [~,errors_bm,ber_bm(realiz)] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),10000); - [~,errors_dcsm,ber_dcsm(realiz)] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),10000); - [~,errors_lvsm,ber_lvsm(realiz)] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),10000); - [~,errors_lvlp,ber_lvlp(realiz)] = calc_ber(d_lvlp.signal(1:end-dbit,:) ,bitpattern(1:end,:),10000); - -% 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))]); + %% PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Digi Demod + d_bm = digimod.demap(Eq_out); + d_dcsm = digimod.demap(yk_dcsm); + d_lvsm = digimod.demap(yk_lvsm); + d_lvlp = digimod.demap(yk_lvlp); + + % 2) BER + + dbit = length(d_bm.signal)-length(bitpattern); + + [~,errors_bm,ber_bm(realiz),loc] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",1); + [~,errors_dcsm,ber_dcsm(realiz)] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); + [~,errors_lvsm,ber_lvsm(realiz)] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); + [~,errors_lvlp,ber_lvlp(realiz)] = calc_ber(d_lvlp.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",0); end @@ -210,15 +222,5 @@ output.ber_dcsm = mean(ber_dcsm); output.ber_lvsm = mean(ber_lvsm); output.ber_lvlp = mean(ber_lvlp); -output.evm_bm = mean(evm_bm); -output.evm_dcsm = mean(evm_dcsm); -output.evm_lsm = mean(evm_lsm); -output.evm_llp = mean(evm_llp); - -output.std_bm = mean(std_bm); -output.std_dcsm = mean(std_dcsm); -output.std_lsm = mean(std_lsm); -output.std_llp = mean(std_llp); - end \ No newline at end of file diff --git a/run_loop.m b/run_loop.m index 6887546..6fe1b56 100644 --- a/run_loop.m +++ b/run_loop.m @@ -1,73 +1,48 @@ -clear +clear -sir_loop = [-22:1.5:-15]; +sir_loop = [-25:2:-13]; +sir_loop = 0; +lw_loop = [1,2,3]; -%dc_tap_loop = [0.001 0.005 0.01 0.05 0.1 0.25]; -bw_loop = [1:10]; - -lw_loop = [0e6:0.2e6:1e6 2e6:2e6:10e6 20e6:20e6:100e6]; - -data = cell(length(sir_loop),length(bw_loop),length(lw_loop)); +data = cell(length(sir_loop),length(lw_loop)); iterations=size(data); -for ix = 1:numel(data) +parfor ix = 1:numel(data) - [u1,u2,u3] = ind2sub(iterations,ix); + [u1,u2] = ind2sub(iterations,ix); - output = imddmodel(sir_loop(u1),bw_loop(u2),lw_loop(u3)); + output = imddmodel(sir_loop(u1),lw_loop(u2)); data{ix} = output; - -% ber(ix) = BER.ber; -% ber_a1(ix) = BER.ber_a1; - - disp(ix) end -%% Plot Winlen Contour -hdfec = 3.8e-3.*ones(size(sir_loop)); - -for dc = 1:size(data,2) - for lw = 1:size(data,3) - for s = 1:size(data,1) - ber_lvlp(s,dc,lw) = data{s,dc,lw}.ber_dcsm; -% ber_dcsm(wl,lw,s) = data{wl,s,lw}.ber_dcsm; - end - - a_bm = InterX([sir_loop;squeeze(ber_lvlp(:,dc,lw))'],[sir_loop;hdfec]); -% a_dcsm = InterX([sirloop;squeeze(ber_dcsm(wl,lw,:))'],[sirloop;hdfec]); - try - thres_a0(dc,lw) = -a_bm(1); - %thres_a1(wl,lw) = -a_dcsm(1); - catch - thres_a0(dc,lw) = NaN; - %thres_a1(wl,lw) = NaN; - end - +for sir = 1:size(data,1) + for lw = 1:size(data,2) + ber_bm(sir,lw) = data{sir,lw}.ber_bm; + ber_dcsm(sir,lw) = data{sir,lw}.ber_dcsm; + ber_lvsm(sir,lw) = data{sir,lw}.ber_lvsm; + ber_lvlp(sir,lw) = data{sir,lw}.ber_lvlp; end end -figure(1) -contour(lw_loop,bw_loop,thres_a0,14:0.2:21,'LineWidth',1.5,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); -a = flip(cbrewer2('seq','Spectral',32)); -a = [a(1:12,:); a(22:end,:)]; -colormap(a); -clim([16 21]); -ylabel("Window Length"); -xlabel("Linewidth in MHz"); -yticks(bw_loop); -yticklabels(bw_loop); +col = [ 0.6510 0.8078 0.8902 + 0.6980 0.8745 0.5412 + 0.9922 0.7490 0.4353]; +figure(22) +for l = 1:numel(lw_loop) + hold on + plot(sir_loop,ber_bm(:,l),'DisplayName',['DC mode:', num2str(lw_loop(l)),', Linewidth= 50 MHz'],'Marker','o','MarkerFaceColor',col(l,:),'Color',col(l,:),'LineWidth',2,'LineStyle','--'); +end set(gca,'yscale','log'); -set(gca,'xscale','log'); -grid minor +yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); +set(gca,'xdir','reverse'); - -%% Plot +%% Plot col = linspecer(8); figure(1) hold on @@ -122,6 +97,43 @@ title("MPI removal - Optimization of Averaging Window Length") +%% Plot Winlen Contour +hdfec = 3.8e-3.*ones(size(sir_loop)); + +for dc = 1:size(data,2) + for lw = 1:size(data,3) + for s = 1:size(data,1) + ber_lvlp(s,dc,lw) = data{s,dc,lw}.ber_dcsm; + % ber_dcsm(wl,lw,s) = data{wl,s,lw}.ber_dcsm; + end + + a_bm = InterX([sir_loop;squeeze(ber_lvlp(:,dc,lw))'],[sir_loop;hdfec]); + % a_dcsm = InterX([sirloop;squeeze(ber_dcsm(wl,lw,:))'],[sirloop;hdfec]); + try + thres_a0(dc,lw) = -a_bm(1); + %thres_a1(wl,lw) = -a_dcsm(1); + catch + thres_a0(dc,lw) = NaN; + %thres_a1(wl,lw) = NaN; + end + + end +end + +figure(1) +contour(lw_loop,bw_loop,thres_a0,14:0.2:21,'LineWidth',1.5,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +a = flip(cbrewer2('seq','Spectral',32)); +a = [a(1:12,:); a(22:end,:)]; +colormap(a); +clim([16 21]); +ylabel("Window Length"); +xlabel("Linewidth in MHz"); +yticks(bw_loop); +yticklabels(bw_loop); +set(gca,'yscale','log'); +set(gca,'xscale','log'); +grid minor + @@ -167,7 +179,7 @@ end xlabel("DC Tap Value"); ylabel("Required SIR to rech FEC in dB"); -%% +%% col = linspecer(7); figure(3) cnt=1; @@ -182,7 +194,7 @@ for lw = [1,2,11] yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); set(gca,'yscale','log'); grid minor - + xlim([15,30]); ylim([1e-4,1e-2]); xlabel("SIR in dB"); diff --git a/setup_simulation.m b/setup_simulation.m index 2218f6f..8e194d0 100644 --- a/setup_simulation.m +++ b/setup_simulation.m @@ -1,26 +1,21 @@ -sir_loop = -22; + lw = [0.1e6 1e6 10e6]; for lp = 1 - dc_mu = 0.05; + rng(9); - sir = sir_loop(1); + %% A) Set Simulation Variables - rng(lp); - - %% Set Simulation Variables - - sir = sir; - - delay = 20; %mpi delay in meter + sir = -18; + delay = 10; %mpi delay in meter fiblen = 0; %main link in km - laser_linewidth =1e6; + laser_linewidth =5e6; - O = 17; %order of prbs + O = 18; %order of prbs N = 2^(O-1); %length of prbs [~,seed] = prbs(O,1); %initialize first seed of prbs @@ -44,7 +39,8 @@ for lp = 1 fsimu = kover * fdac ; - %% CONSTRUCT ALL CLASSES + + %% B) CONSTRUCT ALL CLASSES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% digimod = PAMmapper(M,0); @@ -56,12 +52,12 @@ for lp = 1 u_pi = 3.5; vbias = (0.5*u_pi)-u_pi; - extmodlaser = EML("mode",eml_mode.im_cosinus,"power",15,"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",1550,"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); 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); + reflectionprop = Fiber("fsimu",fdac*kover,"fiber_length",delay/1000,"alpha",0.2,"D",16,"lambda0",1550,"gamma",0); opticatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",-8); @@ -77,101 +73,157 @@ for lp = 1 eq = EQ("K",2,"plottrain",0,"plotfinal",0,... "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,... + "DCmu",0.05,"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",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); + 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); - %% PROCESS + %% C) PROCESS TX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % PRBS Generation + % 1) PRBS Generation for i = 1:log2(M) [bitpattern(:,i),seed] = prbs(O,N,seed); end - % Build Inf. signal class - bits = Informationsignal(bitpattern); + % 2 ) Build Inf. signal class + bits = Informationsignal(bitpattern); - % Digi Mod - mod_out = digimod.map(bits); + % 3) Digi modulation -> PAM-M signal + digimod_out = digimod.map(bits); - % merken für EQ training - reference = mod_out; + % 4) Pulse shaping -> racos + X = pulseform.process(digimod_out); - % shape shape - X = pulseform.process(mod_out); + % 5) AWG (lowpass, quantization, sample and hold) + X = awg.process(X); - % AWG -> ELECTRICAL DOMAIN - X = awg.process(X); + % 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); - - X = lp_laser.process(X); - - X = X.normalize("mode","oneone"); - X.signal = X.signal .* 1.3800; - - % Laser; Modulation -> OPTICAL DOMAIN - [X,extmodlaser] = extmodlaser.process(X); - - % Fiber Propagation - X = fib.process(X); - - %% Reflect with attenuation - R = reflectionpoint.process(X); - - % Propagate back and forth - R = reflectionprop.process(R); + % 7) Normalize signal + X = X.normalize("mode","oneone"); + X.signal = X.signal .* 1.3800; - % Delay - [R,n] = R.delay("delay_meter",delay); + %% D) PROCESS OPTICAL CHANNEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Add together - X.signal = X.signal(n:end); - R.signal = R.signal(n:end); + % 1) Laser; Modulation -> OPTICAL DOMAIN + [X,extmodlaser] = extmodlaser.process(X); - disp(['SIR ',num2str(10*log10(X.power/R.power))]); + % 2) Main fiber propagation + X = fib.process(X); + + if delay ~= 0 + % 3) Reflection + % Reflection is just an attenuation + R = reflectionpoint.process(X); + + % Propagate back and forth (actual fiber propagation) + R = reflectionprop.process(R); + + % Delay the reflected signal + [R,n] = R.delay("delay_meter",delay); + +% % Add together +% col = cbrewer2('qual','Paired',8); +% +% xax = (1:X.length);% / fsimu * (physconst("LightSpeed")/1.4677); +% xax_mtr = (1:X.length) / fsimu * (physconst("LightSpeed")/1.4677); +% thresh = 0.7; +% phaseX = phase(X.signal); +% phaseR = phase(R.signal); +% phasediff = wrapToPi(phaseX-phaseR); +% [pos_high] = find(abs(phasediff)>thresh); +% [pos_low] = find(abs(phasediff)<=thresh); +% +% figure() +% subplot(3,1,1) +% scatter(xax_mtr,wrapToPi(phaseX),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','Phase of original Signal') +% hold on +% scatter(xax_mtr,wrapToPi(phaseR),4,'.','MarkerEdgeColor',col(7,:),'DisplayName','Phase of delayed Signal') +% hold on +% xline(delay,'LineWidth',4,'LineStyle','--','HandleVisibility','off') +% lg = legend; +% lg.Location = "southwest"; +% +% subplot(3,1,2) +% scatter(xax_mtr(pos_low),phasediff(pos_low),4,'.','MarkerEdgeColor',col(6,:),'DisplayName','Phase Difference') +% hold on +% scatter(xax_mtr(pos_high),phasediff(pos_high),4,'.','MarkerEdgeColor',col(4,:),'HandleVisibility','off') +% xline(delay,'LineWidth',4,'LineStyle','--','HandleVisibility','off') +% yline(thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off') +% yline(-thresh,'LineWidth',2,'LineStyle','-','HandleVisibility','off') +% lg = legend; +% lg.Location = "southwest"; + + X = X+R; + +% subplot(3,1,3) +% scatter(xax_mtr(pos_low),abs(X.signal(pos_low).^2),4,'.','MarkerEdgeColor',col(6,:),'DisplayName','Constructive Interference'); +% hold on +% scatter(xax_mtr(pos_high),abs(X.signal(pos_high).^2),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','Destructive Interference'); +% xline(delay,'LineWidth',4,'LineStyle','--','HandleVisibility','off') +% +% lg = legend; +% lg.Location = "southwest"; + + X.signal = X.signal(n:end); + R.signal = R.signal(n:end); + + + disp(['SIR ',num2str(10*log10(X.power/R.power))]); + + % cut reference signal to correct length (nessecary due to MPI delay) + digimod_out.signal = digimod_out.signal(round(n * fsym/fsimu) : end,:); + + bitpattern = bitpattern(round(n * fsym/fsimu):end,:); + + end - X = X+R; +% plot(angle(R.signal)) + + % 4) Attenuation + X = opticatten.process(X); + - %% - X = opticatten.process(X); - % X = edfaamp.process(X); + % 5) Photo Diode -> ELECTRICAL DOMAIN + X = phdiode.process(X); + X = lp_diode.process(X); - % Photo Diode -> ELECTRICAL DOMAIN - X = phdiode.process(X); - X = lp_diode.process(X); - % Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) - X = scp.process(X); + %% E) PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Resample to Symbol Rate - X = X.resample("fs_out",2*fsym,"fs_in",fadc); + % 1) Oscilloscope (Sampling to f_adc; Quantization; Bandwidth Limitation) + X = scp.process(X); - % Normalize - Eq_in = X.normalize("mode","rms"); + % 2) Resample to 2x symbol rate + X = X.resample("fs_out",2*fsym,"fs_in",fadc); - % Equalizer - reference.signal = reference.signal(round(n * fsym/fsimu) : end,:); + % 3) Normalize + Eq_in = X.normalize("mode","rms"); - tracking_speed_bandwidth = eq.FFEmu .* fsym .* 1e-9; + % 4) Equalize - %% 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]); + % MPI reduction DC removal BEFORE EQ + wl = 3000; % symbols + Eq_in.signal = Eq_in.signal - 1/wl .* movsum( Eq_in.signal,[wl/2,wl/2]); - [Eq_out] = eq2.process(Eq_in,reference); + % Equalize Signal + [Eq_out] = eq2.process(Eq_in,digimod_out); - %% MPI reduction DC removal + + + %% A1: MPI reduction DC removal wl = 1000; % symbols yk_dcsm = Eq_out; yk_dcsm.signal = Eq_out.signal - 1/wl .* movsum( Eq_out.signal,[wl/2,wl/2]); - %% MPI reduction Level wise error removal + %% A2: MPI reduction Level wise error removal yk_lvsm = Eq_out; yk_lvlp = Eq_out; pre_decision_level_uni = digimod.decide_pamlevel(Eq_out); @@ -183,7 +235,7 @@ for lp = 1 filtered = lp_mpi.process(e); - wl = 100; % symbols + wl = 30; % symbols smoothed = ( 1/wl .* movsum(e,[wl/2,wl/2]) ); % remove interference @@ -191,14 +243,8 @@ 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]); - - - - %% Calc EVM + + % Calc EVM evm_bm = calc_evm(Eq_out.signal, pre_decision_level_bi); evm_dcsm = calc_evm(yk_dcsm.signal, pre_decision_level_bi); @@ -207,35 +253,35 @@ for lp = 1 % 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); - d_dcsm = digimod.demap(yk_dcsm); - d_lvsm = digimod.demap(yk_lvsm); - d_lvlp = digimod.demap(yk_lvlp); - % BER - bitpattern = bitpattern(round(n * fsym/fsimu):end,:); + + %% PROCESS RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 1) Digi Demod + d_bm = digimod.demap(Eq_out); + d_dcsm = digimod.demap(yk_dcsm); + d_lvsm = digimod.demap(yk_lvsm); + d_lvlp = digimod.demap(yk_lvlp); + + % 2) BER dbit = length(d_bm.signal)-length(bitpattern); - [~,errors_bm,ber_bm,loc] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),"skip",0,"returnErrorLocation",1); + [~,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 + % Display BER disp(['BER benchmark: ', sprintf('%2E',ber_bm), ' ERRORS: ' ,num2str(sum(errors_bm))]); + disp(['BER dc smooth (A1): ', sprintf('%2E',ber_dcsm), ' ERRORS: ' ,num2str(sum(errors_dcsm))]); + disp(['BER lv smooth (A2): ', 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); @@ -244,14 +290,16 @@ for lp = 1 figure(3) sgtitle('') - subplot(1,2,1) + subplot(1,3,1) hold on eq_decision = digimod.decide_pamlevel(Eq_out); - true_symbols = digimod.decide_pamlevel(mod_out); + true_symbols = digimod.decide_pamlevel(digimod_out); xindices = 1:Eq_in.length; errorpos = find(loc~=0)*2; + errorpos(errorpos>length(Eq_in.signal)) = length(Eq_in.signal); correct = find(loc==0)*2; + correct(correct>length(Eq_in.signal)) = length(Eq_in.signal); 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'); @@ -264,7 +312,7 @@ for lp = 1 a = legend; a.Location = "best"; - subplot(1,2,2) + 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'); hold on @@ -279,16 +327,31 @@ for lp = 1 a.Location = "best"; hold off + subplot(1,3,3) + xindices = 1:yk_lvsm.length; + scatter(xindices(loc==0),yk_lvsm.signal(loc==0),4,'.','MarkerEdgeColor',col(4,:),'DisplayName','After A2'); + hold on + scatter(xindices(loc~=0),yk_lvsm.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,:))); + plot((eq2.error_log(1,:)),'DisplayName','First iteration DD mode','Linewidth',0.5); hold on; for i = 1:size(eq2.error_log,1) - plot(diff(eq2.error_log(i,:))); + plot((eq2.error_log(i,:)),'Linewidth',0.5); end loc(loc==0) = NaN; - stem(loc.*mean((eq2.error_log(end,:)))) + stem(loc.*mean((eq2.error_log(end,:))),'DisplayName','Error Positions') hold off end @@ -301,7 +364,7 @@ for lp = 1 subplot(2,1,1) hold on - plot(reference.signal(4150:4175),'DisplayName','Tx','Color',col(1,:),'LineWidth',3); + plot(digimod_out.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); title('Modulated Sequence Zoom'); legend