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

Merge branch 'main' of cau-git.rz.uni-kiel.de:nt/mitarbeiter/silas/imdd_simulation

Browse Source
This commit is contained in:
Silas Oettinghaus
2025-03-10 12:51:36 +01:00
parent 77ac44536d ea8a43b2ae
commit 15fef242a4
24 changed files with 533 additions and 140 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
Classes/00_signals/Signal.m
View File

@@ -698,18 +698,19 @@ classdef Signal
if numel(shifts) > 0
%Cut occurences of ref signal from signal (only positive shifts)
if all(sign(co(pkpos)))
isFlipped = 1;
end
for c = shifts(shifts>=0)
sig = obj.delay(-c,'mode','samples');
sig.signal = sig.signal(1:length(b));
sig.signal = sig.signal(1:length(b)).*-isFlipped;
S{end+1,1} = sig;
end
%
if all(sign(co(pkpos)))
isFlipped = 1;
end
%return/keep the sinal with the highest correlation (only within positive shifts)
[~,idx]=max(pks(shifts>=0));

56
Classes/01_transmit/PAMsource.m
View File

@@ -71,20 +71,17 @@ classdef PAMsource
function [digi_sig,symbols,bits] = process(obj)
%%%%% PRBS Generation in correct shape for Modulation Format %%%%%%
O = obj.order; %order of prbs
N = 2^(O-1); %length of prbs
[~,seed] = prbs(O,1); %initialize first seed of prbs
bitpattern=[];
if obj.useprbs
% for i = 1:log2(obj.M)
% [bitpattern(:,i),seed] = prbs(O,N,seed);
% end
% O = obj.order; %order of prbs
% N = 2^(O-1); %length of prbs
% [~,seed] = prbs(O,1); %initialize first seed of prbs
% % for i = 1:log2(obj.M)
% % [bitpattern(:,i),seed] = prbs(O,N,seed);
% % end
%%%%% MOVE-IT PRMS %%%%
state = struct();
para = struct();
@@ -120,6 +117,7 @@ classdef PAMsource
else
s = RandStream('twister','Seed',obj.randkey);
for i = 1:log2(obj.M)
N = 2^(obj.order-1); %length of prbs
bitpattern(:,i) = randi(s,[0 1], N, 1);
end
end
@@ -155,24 +153,46 @@ classdef PAMsource
sym_max = max(symbols.signal);
end
% symbols.spectrum("fignum",111,"displayname","1) RAW SIGNAL");
pulsf = Moveit_wrapper('pulsef');
pulsf.para.alpharacos = obj.pulseformer.alpha;
pulsf.para.f_sym = obj.fsym;
pulsf.para.fs = obj.fsym;
pulsf.para.pulse = 'racos';
pulsf.process(symbols);
%%%%% Pulse-forming %%%%%%
if obj.applypulseform
%%% MY CODE
digi_sig = obj.pulseformer.process(symbols);
%%% MOVEIT WRAPPER BETA
% symbols.spectrum("fignum",111,"displayname","1) RAW SIGNAL");
%
% pulsf = Moveit_wrapper('pulsef');
% pulsf.para.alpharacos = obj.pulseformer.alpha;
% pulsf.para.f_sym = obj.fsym;
% pulsf.para.fs = obj.fsym;
% pulsf.para.pulse = 2;
%
% digi_sig = pulsf.process(symbols);
% Design raised cosine filter with given order in symbols
% digi_sig = symbols;
% sps = 4;
% nsym = 128;
% rctFilt3 = comm.RaisedCosineTransmitFilter(...
% Shape='Square root', ...
% RolloffFactor=1, ...
% FilterSpanInSymbols=nsym, ...
% OutputSamplesPerSymbol=sps);
%
% digi_sig.signal = rctFilt3([symbols.signal; zeros(nsym/2,1)]);
% % Correct for propagation delay by removing filter transients
% fltDelay = nsym / (2*obj.fsym);
% digi_sig.signal = digi_sig.signal(fltDelay*sps*obj.fsym+1:end);
%
% digi_sig.fs = sps .* symbols.fs;
else
digi_sig = symbols;
end
% digi_sig.spectrum("fignum",112,"displayname","2) after pulseforming");
%%%%% Re-sample to f DAC %%%%%%
n = 10;

17
Classes/01_transmit/Pulseformer.m
View File

@@ -50,7 +50,7 @@ classdef Pulseformer
signalclass_in = signalclass_in.logbookentry(lbdesc);
% write fs to signal
signalclass_in.fs = obj.fdac;
signalclass_in.fs = obj.fdac .* (obj.fdac./obj.fsym);
% write to output
signalclass_out = signalclass_in;
@@ -98,14 +98,14 @@ classdef Pulseformer
%Bau das Filter (hier rc)
racos_len = obj.pulselength*2;
h = rcosdesign(obj.alpha,racos_len,sps,filtertype);
h = h./ max(h);
% h = h./ max(h);
if obj.matched
h = 1./h;
h = conj(fliplr(h));
end
manual_cyclic_convolution = 1;
upfirdn_convolution = 0;
manual_cyclic_convolution = 0;
upfirdn_convolution = 1;
if manual_cyclic_convolution
@@ -131,14 +131,13 @@ classdef Pulseformer
if upfirdn_convolution
%Apply Filter using Matlab build in fctn.
h = rcosdesign(obj.alpha,racos_len,sps,"normal");
h = h./ max(h);
data_out_ = upfirdn(data_in,h,up,dn);
%cut signal, which is longer due to fir filter
st = round(up/dn*racos_len/2)+1; %we need to cut y_out
st = round(up/dn*racos_len/2); %we need to cut y_out
en = round(st + (length(data_in)*up/dn));
data_out_ = data_out_(st:en);
data_out = data_out_(st:en);
end

2
Classes/Moveit_wrapper.m
View File

@@ -110,7 +110,7 @@ classdef Moveit_wrapper < handle
previousLevel = currentLevel;
% Determine input type from obj.comment.type
if
if 1
inputType = obj.comment.type.(fieldName);
else
inputType = 'number'; % Default type

19
Functions/EQ_structures/duobinary_target.m
View File

@@ -9,6 +9,7 @@ arguments
tx_bits
options.precode_mode db_mode
options.showAnalysis = 0;
options.eth_style_symbol_mapping = 0;
end
%Duobinary Targeting
@@ -18,7 +19,7 @@ end
% dir = [1,1];
mlse_sig_sd = mlse_.process(eq_signal);
mlse_sig_hd = PAMmapper(M,0).quantize(mlse_sig_sd);
mlse_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd);
% precoding to mitigate error propagation, most prominently used in
% combination with duobinary signaling to avoid catastrophic error
@@ -37,12 +38,12 @@ end
tx_symbols_precoded = Duobinary().encode(tx_symbols);
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
tx_bits = PAMmapper(M,0).demap(tx_symbols_precoded);
tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded);
case db_mode.db_discard
% normal dsp for precoded sequence == discard/omit/ignore precode
tx_bits = PAMmapper(M,0).demap(tx_symbols);
tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols);
case db_mode.db_encoded
@@ -56,7 +57,7 @@ end
end
% M = numel(unique(tx_symbols.signal));
rx_bits = PAMmapper(M,0).demap(mlse_sig_hd);
rx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd);
[~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
@@ -66,11 +67,15 @@ end
if options.showAnalysis
eq_noise = eq_noise - mean(eq_noise.signal);
rx_signal.spectrum("normalizeTo0dB",1,"fignum",250);
rx_signal.spectrum("normalizeTo0dB",1,"fignum",250,"displayname","Rx Spectrum");
showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise');
Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",250,"displayname","DB encoded reference");
Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",250);
showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise after Equalization');
fprintf('DB tgt BER: %.2e \n',ber);
end

7
Functions/EQ_structures/vnle.m
View File

@@ -74,20 +74,21 @@ function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options)
[~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
[evm_total,evm_lvl] = calc_evm(eq_signal_sd,tx_symbols);
[inf_rate] = calc_air(eq_signal_sd,tx_symbols,"skip_front",10000,"skip_end",10000);
snr(eq_signal_sd.signal,eq_noise.signal)
eq_package.ber_vnle = ber;
eq_package.evm_total = evm_total;
eq_package.evm_lvl = evm_lvl;
eq_package.inf_rate_vnle = inf_rate;
eq_package.snr = snr(eq_signal_sd.signal,eq_noise.signal);
eq_package.signal = eq_signal_sd;
if options.showAnalysis
snr(eq_signal_sd.signal,eq_noise.signal);
fprintf('SNR: %d dB \n',snr(eq_signal_sd.signal,eq_noise.signal));
fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl);
fprintf('VNLE BER: %.2e \n',ber);

11
Functions/EQ_structures/vnle_postfilter_mlse.m
View File

@@ -10,6 +10,7 @@ arguments
tx_bits
options.precode_mode db_mode
options.showAnalysis = 0;
options.eth_style_symbol_mapping = 0;
end
%FFE or VNLE
@@ -23,7 +24,7 @@ end
% [mlse_sig_hd,mlse_sig_sd] = mlse_.process(mlse_sig_sd,tx_symbols);
mlse_sig_sd = mlse_.process(mlse_sig_sd);
mlse_sig_hd = PAMmapper(M,0).quantize(mlse_sig_sd);
mlse_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd);
% precoding to mitigate error propagation, most prominently used in
% combination with duobinary signaling to avoid catastrophic error
@@ -46,12 +47,12 @@ end
tx_symbols_precoded = Duobinary().encode(tx_symbols);
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
tx_bits = PAMmapper(M,0).demap(tx_symbols_precoded);
tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded);
case db_mode.db_discard
% normal dsp for precoded sequence == discard/omit/ignore precode
tx_bits = PAMmapper(M,0).demap(tx_symbols);
tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols);
case db_mode.db_encoded
@@ -68,7 +69,7 @@ end
end
% METRICS OF VNLE %
rx_bits_vnle = PAMmapper(M,0).demap(eq_signal_hd);
rx_bits_vnle = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(eq_signal_hd);
[~,~,ber_vnle,~] = calc_ber(rx_bits_vnle.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% correct TUM implementation of AIR
@@ -76,7 +77,7 @@ end
[evm_vnle_total,evm_vnle_lvl] = calc_evm(eq_signal_sd,tx_symbols);
% METRICS OF MLSE (HD-VITERBI)
rx_bits_mlse = PAMmapper(M,0).demap(mlse_sig_hd);
rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd);
[~,~,ber_mlse,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);

6
projects/IMDD_base_system/imdd_it.m
View File

@@ -7,13 +7,19 @@ if 0
uloops.db_precode = [0];
uloops.bitrate = [300,330,360,390,420,450,480].*1e9; %[300,330,360,390,420,450,480]
% uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4];
uloops.precomp = [0];
uloops.db_precode = [0];
uloops.bitrate = [330].*1e9; %[300,330,360,390,420,450,480]
uloops.laser_wavelength = [1310];
uloops.M = [4,6,8];
uloops.M = [4];
uloops.link_length = [2]; % 1,2,3,5,6,8,10
wh = DataStorage(uloops);
wh.addStorage("ber");
wh = submit_simulations(wh,"parallel",0,"simulation_mode",0);
wh = submit_simulations(wh,"parallel",0,"simulation_mode",1);
end
wh_ana = wh_master;

24
projects/IMDD_base_system/imdd_model.m
View File

@@ -35,7 +35,7 @@ tx_bw_nyquist = 0.8;
link_length = 1;
% RX
rop = -8;
rop = -5;
rx_bw_nyquist = 0.8;
vnle_order1 = 50;
@@ -63,7 +63,7 @@ mu_dfe = 0.0004;
dfe_ = sum(dfe_order)>0;
doub_mode = db_mode.db_precoded;
doub_mode = db_mode.no_db;
%%% change specific parameter if given in varargin
% Parse optional input arguments
@@ -135,9 +135,14 @@ f_nyquist = fsym/2;
if simulation_mode
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
rcalpha = 1;
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha);
db_precode = 0;
db_encode = 0;
apply_pulsef = 1;
[Digi_sig,Symbols,Tx_bits] = PAMsource(...
"fsym",fsym,"M",M,"order",16,"useprbs",0,...
"fsym",fsym,"M",M,"order",18,"useprbs",0,...
"fs_out",fdac,...
"applyclipping",0,"clipfactor",1.5,...
"applypulseform",apply_pulsef,"pulseformer",Pform,...
@@ -145,7 +150,9 @@ if simulation_mode
"db_precode",db_precode,"db_encode",db_encode,...
"mrds_code",0,"mrds_blocklength",512).process();
% Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1);
Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1);
%%%%% AWG
% El_sig = M8199A("kover",kover).process(Digi_sig);
@@ -276,6 +283,13 @@ for occ = 1:proc_occ
Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig);
Scpe_sig = Scpe_sig - mean(Scpe_sig.signal);
%
% Pform = Pulseformer("fsym",Scpe_sig.fs,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha,"matched",0);
%
% Scpe_sig_matched = Pform.process(Scpe_sig);
%
% Scpe_sig.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","scope ");
% Scpe_sig_matched.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","matched");
%%% EQUALIZING
@@ -296,7 +310,7 @@ for occ = 1:proc_occ
end
%%%%% VNLE + PF + MLSE %%%%
if 1
if 0
% len_tr = length(Symbols)-1000;
eq_vnle_ = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);

BIN
projects/Messung_Zürich/Pmod_-0dBm_P_PD_5p8dBm_W03C34-0409E03_128GBd_2PAM__20250217T161301.h5
View File

Binary file not shown.

BIN
projects/Messung_Zürich/Pmod_-0dBm_P_PD_5p8dBm_W03C34-0409E03_32GBd_2PAM__20250217T161634.h5
View File

Binary file not shown.

BIN
projects/Messung_Zürich/Pmod_10p006dBm_P_PD_6p4732dBm_W03C34-0409E03_128GBd_4PAM__20250218T195459.h5
View File

Binary file not shown.

BIN
projects/Messung_Zürich/Pmod_10p01dBm_P_PD_6p4727dBm_W03C34-0409E03_160GBd_2PAM__20250218T195221.h5
View File

Binary file not shown.

BIN
projects/Messung_Zürich/Pmod_7p9772dBm_P_PD_6p477dBm_W03C34-0409E03_200GBd_2PAM__20250218T181515.h5
View File

Binary file not shown.

BIN
projects/Messung_Zürich/Pmod_7p9943dBm_P_PD_6p4526dBm_W03C34-0409E03_160GBd_2PAM__20250218T181637.h5
View File

Binary file not shown.

BIN
projects/Messung_Zürich/Pmod_8p0011dBm_P_PD_6p4762dBm_W03C34-0409E03_64GBd_4PAM__20250218T181158.h5
View File

Binary file not shown.

BIN
projects/Messung_Zürich/Pmod_8p0023dBm_P_PD_6p4776dBm_W03C34-0409E03_64GBd_2PAM__20250218T181323.h5
View File

Binary file not shown.

BIN
projects/Messung_Zürich/Pmod_9p994dBm_P_PD_7p4733dBm_W03C34-0409E03_128GBd_4PAM__20250218T195746.h5
View File

Binary file not shown.

113
projects/Messung_Zürich/dsp_all_files.m Normal file
View File

@@ -0,0 +1,113 @@
M = 8;
if M == 4
file_codes = {"20250221T032043",...
"20250221T032312",...
"20250221T032424",...
"20250221T032529",...
"20250221T032632",...
"20250221T032800",...
"20250221T032936",...
"20250221T033035",...
"20250221T033138",...
"20250221T033138",...
"20250221T033425",...
"20250221T033527",...
"20250221T033642",...
"20250221T033743",...
"20250221T033851",...
"20250221T034314",...
"20250221T034529",...
"20250221T034647",...
"20250221T034756",...
"20250221T034915",...
};
baudrate = [96,120,128,136,144,152,160,168,172,176,180,184,188,192,196,200,204,212,220,224];
pam4.result_vnle = result_vnle;
pam4.result_db = result_db;
elseif M ==8
file_codes = {"20250221T004651",... %96
"20250221T004926",...
"20250221T005844",...
"20250221T010032",...
"20250221T010150",...
"20250221T025505",...%144 ohne balun
"20250221T024629",...%152 ohne balun
"20250221T024315",...%160 ohne balun
"20250221T023534",...%168 ohne balun
"20250221T024929",...%172 ohne balun
"20250221T025305",...%176 ohne balun
};
file_codes = {"20250221T004651"};
baudrate = [96,112,120,128,136,144,152,160,168,172,176];
end
if 1
uloops = struct;
uloops.filecode = file_codes;
uloops.tcorrect = [0];
wh = DataStorage(uloops);
wh.addStorage("ber");
wh = submit_handle(@dsp_ief_file,wh,"parallel",0);
end
a = wh.getStoValue('ber',uloops.filecode, uloops.tcorrect);
% get best results per baudrate
% ber_vnle_values = cellfun(@(a) cellfun(@(y) y.ber_vnle, a.vnle_package, 'UniformOutput', false), a, 'UniformOutput', false);
% best_vnle = cellfun(@(x) min(cell2mat(x)), ber_vnle_values);
ber_vnle_values = cellfun(@(a) cellfun(@(y) y.ber_vnle, a.vnle_pf_package, 'UniformOutput', false), a, 'UniformOutput', false);
best_vnle = cellfun(@(x) min(cell2mat(x)), ber_vnle_values);
ber_mlse_values = cellfun(@(a) cellfun(@(y) y.ber_mlse, a.vnle_pf_package, 'UniformOutput', false), a, 'UniformOutput', false);
best_mlse = cellfun(@(x) min(cell2mat(x)), ber_mlse_values);
%%% PLOT PARAMETER OPTIMIZATION
figure()
cols = cbrewer2('Set1',6);
hold on
title(sprintf('%d km | %d nm | PAM %d',1,1313.5,M));
xax = uloops.tcorrect;
plot(xax,best_vnle,'DisplayName',sprintf('VNLE'),'Color',cols(1,:),'LineStyle','-','HandleVisibility','on');
plot(xax,best_mlse,'DisplayName',sprintf('VNLE + 1 tap post-filter + MLSE'),'Color',cols(2,:),'LineStyle','-','HandleVisibility','on');
% plot(xax,best_db,'DisplayName',sprintf('DB tgt. + MLSE'),'Color',cols(3,:),'LineStyle','-','HandleVisibility','on');
set(gca, 'YScale', 'log');
ylim([5e-5 0.3]);
xlim([min(xax), max(xax) ]);
yline([3.8e-3, 2e-2],'HandleVisibility','off');
legend
beautifyBERplot()
xlabel('Number 1st order taps');
ylabel('BER');
%%% PLOT COMPLETE BAUDRATE
figure(20)
cols = cbrewer2('Set1',6);
hold on
title(sprintf('%d km | %d nm | PAM %d',1,1313.5,M));
plot(baudrate,best_vnle,'DisplayName',sprintf('VNLE'),'Color',cols(1,:),'LineStyle','-','HandleVisibility','on');
plot(baudrate,best_mlse,'DisplayName',sprintf('VNLE + 1 tap post-filter + MLSE'),'Color',cols(2,:),'LineStyle','-','HandleVisibility','on');
plot(baudrate,best_db,'DisplayName',sprintf('DB tgt. + MLSE'),'Color',cols(3,:),'LineStyle','-','HandleVisibility','on');
set(gca, 'YScale', 'log');
ylim([5e-5 0.3]);
xlim([min(baudrate), max(baudrate) ]);
yline([3.8e-3, 2e-2],'HandleVisibility','off');
legend
beautifyBERplot()
xlabel('Bit Rate in Gbps');
ylabel('BER');

300
projects/Messung_Zürich/dsp_ief_file.m
View File

@@ -1,91 +1,215 @@
% Get the script folder
scriptFolder = fileparts(mfilename('fullpath'));
function [results] = dsp_ief_file(varargin)
mu_ffe1 = 0.0001;
mu_ffe2 = 0.0008;
mu_ffe3 = 0.001;
mu_dc = 0.005;
mu_dfe = 0.0004;
vnle_order1 = 50;
vnle_order2 = 3;
vnle_order3 = 3;
tcorrect = 0;
%%% change specific parameter if given in varargin
% Parse optional input arguments
if ~isempty(varargin)
var_s = varargin{1};
if isstruct(var_s)
fields = fieldnames(var_s);
for i = 1:numel(fields)
if isnumeric(fields{i})
eval([fields{i}, ' = ', num2str( var_s.(fields{i}) ), ';']);
fprintf("%s <-- %.2f \n", fields{i}, var_s.(fields{i}));
else
eval([fields{i}, ' = ', 'var_s.(fields{',num2str(i),'})' , ';']);
end
end
else
error('Optional variables should be passed as a struct.');
end
end
ssdPath = "D:\36_IMDD_Kiel\Data";
% code = "20250221T004802";
% code = "20250221T032043";
filename = findFileByCode(ssdPath, filecode{1});
% Extract parameters from the filename using an updated regex
tokens = regexp(filename, 'Pmod_([-0-9p]+)dBm_P_PD_([-0-9p]+)dBm_.*?_(\d+)GBd_(\d+)PAM__\d+T\d+', 'tokens');
if isempty(tokens)
error('Filename format not recognized.');
end
tokens = tokens{1};
% Convert values
config.P_laser = str2double(strrep(tokens{1}, 'p', '.'));
config.P_pd = str2double(strrep(tokens{2}, 'p', '.'));
config.fsym = str2double(tokens{3}) * 1e9; % Convert GBd to Hz
config.M = str2double(tokens{4});
% Display results
fprintf('Loaded file: %s\n', filename);
fprintf('P_laser: %.3f dBm\n', config.P_laser);
fprintf('P_pd: %.3f dBm\n', config.P_pd);
fprintf('fsym: %.1f GBd\n', config.fsym*1e-9);
fprintf('M: %d\n', config.M);
%%% Load Data
filepath = fullfile(filename);
ief_ = h5info(filepath);
config.fs_rx = h5readatt(filepath,'/','fs'); %sampling frequency at Rx
config.fs_tx = h5readatt(filepath, '/','fs_Tx'); %sampling frequency at Tx
config.fsym = h5readatt(filepath, '/','R'); %Baudrate
config.M = h5readatt(filepath, '/','M'); % PAM- 'M'
config.ROF = h5readatt(filepath, '/','ROF');
config.PulseShape =h5readatt(filepath, '/','PulseShape');
yOrg = h5readatt(filepath, ief_.Groups(4).Groups(1).Name, 'YOrg');
yInc = h5readatt(filepath, ief_.Groups(4).Groups(1).Name, 'YInc');
dataRx = double(h5read(filepath, '/Waveforms/Channel 2/Channel 2Data')); % rohdaten des CH4
bitsTx = h5read(filepath, '/Settings/dataTx'); %Binär
bitsTx = reshape(bitsTx,log2(config.M),[])';
%%% Build Tx Signal (Bits, Pam Map, Symbols)
Tx_bits = Informationsignal(bitsTx);
Tx_symbols = PAMmapper(config.M,0,"eth_style",1).map(Tx_bits);
Tx_symbols.fs = config.fsym;
%%% Build Rx Signal (Rx, normalize,remove mean)
dataRx = dataRx*yInc+yOrg;
Rx_Sig = Informationsignal(dataRx,"fs",config.fs_rx);
Rx_Sig.signal = Rx_Sig.signal - mean(Rx_Sig.signal);
Rx_Sig = Rx_Sig.normalize("mode","rms");
% Rx_Sig.spectrum("fignum",3,"displayname",'Rx Signal','normalizeTo0dB',0);
Rx_Sig = Filter('filtdegree',0,"f_cutoff",Tx_symbols.fs.*0.9,"fs",Rx_Sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Rx_Sig);
% Rx_Sig.spectrum("fignum",3,"displayname",'Rx Signal filt','normalizeTo0dB',0);
Rx_Sig = Rx_Sig.delay(tcorrect,"mode","samples");
%%%%%% Sample to 2x fsym %%%%%%
Rx_Sig_resamp = Rx_Sig.resample("fs_out",8*config.fsym);
%%%%%% Sync Rx signal with reference (S is a cell array with all occurences) %%%%%%
[Rx_Sig_sync,S,isFlipped] = Rx_Sig_resamp.tsynch("reference",Tx_symbols,"fs_ref",config.fsym,"debug_plots",0);
% Rx_Sig_sync.eye(fsym,M,"fignum",4,"displayname",'eye diagram');
% timing syncroization ??
% SYMSYNC = comm.SymbolSynchronizer("Modulation","PAM/PSK/QAM","SamplesPerSymbol",2);
%
% synd = SYMSYNC(Rx_Sig_sync.signal);
% Rx_Sig_sync.signal = synd;
%
% [Rx_Sig_sync,S,isFlipped] = Rx_Sig_sync.tsynch("reference",Tx_symbols,"fs_ref",config.fsym,"debug_plots",1);
%
% filter Bw
output = struct();
vnle_package = {};
vnle_pf_package = {};
dbtgt_package = {};
for s = 1:length(S)
Rx_Sig_sync = S{s};
% figure();hold on;plot(Tx_symbols.resample("fs_out",2*config.fsym).normalize("mode","oneone").signal);plot(Rx_Sig_sync.normalize("mode","oneone").signal);
%
Tx_symbols.spectrum("fignum",3,"displayname",'No Matched Filter','normalizeTo0dB',1);
Rx_Sig_sync.spectrum("fignum",3,"displayname",'No Matched Filter','normalizeTo0dB',1);
Nsym = 64;
sampsPerSym = 8;
rcrFilt = comm.RaisedCosineReceiveFilter(...
'Shape', 'Square root', ...
'RolloffFactor', 1, ...
'FilterSpanInSymbols', Nsym, ...
'InputSamplesPerSymbol', sampsPerSym, ...
'DecimationFactor', 1);
yr = rcrFilt([Rx_Sig_sync.signal;zeros(Nsym*sampsPerSym/2, 1)]);
fltDelay = Nsym / (2*config.fsym);
yr = yr(fltDelay*Rx_Sig_sync.fs+1:end);
Rx_matched = Rx_Sig_sync;
Rx_matched.signal = yr;
Rx_matched = Rx_matched.resample("fs_out",2*config.fsym);
symbolSync = comm.SymbolSynchronizer("Modulation","PAM/PSK/QAM","SamplesPerSymbol",2);
Rx_syncd = Rx_matched;
Rx_syncd.signal = real(symbolSync(Rx_syncd.signal));
Rx_syncd.fs = config.fsym;
% figure();hold on;
% stem(Tx_symbols.normalize("mode","oneone").signal);
% stem(Rx_syncd.normalize("mode","oneone").signal);
%
% figure();hold on;
% plot(Rx_Sig_sync.normalize("mode","oneone").signal);
% plot(Rx_matched.normalize("mode","oneone").signal);
Rx_syncd = Rx_syncd.resample("fs_out",2*config.fsym);
Rx_matched.spectrum("fignum",3,"displayname",'No Matched Filter','normalizeTo0dB',1);
mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3];
vnle_order=[vnle_order1,vnle_order2,vnle_order3];
eq_ = EQ("Ne",vnle_order,"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",1,"ideal_dfe",1);
%%%%% VNLE only (or DFE) %%%%
if 1
[result] = vnle(eq_,config.M,Rx_matched,Tx_symbols,Tx_bits,"precode_mode",db_mode.no_db,"showAnalysis",1,'eth_style',1);
vnle_package{s} = result;
end
%%%%% VNLE + PF + MLSE %%%%
if 0
pf_ = Postfilter("ncoeff",1,"useBurg",1);
mlse_ = MLSE_viterbi("duobinary_output",0,'M',config.M,'trellis_states',PAMmapper(config.M,0).levels);
doub_mode = db_mode.no_db;
[result] = vnle_postfilter_mlse(eq_,pf_,mlse_,config.M,Rx_Sig_sync,Tx_symbols,Tx_bits,"precode_mode",doub_mode,'showAnalysis',1,'eth_style_symbol_mapping',1);
vnle_pf_package{s} = result;
end
%%%%% Duobinary Targeting %%%%
if 0
mlse_db = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",config.M,"trellis_states",PAMmapper(config.M,0).levels);
doub_mode = db_mode.db_emulate;
[result] = duobinary_target(eq_, mlse_db, config.M, Rx_Sig_sync, Tx_symbols, Tx_bits, "precode_mode", doub_mode,'showAnalysis',1,'eth_style_symbol_mapping',1);
dbtgt_package{s} = result;
end
end
results.vnle_package = vnle_package;
results.vnle_pf_package = vnle_pf_package;
results.dbtgt_package = dbtgt_package;
% Find the HDF5 file in the folder
files = dir(fullfile(scriptFolder, '*.h5'));
if isempty(files)
error('No HDF5 files found in the script folder.');
end
% Select the first matching file (modify as needed)
filename = files(6).name;
% filename = 'Pmod_8p0023dBm_P_PD_6p4776dBm_W03C34-0409E03_64GBd_2PAM__20250218T181323.h5';
% Extract parameters from the filename using an updated regex
tokens = regexp(filename, 'Pmod_([-0-9p]+)dBm_P_PD_([-0-9p]+)dBm_.*?_(\d+)GBd_(\d+)PAM__\d+T\d+', 'tokens');
if isempty(tokens)
error('Filename format not recognized.');
end
tokens = tokens{1};
% Convert values
P_laser = str2double(strrep(tokens{1}, 'p', '.'));
P_pd = str2double(strrep(tokens{2}, 'p', '.'));
fsym = str2double(tokens{3}) * 1e9; % Convert GBd to Hz
M = str2double(tokens{4});
% Display results
fprintf('Loaded file: %s\n', filename);
fprintf('P_laser: %.3f dBm\n', P_laser);
fprintf('P_pd: %.3f dBm\n', P_pd);
fprintf('fsym: %.3f Hz\n', fsym);
fprintf('M: %d\n', M);
%%% Load Data
folder = "/Users/silasoettinghaus/Documents/MATLAB/imdd_simulation/projects/Messung_Zürich";
filepath = fullfile(folder,filename);
ief_ = h5info(filepath);
fs_rx = h5readatt(filepath,'/','fs'); %sampling frequency at Rx
fs_tx = h5readatt(filepath, '/','fs_Tx'); %sampling frequency at Tx
fsym = h5readatt(filepath, '/','R'); %Baudrate
M = h5readatt(filepath, '/','M'); % PAM- 'M'
ROF = h5readatt(filepath, '/','ROF');
PulseShape =h5readatt(filepath, '/','PulseShape');
yOrg = h5readatt(filepath, ief_.Groups(4).Groups(1).Name, 'YOrg');
yInc = h5readatt(filepath, ief_.Groups(4).Groups(1).Name, 'YInc');
plotDisplayName = sprintf('P_{cw} = %.3f dBm, P_{PD} = %.3f dBm, f_{sym} = %.3f GBd, M = %d', ...
P_laser, P_pd, fsym / 1e9, M);
dataRx = double(h5read(filepath, '/Waveforms/Channel 2/Channel 2Data')); % rohdaten des CH4
bitsTx = h5read(filepath, '/Settings/dataTx'); %Binär
bitsTx = reshape(bitsTx,log2(M),[])';
%%% Build Tx Signal (Bits, Pam Map, Symbols)
Tx_bits = Informationsignal(bitsTx);
Tx_symbols = PAMmapper(M,0,"eth_style",1).map(Tx_bits);
Tx_symbols.fs = fsym;
% Tx_symbols.plot("fignum",1,"displayname",'Rx Signal');
% Tx_symbols.move_it_spectrum("fignum",2,"displayname",'tx spectrum');
% Tx_symbols.spectrum("fignum",3,"displayname",'Rx Signal','normalizeTo0dB',0);
%%% Build Rx Signal (Rx, normalize,remove mean)
dataRx = dataRx*yInc+yOrg;
Rx_Sig = Informationsignal(dataRx,"fs",fs_rx);
Rx_Sig = Rx_Sig.normalize("mode","rms");
% Rx_Sig.power();
% Rx_Sig.eye(fsym,M,"fignum",4,"displayname",'eye diagram');
% Rx_Sig.plot("fignum",5,"displayname",'Rx Signal');
% Rx_Sig.move_it_spectrum("fignum",6,"displayname",'tx spectrum');
Rx_Sig.spectrum("fignum",3,"displayname",'Rx Signal','normalizeTo0dB',0);
%%%%%% Sample to 2x fsym %%%%%%
Rx_Sig_resamp = Rx_Sig.resample("fs_out",2*fsym);
%%%%%% Sync Rx signal with reference (S is a cell array with all occurences) %%%%%%
[Rx_Sig_sync,S,isFlipped] = Rx_Sig_resamp.tsynch("reference",Tx_symbols,"fs_ref",fsym,"debug_plots",0);
eq_vnle_dfe = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.05,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",1,"ideal_dfe",0);
[result] = vnle(eq_vnle_dfe,M,Rx_Sig_sync,Tx_symbols,Tx_bits,"precode_mode",db_mode.no_db,"showAnalysis",1,'eth_style',1);

40
projects/Messung_Zürich/findFileByCode.m Normal file
View File

@@ -0,0 +1,40 @@
function fileName = findFileByCode(folderPath, code)
% FINDFILEBYCODE Searches for a file in a folder structure by a given code.
% fileName = findFileByCode(folderPath, code) searches recursively in
% folderPath for a file containing 'code' in its name and returns the
% full file name if found.
%
% Inputs:
% folderPath - The root directory to search in
% code - The unique code to search for in file names
%
% Output:
% fileName - The full file name if found, empty if not found
% Initialize output
fileName = '';
% Get list of all files and folders in the folderPath
files = dir(folderPath);
% Iterate through the list
for i = 1:length(files)
% Skip '.' and '..'
if files(i).isdir
if ~startsWith(files(i).name, '.') % Avoid hidden folders
% Recursive search in subdirectories
subFolder = fullfile(folderPath, files(i).name);
fileName = findFileByCode(subFolder, code);
if ~isempty(fileName)
return; % Stop searching if found
end
end
else
% Check if the file name contains the code
if contains(files(i).name, code)
fileName = fullfile(folderPath, files(i).name);
return;
end
end
end
end

BIN
projects/Messung_Zürich/pam4.mat Normal file
View File

Binary file not shown.

BIN
projects/Messung_Zürich/pam4_eq.fig Normal file
View File

Binary file not shown.

69
projects/Messung_Zürich/submit_handle.m Normal file
View File

@@ -0,0 +1,69 @@
function wh = submit_handle(funcHandle, wh, options)
arguments
funcHandle
wh
options.parallel = 1;
end
%%% 2) SUBMIT SIMULATION
% Initialize job results
if options.parallel
curpool = gcp('nocreate');
if isempty(curpool)
parpool;
else
% stop all forgotten or unfetched jobs from queue
if ~isempty(curpool.FevalQueue.QueuedFutures) || ~isempty(curpool.FevalQueue.RunningFutures)
oldq = length(curpool.FevalQueue.QueuedFutures) + length(curpool.FevalQueue.RunningFutures);
curpool.FevalQueue.cancelAll
fprintf('Canceled %d unfetched jobs from old queue.', oldq);
end
end
results = parallel.FevalFuture.empty();
else
results = [];
end
fprintf('Requested %d loops\n', wh.getLastLinIndice);
for lin_idx = 1:wh.getLastLinIndice
optionalVars = struct();
if ~isempty(wh.getDimension)
% Build the optionalVars struct
[parametervalues, parameternames] = wh.getPhysIndicesByLinIndex(lin_idx);
for pidx = 1:numel(parameternames)
optionalVars.(parameternames{pidx}) = parametervalues{pidx};
end
end
%%% SIMULATION HERE
if options.parallel
numOutputs = 1;
results(lin_idx) = parfeval(funcHandle, numOutputs, optionalVars);
else
finalresults{lin_idx} = feval(funcHandle, optionalVars);
wh.addValueToStorageByLinIdx(finalresults{lin_idx}, 'ber', lin_idx);
end
end
if options.parallel
%%% 4) Setup waitbar
h = waitbar(0, 'Processing Simulations...');
updateWaitbar = @(~) waitbar(mean(arrayfun(@(f) strcmp(f.State, 'finished'), results)), h);
fprintf('Fetching results... \n');
updateWaitbarFutures = afterEach(results, updateWaitbar, 0);
afterAll(updateWaitbarFutures, @(~) delete(h), 0);
%%% 7) Fetch final results after all computations
fetchOutputs(results);
for ridx = 1:length(results)
wh.addValueToStorageByLinIdx(results(ridx).OutputArguments{1}, 'ber', ridx);
end
end
end