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

EQ Structures for MPI mitigation

Browse Source
This commit is contained in:
Silas
2024-08-15 09:40:59 +02:00
parent 34f9149346
commit 23f386da77
13 changed files with 880 additions and 135 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
Classes/01_transmit/PAMsource.m
View File

@@ -51,7 +51,7 @@ classdef PAMsource
warning('No Pulseformer given. Proceeding with RRC and alpha 0.05'); warning('No Pulseformer given. Proceeding with RRC and alpha 0.05');
options.pulseformer = Pulseformer("fsym",obj.fsym,"fdac",obj.fs_out,"pulse","rrc","pulselength",16,"rrcalpha",0.05); options.pulseformer = Pulseformer("fsym",obj.fsym,"fdac",obj.fs_out,"pulse","rrc","pulselength",16,"rrcalpha",0.05);
end end
end end
@@ -67,12 +67,12 @@ classdef PAMsource
if obj.useprbs if obj.useprbs
for i = 1:log2(obj.M) for i = 1:log2(obj.M)
[bitpattern(:,i),seed] = prbs(O,N,seed); [bitpattern(:,i),seed] = prbs(O,N,seed);
end end
else else
s = RandStream('twister','Seed',obj.randkey); s = RandStream('twister','Seed',obj.randkey);
for i = 1:log2(obj.M) for i = 1:log2(obj.M)
bitpattern(:,i) = randi(s,[0 1], N, 1); bitpattern(:,i) = randi(s,[0 1], N, 1);
end end
end end
if obj.M == 6 if obj.M == 6
@@ -102,7 +102,12 @@ classdef PAMsource
%%%%% Hard clip digital signal to PAM range before DAC %%%%%% %%%%% Hard clip digital signal to PAM range before DAC %%%%%%
if obj.applyclipping if obj.applyclipping
digi_sig.signal = clip(digi_sig.signal , sym_min * obj.clipfactor , sym_max * obj.clipfactor); try
digi_sig.signal = clip(digi_sig.signal , sym_min * obj.clipfactor , sym_max * obj.clipfactor);
catch
digi_sig.signal(digi_sig.signal > sym_max * obj.clipfactor) = sym_max * obj.clipfactor;
digi_sig.signal(digi_sig.signal < sym_min * obj.clipfactor) = sym_min * obj.clipfactor;
end
end end
% append to logbook % append to logbook

1
Classes/03_receive/Photodiode.m
View File

@@ -72,6 +72,7 @@ classdef Photodiode
% Thermal Noise % Thermal Noise
therm_current_psd = (2 * k * T / R ) ; %squared therm_current_psd = (2 * k * T / R ) ; %squared
therm_current_psd = 1.8e-11^2;
Bw = obj.fsimu; Bw = obj.fsimu;
therm_noise_pow = therm_current_psd * Bw; %squared therm_noise_pow = therm_current_psd * Bw; %squared

179
Classes/04_DSP/FFE_DCremoval.m Normal file
View File

@@ -0,0 +1,179 @@
classdef FFE_DCremoval < handle
% Implementation of plain and simple FFE.
% 1) Training mode (stable performance when you use NLMS)
% 2) Decision directed mode
% Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0);
properties
sps % usually 2
order
e
error
len_tr
mu_tr
epochs_tr
mu_dd
epochs_dd
mu_dc
constellation
decide
end
methods
function obj = FFE_DCremoval(options)
arguments(Input)
options.sps = 2;
options.order = 15;
options.len_tr = 4096;
options.mu_tr = 0;
options.epochs_tr = 5;
options.mu_dd = 1e-5;
options.epochs_dd = 5;
options.mu_dc = 0.05;
options.decide = false;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
obj.e = zeros(obj.order,1);
obj.error = 0;
end
function [X] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
X = X.normalize("mode","rms");
obj.constellation = unique(D.signal);
% Training Mode
training = 1;
showviz = 0;
obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz);
% Decision Directed Mode
N = X.length;
training = 0;
showviz = 0;
[signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz);
% Output Signal
if obj.decide
X.signal = decision;
else
X.signal = signal;
end
X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym
lbdesc = [num2str(obj.order),' tap FFE'];
X = X.logbookentry(lbdesc); % append to logbook
end
function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz)
arguments
obj
x
d
mio
epochs
N
training
showviz
end
x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)];
if showviz
f = figure(111);
subplot(2,2,1:2);
hold on
a = scatter(1:numel(x),x,1,'.');
a2 = scatter(1,1,1,'.');
a3 = scatter(1,1,2,'.');
a4 = xline(1);
ylim([-3 3])
xlim([0 length(x)]);
% subplot(2,2,3)
% dplot = x(1:1+500);
% b = scatter(1:numel(dplot),dplot,5,'x');
% xline(1)
% xline(obj.order)
% ylim([-3 3])
% xlim([0 500]);
subplot(2,2,3:4)
c = stem(obj.e);
ylim([-1 1])
drawnow
end
err = 0;
e_dc = 0;
for epoch = 1 : epochs
symbol = 0;
for sample = 1 : obj.sps : N
symbol = symbol+1;
U = x(obj.order+sample-1:-1:sample);
y(symbol,1) = e_dc + obj.e.' * U; % Calculating output of LMS __ * |
if training
d_hat(symbol,1) = d(symbol);
else
[~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point
d_hat(symbol,1) = obj.constellation(symbol_idx);
end
err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error
if mio ~= 0
obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS
else
normalizationfactor = (U.' * U);
obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS
end
if mod(sample,100) == 1 && showviz
a2.XData = 1:2*numel(y);
a2.YData = repelem(y, 2);
a3.XData = 1:2*numel(d_hat);
a3.YData = repelem(d_hat, 2);
a4.Value = sample;
% b.YData = x(symbol:symbol+500);
c.YData = obj.e;
drawnow;
end
e_dc = e_dc - obj.mu_dc * err(symbol);
obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error
end
end
end
end
end

193
Classes/04_DSP/FFE_FFDCAVG.m Normal file
View File

@@ -0,0 +1,193 @@
classdef FFE_FFDCAVG < handle
% Implementation of plain and simple FFE.
% 1) Training mode (stable performance when you use NLMS)
% 2) Decision directed mode
% Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0);
properties
sps % usually 2
order
e
error
len_tr
mu_tr
epochs_tr
mu_dd
epochs_dd
mu_buff
constellation
decide
end
methods
function obj = FFE_FFDCAVG(options)
arguments(Input)
options.sps = 2;
options.order = 15;
options.len_tr = 4096;
options.mu_tr = 0;
options.epochs_tr = 5;
options.mu_dd = 1e-5;
options.epochs_dd = 5;
options.mu_buff = 0;
options.decide = false;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
obj.e = zeros(obj.order,1);
obj.error = 0;
end
function [X] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
X = X.normalize("mode","rms");
obj.constellation = unique(D.signal);
% Training Mode
training = 1;
showviz = 0;
obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz);
% Decision Directed Mode
N = X.length;
training = 0;
showviz = 0;
[signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz);
% Output Signal
if obj.decide
X.signal = decision;
else
X.signal = signal;
end
X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym
lbdesc = [num2str(obj.order),' tap FFE'];
X = X.logbookentry(lbdesc); % append to logbook
end
function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz)
arguments
obj
x
d
mio
epochs
N
training
showviz
end
x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)];
if showviz
f = figure(111);
subplot(2,2,1:2);
hold on
a = scatter(1:numel(x),x,1,'.');
a2 = scatter(1,1,1,'.');
a3 = scatter(1,1,2,'.');
a4 = xline(1);
ylim([-3 3])
xlim([0 length(x)]);
% subplot(2,2,3)
% dplot = x(1:1+500);
% b = scatter(1:numel(dplot),dplot,5,'x');
% xline(1)
% xline(obj.order)
% ylim([-3 3])
% xlim([0 500]);
subplot(2,2,3:4)
c = stem(obj.e);
ylim([-1 1])
drawnow
end
for epoch = 1 : epochs
symbol = 0;
err_buffer = zeros(numel(obj.constellation),50);
for sample = 1 : obj.sps : N
symbol = symbol+1;
U = x(obj.order+sample-1:-1:sample);
y(symbol,1) = obj.e.' * U; % Calculating output of LMS __ * |
if training
[~,symbol_idx] = min(abs(d(symbol) - obj.constellation)); % decision for closest constellation point
d_hat(symbol,1) = d(symbol);
else
[~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point
d_hat(symbol,1) = obj.constellation(symbol_idx);
end
err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error
err_buffer(symbol_idx,1) = err(symbol);
err_buffer(symbol_idx,:) = circshift(err_buffer(symbol_idx,:),1);
y(symbol) = y(symbol) - obj.mu_buff*mean(err_buffer(symbol_idx,:));
if training
[~,symbol_idx] = min(abs(d(symbol) - obj.constellation)); % decision for closest constellation point
d_hat(symbol,1) = d(symbol);
else
[~,symbol_idx] = min(abs(y(symbol) - obj.constellation)); % decision for closest constellation point
d_hat(symbol,1) = obj.constellation(symbol_idx);
end
err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error
if mio ~= 0
obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS
else
normalizationfactor = (U.' * U);
obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS
end
if mod(sample,100) == 1 && showviz
a2.XData = 1:2*numel(y);
a2.YData = repelem(y, 2);
a3.XData = 1:2*numel(d_hat);
a3.YData = repelem(d_hat, 2);
a4.Value = sample;
% b.YData = x(symbol:symbol+500);
c.YData = obj.e;
drawnow;
end
obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error
end
end
end
end
end

149
Classes/04_DSP/FFE_adaptive_decision.m Normal file
View File

@@ -0,0 +1,149 @@
classdef FFE_adaptive_decision < handle
% Implementation of plain and simple FFE.
% 1) Training mode (stable performance when you use NLMS)
% 2) Decision directed mode
% Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0);
properties
sps % usually 2
order
e
error
len_tr
mu_tr
epochs_tr
mu_dd
epochs_dd
buffer_length
constellation
decide
end
methods
function obj = FFE_adaptive_decision(options)
arguments(Input)
options.sps = 2;
options.order = 15;
options.len_tr = 4096;
options.mu_tr = 0;
options.epochs_tr = 5;
options.mu_dd = 1e-5;
options.epochs_dd = 5;
options.buffer_length = 100;
options.decide = false;
end
fn = fieldnames(options);
for n = 1:numel(fn)
obj.(fn{n}) = options.(fn{n});
end
obj.e = zeros(obj.order,1);
obj.error = 0;
end
function [X] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS
X = X.normalize("mode","rms");
obj.constellation = unique(D.signal);
% Training Mode
training = 1;
showviz = 0;
obj.equalize(X.signal, D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz);
% Decision Directed Mode
N = X.length;
training = 0;
showviz = 0;
[signal,decision]=obj.equalize(X.signal, D.signal,obj.mu_dd,obj.epochs_dd,N,training,showviz);
% Output Signal
if obj.decide
X.signal = decision;
else
X.signal = signal;
end
X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym
lbdesc = [num2str(obj.order),' tap FFE'];
X = X.logbookentry(lbdesc); % append to logbook
end
function [y,d_hat] = equalize(obj,x,d,mio,epochs,N,training,showviz)
arguments
obj
x
d
mio
epochs
N
training
showviz
end
x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)];
for epoch = 1 : epochs
symbol = 0;
% y_buffer = zeros(numel(obj.constellation),500);
y_buffer = repmat(obj.constellation,1,obj.buffer_length);
for sample = 1 : obj.sps : N
symbol = symbol+1;
U = x(obj.order+sample-1:-1:sample);
y(symbol,1) = obj.e.' * U; % Calculating output of LMS __ * |
if training
[~,symbol_idx] = min(abs(d(symbol) - obj.constellation)); % decision for closest constellation point
d_hat(symbol,1) = d(symbol);
else
y_buffer(y_buffer==0) = NaN;
adap_constellation = mean(y_buffer,2,"omitnan");
[~,symbol_idx] = min(abs(y(symbol) - adap_constellation)); % decision for closest constellation point
d_hat(symbol,1) = obj.constellation(symbol_idx);
end
y_buffer(symbol_idx,1) = y(symbol);
y_buffer(symbol_idx,:) = circshift(y_buffer(symbol_idx,:),1);
err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error
if mio ~= 0
obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS
else
normalizationfactor = (U.' * U);
obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS
end
obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error
end
end
end
end
end

42
projects/MPI_August/EQ_development.m Normal file
View File

@@ -0,0 +1,42 @@
Scope_sig = load("projects/MPI_August/scpe_sig.mat","Scpe_sig");
Scope_sig = Scope_sig.Scpe_sig;
Symbols = load("projects/MPI_August/symbols.mat","Symbols");
Symbols = Symbols.Symbols;
Bits = load("projects/MPI_August/bits.mat","Bits");
Bits = Bits.Bits;
mdc = [1e-3,1e-2,2e-2,3e-2,4e-2,5e-2,1e-1];
for m = 1:numel(mdc)
Eq = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"mu_dc",mdc(m));
[EQ_sig] = Eq.process(Scope_sig,Symbols);
%%%%% DEMAP %%%%%%
Rx_bits = PAMmapper(4,0).demap(EQ_sig);
% BER
[~,errors_bm,ber(m),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
end
figure()
scatter(mdc,ber,20);
yline(3.8e-3,'DisplayName','HD-FEC');
xlabel('Mu DC');
ylabel('Bit Error Rate (BER)');
title('Bit Error Rate vs. Mu DC');
set(gca,'yscale','log');
set(gca,'xscale','log');
grid on;
legend
%%%% Look at Pam levels %%%%%
if 1
a = PAMmapper(4,0).separate_pamlevels(EQ_sig);
figure(14);hold on;scatter(1:EQ_sig.length,a,1,'.');
end

56
projects/MPI_August/EQ_development_ffe_adap_decision.m Normal file
View File

@@ -0,0 +1,56 @@
Scope_sig = load("projects/MPI_August/scpe_sig.mat","Scpe_sig");
Scope_sig = Scope_sig.Scpe_sig;
Symbols = load("projects/MPI_August/symbols.mat","Symbols");
Symbols = Symbols.Symbols;
Bits = load("projects/MPI_August/bits.mat","Bits");
Bits = Bits.Bits;
lenbuff = [80];
for m = 1:numel(lenbuff)
Eq = FFE_adaptive_decision("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"buffer_length",lenbuff(m));
[EQ_sig] = Eq.process(Scope_sig,Symbols);
%%%%% DEMAP %%%%%%
Rx_bits = PAMmapper(4,0).demap(EQ_sig);
% BER
[~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
disp(['BER: ',sprintf('%.1E',ber)]);
end
% fs = 2*112e9;
% N = 2^(nextpow2(length(Eq.error(1,:)))-6);
% [p_lin,w] = pwelch(Eq.error(1,:),hanning(N),N/2,N,fs,"centered","power","mean");
% p_dbm = 10*log10(p_lin)+30; %dB to dBm in case of "power"
% figure(123); % If figure does not exist, create new figure
% hold on
% plot(w.*1e-9,p_dbm,'DisplayName',['bla'],'LineWidth',1);
% xlabel("Frequency in GHz");
% %ylabel("Power/frequency (dB/Hz)");
% ylabel("Power (dBm)");
% xlim([-fs/2 fs/2].*1e-9)
% edgetick = 2^(nextpow2(fs*1e-9));
% xticks([-edgetick:16:edgetick]);
% xlim([-244, 244])
% ylim([-120,-0]);
% yticks([-200:10:10]);
% legend
%%%% Look at Pam levels %%%%%
if 1
a = PAMmapper(4,0).separate_pamlevels(EQ_sig);
figure(15);hold on;scatter(1:EQ_sig.length,a,1,'.');
end

56
projects/MPI_August/EQ_development_ffe_dc_avg.m Normal file
View File

@@ -0,0 +1,56 @@
Scope_sig = load("projects/MPI_August/scpe_sig.mat","Scpe_sig");
Scope_sig = Scope_sig.Scpe_sig;
Symbols = load("projects/MPI_August/symbols.mat","Symbols");
Symbols = Symbols.Symbols;
Bits = load("projects/MPI_August/bits.mat","Bits");
Bits = Bits.Bits;
mubuff = [0.7];
for m = 1:numel(mubuff)
Eq = FFE_FFDCAVG("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"mu_buff",mubuff(m));
[EQ_sig] = Eq.process(Scope_sig,Symbols);
%%%%% DEMAP %%%%%%
Rx_bits = PAMmapper(4,0).demap(EQ_sig);
% BER
[~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
disp(['BER: ',sprintf('%.1E',ber)]);
end
% fs = 2*112e9;
% N = 2^(nextpow2(length(Eq.error(1,:)))-6);
% [p_lin,w] = pwelch(Eq.error(1,:),hanning(N),N/2,N,fs,"centered","power","mean");
% p_dbm = 10*log10(p_lin)+30; %dB to dBm in case of "power"
% figure(123); % If figure does not exist, create new figure
% hold on
% plot(w.*1e-9,p_dbm,'DisplayName',['bla'],'LineWidth',1);
% xlabel("Frequency in GHz");
% %ylabel("Power/frequency (dB/Hz)");
% ylabel("Power (dBm)");
% xlim([-fs/2 fs/2].*1e-9)
% edgetick = 2^(nextpow2(fs*1e-9));
% xticks([-edgetick:16:edgetick]);
% xlim([-244, 244])
% ylim([-120,-0]);
% yticks([-200:10:10]);
% legend
%%%% Look at Pam levels %%%%%
if 1
a = PAMmapper(4,0).separate_pamlevels(EQ_sig);
figure(15);hold on;scatter(1:EQ_sig.length,a,1,'.');
end

35
projects/MPI_August/auswertung/rop_curve.m Normal file
View File

@@ -0,0 +1,35 @@
M = wh.parameter.M.values(1);
datarate = wh.parameter.datarate.values(1);
sir = wh.parameter.sir.values(1);
laser_linewidth = wh.parameter.laser_linewidth.values(1);
pn_key = wh.parameter.pn_key.values;
rop = wh.parameter.rop.values;
cols = linspecer(4);
cnt = 0;
for pnk = pn_key
cnt = cnt+1;
ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,laser_linewidth,pnk,rop);
ber_dcavg = wh.getStoValue('ber_dcavg',M,datarate,sir,laser_linewidth,pnk,rop);
ber_adapt = wh.getStoValue('ber_adapt',M,datarate,sir,laser_linewidth,pnk,rop);
ber_derem = wh.getStoValue('ber_dcrem',M,datarate,sir,laser_linewidth,pnk,rop);
% Create the initial plot
figure(43);
hold on; % Retain the plot so new points can be added without complete redraw
plot(rop,ber_ffe',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber ffe",'Color',cols(1,:));
plot(rop,ber_dcavg',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber dcavg",'Color',cols(2,:));
plot(rop,ber_adapt',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber adapt",'Color',cols(3,:));
plot(rop,ber_derem',"LineWidth",1,"LineStyle","-","Marker",".","MarkerSize",10,"DisplayName","ber derem",'Color',cols(4,:));
end
yline(3.8e-3,'DisplayName','HD-FEC');
xlabel('Signal to Interference Ratio (dB)');
ylabel('Bit Error Rate (BER)');
title('Bit Error Rate vs. ROP');
set(gca,'yscale','log');
grid on;
legend

BIN
projects/MPI_August/bits.mat Normal file
View File

Binary file not shown.

291
projects/MPI_August/mpi_uncorrelated.m
View File

@@ -2,25 +2,21 @@
%% Parameter to simulate and save %% Parameter to simulate and save
params = struct; params = struct;
params.M = [4]; params.M = [8];
params.datarate = [224]; params.datarate = [224];
params.sir = [24]; %decibel = attenuation of interference path params.sir = [30]; %decibel = attenuation of interference path
params.laser_linewidth = [1e6]; params.laser_linewidth = [1e6];
params.pn_key = [11]; params.pn_key = [1];
params.vbias_rel = [0.5]; params.rop = [-12:0];
params.rop = -5;
params.clipfactor = [1.5];
params.rrcalpha = [0.1];
name = ['wh_',strrep(num2str(now),'.','')]; name = ['wh_',strrep(num2str(now),'.','')];
wh = DataStorage(params); wh = DataStorage(params);
wh.addStorage("ber"); wh.addStorage("ber_ffe");
wh.addStorage("rop"); wh.addStorage("ber_dcavg");
wh.addStorage("txpapr"); wh.addStorage("ber_adapt");
wh.addStorage("er"); wh.addStorage("ber_dcrem");
%% Init Params %% Init Params
link_length = 10000; %meter link_length = 10000; %meter
@@ -32,7 +28,9 @@ disp(['Start Simulation of ',num2str(endcnt),' loops...'])
tic tic
for M = wh.parameter.M.values for M = wh.parameter.M.values
for datarate = wh.parameter.datarate.values for datarate = wh.parameter.datarate.values
for rrcalpha = wh.parameter.rrcalpha.values for pn_key = wh.parameter.pn_key.values
%% SETUP HERE: %% %% SETUP HERE: %%
kover = 8; kover = 8;
@@ -44,157 +42,150 @@ for M = wh.parameter.M.values
% MAIN SIGNAL % MAIN SIGNAL
%%%%% Symbol Generation %%%%%% %%%%% Symbol Generation %%%%%%
[Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",2).process(); [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key).process();
%%%%% AWG %%%%%% %%%%% AWG %%%%%%
El_sig = M8199.process(Digi_sig); El_sig = M8199.process(Digi_sig);
El_sig = El_sig.*0.7222; El_sig = El_sig.*0.7222;
El_sig.signal = awgn(El_sig.signal,20,'measured',1); El_sig.signal = awgn(El_sig.signal,20,'measured',pn_key);
%%%%% Lowpass before Modulator %%%%%% %%%%% Lowpass before Modulator %%%%%%
El_sig = Filter('filtdegree',2,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); El_sig = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig);
El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",15).process(El_sig); El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",10).process(El_sig);
% INTERFERENCE SIGNAL % INTERFERENCE SIGNAL
%%%%% Symbol Generation %%%%%% %%%%% Symbol Generation %%%%%%
[Digi_sig_i,Symbols_i,Bits_i] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",1).process(); [Digi_sig_i,Symbols_i,Bits_i] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",0,"fs_out",M8199.fdac,"applyclipping",1,"clipfactor",1.3,"applypulseform",1,"pulseformer",Pform,"randkey",pn_key*2).process();
%%%%% AWG %%%%%% %%%%% AWG %%%%%%
El_sig_i = M8199.process(Digi_sig_i); El_sig_i = M8199.process(Digi_sig_i);
El_sig_i = El_sig_i.*0.7222; El_sig_i = El_sig_i.*0.7222;
El_sig_i.signal = awgn(El_sig_i.signal,20,'measured',2); El_sig_i.signal = awgn(El_sig_i.signal,20,'measured',pn_key*2);
%%%%% Lowpass before Modulator %%%%%% %%%%% Lowpass before Modulator %%%%%%
El_sig_i = Filter('filtdegree',2,"f_cutoff",100e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig_i); El_sig_i = Filter('filtdegree',2,"f_cutoff",60e9,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig_i);
El_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",15).process(El_sig_i); El_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",10).process(El_sig_i);
for laser_linewidth = wh.parameter.laser_linewidth.values for laser_linewidth = wh.parameter.laser_linewidth.values
for pn_key = wh.parameter.pn_key.values
for vbias_rel = wh.parameter.vbias_rel.values % MAIN SIGNAL
%%%%% MODULATE E/O CONVERSION %%%%%%
vbias_rel = 0.5;
u_pi = 2.9;
vbias = -vbias_rel*u_pi;
[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key).process(El_sig);
Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true);
Opt_sig = Optfilter.process(Opt_sig);
er = Opt_sig.extinctionratio(fsym,M);
% INTERFERENCE SIGNAL
%%%%% MODULATE E/O CONVERSION %%%%%%
u_pi = 2.9;
vbias = -vbias_rel*u_pi;
[Opt_sig_i] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key+1).process(El_sig_i);
Opt_sig_i = Optfilter.process(Opt_sig_i);
% MAIN SIGNAL j_ = wh.parameter.sir.length;
%%%%% MODULATE E/O CONVERSION %%%%%% i_ = wh.parameter.rop.length;
u_pi = 2.9; ber=zeros(j_,i_);
vbias = -vbias_rel*u_pi; patten=zeros(j_,i_);
[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key).process(El_sig);
Optfilter = Filter('filtdegree',6,"f_cutoff",fsym.*0.7,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true); for j = 1:j_
Opt_sig = Optfilter.process(Opt_sig); sir = wh.parameter.sir.values(j);
% INTERFERENCE SIGNAL %%%%% Interference Signal Fiber Prop %%%%%%
%%%%% MODULATE E/O CONVERSION %%%%%% Opt_sig_i = Fiber("fsimu",Opt_sig_i.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig_i);
u_pi = 2.9;
vbias = -vbias_rel*u_pi; Opt_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",Opt_sig.power-sir).process(Opt_sig_i);
[Opt_sig_i] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",1290,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",pn_key+1).process(El_sig_i);
Opt_sig_i = Optfilter.process(Opt_sig_i); %%%%% ADD Interference and Main Signal %%%%%%
Opt_sig = Opt_sig_i + Opt_sig;
%%%%% Interference Signal Fiber Prop %%%%%%
Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig);
% % MPI Channel
% Opt = channel_model_mpi(Opt_sig,link_length,mpi_path,sir);
j_ = wh.parameter.sir.length; % Receiver ROP curve
i_ = wh.parameter.rop.length; parfor i = 1:i_
ber=zeros(j_,i_); rop=wh.parameter.rop.values(i);
patten=zeros(j_,i_);
for j = 1:j_ % Set ROP
sir = wh.parameter.sir.values(j); Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig);
patten(j,i) = Rx_sig.power;
%%%%%% Square Law %%%%%%
Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20).process(Rx_sig);
%%%%%% Lowpass PhDiode %%%%%%
Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig);
% Rx_sig.spectrum("displayname","Received Signal after PhD","fignum",201);
%%%%% Interference Signal Fiber Prop %%%%%% %%%%%% Scope %%%%%%
Opt_sig_i = Fiber("fsimu",Opt_sig_i.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig_i); fadc = 256e9;
Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,...
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
"adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig);
%%%%%% Sample to 2x fsym %%%%%%
Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym);
Opt_sig_i = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",Opt_sig.power-sir).process(Opt_sig_i); %%%%%% Sync Rx signal with reference %%%%%%
[Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym);
%%%%% ADD Interference and Main Signal %%%%%% %%%%% EQUALIZE %%%%%%
Opt_sig = Opt_sig_i + Opt_sig; Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1);
[EQ_sig] = Eq.process(Scpe_sig,Symbols);
Rx_bits = PAMmapper(M,0).demap(EQ_sig);
[~,errors_bm,ber_ffe(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
disp(['BER: ',sprintf('%.1E',ber_ffe(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']);
%%%%% Interference Signal Fiber Prop %%%%%% Eq = FFE_FFDCAVG("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"mu_buff",0.7);
Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); [EQ_sig] = Eq.process(Scpe_sig,Symbols);
Rx_bits = PAMmapper(M,0).demap(EQ_sig);
[~,errors_bm,ber_dcavg(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
disp(['BER: ',sprintf('%.1E',ber_dcavg(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']);
Eq = FFE_adaptive_decision("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"buffer_length",85); [EQ_sig] = Eq.process(Scpe_sig,Symbols);
Rx_bits = PAMmapper(M,0).demap(EQ_sig);
[~,errors_bm,ber_adapt(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
disp(['BER: ',sprintf('%.1E',ber_adapt(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']);
% % MPI Channel Eq = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",1,"mu_dc",0.07);
% Opt = channel_model_mpi(Opt_sig,link_length,mpi_path,sir); [EQ_sig] = Eq.process(Scpe_sig,Symbols);
Rx_bits = PAMmapper(M,0).demap(EQ_sig);
% Receiver ROP curve [~,errors_bm,ber_dcrem(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
for i = 1:i_ disp(['BER: ',sprintf('%.1E',ber_dcrem(j,i)),' - - ROP: ',num2str(patten(j,i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']);
rop=wh.parameter.rop.values(i);
% Set ROP
Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig);
patten(j,i) = Rx_sig.power;
%%%%%% Square Law %%%%%%
Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20).process(Rx_sig);
%%%%%% Lowpass PhDiode %%%%%%
Rx_sig = Filter('filtdegree',2,"f_cutoff",70e9,"fs",fdac*kover,"filterType",filtertypes.gaussian,"active",true).process(Rx_sig);
% Rx_sig.spectrum("displayname","Received Signal after PhD","fignum",201);
%%%%%% Scope %%%%%%
fadc = 256e9;
Lp_scpe = Filter('filtdegree',4,"f_cutoff",63e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,...
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,...
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,...
"adcresolution",16,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig);
Scpe_sig.plot("displayname","SIgnal after Scope","fignum",999);
%%%%%% Sample to 2x fsym %%%%%%
Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym);
%%%%%% Sync Rx signal with reference %%%%%%
[Scpe_sig,D,cuts] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym);
%%%%% EQUALIZE %%%%%%
% Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0);
% Eq = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",25,"dfe_order",2,"sps",2,"decide",1);
% Eq = EQ("Ne",[25,3,3],"Nb",[0,0,0],"training_length",4096,"training_loops",4,"dd_loops",4,"K",2,"DCmu",0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",1);
Eq = VNLE("epochs_tr",7,"epochs_dd",7,"len_tr",4096,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[25,3,3],"sps",2,"decide",0);
[EQ_sig] = Eq.process(Scpe_sig,Symbols);
% EQ_sig.normalize("mode","rms").plot('fignum',23,'displayname','before eq')
%%%%% DEMAP %%%%%%
Rx_bits = PAMmapper(M,0).demap(EQ_sig);
%%%% Look at Pam levels %%%%%
if 1
a = PAMmapper(M,0).separate_pamlevels(EQ_sig);
figure(14);hold on;scatter(1:EQ_sig.length,a,1,'.');
end
% BER
[~,errors_bm,ber(j,i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
cnt = cnt+1;
disp(['BER: ',sprintf('%.1E',ber(j,i)),' - - ROP: ',num2str(Rx_sig.power),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']);
end
end
for j = 1:j_
sir = wh.parameter.sir.values(j);
for i = 1:i_
rop=wh.parameter.rop.values(i);
wh.addValueToStorage(ber(j,i),'ber',M,datarate,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha);
wh.addValueToStorage(patten(j,i),'rop',M,datarate,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha);
wh.addValueToStorage(er,'er',M,datarate,sir,laser_linewidth,pn_key,vbias_rel,rop,clipfactor,rrcalpha);
end
end
end end
toc
disp(['Simulated: ',num2str(cnt/endcnt*100),' %']);
save(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\',name,'.mat'],"wh");
end end
for j = 1:j_
sir = wh.parameter.sir.values(j);
for i = 1:i_
rop=wh.parameter.rop.values(i);
wh.addValueToStorage(ber_ffe(j,i) ,'ber_ffe',M,datarate,sir,laser_linewidth,pn_key,rop);
wh.addValueToStorage(ber_dcavg(j,i),'ber_dcavg',M,datarate,sir,laser_linewidth,pn_key,rop);
wh.addValueToStorage(ber_adapt(j,i),'ber_adapt',M,datarate,sir,laser_linewidth,pn_key,rop);
wh.addValueToStorage(ber_dcrem(j,i),'ber_dcrem',M,datarate,sir,laser_linewidth,pn_key,rop);
end
end
toc
disp(['Simulated: ',num2str(cnt/endcnt*100),' %']);
save(['C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_Juni\',name,'.mat'],"wh");
end end
end end
end end
@@ -202,5 +193,43 @@ end
M = wh.parameter.M.values(1);
datarate = wh.parameter.datarate.values(1);
sir = wh.parameter.sir.values(1);
laser_linewidth = wh.parameter.laser_linewidth.values(1);
pn_key = wh.parameter.pn_key.values;
rop = wh.parameter.rop.values;
cols = linspecer(4);
cnt = 0;
for pnk = pn_key
cnt = cnt+1;
ber_ffe = wh.getStoValue('ber_ffe',M,datarate,sir,laser_linewidth,pnk,rop);
ber_dcavg = wh.getStoValue('ber_dcavg',M,datarate,sir,laser_linewidth,pnk,rop);
ber_adapt = wh.getStoValue('ber_adapt',M,datarate,sir,laser_linewidth,pnk,rop);
ber_derem = wh.getStoValue('ber_dcrem',M,datarate,sir,laser_linewidth,pnk,rop);
% Create the initial plot
figure(43);
hold on; % Retain the plot so new points can be added without complete redraw
plot(rop,ber_ffe',"LineWidth",1,"LineStyle",":","Marker",".","MarkerSize",10,"DisplayName","ber ffe",'Color',cols(1,:));
plot(rop,ber_dcavg',"LineWidth",1,"LineStyle",":","Marker",".","MarkerSize",10,"DisplayName","ber dcavg",'Color',cols(2,:));
plot(rop,ber_adapt',"LineWidth",1,"LineStyle",":","Marker",".","MarkerSize",10,"DisplayName","ber adapt",'Color',cols(3,:));
plot(rop,ber_derem',"LineWidth",1,"LineStyle",":","Marker",".","MarkerSize",10,"DisplayName","ber derem",'Color',cols(4,:));
end
yline(3.8e-3,'DisplayName','HD-FEC');
xlabel('Signal to Interference Ratio (dB)');
ylabel('Bit Error Rate (BER)');
title('Bit Error Rate vs. SIR (MPI)');
set(gca,'yscale','log');
grid on;
legend

BIN
projects/MPI_August/scpe_sig.mat Normal file
View File

Binary file not shown.

BIN
projects/MPI_August/symbols.mat Normal file
View File

Binary file not shown.