silas/imdd_silas
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

version where the signal is flipped and added together

Browse Source
This commit is contained in:
Silas Oettinghaus
2023-07-10 15:08:25 +02:00
parent 93dd63bb18
commit 3c11e4b2e5
10 changed files with 778 additions and 421 deletions
Download Patch File Download Diff File Expand all files Collapse all files

38
Classes/00_signals/Signal.m
View File

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

12
Classes/01_transmit/PAMmapper.m
View File

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

17
Classes/02_etc/Filter.m
View File

@@ -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);
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.'];
signalclass_out = signalclass_out.logbookentry(filterdesc);
signal_out = signal_out.logbookentry(filterdesc);
elseif isa(signal_in,"double")
signal_out = obj.process_(signal_in);
end
end
@@ -179,7 +186,7 @@ classdef Filter
% xlabel('Freq in GHz')
% ylabel('Magnitude (dB)')
% figure(30)
% freqz(B, A, 2048, obj.fsamp);
% hfvt = fvtool(B,A);
end

17
Classes/02_optical/EML.m
View File

@@ -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)

32
Classes/04_DSP/EQ.m
View File

@@ -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)

9
Functions/calc_ber.m
View File

@@ -2,20 +2,19 @@ function [bits,errors,BER] = calc_ber(data_in,data_ref,skip)
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

19
Functions/calc_evm.m Normal file
View File

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

109
imddmodel.m
View File

@@ -1,19 +1,20 @@
function BER = imddmodel(dc_mu,sir,linewidth)
function output = imddmodel(sir,winlen,linewidth)
rng(2023);
for realiz = 1:3
rng(realiz);
%% Set Simulation Variables
sir = sir;
delay = 300; %mpi delay in meter
delay = 50; %mpi delay in meter
fiblen = 0.5; %main link in km
fiblen = 0; %main link in km
laser_linewidth = linewidth;
O = 18; %order of prbs
O = 17; %order of prbs
N = 2^(O-1); %length of prbs
[~,seed] = prbs(O,1); %initialize first seed of prbs
@@ -57,11 +58,11 @@ reflectionpoint = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amp
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);
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",0,"responsivity",1,"temperature",20);
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);
@@ -73,7 +74,7 @@ scp = Scope("fsimu",fdac*kover,"fadc",fadc,...
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,...
"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);
@@ -84,7 +85,7 @@ eq = EQ("K",2,"plottrain",0,"plotfinal",0,...
for i = 1:log2(M)
[bitpattern(:,i),seed] = prbs(O,N,seed);
end
bitpattern = [bitpattern ; flip(bitpattern)];
%bitpattern = prms_out';
% Build Inf. signal class
@@ -98,7 +99,6 @@ reference = mod_out;
% shape shape
X = pulseform.process(mod_out);
test = applyPulseShaping(reference.signal,fsym,fdac);
% AWG -> ELECTRICAL DOMAIN
awg_out = awg.process(X);
@@ -121,7 +121,7 @@ R = reflectionpoint.process(X);
% Propagate
R = reflectionprop.process(R);
disp(['SIR ',num2str(10*log10(X.power/R.power))]);
% disp(['SIR ',num2str(10*log10(X.power/R.power))]);
% Delay
R = R.delay("delay_meter",delay);
@@ -131,7 +131,7 @@ X = X+R;
%%
X = opticatten.process(X);
X = edfaamp.process(X);
% X = edfaamp.process(X);
% Photo Diode -> ELECTRICAL DOMAIN
X = phdiode.process(X);
@@ -145,29 +145,80 @@ X = scp.process(X);
X = X.resample("fs_out",2*fsym,"fs_in",fadc);
% Normalize
X = X.normalize("mode","rms");
Eq_in = X.normalize("mode","rms");
% Equalizer
X = eq.process(X,reference);
Eq_out = eq.process(Eq_in,reference);
%% MPI reduction
wl = 10000;
mpi_red = X;
mpi_red.signal = X.signal - 1/wl .* movsum( X.signal,[wl/2,wl/2]);
%% 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);
% Digi Demod
demap_out = digimod.demap(X);
demap_out_a1 = digimod.demap(mpi_red);
e = Eq_out.signal - pre_decision_level_bi;
% 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);
lp_mpi = Filter('filtdegree',1,"f_cutoff",max(10*linewidth,1e6),"fsamp",fsym,"filterType",filtertypes.bessel_inp);
disp(['BER: ', sprintf('%2E',ber), ' ERRORS: ' ,num2str(sum(errors))]);
disp(['BER: ', sprintf('%2E',ber_a1), ' ERRORS: ' ,num2str(sum(errors))]);
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))]);
end
output.ber_bm = mean(ber_bm);
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);
BER.ber = ber;
BER.ber_a1 = ber_a1;
end

209
run_loop.m
View File

@@ -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
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
%% 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 = InterX([sirloop;ber(dc,:,lw)],[sirloop;hdfec]);
thres(dc,lw) = -a(1);
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
contour(thres,'LineWidth',2);
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")

183
setup_simulation.m
View File

@@ -1,16 +1,25 @@
sir_loop = -30:5:-10;
sir_loop = -16;
for lp = 1:5
for lp = 1
dc_mu = 0.05;
sir = sir_loop(lp);
rkey = 13;
rng(rkey);
rng(10);
%% Set Simulation Variables
O = 18; %order of prbs
sir = sir;
delay = 40; %mpi delay in meter
fiblen = 0; %main link in km
laser_linewidth = 1e6;
O = 17; %order of prbs
N = 2^(O-1); %length of prbs
[~,seed] = prbs(O,1); %initialize first seed of prbs
@@ -38,7 +47,7 @@ fsimu = kover * fdac ;
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);
@@ -46,22 +55,16 @@ lp_laser = Filter('filtdegree',1,"f_cutoff",50e9,"fsamp",fdac*kover,"filterTy
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);
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);
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);
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);
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);
@@ -71,22 +74,20 @@ scp = Scope("fsimu",fdac*kover,"fadc",fadc,...
"adcresolution",6,"quantbuffer",0.1,'block_dc',1);
eq = EQ("K",2,"plottrain",0,"plotfinal",0,...
"training_length",4096,"training_loops",5,...
"training_length",2048,"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,...
"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
% 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
@@ -108,10 +109,11 @@ 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)));
X.signal = X.signal .* 1.3800;
% Laser; Modulation -> OPTICAL DOMAIN
X = extmodlaser.process(X);
[X,extmodlaser] = extmodlaser.process(X);
% Fiber Propagation
X = fib.process(X);
@@ -122,17 +124,17 @@ R = reflectionpoint.process(X);
% Propagate
R = reflectionprop.process(R);
% disp(['SIR ',num2str(10*log10(X.power/R.power))]);
% Delay
R = R.delay("delay_meter",200);
[R,n] = R.delay("delay_meter",delay);
% Add together
X = X+R;
X.spectrum(fsimu,"displayname",['SIR: ',num2str(sir)],"figurename",'spectrum')
%%
X = opticatten.process(X);
X = edfaamp.process(X);
% X = edfaamp.process(X);
% Photo Diode -> ELECTRICAL DOMAIN
X = phdiode.process(X);
@@ -146,19 +148,105 @@ X = scp.process(X);
X = X.resample("fs_out",2*fsym,"fs_in",fadc);
% Normalize
X = X.normalize("mode","rms");
Eq_in = X.normalize("mode","rms");
% Equalizer
eq_out = eq.process(X,reference);
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
demap_out = digimod.demap(eq_out);
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
[bits,errors,BER] = calc_ber(demap_out.signal(1:end-1,:),bitpattern(1:end,:),1000);
disp(['BER: ', sprintf('%2E',BER), ' ERRORS: ' ,num2str(sum(errors))]);
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
@@ -176,17 +264,17 @@ if 0
hold off
subplot(2,1,2)
hold on
stem(demap_out.signal(4150:4175,1),'DisplayName','Tx','Color',col(1,:),'LineStyle','-','LineWidth',5)
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 1
if 0
col = cbrewer2('Paired',8);
xax = 1:eq_out.length;
xax = 1:Eq_out.length;
figure(3)
sgtitle('')
@@ -197,9 +285,24 @@ if 1
% 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,5,lp)
scatter(xax,eq_out.signal,4,'.','MarkerEdgeColor',col(6,:),'DisplayName','After EQ');
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')
@@ -208,3 +311,5 @@ if 1
end
end