Files
imdd_silas/projects/FSO_transmission/first_analysis.m
Silas Oettinghaus 1bf17ec6be FSO:
- try to analyze spectrum of FSO and use the coeffs in the ML analysis. simple channel coeff + awgn example
2026-01-20 13:48:15 +01:00

252 lines
9.4 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.
base = "C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC";
mode = 0; %0 oder 1
M = 2;
all_files = dir(fullfile(base, "**/*.mat"));
if M == 2
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");
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
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 = 1;
k = 1;
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();
%%
coefficients = arburg(Rx_matched.signal,25);
figure()
[h,w] = freqz(1,coefficients,Rx_matched.length,"whole",Rx_matched.fs);
h = h/max(abs(h));
hold on
w_ = (w - Rx_matched.fs/2);
plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(coefficients,2)), ' '],'LineWidth',2);
%% 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
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 ');