Files
magf 005e821131 Added:
- Class folder 'Timing Recovery' with different timing recoveries
- Minimal example for the timing recovery on the FSO data
- New evaluation scripts in the FSO project folder
2026-02-02 10:56:12 +01:00

293 lines
11 KiB
Matlab
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<<<<<<< HEAD
base = "C:\Users\magf\Desktop\Desktop\MATLAB-Zeugs\FSO Equalizer\FSO_FP_QCL_60umUTC\";
=======
base = "C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC";
>>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0
mode = 0; %0 oder 1
M = 2;
all_files = dir(fullfile(base, "**/*.mat"));
if M == 2
<<<<<<< HEAD
tx_data_path = fullfile(base, "14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat");
filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat");
elseif M == 4
tx_data_path = fullfile(base, "6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat");
=======
tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat");
filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat");
elseif M == 4
tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat");
>>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0
filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat");
end
if mode == 1
[f, p] = uigetfile(fullfile(base, "**/*.mat"));
if f~=0
filename = fullfile(p,f);
end
end
tx_data = load(tx_data_path);
datas = load(filename);
%%
str = filename;
M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once'));
assert(M==M_);
fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once'));
fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once'));
I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f');
rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f');
L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f');
pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once'));
rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once'));
mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once'));
%%
% Tx data
Bits = Informationsignal(tx_data.tx_data,"fs",fsym);
Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym);
mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there
PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme
Symbols_ = PM.map(Bits) .* PM.scaling;
assert(isequal(Symbols.signal,Symbols_.signal));
Bits_ = PM.demap(Symbols);
[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal);
assert(ber == 0);
%% For comparison, apply pulsef on Tx Symbols
Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff);
Digi_sig_tx_compare = Pform.process(Symbols);
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff);
Digi_sync = Pform.process(Symbols);
MF = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff);
Rx_sig_compare = MF.process(Digi_sig_tx_compare);
%%
% Rx Data
traceData = datas.tr.lastData(2).trace.ch3;
%FYI: Voltage=(RawDataYReference)×YIncrement+YOrigin
scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin;
demystified = isequal(traceData.YData,scoperead_volts);
assert(demystified);
timesig_compare = [0:1:datas.tr.lastData(1).trace.ch3.Points-1] ./ fs;
timesig = datas.tr.lastData(1).trace.ch3.XData;
assert(isequal(traceData.YData,scoperead_volts));
Scope_sig = Electricalsignal(traceData.YData,"fs",fs);
%%
% 1) matched filter
% pulse is symmetric, hence we can use pulsef directly as matched filter.
% It feels off (bit I think correct) that the fsym is now the output freq.!!
% -> output 2 sps to omit timing recovery!?
apply_matched_filter = 0;
k = 4;
if apply_matched_filter
Pform = Pulseformer("fsym",fsym,"fdac",k*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1);
Rx_matched = Pform.process(Scope_sig);
else
Rx_matched = Filter('filtdegree',4,"f_cutoff",fsym*0.5,"fs",Scope_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scope_sig);
Rx_matched = Rx_matched.resample("fs_out",k*fsym);
end
Rx_matched.spectrum();
<<<<<<< HEAD
%%
sys = comm.SymbolSynchronizer('TimingErrorDetector', 'Gardner (non-data-aided)', ...
'SamplesPerSymbol', 2, ...
'DampingFactor', 0.7, ...
'NormalizedLoopBandwidth', 0.01);
Rx_symbolsync = Rx_matched;
[Rx_symbolsync.signal, timing_error] = sys(Rx_matched.signal);
plot(timing_error); % If this is a ramp, you have drift!
%% timing sync -> at this point we still have no symbol timing recovery, we
% % try to do this with 2sps EQ!
[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 1);
% Rx_matched_1 = Rx_synced_cell{1};
%
% % Rx_Time_Rec = Rx_matched;
% [Rx_Time_Rec, Timing_Error] = Timing_Recovery("modulation", 'PAM/PSK/QAM', "timing_error_detector",'Gardner (non-data-aided)','sps',2,'damping_factor',1,'normalized_loop_bandwidth',0.01,'detector_gain',2.7).process(Rx_matched_1);
% figure;plot(Timing_Error);
% % Rx_Time_Rec.signal = resample(Rx_Time_Rec.signal, 28e9, 14e9);
%% not working..
% Rx_synced = Rx_Time_Rec;
Rx_synced = Rx_synced_cell{1};
len_tr = 4096*2;
mu_ffe1 = 0.0001;
mu_ffe2 = 0.0008;
mu_ffe3 = 0.001;
mu_dc = 0.005;
mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3];
mu_dfe = 0.0004;
duob_mode = db_mode.no_db;
Rx_synced.plot("displayname",'RX: Matched+Sync+2sps','fignum',2);
Digi_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped','fignum',1,'normalizeTo0dB',0);
Rx_sig_compare.normalize("mode","rms").spectrum("displayname",'Tx: RC-shaped + matched filtered ','fignum',1,'normalizeTo0dB',0);
Rx_synced.normalize("mode","rms").spectrum("displayname",'RX: matched filtered + synced','fignum',1,'normalizeTo0dB',0);
=======
%% Timing Rec
apply_timing_rec = 1;
if apply_timing_rec
[Rx_symbolsync, timing_error] = Timing_Recovery("timing_error_detector",'Gardner (non-data-aided)','sps',k,'damping_factor',0.1,'normalized_loop_bandwidth',0.1,'detector_gain',2.7).process(Rx_matched);
figure();plot(timing_error);
Rx_symbolsync.fs = fsym;
% Tsynch
[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0);
Rx_synced = Rx_synced_cell{1};
sps = 1;
else
% Tsynch
[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0);
Rx_synced = Rx_synced_cell{1};
Rx_synced = Rx_synced.resample("fs_out",2*fsym);
sps = 2;
end
>>>>>>> 719e5508e776c18e3dcf72a954c71bdedda179e0
if M == 2
ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature
elseif M == 4
ber_in_paper = 10^(-2.5);
end
%%
Rx_synced = Rx_synced_cell{1};
% -------------------- FFE --------------------
% requires some more digging what is going on :-)
eq_ffe = EQ("Ne",[50, 1, 1],"Nb",[2,0,0], ...
"training_length",512,"training_loops",5,"dd_loops",5, ...
"K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
vars = logspace(-4,-3,36);
parfor i = 1:numel(vars)
len_tr = 4096;
mu_ffe1 = 0.01;% mus(i);%0.0001;
mu_ffe2 = 0.0008;
mu_ffe3 = 0.001;
mu_dc = 0.005;
mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3];
mu_dfe = vars(i);
duob_mode = db_mode.no_db;
% requires some more digging what is going on :-)
eq_ffe_1 = EQ("Ne",[150, 1, 0],"Nb",[50,0,0], ...
"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
"K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
"FFEmu",0,"plotfinal",0,"ideal_dfe",0);
eq_ffe_2 = FFE("epochs_tr",1,"epochs_dd",vars(i),"len_tr",4096,"mu_dd",vars(i),"mu_tr",vars(i),"order",999,"sps",1,"decide",0, "adaption",adaption_method.nlms,"dd_mode",0);
% eq_ffe_2 = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-5,"dfe_mu_dd",mus(i),"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",50,"dfe_order",10,"sps",1,"decide",1);
ffe_results = ffe(eq_ffe_1,M,Rx_synced,Symbols,Bits, ...
"precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ...
"eth_style_symbol_mapping",mapping_style);
ffe_results.metrics.BER
bers(i) = ffe_results.metrics.BER;
end
figure();
plot(vars,bers);
yline(ber_in_paper)
beautifyBERplot();
%
fprintf('Paper: %.1e \n \n',ber_in_paper);
ffe_results.metrics.print("description",'FFE');
fprintf('FFE: %.1e \n',ffe_results.metrics.BER);
%% -------------------- VNLE + MLSE --------------------
len_tr = 4096;
mu_ffe1 = 0.0001;% mus(i);%0.0001;
mu_ffe2 = 0.0008;
mu_ffe3 = 0.001;
mu_dc = 0.005;
mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3];
mu_dfe = 0.0004;
duob_mode = db_mode.no_db;
pf_ncoeffs = 4;
vars = 1:7;
bers = zeros(size(vars));
parfor i = 1:numel(vars)
eqv = EQ("Ne",[200, 1, 0],"Nb",[2, 0, 0], ...
"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
"K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
pf_ncoeffs = vars(i);
pf = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels);
[vnle_results, mlse_results] = vnle_postfilter_mlse(eqv, pf, mlse_, M, Rx_synced, Symbols, Bits, ...
"precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style);
fprintf('Paper: %.1e \n \n',ber_in_paper);
vnle_results.metrics.print("description",'VNLE');
mlse_results.metrics.print("description",'MLSE');
bers(i) = mlse_results.metrics.BER;
end
%%
figure();hold on
plot(vars,bers_ffe,'DisplayName','FFE [200,0,0] + PF + MLSE');
plot(vars,bers_vnle,'DisplayName','VNLE [200,1,0] + PF + MLSE');
plot(vars,bers_vnledfe,'DisplayName','VNLE [200,1,0] + DFE [2] + PF + MLSE');
plot(vars,bers_vnledfe_ideal,'DisplayName','VNLE [200,1,0] + ideal DFE [2] + PF + MLSE');
yline(ber_in_paper);
yline([2e-2, 4.85e-3, 3.8e-3, 2,2e-4],'LineWidth',2,'Color',[0.8,0.8,0.8],'LineStyle',':','HandleVisibility','off');
ylim([1e-5,0.1]);
beautifyBERplot();
%% -------------------- DB target --------------------
mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels);
eq_ = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
"K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ...
"precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style);
fprintf('Paper: %.1e \n \n',ber_in_paper);
dbt_results.metrics.print("description",'Duobinary');
%%
%ML-based MLSE (L=2)
mu_ml = 0.01; training_epochs = 100;
ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ...
"len_tr",len_tr,"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",1, ...
"traceback_depth",128,"L",3,"delta",4,"adaptive_mu",0);
[ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode);
ml_mlse_results.metrics.print("description",'ML ');