From 3c11e4b2e54234b93c987d6fa3e50f162aa7fc2c Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Mon, 10 Jul 2023 15:08:25 +0200 Subject: [PATCH] version where the signal is flipped and added together --- Classes/00_signals/Signal.m | 38 +-- Classes/01_transmit/PAMmapper.m | 12 + Classes/02_etc/Filter.m | 25 +- Classes/02_optical/EML.m | 17 +- Classes/04_DSP/EQ.m | 42 ++- Functions/calc_ber.m | 13 +- Functions/calc_evm.m | 19 ++ imddmodel.m | 373 ++++++++++++++------------ run_loop.m | 213 ++++++++++++--- setup_simulation.m | 447 ++++++++++++++++++++------------ 10 files changed, 778 insertions(+), 421 deletions(-) create mode 100644 Functions/calc_evm.m diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index f6ea035..051c6ca 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -280,7 +280,7 @@ classdef Signal end %% - function obj = delay(obj,options) + function [obj,delay_n] = delay(obj,options) arguments obj Opticalsignal @@ -291,27 +291,29 @@ classdef Signal delay_n = round(delay_t .* obj.fs); - % build "long" hann window to fade the signal in and out - % -> prevent hard step in the signal! - hann_wind = hann(200); - ones_wind = ones(size(obj.signal)); - ones_wind(1:100) = hann_wind(1:100); - ones_wind(end-100:end) = hann_wind(end-100:end); - - % subtract average - mu = mean(obj.signal,"all"); - - obj.signal = obj.signal - mu; - - %apply hann - obj.signal = obj.signal .* ones_wind; - - %add average again - obj.signal = obj.signal + mu; +% % build "long" hann window to fade the signal in and out +% % -> prevent hard step in the signal! +% hann_wind = hann(200); +% ones_wind = ones(size(obj.signal)); +% ones_wind(1:100) = hann_wind(1:100); +% ones_wind(end-100:end) = hann_wind(end-100:end); +% +% % subtract average +% mu = mean(obj.signal,"all"); +% +% obj.signal = obj.signal - mu; +% +% %apply hann +% obj.signal = obj.signal .* ones_wind; +% +% %add average again +% obj.signal = obj.signal + mu; % finally circshift the signal obj.signal=circshift(obj.signal,delay_n); +% obj.signal=[obj.signal(delay_n:end); zeros(delay_n-1,1)]; + end end end diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index 89bdb54..e7356d0 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -180,6 +180,18 @@ classdef PAMmapper end + function [data_out] = decide_pamlevel(obj,data_in) + a = squeeze(repmat(real(data_in.signal),[1 1 length(obj.thresholds)])); %Eingangssignal in 3 spalten + b = squeeze(repmat(reshape(obj.thresholds(:).',[1 1 length(obj.thresholds)]),[1 length(data_in.signal) 1])); %Threshold in 3 Spalten + comp_real = a > b; %check for each symbol/ sampling if it exeeds the obj.thresholdseshold 1, 2 or 3 + + data_out = sum(comp_real,2); + + %data_out = (data_out*2)-3; + + %data_out = pam_level_decision .* 1/sqrt(5); + end + end end diff --git a/Classes/02_etc/Filter.m b/Classes/02_etc/Filter.m index 2a0e8b7..a38cab9 100644 --- a/Classes/02_etc/Filter.m +++ b/Classes/02_etc/Filter.m @@ -40,16 +40,23 @@ classdef Filter end - function signalclass_out = process(obj,signalclass_in) + function signal_out = process(obj,signal_in) + + if isa(signal_in,'Signal') - % actual processing of the signal - signalclass_out = signalclass_in; - signalclass_out.signal = obj.process_(signalclass_in.signal); + % actual processing of the signal + signal_out = signal_in; + signal_out.signal = obj.process_(signal_in.signal); + + % append to logbook + filterdesc = [num2str(obj.filtdegree),'. order ',char(obj.filterType),' filter with f_cutoff at ', num2str(obj.f_cutoff*1e-9), ' GHz.']; + signal_out = signal_out.logbookentry(filterdesc); - % append to logbook - filterdesc = [num2str(obj.filtdegree),'. order ',char(obj.filterType),' filter with f_cutoff at ', num2str(obj.f_cutoff*1e-9), ' GHz.']; - signalclass_out = signalclass_out.logbookentry(filterdesc); + elseif isa(signal_in,"double") + signal_out = obj.process_(signal_in); + + end end @@ -179,8 +186,8 @@ classdef Filter % xlabel('Freq in GHz') % ylabel('Magnitude (dB)') - -% freqz(B, A, 2048, obj.fsamp); +% figure(30) +% freqz(B, A, 2048, obj.fsamp); % hfvt = fvtool(B,A); end diff --git a/Classes/02_optical/EML.m b/Classes/02_optical/EML.m index 9eda123..d054822 100644 --- a/Classes/02_optical/EML.m +++ b/Classes/02_optical/EML.m @@ -61,10 +61,10 @@ classdef EML end - function signalclass_out = process(obj,signalclass_in) + function [signalclass_out,obj] = process(obj,signalclass_in) % actual processing of the signal (steps 1. - 3.) - signalclass_in.signal = obj.process_(signalclass_in.signal); + [signalclass_in.signal,obj] = obj.process_(signalclass_in.signal); % cast the inform. signal to electrical signal signalclass_in = Opticalsignal(signalclass_in,"fs",obj.fsimu,"logbook",signalclass_in.logbook,"lambda",obj.lambda*1e-9,"nase",0); @@ -92,7 +92,7 @@ classdef EML ph_noi = obj.createPhaseNoise; laserfield = obj.field.*exp(1i*ph_noi); %remember phase (! you need to receive the altered eml object in you sim program !) - obj.phase = ph_noi(end); + obj.phase = ph_noi; else laserfield = obj.field; end @@ -106,13 +106,18 @@ classdef EML function noi = createPhaseNoise(obj) %create random vector - noi = randn(obj.signal_len,1); + noi = randn(obj.signal_len/2,1); + %scale with noisefactor noi = noi * obj.noisefactor; + + %cumsum to accumulate noise over time vector noi = cumsum(noi); - %add phase from previous block/ loop of simulation - noi = noi + obj.phase; + + noi = [noi ; flip(noi)]; + + end function modulated_laserfield = externalmodulation(obj,laserfield,electrical_in) diff --git a/Classes/04_DSP/EQ.m b/Classes/04_DSP/EQ.m index 8784ed2..e34d62a 100644 --- a/Classes/04_DSP/EQ.m +++ b/Classes/04_DSP/EQ.m @@ -104,15 +104,15 @@ classdef EQ function signalclass_out = process(obj,signalclass_in, reference_signalclass_in) - % actual processing of the signal (steps 1. - 3.) - signalclass_in.signal = obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); + % actual processing of the signal (steps 1. - 3.) + signalclass_in.signal = obj.process_(signalclass_in.signal', reference_signalclass_in.signal'); signalclass_in.signal = signalclass_in.signal'; - % append to logbook - lbdesc = ['EQ ']; + % append to logbook + lbdesc = ['EQ ']; signalclass_in = signalclass_in.logbookentry(lbdesc); - % write to output + % write to output signalclass_out = signalclass_in; end @@ -218,14 +218,20 @@ classdef EQ e_ = zeros(obj.Ne(1)+N2+N3,1); % initialization of filter coefficients % e(ceil(obj.Ne(1)/2)) = 1; % set central tap to 1 (better starting point since it's closer to the expected solution) b_ = zeros(obj.Nb(1)+Nb2+Nb3,1); + e_dc = mean(data_in); % initilaization of the dc tap with the mean value of the data % e_save = NaN(361,8.6e5); % save_ind = 1; for trainloops = 1:obj.training_loops m = obj.k0+1; % starting symbol index at the delay compared to the training sequence + + error_log = []; for n = obj.K*obj.k0+1:obj.K:obj.K*obj.training_length m = m+1; X_1 = data(obj.Ne(1)+n+(obj.K-1):-1:n+obj.K).'; + +% X_1(X_1~=0) = X_1(X_1~=0) - mean(X_1(X_1~=0)); + [X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx); D_1 = ref(obj.Nb(1)-obj.k0+m-2:-1:m-obj.k0-1).'; @@ -255,14 +261,22 @@ classdef EQ % e_save(:,save_ind) = e; % save_ind = save_ind+1; e_dc = e_dc - obj.DCmu*error; + error_log(end+1) = e_dc; if obj.Nb(1) > 0 b_ = b_ + obj.DFEmu*error*reference_vec; % Seems like normalized DFE has worse performance end - end - end - %% + end + + + + + + end + + + %% % Plot the intermediate coefficients after training mode obj.b = b_(1:obj.Nb(1)); obj.b2 = b_(obj.Nb(1)+1:obj.Nb(1)+Nb2); @@ -273,7 +287,8 @@ classdef EQ if obj.plottrain figure(8052) - subplot(2,3,1); stem(abs(obj.e),'Markersize',2); + sgtitle('Training Coeff') + subplot(2,3,1); stem((obj.e),'Markersize',2); title('FFE coeff linear') xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) subplot(2,3,2); stem(obj.e2,'Markersize',2); @@ -291,9 +306,10 @@ classdef EQ subplot(2,3,6);stem(obj.b3,'Markersize',2); title('DFE coeff nl 3rd') xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) - set(gcf,'Position',[200 500 700 400]) + %set(gcf,'Position',[200 500 700 400]) end + if obj.l1act neg_lin = find(abs(obj.e) < obj.thres(1)); neg_2nd = find(abs(obj.e2) < obj.thres(2)); @@ -369,6 +385,9 @@ classdef EQ for k = 1:obj.K:length(data_in) m=m+1; % Symbol index X_1 = data(obj.Ne(1)+k-1:-1:k).'; + + + [X_2,X_3] = obj.calc_nl_vecs(X_1,ind_mat_2nd,ind_mat_3rd,norm_fac2,norm_fac3,delta_2,delta_3,cplx); if obj.l1act @@ -405,6 +424,7 @@ classdef EQ if mu_mat ~= 0 e_dc = e_dc - obj.DCmu*error; + error_log(end+1) = e_dc; end end @@ -434,6 +454,7 @@ classdef EQ if obj.plotfinal figure(8054) if obj.l1act + sgtitle('Final Coeff') subplot(2,3,1); stem(rel_lin,obj.e,'Markersize',2); title('FFE coeff linear') xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) @@ -453,6 +474,7 @@ classdef EQ title('DFE coeff nl 3rd') xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) else + sgtitle('Final Coeff') subplot(2,3,1); stem(obj.e/max(e_),'Markersize',2); title('FFE coeff linear') xlabel('coefficient index'); ylabel('value'); set(gca,'Fontsize',12) diff --git a/Functions/calc_ber.m b/Functions/calc_ber.m index 5bc9686..2d79f38 100644 --- a/Functions/calc_ber.m +++ b/Functions/calc_ber.m @@ -1,21 +1,20 @@ function [bits,errors,BER] = calc_ber(data_in,data_ref,skip) -data_ref=logical(data_ref)'; -data_in = logical(data_in)'; +data_ref = logical(data_ref)'; +data_in = logical(data_in)'; + bits = 0; - -% Determine BER - -bits = bits+size(data_in,2)-skip; +bits = bits+numel(data_in(:,skip+1:end)); try - errors = sum( data_in(:,skip+1:end,:) ~= data_ref(:,skip+1:end,:),2 ); + errors = sum( data_in(:,skip+1:end) ~= data_ref(:,skip+1:end),"all" ); catch %warning('BER calculation not optimal: Arrays have incompatible sizes for this operation.') errors = NaN; end +% Determine BER BER = sum(errors)/sum(bits); end \ No newline at end of file diff --git a/Functions/calc_evm.m b/Functions/calc_evm.m new file mode 100644 index 0000000..9dfe591 --- /dev/null +++ b/Functions/calc_evm.m @@ -0,0 +1,19 @@ +function [evm,stdev] = calc_evm(vector_received,vector_ideal) + + error_vector = (vector_received-vector_ideal); + + k = unique(vector_ideal); + + error = repmat(zeros(size(vector_ideal)),1,length(k)); + + for lvl = 1:length(k) + error(vector_ideal==k(lvl),lvl) = error_vector(vector_ideal==k(lvl)); + end + + error(error==0) = NaN; + + stdev = std(error,"omitnan"); + + error = sqrt(error.^2); + evm = sqrt( 1/length(error) .* sum(error.^2,1,'omitnan') ) ; +end \ No newline at end of file diff --git a/imddmodel.m b/imddmodel.m index cb5f9bf..3743fc8 100644 --- a/imddmodel.m +++ b/imddmodel.m @@ -1,173 +1,224 @@ -function BER = imddmodel(dc_mu,sir,linewidth) +function output = imddmodel(sir,winlen,linewidth) -rng(2023); +for realiz = 1:3 -%% Set Simulation Variables - -sir = sir; - -delay = 300; %mpi delay in meter - -fiblen = 0.5; %main link in km - -laser_linewidth = linewidth; + rng(realiz); + + %% Set Simulation Variables + + sir = sir; + + delay = 50; %mpi delay in meter + + fiblen = 0; %main link in km + + laser_linewidth = linewidth; + + O = 17; %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 + + digimod = PAMmapper(M,0); + + pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.1); + + 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); + + 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); + + 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 + 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 + R = reflectionpoint.process(X); + + % Propagate + R = reflectionprop.process(R); + + % disp(['SIR ',num2str(10*log10(X.power/R.power))]); + + % Delay + R = R.delay("delay_meter",delay); + + % Add together + X = X+R; + + %% + 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 + 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 + 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); + + filtered = lp_mpi.process(e); + + wl = winlen; % 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))]); -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 - -digimod = PAMmapper(M,0); - -pulseform = Pulseformer("pulseform","rrc","fdac",fdac,"fsym",fsym,"pulselength",32,"rrcalpha",0.1); - -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); - -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",-5); - -edfaamp = Amplifier("amp_mode","edfa_increase_nase","gain_mode","output_power","amplification_db",0,"nase_mode","generate_ase","noifig",10.5); - -phdiode = Photodiode("fsimu",fdac*kover,"dark_current",0,"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",dc_mu,"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 -for i = 1:log2(M) - [bitpattern(:,i),seed] = prbs(O,N,seed); end -%bitpattern = prms_out'; +output.ber_bm = mean(ber_bm); +output.ber_dcsm = mean(ber_dcsm); +output.ber_lvsm = mean(ber_lvsm); +output.ber_lvlp = mean(ber_lvlp); -% Build Inf. signal class -bits = Informationsignal(bitpattern); +output.evm_bm = mean(evm_bm); +output.evm_dcsm = mean(evm_dcsm); +output.evm_lsm = mean(evm_lsm); +output.evm_llp = mean(evm_llp); -% Digi Mod -mod_out = digimod.map(bits); - -% merken für EQ training -reference = mod_out; - -% shape shape -X = pulseform.process(mod_out); -test = applyPulseShaping(reference.signal,fsym,fdac); - -% 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 -R = reflectionpoint.process(X); - -% Propagate -R = reflectionprop.process(R); - -disp(['SIR ',num2str(10*log10(X.power/R.power))]); - -% Delay -R = R.delay("delay_meter",delay); - -% Add together -X = X+R; - -%% -X = opticatten.process(X); -X = edfaamp.process(X); - -% Photo Diode -> ELECTRICAL DOMAIN -X = phdiode.process(X); -X = lp_diode.process(X); +output.std_bm = mean(std_bm); +output.std_dcsm = mean(std_dcsm); +output.std_lsm = mean(std_lsm); +output.std_llp = mean(std_llp); -% 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 -X = X.normalize("mode","rms"); - -% Equalizer -X = eq.process(X,reference); - -%% MPI reduction -wl = 10000; -mpi_red = X; -mpi_red.signal = X.signal - 1/wl .* movsum( X.signal,[wl/2,wl/2]); - - -% Digi Demod -demap_out = digimod.demap(X); -demap_out_a1 = digimod.demap(mpi_red); - -% BER -[bits,errors,ber] = calc_ber(demap_out.signal(1:end-1,:),bitpattern(1:end,:),0); -[bits,errors,ber_a1] = calc_ber(demap_out_a1.signal(1:end-1,:),bitpattern(1:end,:),0); - -disp(['BER: ', sprintf('%2E',ber), ' ERRORS: ' ,num2str(sum(errors))]); -disp(['BER: ', sprintf('%2E',ber_a1), ' ERRORS: ' ,num2str(sum(errors))]); - - -BER.ber = ber; -BER.ber_a1 = ber_a1; end \ No newline at end of file diff --git a/run_loop.m b/run_loop.m index fe35284..6887546 100644 --- a/run_loop.m +++ b/run_loop.m @@ -1,61 +1,196 @@ -dc_mu = [0 0.001 0.1]; -sirloop = -25:1:0; -lw_loop = [0e6:1e6:10e6]; +clear -ber = zeros(length(dc_mu),length(sirloop),length(lw_loop)); -ber_a1 = zeros(length(dc_mu),length(sirloop),length(lw_loop)); +sir_loop = [-22:1.5:-15]; -iterations=size(ber); +%dc_tap_loop = [0.001 0.005 0.01 0.05 0.1 0.25]; +bw_loop = [1:10]; -parfor ix = 1:numel(ber) +lw_loop = [0e6:0.2e6:1e6 2e6:2e6:10e6 20e6:20e6:100e6]; + +data = cell(length(sir_loop),length(bw_loop),length(lw_loop)); + +iterations=size(data); + +for ix = 1:numel(data) [u1,u2,u3] = ind2sub(iterations,ix); - BER = imddmodel(dc_mu(u1),sirloop(u2),lw_loop(u3)); + output = imddmodel(sir_loop(u1),bw_loop(u2),lw_loop(u3)); - ber(ix) = BER.ber; - ber_a1(ix) = BER.ber_a1; + 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)); -hdfec = 3.8e-3.*ones(size(sirloop)); +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 dc = 1:size(ber,1) - for lw = 1:size(ber,3) - a = InterX([sirloop;ber(dc,:,lw)],[sirloop;hdfec]); - thres(dc,lw) = -a(1); end end -contour(thres,'LineWidth',2); +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 + + + + +%% Plot +col = linspecer(8); +figure(1) +hold on +m = ["x","o","pentagram","hexagram","*","+"]; +cnt = 1; +for d = 1%:size(data,2) + for dc = 1:size(data,1) + for lw = 1:size(data,3) + ber_bm(dc,d,lw) = (data{dc,d,lw}.ber_bm); + ber_dcsm(dc,d,lw) = data{dc,d,lw}.ber_dcsm; + ber_lvsm(dc,d,lw) = data{dc,d,lw}.ber_lvsm; + ber_lvlp(dc,d,lw) = (data{dc,d,lw}.ber_lvlp); + end + end + i = 1; + comm_dn = [];% ['EQ DC Tap: ',num2str(dc_tap_loop(d)),' m']; + title("Dependency on Laser Linewidth; B2B; Delay : 2*50m") + plot(lw_loop*1e-6,mean(squeeze(ber_bm(1:end,d,:)),1),"LineWidth",1.2,"Marker",m(1),"MarkerSize",5,'Color',col(cnt,:),'DisplayName',[' ',comm_dn]); + plot(lw_loop*1e-6,mean(squeeze(ber_dcsm(1:end,d,:)),1),"LineWidth",1,"Marker",m(1+4),"MarkerSize",5,'LineStyle','--','Color',col(cnt,:),'DisplayName',['DC smoothing ',comm_dn],'HandleVisibility','on'); + plot(lw_loop*1e-6,mean(squeeze(ber_lvsm(1:end,d,:)),1),"LineWidth",1.2,"Marker",m(1+1),"MarkerSize",5,'LineStyle',':','Color',col(cnt,:),'DisplayName',['Lvl Smoothing ',comm_dn],'HandleVisibility','on'); + plot(lw_loop*1e-6,mean(squeeze(ber_lvlp(1:end,d,:)),1),"LineWidth",1,"Marker",m(1+3),"MarkerSize",5,'LineStyle','-.','Color',col(cnt,:),'DisplayName',['Lvl Lowpass ',comm_dn],'HandleVisibility','on'); + set(gca,'yscale','log'); + set(gca,'xscale','log'); + %xticklabels([10 100 1000 10000]); + grid minor + yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); + ylim([1e-3 4e-2]); + ylabel("BER"); + xlabel("Linewidth in MHz") + legend + cnt = cnt+1; +end + + +figure(1) +contour(lw_loop,sir_loop,thres_a0,16:0.4:21,'LineWidth',2,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +clim([16 21]); +ylabel("Bandwidth in Multiples of Linewidth"); +xlabel("Linewidth in MHz"); +% set(gca,'yscale','log'); +grid minor +title("MPI removal - Optimization of Lowpass Filter Bandwidth") + +figure(2) +contour(lw_loop,sir_loop([1:10,12:end]),thres_a0([1:10,12:end],:),16:0.3:21,'LineWidth',2,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +clim([16 21]); +ylabel("Window Length"); +xlabel("Linewidth in MHz"); +set(gca,'yscale','log'); +grid minor +title("MPI removal - Optimization of Averaging Window Length") + + + + + + +%% Plot DC Tap Contour +hdfec = 3.8e-3.*ones(size(bw_loop)); +thres_a0 = zeros(size(ber,1),size(ber,3)); +thres_a1 = zeros(size(ber,1),size(ber,3)); +for dc = 1:size(ber,1) + for lw = 1:size(ber,3) + a_lvsm = InterX([bw_loop;ber(dc,:,lw)],[bw_loop;hdfec]); + thres_a0(dc,lw) = -a_lvsm(1); + + a1 = InterX([bw_loop;ber_a1(dc,:,lw)],[bw_loop;hdfec]); + thres_a1(dc,lw) = -a1(1); + end +end + +figure(1) +subplot(2,1,1) +contour(lw_loop,sir_loop,thres_a0,16:0.5:21,'LineWidth',3,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +clim([16 21]); ylabel("DC Tap"); xlabel("Linewidth in MHz"); +set(gca,'yscale','log'); +grid minor +subplot(2,1,2) +contour(lw_loop,sir_loop,thres_a1,16:0.5:21,'LineWidth',3,'FaceAlpha',0.3,'ShowText','on',"LabelFormat","%0.1f dB"); +clim([16 21]); +ylabel("DC Tap"); +xlabel("Linewidth in MHz"); +set(gca,'yscale','log'); +grid minor -yticks(1:length(dc_mu)) -yticklabels(dc_mu); +%% Plot Curves of required SIR to see the minimum a bit better +col = flip(cbrewer2('seq','Spectral',16)); +col = col([1:4, 10:end],:); +figure(2) +hold on +for lw = [size(ber,3):-2:2 2 1] + plot(sir_loop,thres_a0(:,lw),'DisplayName',[' Linewidth: ',num2str(lw_loop(lw)*1e-6), ' MHz'],'Color',col(lw,:),'Marker','o','MarkerFaceColor',col(lw,:),'LineWidth',2); + set(gca,'xscale','log'); +end +xlabel("DC Tap Value"); +ylabel("Required SIR to rech FEC in dB"); -xticks(1:length(lw_loop)); -xticklabels(lw_loop*1e-6); +%% +col = linspecer(7); +figure(3) +cnt=1; +for lw = [1,2,11] -grid on + for i = 1:length(sir_loop)-1 + + subplot(1,3,cnt) + hold on + plot(-1.*bw_loop,ber(i,:,lw),"LineWidth",2,"Marker","o","MarkerSize",5,'Color',col(i,:),'DisplayName',['DC tap ',num2str(sir_loop(i))]); + plot(-1.*bw_loop,ber_a1(i,:,lw),"LineWidth",2,"Marker","x","MarkerSize",5,'LineStyle','--','Color',col(i,:),'DisplayName',['A1. DC tap ',num2str(sir_loop(i))]); + 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"); + ylabel("BER"); + title(['BER for different SIR;',' Linewidth: ',num2str(lw_loop(lw)*1e-6), ' MHz']) + text(25,4.2e-3,"FEC $3.8 e^{-3}$"); + end + cnt = cnt+1; +end +legend -% -% col = cbrewer2('Paired',12); -% scatter(-1.*thres(1),thres(2),'MarkerEdgeColor',col(2*i,:),'LineWidth',6,'DisplayName',['DC tap ',num2str(dc_mu(i))]); -% -% figure(1) -% hold on -% for i = 1:length(dc_mu) -% plot(-1.*sirloop,ber(i,:),"LineWidth",2,"Marker","o","MarkerSize",3,'Color',col(2*i,:),'DisplayName',['DC tap ',num2str(dc_mu(i))]); -% plot(-1.*sirloop,ber_a1(i,:),"LineWidth",2,"Marker","o","MarkerSize",3,'Color',col(2*i-1,:),'DisplayName',['A1. Linewidth in MHz ',num2str(dc_mu(i)*1e-6)]); -% end -% -% yline(3.8e-3,'LineWidth',2,'LineStyle','--','HandleVisibility','off'); -% set(gca,'yscale','log'); -% grid minor -% legend -% xlabel("SIR"); -% ylabel("BER"); -% title("BER for different SIR") \ No newline at end of file diff --git a/setup_simulation.m b/setup_simulation.m index 14f779c..a732aa9 100644 --- a/setup_simulation.m +++ b/setup_simulation.m @@ -1,210 +1,315 @@ -sir_loop = -30:5:-10; +sir_loop = -16; -for lp = 1:5 +for lp = 1 -sir = sir_loop(lp); + dc_mu = 0.05; -rkey = 13; -rng(rkey); + sir = sir_loop(lp); -%% Set Simulation Variables + rng(10); -O = 18; %order of prbs -N = 2^(O-1); %length of prbs -[~,seed] = prbs(O,1); %initialize first seed of prbs + %% Set Simulation Variables -% Modulation -M = 4; %PAM-M -bitpattern = zeros(N,log2(M)); + sir = sir; -% Symbol Rate -fsym = 112e9; + delay = 40; %mpi delay in meter -% DAC Rate -fdac = 120e9; + fiblen = 0; %main link in km -% Simulation oversampling rate "k"; -kover = 16; + laser_linewidth = 1e6; -% ADC Rate -fadc = 256e9; + O = 17; %order of prbs + N = 2^(O-1); %length of prbs + [~,seed] = prbs(O,1); %initialize first seed of prbs -% Simulation frequency in "analog domain" -fsimu = kover * fdac ; + % 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 + %% CONSTRUCT ALL CLASSES -digimod = PAMmapper(M,0); + 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); + 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); + 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",5,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1e6); + 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); -fib = Fiber("fsimu",fdac*kover,"fiber_length",0,"alpha",0.2,"D",16,"lambda0",thz2nm(193.1),"gamma",0); + 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",1,"alpha",0.2,"D",16,"lambda0",1550,"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); + + 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",2048,"training_loops",5,... + "Ne",[50,0,0],"Nb",[0,0,0],... + "DCmu",dc_mu,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"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",1,"rho",0.0005,"ideal_dfe",0,"DB_aim",0); + %% PROCESS -opticatten = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",-5); + % PRBS Generation + for i = 1:log2(M) + [bitpattern(:,i),seed] = prbs(O,N,seed); + end + bitpattern = [bitpattern ; flip(bitpattern)]; + %bitpattern = prms_out'; -edfaamp = Amplifier("amp_mode","edfa_increase_nase","gain_mode","output_power","amplification_db",0,"nase_mode","generate_ase","noifig",10.5); + % 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] = extmodlaser.process(X); + % Fiber Propagation + X = fib.process(X); + + %% Reflect with attenuation + R = reflectionpoint.process(X); -phdiode = Photodiode("fsimu",fdac*kover,"dark_current",0,"responsivity",1,"temperature",20); + % Propagate + R = reflectionprop.process(R); -lp_diode = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterType",filtertypes.bessel_inp); + % disp(['SIR ',num2str(10*log10(X.power/R.power))]); -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); + % Delay + [R,n] = R.delay("delay_meter",delay); -eq = EQ("K",2,"plottrain",0,"plotfinal",0,... - "training_length",4096,"training_loops",5,... - "Ne",[50,0,0],"Nb",[0,0,0],... - "DCmu",0.005,"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); + % Add together + X = X+R; + + %% + X = opticatten.process(X); + % X = edfaamp.process(X); + + % Photo Diode -> ELECTRICAL DOMAIN + X = phdiode.process(X); + X = lp_diode.process(X); -%% PROCESS + % 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 = 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 + 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",2e6,"fsamp",fsym,"filterType",filtertypes.bessel_inp); + + filtered = lp_mpi.process(e); + + wl = 100; % 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 + %% MPI Reduction Tunable Notch Filter +% yk = Eq_in.normalize("mode","oneone"); +% yk = yk.signal; +% dk = Eq_out.resample("fs_in",fsym,"fs_out",2*fsym); +% dk = ( digimod.decide_pamlevel(dk) *2-3 ) .* 1/sqrt(5); +% dk = dk ./ max(abs(dk)); +% +% xk = yk(4:end)-dk(1:end-3); +% +% ya = hilbert(xk); +% figure;periodogram(abs(ya),[],length(ya),2*fsym,'centered') +% plot(abs(ya)); +% grid on; + + %% Calc EVM + + evm_bm = calc_evm(Eq_out.signal, pre_decision_level_bi); + evm_dcsm = calc_evm(yk_dcsm.signal, pre_decision_level_bi); + evm_lsm = calc_evm(yk_lvsm.signal, pre_decision_level_bi); + evm_llp = calc_evm(yk_lvlp.signal, pre_decision_level_bi); + + %figure(1);bar([evm_bm' evm_dcsm' evm_llp' evm_lsm']);ylim([0.01 0.1]);set(gca,'yscale','log'); + + %% Plot stuff +% figure(200) +% hold on +% for i = 1:4 +% subplot(4,1,i) +% scatter(1:length(err_),err_(:,i),4,'Marker','.','MarkerEdgeColor',col(i,:)); +% end +% +% figure(12); +% hold on; +% scatter(1:length(yk_sm.signal),yk_sm.signal,4,'Marker','.','MarkerEdgeColor',col(6,:)); +% for i = 1:4 +% x = err_(:,i); +% nanx = isnan(x); +% t = 1:numel(x); +% x(nanx) = interp1(t(~nanx), x(~nanx), t(nanx)); +% +% X(i,:) = x; +% +% plot(x(2000:end-2000,1),'LineWidth',2) +% end + %% MPI reduction A3 + + + % 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] = calc_ber(d_bm.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); + [~,errors_dcsm,ber_dcsm] = calc_ber(d_dcsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); + [~,errors_lvsm,ber_lvsm] = calc_ber(d_lvsm.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); + [~,errors_lvlp,ber_lvlp] = calc_ber(d_lvlp.signal(1:end-dbit,:) ,bitpattern(1:end,:),0); + + % Display + 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))]); + + %% Generate some Plots + if 0 + col = cbrewer2('Paired',8); + + figure(11) + clf + + subplot(2,1,1) + hold on + plot(reference.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 + hold off + subplot(2,1,2) + hold on + stem(d_bm.signal(4150:4175,1),'DisplayName','Tx','Color',col(1,:),'LineStyle','-','LineWidth',5) + stem(bitpattern(4150:4175,1)','DisplayName','Rx','Color',col(6,:),'LineStyle','--','LineWidth',2) + title('Bitpattern Tx - Rx'); + legend + hold off + + end + if 0 + col = cbrewer2('Paired',8); + + xax = 1:Eq_out.length; + + figure(3) + sgtitle('') + % subplot(1,4,1:2) + % scatter(1:4:X.length,X.signal(1:4:end),4,'.','MarkerEdgeColor',col(6,:),'DisplayName','Before EQ'); + % xlim([1, xax(end)]); + % %ylim([-2 2]); + % xlabel('Sampling Index') + % ylabel('Amplitude') + % legend + subplot(1,3,1) + scatter(1:Eq_in.length,Eq_in.signal,4,'.','MarkerEdgeColor',col(6,:),'DisplayName','After EQ'); + xlim([1, xax(end)]); + ylim([-2 2]); + xlabel('Sampling Index') + ylabel('Amplitude') + legend + + subplot(1,3,2) + scatter(xax,Eq_out.signal,4,'.','MarkerEdgeColor',col(6,:),'DisplayName','After EQ'); + xlim([1, xax(end)]); + ylim([-2 2]); + xlabel('Sampling Index') + ylabel('Amplitude') + legend + + subplot(1,3,3) + scatter(xax,yk_lvsm.signal,4,'.','MarkerEdgeColor',col(6,:),'DisplayName','After EQ'); + xlim([1, xax(end)]); + ylim([-2 2]); + xlabel('Sampling Index') + ylabel('Amplitude') + legend + end -% PRBS Generation -for i = 1:log2(M) - [bitpattern(:,i),seed] = prbs(O,N,seed); end -%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 .* 0.6*(u_pi/2-abs(((0.5*u_pi)-u_pi/2))); - -% Laser; Modulation -> OPTICAL DOMAIN -X = extmodlaser.process(X); - -% Fiber Propagation -X = fib.process(X); - -%% Reflect with attenuation -R = reflectionpoint.process(X); - -% Propagate -R = reflectionprop.process(R); - -% Delay -R = R.delay("delay_meter",200); - -% Add together -X = X+R; - -X.spectrum(fsimu,"displayname",['SIR: ',num2str(sir)],"figurename",'spectrum') - -%% -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 -X = X.normalize("mode","rms"); - -% Equalizer -eq_out = eq.process(X,reference); - -% Digi Demod -demap_out = digimod.demap(eq_out); - -% BER -[bits,errors,BER] = calc_ber(demap_out.signal(1:end-1,:),bitpattern(1:end,:),1000); - -disp(['BER: ', sprintf('%2E',BER), ' ERRORS: ' ,num2str(sum(errors))]); - - -%% Generate some Plots -if 0 - col = cbrewer2('Paired',8); - - figure(11) - clf - - subplot(2,1,1) - hold on - plot(reference.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 - hold off - subplot(2,1,2) - hold on - stem(demap_out.signal(4150:4175,1),'DisplayName','Tx','Color',col(1,:),'LineStyle','-','LineWidth',5) - stem(bitpattern(4150:4175,1)','DisplayName','Rx','Color',col(6,:),'LineStyle','--','LineWidth',2) - title('Bitpattern Tx - Rx'); - legend - hold off - -end -if 1 - col = cbrewer2('Paired',8); - - xax = 1:eq_out.length; - - figure(3) - sgtitle('') -% subplot(1,4,1:2) -% scatter(1:4:X.length,X.signal(1:4:end),4,'.','MarkerEdgeColor',col(6,:),'DisplayName','Before EQ'); -% xlim([1, xax(end)]); -% %ylim([-2 2]); -% xlabel('Sampling Index') -% ylabel('Amplitude') -% legend - - subplot(1,5,lp) - scatter(xax,eq_out.signal,4,'.','MarkerEdgeColor',col(6,:),'DisplayName','After EQ'); - xlim([1, xax(end)]); - ylim([-2 2]); - xlabel('Sampling Index') - ylabel('Amplitude') - legend -end - -end