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

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

Browse Source
This commit is contained in:
Silas Labor Zizou
2025-03-03 09:13:51 +01:00
parent 6570746d17 d099efea03
commit 4fedd806b5
295 changed files with 6849 additions and 2695 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
Classes/00_signals/Frame 1.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 101 KiB

219
Classes/00_signals/Signal.m
View File

@@ -133,6 +133,7 @@ classdef Signal
end end
%%
function plot(obj, options) function plot(obj, options)
% signal to plot: obj.signal % signal to plot: obj.signal
% fsamp : obj.fs (e.g. 92e9 => 92 GHz) % fsamp : obj.fs (e.g. 92e9 => 92 GHz)
@@ -144,6 +145,7 @@ classdef Signal
options.displayname = []; options.displayname = [];
options.timeframe = 0; options.timeframe = 0;
options.clear = 0; options.clear = 0;
options.color = [];
end end
figure(options.fignum); % If figure does not exist, create new figure figure(options.fignum); % If figure does not exist, create new figure
@@ -169,8 +171,11 @@ classdef Signal
end end
hold on; hold on;
if isempty(options.color)
plot(t, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1); plot(t, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
else
plot(t, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color);
end
% 2 c) % 2 c)
% - xlabel if not already here: time in readable format (1 ms and not 1e-3 s) % - xlabel if not already here: time in readable format (1 ms and not 1e-3 s)
% - ylabel amplitude % - ylabel amplitude
@@ -254,7 +259,12 @@ classdef Signal
CallingModifier = evalin('caller','obj'); CallingModifier = evalin('caller','obj');
end end
if isa(CallingModifier,"Signal")
CallingModifierStruct = obj.objToStructFilteredRecursive(CallingModifier); CallingModifierStruct = obj.objToStructFilteredRecursive(CallingModifier);
ModifierCopy = {CallingModifierStruct};
end
ModifierCopy = {CallingModifier};
SignalType = [string(class(obj))]; SignalType = [string(class(obj))];
TimeStamp = [(datetime('now','TimeZone','local','Format','HH:mm:ss'))]; TimeStamp = [(datetime('now','TimeZone','local','Format','HH:mm:ss'))];
@@ -263,7 +273,7 @@ classdef Signal
Nase = [0]; Nase = [0];
SignalCopy = obj.signal; SignalCopy = obj.signal;
ModifierName = class(CallingModifier); ModifierName = class(CallingModifier);
ModifierCopy = {CallingModifierStruct};
cell = {SignalType , TimeStamp , Length , SignalPower(1) , Nase, SignalCopy, ModifierName, ModifierCopy, Description}; cell = {SignalType , TimeStamp , Length , SignalPower(1) , Nase, SignalCopy, ModifierName, ModifierCopy, Description};
@@ -340,59 +350,193 @@ classdef Signal
options.color = []; options.color = [];
options.normalizeToNyquist = 0; options.normalizeToNyquist = 0;
options.normalizeTo0dB = 0; options.normalizeTo0dB = 0;
options.max_num_lines = []; % Leave empty or omit to disable line rotation
options.fft_length = [];
end end
% spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname); if isempty(options.fft_length)
options.fft_length = 2^(nextpow2(length(obj.signal))-7);
N = 2^(nextpow2(length(obj.signal))-10); end
if options.normalizeToNyquist == 0 if options.normalizeToNyquist == 0
[p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,obj.fs,"centered","psd","mean");
w = w.*1e-9; w = w.*1e-9;
else else
[p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,"centered","power","mean"); [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,"centered","psd","mean");
p_lin = smooth(p_lin,0.05,'rloess');
end end
if options.normalizeTo0dB if options.normalizeTo0dB
p_lin = p_lin./ max(p_lin); p_lin = p_lin./ max(p_lin);
p_dbm = 10*log10(p_lin); %dB to dBm in case of "power" p_dbm = 10*log10(p_lin); % normalized to 0 dB
ylab = "normalized to 0 dB"; ylab = "normalized to 0 dB";
else else
p_dbm = 10*log10(p_lin)+30; %dB to dBm in case of "power" p_dbm = 10*log10(p_lin);
ylab = "Power (dBm)"; ylab = "Power (dB/Hz)";
end end
figure(options.fignum); % If figure does not exist, create new figure figure(options.fignum);
ax = gca; ax = gca;
hold on hold on
if isempty(options.color) if isempty(options.color)
plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1); hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1);
else else
plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color); hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color);
end
% If user wants to limit the number of lines, check and remove old lines
if ~isempty(options.max_num_lines) && options.max_num_lines > 0
allLines = findall(ax, 'Type', 'Line');
if length(allLines) > options.max_num_lines
% Sort lines by creation order. Usually, the oldest lines appear first in allLines.
% If needed, you can sort by UserData or other criteria.
numToRemove = length(allLines) - options.max_num_lines;
delete(allLines(1:numToRemove));
end
end end
if options.normalizeToNyquist == 0 if options.normalizeToNyquist == 0
xlabel("Frequency in GHz"); xlabel("Frequency in GHz");
%xlim([-obj.fs/2 obj.fs/2].*1e-9)
edgetick = 2^(nextpow2(obj.fs*1e-9)); edgetick = 2^(nextpow2(obj.fs*1e-9));
xticks([-edgetick:16:edgetick]); xticks(-edgetick:16:edgetick);
xlim([100*round(min(w)/100,1)-10, 100*round(max(w)/100,1)+10]) xlim([100*round(min(w)/100,1)-10, 100*round(max(w)/100,1)+10])
xlim([-128 128]);%256GSa/s
else else
xlabel("Normalized Frequency"); xlabel("Normalized Frequency");
xlim([-pi, pi]); xlim([-pi, pi]);
end end
ylabel("Power/frequency (dB/Hz)");
ylabel(ylab); ylabel(ylab);
ylim([min(floor( min(p_dbm))-3 , ax.YLim(1)), max(ceil( max(p_dbm) )+(3), ax.YLim(2))]); try
yticks([-200:10:10]); ylim([max(min(floor(min(p_dbm))-3, ax.YLim(1)),-40), min(max(ceil(max(p_dbm))+3, ax.YLim(2)),10)]);
grid on catch
grid minor ylim([min(floor(min(p_dbm))-3, ax.YLim(1)), max(ceil(max(p_dbm))+3, ax.YLim(2))]);
legend('Interpreter','none'); end
yticks(-200:10:10);
grid on; grid minor;
legend
% legend('Interpreter','none');
end
function move_it_spectrum(obj,options)
arguments
obj
options.fignum
options.displayname = "";
options.color = [];
options.normalizeToNyquist = 0;
options.normalizeTo0dB = 0;
end
data_in = obj.signal;
if size(data_in,1) > size(data_in,2)
data_in = data_in';
end
for pol = 1:size(data_in,1)
%compute FFT of input
Data_in = fft( data_in(pol,:) );
%psd = Data_in.*conj(Data_in);
psd = Data_in;
%Use only magnitude of FFT (which was complex)
psd = abs(psd);
%Shift the spectrum to yield
psd = fftshift(psd);
%divide by N
psd = psd/length(data_in(pol,:));
psd_plot = 20*log10(psd);
% psd_plot = psd_plot - max(psd_plot);
%smoothing
% psd_smoothed = smooth(psd,1000);
%
% psd_smoothed = 10*log10(psd_smoothed);
% psd_smoothed = psd_smoothed - max(psd_smoothed);
carrier_power_time_dbm = 20*log10( mean(abs(data_in)) .^2 )+30; % dB -> +30 -> dBm
carrier_power_freq_dbm = max(psd_plot);
% psd_plot = psd_plot - max(psd_plot);
%% cspr
c = mean(data_in).^2;
s = mean(data_in.^2);
cspr = 10*log10(c / s);
testParseval = 1;
if testParseval == 1
E_FreqDomain =1/length(psd) * sum((psd.*length(psd)).^2);
%test parseval
E_TimeDomain = sum( (data_in(pol,:).^2) );
if isequal(round(E_FreqDomain,1),round(E_TimeDomain,1))
disp('Parseval is right!');
else
% disp('Something is wrong here?!');
end
end
figure(options.fignum); % If figure does not exist, create new figure
if 1
%Frequency Axis
freq_vec = linspace(-obj.fs/2,obj.fs/2,length(psd));
freq_vec = reshape(freq_vec,size(psd_plot));
if nargin == 4
p = plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5,'DisplayName',options.displayname);
% plot(freq_vec*1e-9,psd_smoothed','Linewidth',1,'Color',[0 0 0],'DisplayName',[char(varargin{2}),' smoothed']);
else
p = plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5);
% plot(freq_vec*1e-9,psd_smoothed','Linewidth',1,'Color',[1 1 1],'LineStyle',':');
end
%xlim([freq_vec(1)/1e9-2 freq_vec(end)/1e9+2])
xlabel('frequency [GHz]')
else
%Wavelength Axis
freq_vec = physconst('LightSpeed')*linspace(-obj.fs/2,obj.fs/2,length(psd))./((physconst('LightSpeed')/1310e-9)^2)*1e9;
freq_vec = freq_vec+1310;
if nargin == 4
plot(freq_vec',psd_plot,'Linewidth',0.5,'DisplayName',options.displayname)
else
plot(freq_vec,psd_plot','Linewidth',0.5);
end
%xlim([freq_vec(1)/1e9-2 freq_vec(end)/1e9+2])
xlabel('wavelength [nm]')
end
hold on
end
% xlim([-150 150])
% ylim([-100,0]);
ylabel('magnitude [dBm]')
legend
grid minor;
end end
@@ -500,6 +644,7 @@ classdef Signal
if options.mode == delay_mode.samples if options.mode == delay_mode.samples
obj.signal=delayseq(obj.signal,delay); obj.signal=delayseq(obj.signal,delay);
% obj.signal=circshift(obj.signal,delay);
elseif options.mode == delay_mode.time elseif options.mode == delay_mode.time
@@ -534,25 +679,32 @@ classdef Signal
%estimate start pos of signal %estimate start pos of signal
maxpeaknum = floor(length(a)/length(b)); maxpeaknum = floor(length(a)/length(b));
[pks,pkpos] = findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum); [pks,pkpos] = findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend');
shifts = lags(pkpos); shifts = lags(pkpos);
%Cut occurences of ref signal from signal (only positive shifts) shifts = shifts(shifts>=0);
S = {}; S = {};
for c = shifts(shifts>0) isFlipped=0;
if numel(shifts) > 0
%Cut occurences of ref signal from signal (only positive shifts)
for c = shifts(shifts>=0)
sig = obj.delay(-c,'mode','samples'); sig = obj.delay(-c,'mode','samples');
sig.signal = sig.signal(1:length(b)); sig.signal = sig.signal(1:length(b));
S{end+1,1} = sig; S{end+1,1} = sig;
end end
% %
isFlipped=0;
if all(sign(co(pkpos))) if all(sign(co(pkpos)))
isFlipped = 1; isFlipped = 1;
end end
%return/keep the sinal with the highest correlation (only within positive shifts) %return/keep the sinal with the highest correlation (only within positive shifts)
[~,idx]=max(pks(shifts>0)); [~,idx]=max(pks(shifts>=0));
obj.signal = S{idx}.signal; obj.signal = S{idx}.signal;
%put signal with highest corr. to first index in S array %put signal with highest corr. to first index in S array
swap = S{1}; swap = S{1};
@@ -563,6 +715,12 @@ classdef Signal
S{c}.logbook = []; S{c}.logbook = [];
end end
else
%do nothing when shifts are negative or there are none...
end
%plot all synced signals and the ref signal %plot all synced signals and the ref signal
debug = 0; debug = 0;
if debug if debug
@@ -577,7 +735,7 @@ classdef Signal
end end
%%
function obj = filter(obj,a,b) function obj = filter(obj,a,b)
lbdesc = ['Filtering signal with H = a: ',num2str(a),' / b: ',num2str(b)]; lbdesc = ['Filtering signal with H = a: ',num2str(a),' / b: ',num2str(b)];
@@ -586,6 +744,7 @@ classdef Signal
obj.signal = filter(a,b,obj.signal); obj.signal = filter(a,b,obj.signal);
end end
%%
function er = extinctionratio(obj,fsym,M) function er = extinctionratio(obj,fsym,M)
histpoints = 1024; %% verticale resolution histpoints = 1024; %% verticale resolution
histpoints = floor(histpoints/2)*2+1; %% to have the eye digram centered around one point make the vertical resolution uneven histpoints = floor(histpoints/2)*2+1; %% to have the eye digram centered around one point make the vertical resolution uneven
@@ -665,7 +824,7 @@ classdef Signal
end end
%%
function eye(obj,fsym,M,options) function eye(obj,fsym,M,options)
arguments arguments

10
Classes/01_transmit/AWG.m
View File

@@ -8,6 +8,7 @@ classdef AWG < handle
upsampling_method upsampling_method
repetitions %repeat the signal to generate a longer sequence? repetitions %repeat the signal to generate a longer sequence?
fdac %needed fdac %needed
precomp_sinc_rolloff = 1;
normalize2dac %want to normalize at first? either 0 or 1 normalize2dac %want to normalize at first? either 0 or 1
bit_resolution %bit res. of quantizer (e.g. 5 bit) bit_resolution %bit res. of quantizer (e.g. 5 bit)
dac_min dac_min
@@ -33,7 +34,8 @@ classdef AWG < handle
arguments arguments
options.kover = 16; options.kover = 16;
options.upsampling_method upsampling_mode options.upsampling_method upsampling_mode = upsampling_mode.samplehold;
options.precomp_sinc_rolloff = 1;
options.repetitions = 1; options.repetitions = 1;
options.normalize2dac = 1; options.normalize2dac = 1;
options.fdac = 92e9; options.fdac = 92e9;
@@ -43,7 +45,7 @@ classdef AWG < handle
options.skew_active = 0; options.skew_active = 0;
options.awg_skew = 0; options.awg_skew = 0;
options.lpf_active = 0; options.lpf_active = 0;
options.lpf_type = 0; options.lpf_type = filtertypes.butterworth;
options.f_cutoff = 32e9; options.f_cutoff = 32e9;
options.H_lpf Filter options.H_lpf Filter
@@ -95,7 +97,7 @@ classdef AWG < handle
signalclass_in = obj.H_lpf.process(signalclass_in); signalclass_in = obj.H_lpf.process(signalclass_in);
else else
%4.B) just use a standard filter %4.B) just use a standard filter
lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",obj.lpf_type); lpf = Filter('filtdegree',5,"f_cutoff",obj.f_cutoff,"fs",obj.kover*obj.fdac,"filterType",obj.lpf_type);
signalclass_in = lpf.process(signalclass_in); signalclass_in = lpf.process(signalclass_in);
end end
end end
@@ -133,7 +135,7 @@ classdef AWG < handle
obj.signal_length = length(data_in); obj.signal_length = length(data_in);
%%%%%%%%% PRECOMP SINC ROLLOFF %%%%%%%%% %%%%%%%%% PRECOMP SINC ROLLOFF %%%%%%%%%
if 1 if obj.precomp_sinc_rolloff
% X: design FIR filter for sinc precomp % X: design FIR filter for sinc precomp
% https://www.dsprelated.com/showarticle/1191.php % https://www.dsprelated.com/showarticle/1191.php
ntaps = 13; ntaps = 13;

26
Classes/01_transmit/ChannelFreqResp.m
View File

@@ -202,19 +202,19 @@ classdef ChannelFreqResp < handle
function plot(obj) function plot(obj)
figure(55); figure(55);
%clf; clf;
Havg = obj.H; Havg = obj.H;
%1) %1)
subplot(2,1,1);hold on;box on;title('Magnitude Freq. Response'); subplot(2,1,1);hold all;box on;title('Magnitude Freq. Response');
plot(obj.faxis/1e9, 20*log10(abs(obj.H_all)),'linewidth',0.1,'LineStyle','-','Color','#808080') ; plot(obj.faxis/1e9, 20*log10(abs(obj.H_all)),'linewidth',0.1,'LineStyle','-','Color','#808080') ;
xlim([0.2 .5*max(obj.faxis)*1e-9]); xlim([0.2 .5*max(obj.faxis)*1e-9]);
plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2); plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2);
grid on; grid on;
%2) %2)
subplot(2,1,2); hold on; box on; title('Phase Freq. Response'); subplot(2,1,2); hold all; box on; title('Phase Freq. Response');
plot(obj.faxis/1e9, angle(obj.H_all),'linewidth',0.1,'LineStyle','-','Color','#808080') ; plot(obj.faxis/1e9, angle(obj.H_all),'linewidth',0.1,'LineStyle','-','Color','#808080') ;
plot(obj.faxis/1e9, unwrap(angle(Havg)),'LineWidth',2) ; plot(obj.faxis/1e9, unwrap(angle(Havg)),'LineWidth',2) ;
xlim([0.2 .5*max(obj.faxis)*1e-9]); xlim([0.2 .5*max(obj.faxis)*1e-9]);
@@ -228,7 +228,7 @@ classdef ChannelFreqResp < handle
Havg = Havg./mean(Havg(2:10)); Havg = Havg./mean(Havg(2:10));
%3) %3)
subplot(2,1,1); hold on; box on; title('Inverse Magnitude Freq. Response'); subplot(2,1,1); hold all; box on; title('Inverse Magnitude Freq. Response');
plot(obj.faxis/1e9, 20*log10(abs(1./Havg)),"LineWidth",2,"Color",[0.3467 0.5360 0.6907]) ; plot(obj.faxis/1e9, 20*log10(abs(1./Havg)),"LineWidth",2,"Color",[0.3467 0.5360 0.6907]) ;
xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on; xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on;
ylim([-1 15]); ylim([-1 15]);
@@ -236,17 +236,19 @@ classdef ChannelFreqResp < handle
yline(3,'LineWidth',2,'LineStyle','--'); yline(3,'LineWidth',2,'LineStyle','--');
%4) %4)
subplot(2,1,2); hold on; box on; title('Inverse Phase Freq. Response'); subplot(2,1,2); hold all; box on; title('Inverse Phase Freq. Response');
plot(obj.faxis/1e9, unwrap(angle(1./Havg)),"LineWidth",2,"Color",[0.3467 0.5360 0.6907]) ; plot(obj.faxis/1e9, unwrap(angle(1./Havg)),"LineWidth",2,"Color",[0.3467 0.5360 0.6907]) ;
xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on; xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on;
%%% plot for publication %%% plot for publication
figure(30);hold on;box on;title('Magnitude Freq. Response'); figure(1234);hold all;box on;title('Magnitude Freq. Response');
% xlim([0 max(obj.faxis)*1e-9]); %xlim([0.2 .5*max(obj.faxis)*1e-9]);
% ylim([-20, 10]); %ylim([-40, 2]);
fax = obj.faxis - obj.f_ref/2; Havg_smooth = smooth(Havg,50);
Havg = Havg ./ max(abs(Havg)); symaxis = (obj.faxis-(obj.f_ref/2))/1e9;
plot(fax/1e9, 20*log10(abs(fftshift(Havg)))+7,'LineWidth',2); Havg = fftshift(Havg);
%Havg = smooth(Havg);
plot(symaxis, 20*log10(abs(Havg)),'LineWidth',0.5);
grid on; grid on;
@@ -362,7 +364,7 @@ classdef ChannelFreqResp < handle
end end
% fprintf('Frequency response information successfully loaded from %s\n', fullFileName); fprintf('Frequency response information successfully loaded from %s\n', fullFileName);
end end
end end

55
Classes/01_transmit/PAMmapper.m
View File

@@ -23,7 +23,6 @@ classdef PAMmapper
obj.scaling = rms(obj.get_levels()); obj.scaling = rms(obj.get_levels());
end end
function out = map(obj,signal_in) function out = map(obj,signal_in)
@@ -41,11 +40,25 @@ classdef PAMmapper
end end
function signalclass_out = demap(obj,signalclass_in) function signal_out = demap(obj,signal_in)
signalclass_in.signal = obj.demap_(signalclass_in.signal); issignalclass = 0;
if isa(signal_in,'Signal')
signalclass = signal_in;
signal_in = signal_in.signal;
issignalclass = 1;
end
signal_out = obj.demap_(signal_in);
if issignalclass
lbdesc = ['Demap PAM ',num2str(obj.M),' symbols to bit stream']; lbdesc = ['Demap PAM ',num2str(obj.M),' symbols to bit stream'];
signalclass_in = signalclass_in.logbookentry(lbdesc,obj); signal_in = signalclass;
signalclass_out = signalclass_in; signal_in = signal_in.logbookentry(lbdesc,obj);
signal_in.signal = signal_out;
signal_out = signal_in;
end
end end
function pam_sig = map_(obj,bitpattern) function pam_sig = map_(obj,bitpattern)
@@ -156,6 +169,8 @@ classdef PAMmapper
case 6 %PAM 6 case 6 %PAM 6
thres = [-3 5;-1 5;-3 -5;-1 -5;-5 3;-5 1;-5 -3;-5 -1;-1 3;-1 1;-1 -3;-1 -1;-3 3;-3 1;-3 -3;-3 -1;3 5;1 5;3 -5;1 -5;5 3;5 1;5 -3;5 -1;1 3;1 1;1 -3;1 -1;3 3;3 1;3 -3;3 -1]; thres = [-3 5;-1 5;-3 -5;-1 -5;-5 3;-5 1;-5 -3;-5 -1;-1 3;-1 1;-1 -3;-1 -1;-3 3;-3 1;-3 -3;-3 -1;3 5;1 5;3 -5;1 -5;5 3;5 1;5 -3;5 -1;1 3;1 1;1 -3;1 -1;3 3;3 1;3 -3;3 -1];
% thres = [-4 -2 0 2 4];
% thres = thres ./ sqrt(10);
case 8 case 8
% 8-ASK % 8-ASK
@@ -194,7 +209,6 @@ classdef PAMmapper
end end
end end
function [data_out] = demap_(obj,data_in) function [data_out] = demap_(obj,data_in)
data_in= data_in'; data_in= data_in';
if obj.M ~= 6 if obj.M ~= 6
@@ -309,14 +323,35 @@ classdef PAMmapper
end end
function [Signal_out] = quantize(obj,Signal_in) function [Signal_out] = quantize(obj,Signal_in)
constellation = obj.get_levels(); constellation = obj.get_levels();
constellation = constellation ./ rms(constellation); constellation = constellation ./ rms(constellation);
Signal_out = Signal_in;
dist = abs(Signal_in.signal - constellation); issignalclass = 0;
if isa(Signal_in,'Signal')
issignalclass = 1;
Sig_class = Signal_in;
Signal_in = Signal_in.signal;
end
[~,high_dim_sig] = max(size(Signal_in));
[~,high_dim_const] = max(size(constellation));
if high_dim_sig == high_dim_const
Signal_in = Signal_in';
end
dist = abs(Signal_in - constellation);
[~,symbol_idx] = min(dist,[],2); % decision for closest constellation point [~,symbol_idx] = min(dist,[],2); % decision for closest constellation point
Signal_out.signal = constellation(symbol_idx); Signal_out = constellation(symbol_idx);
Signal_out = reshape(Signal_out,size(Signal_in));
if issignalclass
Sig_class.signal = Signal_out;
Signal_out = Sig_class;
end
Signal_out.signal = reshape(Signal_out.signal,size(Signal_in.signal));
end end

6
Classes/01_transmit/PAMsource.m
View File

@@ -151,6 +151,12 @@ classdef PAMsource
sym_max = max(symbols.signal); sym_max = max(symbols.signal);
end end
% symbols.move_it_spectrum("fignum",222,"displayname","Symbols only");
% symbols.spectrum("fignum",222,"displayname","Symbols only","normalizeTo0dB",1);
%%%%% Pulse-forming %%%%%% %%%%% Pulse-forming %%%%%%
if obj.applypulseform if obj.applypulseform
digi_sig = obj.pulseformer.process(symbols); digi_sig = obj.pulseformer.process(symbols);

2
Classes/01_transmit/Pulseformer.m
View File

@@ -4,6 +4,8 @@ classdef Pulseformer
properties(Access=public) properties(Access=public)
fdac fdac
end
properties(Access=private)
fsym fsym
pulse pulse
pulselength pulselength

52
Classes/04_DSP/Coding/Duobinary.m
View File

@@ -4,7 +4,6 @@ classdef Duobinary
properties(Access=public) properties(Access=public)
end end
methods (Access=public) methods (Access=public)
@@ -12,6 +11,7 @@ classdef Duobinary
function obj = Duobinary() function obj = Duobinary()
%NAME Construct an instance of this class %NAME Construct an instance of this class
end end
function signal = precode(~,signal) function signal = precode(~,signal)
@@ -78,7 +78,13 @@ classdef Duobinary
end end
function signal = encode(~,signal) function signal = encode(~,signal,options)
arguments
~
signal
options.M = [];
end
if isa(signal,'Signal') if isa(signal,'Signal')
data = signal.signal; data = signal.signal;
@@ -88,17 +94,19 @@ classdef Duobinary
data = data./rms(data); data = data./rms(data);
if isempty(options.M)
u = unique(data); u = unique(data);
M = numel(u); options.M = numel(u);
end
%make unipolar %make unipolar
if M == 4 if options.M == 4
data = data .* sqrt(5); data = data .* sqrt(5);
elseif M == 6 elseif options.M == 6
data = data .* sqrt(10); data = data .* sqrt(10);
elseif M == 8 elseif options.M == 8
data = data .* sqrt(21); data = data .* sqrt(21);
elseif M == 16 elseif options.M == 16
data = data .* sqrt(85); data = data .* sqrt(85);
warning('Check if PAM16 implementation, mapping and scaling is correct!') warning('Check if PAM16 implementation, mapping and scaling is correct!')
end end
@@ -108,7 +116,7 @@ classdef Duobinary
data = data - b; data = data - b;
data = data ./ 2; data = data ./ 2;
% assert(isequal((0:M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar % assert(isequal((0:options.M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar
% duobinary coding (1+D) % duobinary coding (1+D)
% coeff = [1,1]; % coeff = [1,1];
@@ -132,12 +140,14 @@ classdef Duobinary
mean_power = sum((unique_points .^ 2) .* probabilities); mean_power = sum((unique_points .^ 2) .* probabilities);
scaling_factor = sqrt(mean_power); scaling_factor = sqrt(mean_power);
if M == 4 if options.M == 4
data = data ./ sqrt(2.5); % 7-level constellation weighted with probability after DB code i.e. mean([-3 3 -2 -2 2 2 -1 -1 -1 1 1 1 0 0 0 0].^2) = 2.5 --> sqrt(2.5) == rms(constellation) data = data ./ sqrt(2.5); % 7-level constellation weighted with probability after DB code i.e. mean([-3 3 -2 -2 2 2 -1 -1 -1 1 1 1 0 0 0 0].^2) = 2.5 --> sqrt(2.5) == rms(constellation)
elseif M == 6 elseif options.M == 6
data = data ./ sqrt(5.8); data = data ./ sqrt(5.8);
elseif M == 8 elseif options.M == 8
data = data ./ sqrt(10.5); % 15-level constellation weighted with probability after DB code i.e. data = data ./ sqrt(10.5); % 15-level constellation weighted with probability after DB code i.e.
else
error("Error in: Duobinary encode > scale unipolar to bipolar > The data is not PAM4, PAM6 or PAM 8? ")
end end
if isa(signal,'Signal') if isa(signal,'Signal')
@@ -158,24 +168,30 @@ classdef Duobinary
end end
function signalclass = decode(~,signalclass) function signalclass = decode(~,signalclass,options)
arguments
~
signalclass
options.M = [];
end
data = signalclass.signal; data = signalclass.signal;
u = unique(data); u = unique(data);
I = numel(u); %number of duobinary coded const. points I = numel(u); %number of duobinary coded const. points
if isempty(options.M)
M = (I+1)/2; %PAM-M order M = (I+1)/2; %PAM-M order
else
M = options.M; %PAM-M order
end
%make unipolar %make unipolar
if I == 7 if I == 7 || I == 6
data = data .* sqrt(2.5); data = data .* sqrt(2.5);
elseif I == 11 elseif I == 11
%todo
data = data .* sqrt(5.8); data = data .* sqrt(5.8);
warning('Check db decode implementation, mapping and scaling is correct!')
elseif I == 15 elseif I == 15
data = data .* sqrt(10.5); data = data .* sqrt(10.5);
elseif I == 16
warning('Check db decode implementation, mapping and scaling is correct!')
end end
data = round(data); data = round(data);
@@ -196,6 +212,8 @@ classdef Duobinary
data = data ./ sqrt(10); data = data ./ sqrt(10);
elseif M == 8 elseif M == 8
data = data ./ sqrt(21); data = data ./ sqrt(21);
else
errormsg("Error in: Duobinary decode > scale unipolar to bipolar > The data is not PAM4, PAM6 or PAM 8? ")
end end
signalclass.signal = data; signalclass.signal = data;

10
Classes/04_DSP/Equalizer/EQ.m
View File

@@ -1,4 +1,4 @@
classdef EQ %< handle classdef EQ < handle
%EQ Summary of this class goes here %EQ Summary of this class goes here
% Detailed explanation goes here % Detailed explanation goes here
@@ -289,9 +289,9 @@ classdef EQ %< handle
e_dc = e_dc - obj.DCmu*error; e_dc = e_dc - obj.DCmu*error;
obj.error_log.e_ffe(cnt,trainloops) = e_ffe; % obj.error_log.e_ffe(cnt,trainloops) = e_ffe;
obj.error_log.e_dfe(cnt,trainloops) = e_dfe; % obj.error_log.e_dfe(cnt,trainloops) = e_dfe;
obj.error_log.e_(cnt,trainloops) = error; % obj.error_log.e_(cnt,trainloops) = error;
cnt = cnt+1; cnt = cnt+1;
if obj.Nb(1) > 0 if obj.Nb(1) > 0
@@ -460,7 +460,7 @@ classdef EQ %< handle
if 1%mu_mat ~= 0 if 1%mu_mat ~= 0
e_dc = e_dc - obj.DCmu*error; e_dc = e_dc - obj.DCmu*error;
error_log(cnt,dd_loop) = e_dc; % error_log(cnt,dd_loop) = e_dc;
cnt = cnt+1; cnt = cnt+1;
end end

98
Classes/04_DSP/Equalizer/Postfilter.m Normal file
View File

@@ -0,0 +1,98 @@
classdef Postfilter < handle
%NAME Summary of this class goes here
% Detailed explanation goes here
properties(Access=public)
ncoeff = 2;
coefficients = [];
useBurg = NaN
end
methods (Access=public)
function obj = Postfilter(options)
%NAME Construct an instance of this class
% Detailed explanation goes here
arguments
options.useBurg = NaN
options.coefficients = []
options.ncoeff double = 2
end
%
fn = fieldnames(options);
for n = 1:numel(fn)
try
obj.(fn{n}) = options.(fn{n});
end
end
% do more stuff
end
function signalclass_out = process(obj,signalclass_in,noiseclass_in,options)
arguments
obj
signalclass_in
noiseclass_in
options.useBurg = obj.useBurg;
options.coefficients = obj.coefficients;
end
if ~isnan(options.useBurg) && options.useBurg
disp('using burg alg')
obj.coefficients = arburg(noiseclass_in.signal,obj.ncoeff);
elseif ~isempty(options.coefficients)
disp('using given taps')
obj.coefficients = options.coefficients;
obj.useBurg = 0;
else
end
signalclass_in = signalclass_in.filter(obj.coefficients,1);
% append to logbook
lbdesc = ['Postfilter'];
signalclass_in = signalclass_in.logbookentry(lbdesc);
% write to output
signalclass_out = signalclass_in;
end
function showFilter(obj,noiseclass_in,options)
arguments
obj
noiseclass_in
options.fignum = 121
options.color = []
end
% noiseclass_in.spectrum('displayname','Noise PSD shifted to 0dBm','fignum',options.fignum,'normalizeTo0dB',1);
figure(options.fignum)
[h,w] = freqz(1,obj.coefficients,length(noiseclass_in),"whole",noiseclass_in.fs);
h = h/max(abs(h));
hold on
w_ = (w - noiseclass_in.fs/2);
if isempty(options.color)
plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(obj.coefficients,2)), ' '],'LineWidth',2);
else
plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(obj.coefficients,2)), ' '],'Color',options.color,'LineWidth',2);
end
ylim([-30,2]);
end
end
end

6
Classes/04_DSP/Equalizer/VNLE.m
View File

@@ -57,7 +57,7 @@ classdef VNLE < handle
end end
function [X] = process(obj, X, D) function [X,N] = process(obj, X, D)
% actual processing of the signal (steps 1. - 3.) % actual processing of the signal (steps 1. - 3.)
% 1 normalize RMS % 1 normalize RMS
@@ -91,6 +91,9 @@ classdef VNLE < handle
lbdesc = [num2str(obj.order),' tap FFE']; lbdesc = [num2str(obj.order),' tap FFE'];
X = X.logbookentry(lbdesc); % append to logbook X = X.logbookentry(lbdesc); % append to logbook
N = X;
N = X - D;
end end
@@ -156,7 +159,6 @@ classdef VNLE < handle
end end
if ~all(mu==0,'all') %mu has not only zeros if ~all(mu==0,'all') %mu has not only zeros
% obj.e = obj.e - (mu * err * x_in) ; % Weight update rule of LMS
obj.e = obj.e - ( (mu * x_in) * err ) ; % Weight update rule of LMS obj.e = obj.e - ( (mu * x_in) * err ) ; % Weight update rule of LMS
else else
normalizationfactor = (x_in.' * x_in); normalizationfactor = (x_in.' * x_in);

400
Classes/04_DSP/Sequence Detection/MLSE.m
View File

@@ -1,4 +1,4 @@
classdef MLSE classdef MLSE < handle
%MLSE calculates the most probable sequence for an input signal with given/ known channel impulse response of any length %MLSE calculates the most probable sequence for an input signal with given/ known channel impulse response of any length
properties(Access=public) properties(Access=public)
@@ -34,17 +34,22 @@ classdef MLSE
end end
function signalclass = process(obj,signalclass) function [signalclass_hd,signalclass_sd] = process(obj,signalclass,ref_symbolclass)
data_in = signalclass.signal; data_in = signalclass.signal;
data_ref = ref_symbolclass.signal;
data_out = obj.process_(data_in); [data_out_hd,data_out_sd] = obj.process_(data_in,data_ref);
signalclass.signal = data_out; signalclass_hd = signalclass;
signalclass_hd.signal = data_out_hd;
signalclass_sd = signalclass;
signalclass_sd.signal = data_out_sd;
end end
function data_out = process_(obj,data_in) function [VITERBI_ESTIMATION_SYMBOLS,soft_decisions] = process_(obj,data_in,data_ref)
% remove unnecessary zeros at start of impulse response to keep % remove unnecessary zeros at start of impulse response to keep
@@ -58,6 +63,17 @@ classdef MLSE
obj.DIR = [0 obj.DIR]; obj.DIR = [0 obj.DIR];
end end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% PREPARATIONS %%%%%%%%
%%%% Separate the equalized signal into the respective levels based on the actually transmitted level
constellation = unique(data_ref);
decisionLevels = (constellation(1:end-1) + constellation(2:end)) / 2;
tx_bits = PAMmapper(numel(constellation),0).demap(data_ref);
% impulse respnse i.e. [0.5, 1.0000]
obj.DIR = flip(obj.DIR); obj.DIR = flip(obj.DIR);
% RMS normalization of input data % RMS normalization of input data
@@ -77,6 +93,7 @@ classdef MLSE
% Calculate all possible input symbols for the desired impulse % Calculate all possible input symbols for the desired impulse
% response. Row number is the index of the previous state, % response. Row number is the index of the previous state,
% column number is the index of the next state % column number is the index of the next state
% noise free received == branch metrics
noise_free_received = zeros(length(states),length(states)); noise_free_received = zeros(length(states),length(states));
count_row = 1; count_row = 1;
count_col = 1; count_col = 1;
@@ -101,127 +118,302 @@ classdef MLSE
end end
end end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% FORWARD PASS %%%%%
% Initialize the output vector
pm = zeros(length(states),length(states));
bm_fw = zeros(length(states),length(states),length(data_in));
% Forward Recursion (FSM Computation)
for n = 1:length(data_in)
% bm = abs(data_in(n) - noise_free_received).^2;
% %% Optimized pm = pm + bm;
% % Preallocate and initialize variables [pm_survivor_fw(:,n),pm_survivor_fw_idx(:,n)] = min(pm,[],2); % choose lowest path metric as new state
% data_out = NaN(size(data_in)); % output vector pm = repmat(pm_survivor_fw(:,n).',length(states),1); % update pm (chosen state to 2nd dimension -> FROM state)
% sum_path_metrics_res = zeros(length(states), length(data_in));
% path_idx = zeros(length(states), length(data_in));
%
% % Precompute repmat size
% num_states = length(states);
% num_signals = numel(noise_free_received);
%
% % First trellis path (initialize)
% sum_path_metrics = zeros(num_states, num_states);
% path_metrics = (abs(data_in(1) - noise_free_received)).^2; % Use broadcasting instead of repmat
% sum_path_metrics = sum_path_metrics + path_metrics; % Compute initial path metrics
%
% [sum_path_metrics_res(:,1), path_idx(:,1)] = min(sum_path_metrics, [], 2); % Best path for first step
%
% % Preallocate path_metrics and sum_path_metrics to avoid reallocating in each loop
% path_metrics = zeros(num_states, num_signals);
%
% % Loop over remaining trellis paths
% for n = 2:length(data_in)
% % Avoid reallocation of sum_path_metrics, reuse the same matrix and update
% previous_sum_path_metrics = sum_path_metrics_res(:,n-1).'; % Transpose once for broadcasting
% sum_path_metrics = repmat(previous_sum_path_metrics, num_states, 1); % Avoid dynamic resizing
%
% % Calculate path metrics using broadcasting
% path_metrics = (abs(data_in(n) - noise_free_received)).^2; % Avoid repmat
%
% % Update sum path metrics
% sum_path_metrics = sum_path_metrics + path_metrics;
%
% % Find the best path for each state and store results
% [sum_path_metrics_res(:,n), path_idx(:,n)] = min(sum_path_metrics, [], 2);
% end
%
% %% Traceback
% ideal_path = NaN(1, length(data_in)+1); % Preallocate ideal path
% [~, ideal_path(length(data_in)+1)] = min(sum_path_metrics_res(:,length(data_in))); % Start from final state
%
% % Trace back through trellis
% for h = length(data_in):-1:1
% ideal_path(h) = path_idx(ideal_path(h+1),h); % Follow the best path back
% end
%
% % Extract the output indices
% idx_out = ideal_path(2:end);
bm_fw(:,:,n) = bm;
% OLD
% initilaize the output vector
data_out = NaN(size(data_in));
sum_path_metrics = zeros(length(states),length(states));
% first trellis path
% euclidian distance as path metric
path_metrics = (abs(repmat(data_in(1),size(noise_free_received)) - noise_free_received)).^2;
sum_path_metrics = sum_path_metrics + path_metrics; % calculation of all possible sum path metrics
[sum_path_metrics_res(:,1),path_idx(:,1)] = min(sum_path_metrics,[],2); % find the best path to each state, store sum path metric and predecessor for each state
% remaining trellis paths
for n = 2:length(data_in)
sum_path_metrics = repmat(sum_path_metrics_res(:,n-1).',length(states),1);
% path_metrics = (repmat(data_in(n),size(noise_free_received)) - noise_free_received).^2;%.*prob_mat;
path_metrics = (abs(repmat(data_in(n),size(noise_free_received)) - noise_free_received)).^2;
sum_path_metrics = sum_path_metrics + path_metrics;
[sum_path_metrics_res(:,n),path_idx(:,n)] = min(sum_path_metrics,[],2);
end end
%% trace back % we can now get the best path as min
ideal_path = NaN(1,length(data_in)+1); best_fw_path = NaN(1,length(data_in)+1);
% find ideal trellis path by going through the trellis
% backwards % find ideal trellis path by going through the trellis backwards
[~,ideal_path(length(data_in)+1)] = min(sum_path_metrics_res(:,length(data_in))); [~,best_fw_path(length(data_in)+1)] = min(pm_survivor_fw(:,length(data_in)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% BACKWARD PASS %%%%%
% Initialize the output vector
pm = zeros(length(states),length(states));
pm_survivor_bw = zeros(length(states),length(data_in)+1);
pm_survivor_bw_idx = zeros(length(states),length(data_in)+1);
bm_bw = zeros(length(states),length(states),length(data_in)+1);
% starting with the state that has the lowest sum path % starting with the state that has the lowest sum path
% metric, follow the stored information about the % metric, follow the stored information about the
% predecessor % predecessor
for h = length(data_in):-1:1 for h = length(data_in):-1:1
ideal_path(h) = path_idx(ideal_path(h+1),h);
end
idx_out = ideal_path(2:length(data_in)+1); best_fw_path(h) = pm_survivor_fw_idx(best_fw_path(h+1),h);
bm = abs(data_in(h) - noise_free_received).^2;
pm = pm + bm.';
[pm_survivor_bw(:,h),pm_survivor_bw_idx(:,h)] = min(pm,[],2); % choose lowest path metric as new state
pm = repmat(pm_survivor_bw(:,h).',length(states),1); % update pm (chosen state to 2nd dimension -> FROM state)
bm_bw(:,:,h) = bm;
if obj.duobinary_output
%use duobinary encoder, output is already scaled inside this
%one
data_out = Duobinary().encode(first_sym(idx_out));
else
%
data_out(1:length(data_in)) = first_sym(idx_out);
% scale to rms = 1 using the standard sqrt() expressions
if obj.M == 4
data_out = data_out./sqrt(5);
elseif obj.M == 6
data_out = data_out./sqrt(10);
elseif obj.M == 8
data_out = data_out./sqrt(21);
end
end end
VITERBI_ESTIMATION_IDX(1:length(data_in)) = first_sym(best_fw_path(2:end));
VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = constellation(best_fw_path(2:end));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% FORWARD PASS %%%%%
%calc the log probabilities (llp's)
for k = 2:length(data_in)
fsm = repmat(pm_survivor_fw(:,k-1)',[length(states),1]); %pm_survivor_fw size: length(states)xlength(sequence)
bm = bm_fw(:,:,k); %bm_fw size: length(states)xlength(states)xlength(sequence)
bsm = repmat(pm_survivor_bw(:,k)',[length(states),1]); %pm_survivor_bw size: length(states)xlength(sequence)+1
llp(:,k) = min(fsm+bm+bsm,[],2);
end
% directly decide based on lowest LLP index
[~,llp_based_state_seq]=min(llp);
LLP_EST(1:length(data_in)) = constellation(llp_based_state_seq);
rx_bits = PAMmapper(numel(constellation),0).demap(LLP_EST');
[~,~,ber_llp,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
fprintf('LLP BER : %.2e \n',ber_llp);
% [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% fprintf('LLP BER +1: %.2e \n',ber_llp);
% [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% fprintf('LLP BER -1: %.2e \n',ber_llp);
%%%% DECIDE based on Viterbi traceback
rx_bits = PAMmapper(numel(constellation),0).demap(VITERBI_ESTIMATION_SYMBOLS');
[~,~,ber_viterbi,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
fprintf('Viterbi BER: %.2e \n',ber_viterbi);
% [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% fprintf('Viterbi BER: %.2e \n',ber_viterbi);
% directly decide based on the FW path metrics
[~,fw_direct_state_seq]=min(pm_survivor_fw);
FW_EST(1:length(data_in)) = constellation(fw_direct_state_seq);
rx_bits = PAMmapper(numel(constellation),0).demap(FW_EST');
[~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
fprintf('FW BER: %.2e \n',ber_fw);
% directly decide based on the BW path metrics
[~,bw_direct_state_seq]=min(pm_survivor_bw);
BW_EST(1:length(data_in)) = constellation(bw_direct_state_seq(2:end));
rx_bits = PAMmapper(numel(constellation),0).demap(BW_EST');
[~,~,ber_bw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
fprintf('BW BER: %.2e \n',ber_bw);
% [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% fprintf('BW BER: %.2e \n',ber_viterbi);
% [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% fprintf('BW BER: %.2e \n',ber_viterbi);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tx_symbolpos = zeros(numel(constellation),length(data_ref));
est_symbolpos = zeros(numel(constellation),length(data_ref));
for lvl = 1:numel(constellation)
tx_symbolpos(lvl,data_ref==constellation(lvl)) = 1;
est_symbolpos(lvl,VITERBI_ESTIMATION_IDX==first_sym(lvl)) = 1;
est_symbolpos_llm(lvl,llp_based_state_seq==lvl) = 1;
end
est_symbolpos= logical(est_symbolpos);
tx_symbolpos = logical(tx_symbolpos);
est_symbolpos_llm = logical(est_symbolpos_llm);
figure(299);clf;hold on;
llp_filt = NaN(length(states),length(data_in));
llp_ = llp - min(llp,[],1);
for c = 2%:numel(constellation)
for c2 = 1:3%numel(constellation)
idx = est_symbolpos_llm(c,:);
llp_filt(c2,idx) = (llp(c2,idx));
plotidx = 2:200000;
scatter(plotidx,llp_filt(c2,plotidx),1,'o','LineWidth',1,'DisplayName',sprintf('LLP of Symbol %d | %d was transmitted',c2,c));
end
end
% Assume llp is a 4 x N matrix (4 PAM-4 symbols, N time steps)
[numSymbols, numTimeSteps] = size(llp);
P_symbols = zeros(numSymbols, numTimeSteps); % To hold the soft output probabilities
for k = 1:numTimeSteps
% Compute the exponentials. (Use -llp_shifted if llp's are costs.)
exponents = -llp_(:, k);
% Normalize to form a probability vector.
P_symbols(:, k) = exponents / sum(exponents);
end
% Optionally, if you prefer a single soft output value per symbol (an expectation),
% define your PAM-4 constellation levels, e.g.:
soft_output = sum(P_symbols .* constellation, 1); % 1 x N vector of soft outputs
%%%% BITWISE LLR's ??? %%%%%
figure(100);
clf
hold on
title('LLR between inner and outer bits')
bitmapping = PAMmapper(numel(constellation),0).demap(constellation);
for bp = 1:size(bitmapping,2)
pos_bitone = find(bitmapping(:,bp)==1);
pos_bitzero = find(bitmapping(:,bp)~=1);
llr_bits(bp,:) = min(llp(pos_bitone,:),[],1) - min(llp(pos_bitzero,:),[],1);
subplot(size(bitmapping,2),1,bp)
scatter(1:length(data_ref),llr_bits(bp,:),1,'.');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% From: Log-Likelihood Probabilities LLP --> To: Log-Likelihood Ratios LLR
for s = 1:length(states)-1
llr(s,:) = llp_(s,:) - llp_(s+1,:); % subtract llp's of successive const. points to get llr's
end
% Build Sets of LLR's that contain only those values at
% timepoint k, where the symbols:
% a) were actually transmitted: set "S"
% b) were decoded by viterbi: set "SC"
S = NaN(length(states)-1,length(data_in));
SC = NaN(length(states)-1,length(data_in));
llp_inf = llp_;
llp_inf(llp_==0) = Inf;
[~,scnd_idx] = min(llp_inf,[],1) ;
for s = 1:length(states)-1
%%% SET "S"
% find all time indices k, where const point s was actually transmitted
indice = tx_symbolpos(s,:)==1;
S(s,indice) = llr(s,indice);
%calc mean of all llr's where symbol s was transmitted
K(s,1) = mean(S(s,indice),'omitnan');
% find all time indices k, where const. point s+1 was actually transmitted
indice_plusone = tx_symbolpos(s+1,:)==1;
S(s,indice_plusone) = llr(s,indice_plusone);
K(s,2) = mean(S(s,indice_plusone),'omitnan');
%%% SET "SC"
% find all time indices k, where symbol s was decoded
% idx: find positions in time where a symbol s was decoded
idx = est_symbolpos_llm(s,:);
% idx 2: find positions where the strongest competitor is
% s+1, i.e. the second best llp is at s+1
idx2 = scnd_idx == s+1;
SC(s,idx&idx2) = llr(s,idx&idx2);
KC(s,1) = mean(SC(s,idx&idx2),'omitnan');
% find all time indices k, where symbol s+1 was decoded
% idx: find positions in time where a symbol s+1 was decoded
idx = est_symbolpos_llm(s+1,:);
idx2 = scnd_idx == s;
% idx 2: find positions where the strongest competitor is
% s, i.e. the second best llp is at s
SC(s,idx&idx2) = llr(s,idx&idx2);
KC(s,2) = mean(SC(s,idx&idx2),'omitnan');
% scale the sets, using the average (?) of the
llrcn(s,:) = SC(s,:) * ( constellation(s+1)-constellation(s) ) ./ ( K(s,2)-K(s,1)) + decisionLevels(s);
end
figure(111)
title("Bla")
clf
for s = 1:length(states)-1
figure(1111)
hold on
title(sprintf('llrcn = llp %d - llp %d',s, s+1,s, s+1));
scatter(1:length(llrcn),llrcn(s,:),1,'.');
% STUFENARTIGE LLR'S UND LLP'S %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(111)
subplot(1,length(states),s)
hold on
title(sprintf('llr = llp %d - llp %d \n SC=llr(tx sym = %d or %d)',s, s+1,s, s+1));
scatter(1:length(llr),llr(s,:),1,'.');
scatter(1:length(SC),SC(s,:),1,'.');
yline([K(s,1),K(s,2)]);
low = min(min(llr));
hi = max(max(llr));
ylim([low,hi]);
subplot(1,length(states),length(states))
hold on
scatter(1:length(llr),llr(s,:),1,'.');
ylim([low,hi]);
% HISTOGRAMME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(112)
subplot(length(states),1,s)
hold on
title(sprintf('histogram of llr; between const point %d - %d',s, s+1));
histogram(llrcn(s,:),10000,'EdgeAlpha',0)
% histogram(SC(s,:),10000,'EdgeAlpha',0)
xlim([-20, 20]);
subplot(length(states),1,length(states))
hold on
histogram(llrcn(s,:),10000,'EdgeAlpha',0)
% histogram(SC(s,:),10000,'EdgeAlpha',0)
% xlim([-20, 20]);
end
softdecisions = mean(llrcn,1,'omitnan');
distance = abs(VITERBI_ESTIMATION_SYMBOLS-llrcn);
[win_cost,win_idx] = min(distance,[],1);
for i = 1:length(data_in)
soft_decisions(i) = llrcn(win_idx(i),i);
end
soft_decisions = max(llrcn,[],1);
showLevelHistogram(soft_decisions,data_ref)
figure()
% scatter(1:length(data_in),VITERBI_ESTIMATION_SYMBOLS,1,'.');
scatter(1:length(data_in),soft_decisions,1,'.');
MLM_ESTIMATION(1:length(data_in)) = PAMmapper(numel(constellation),0).quantize(soft_decisions)';
rx_bits = PAMmapper(numel(constellation),0).demap(MLM_ESTIMATION');
[~,~,ber_mlm,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
fprintf('MLM BER: %.2e \n',ber_mlm);
end end
end end
methods (Access=private) methods (Access=private)
% Cant be seen from outside! So put all your functions here that can/
% shall not be called from outside
end end
end end

177
Classes/04_DSP/Sequence Detection/MLSE_viterbi.m Normal file
View File

@@ -0,0 +1,177 @@
classdef MLSE_viterbi < handle
%MLSE calculates the most probable sequence for an input signal with given/ known channel impulse response of any length
properties(Access=public)
M %PAM-M
DIR
trellis_states
duobinary_output
end
methods (Access=public)
function obj = MLSE_viterbi(options)
%NAME Construct an instance of this class
% Detailed explanation goes here
arguments
options.M double = 4;
options.DIR double = [1];
options.trellis_states double = [-3 -1 1 3];
options.duobinary_output logical = false;
end
%
fn = fieldnames(options);
for n = 1:numel(fn)
try
obj.(fn{n}) = options.(fn{n});
end
end
% do more stuff
end
function signalclass = process(obj,signalclass)
data_in = signalclass.signal;
data_out = obj.process_(data_in);
signalclass.signal = data_out;
end
function data_out = process_(obj,data_in)
% remove unnecessary zeros at start of impulse response to keep
% number of trellis states minimal
DIR_nonzero = find(obj.DIR ~= 0);
if DIR_nonzero(1) > 1
obj.DIR(1:DIR_nonzero(1)-1) = [];
end
if isscalar(obj.DIR)
obj.DIR = [0 obj.DIR];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%% WORKING
% impulse respnse i.e. [0.5, 1.0000]
obj.DIR = flip(obj.DIR);
% RMS normalization of input data
data_in = data_in ./ rms(data_in);
% seems to be the only way to use combvec for a flexible amount
% of vectors. 'combs' contains all trellis states
pre_comb_mat = repmat(obj.trellis_states,length(obj.DIR)-1,1);
pre_comb_cell = mat2cell(pre_comb_mat,ones(1,size(pre_comb_mat,1)),size(pre_comb_mat,2));
combs = fliplr(combvec(pre_comb_cell{:}).');
% Save first and last symbol of each state
first_sym = combs(:,1);
last_sym = combs(:,end);
states = sum(combs,2);
% Calculate all possible input symbols for the desired impulse
% response. Row number is the index of the previous state,
% column number is the index of the next state
% noise free received == branch metrics
noise_free_received = zeros(length(states),length(states));
count_row = 1;
count_col = 1;
for l1 = 1:length(states)
for l2 = 1:length(states)
if sum(combs(l2,2:end) == combs(l1,1:end-1)) == size(combs,2)-1
noise_free_received(count_row,count_col) = sum(combs(l2,:).*obj.DIR(end:-1:2)) + last_sym(l1)*obj.DIR(1);
else
noise_free_received(count_row,count_col) = inf;
end
count_row = count_row + 1;
end
count_col = count_col + 1;
count_row = 1;
end
% match amplitude levels of input signal to those of the calculated ideal symbols
% i.e. match the rms values of data_in to noise_free_received
if isreal(data_in)
if obj.M == round(obj.M)
data_in = data_in * rms(noise_free_received(noise_free_received ~= inf),'all','omitnan');
end
end
% OLD
% initilaize the output vector
data_out = NaN(size(data_in));
sum_path_metrics = zeros(length(states),length(states));
% first trellis path
% euclidian distance as path metric
path_metrics = (abs(repmat(data_in(1),size(noise_free_received)) - noise_free_received)).^2;
sum_path_metrics = sum_path_metrics + path_metrics; % calculation of all possible sum path metrics
[sum_path_metrics_res(:,1),path_idx(:,1)] = min(sum_path_metrics,[],2); % find the best path to each state, store sum path metric and predecessor for each state
% remaining trellis paths
for n = 2:length(data_in)
sum_path_metrics = repmat(sum_path_metrics_res(:,n-1).',length(states),1);
% path_metrics = (repmat(data_in(n),size(noise_free_received)) - noise_free_received).^2;%.*prob_mat;
path_metrics = (abs(repmat(data_in(n),size(noise_free_received)) - noise_free_received)).^2;
sum_path_metrics = sum_path_metrics + path_metrics;
[sum_path_metrics_res(:,n),path_idx(:,n)] = min(sum_path_metrics,[],2);
end
%% trace back
ideal_path = NaN(1,length(data_in)+1);
% find ideal trellis path by going through the trellis
% backwards
[~,ideal_path(length(data_in)+1)] = min(sum_path_metrics_res(:,length(data_in)));
% starting with the state that has the lowest sum path
% metric, follow the stored information about the
% predecessor
for h = length(data_in):-1:1
ideal_path(h) = path_idx(ideal_path(h+1),h);
end
idx_out = ideal_path(2:length(data_in)+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%% WORKING
if obj.duobinary_output
%use duobinary encoder, output is already scaled inside this
%one
data_out = Duobinary().encode(first_sym(idx_out));
else
%
data_out(1:length(data_in)) = first_sym(idx_out);
% scale to rms = 1 using the standard sqrt() expressions
if obj.M == 4
data_out = data_out./sqrt(5);
elseif obj.M == 6
data_out = data_out./sqrt(10);
elseif obj.M == 8
data_out = data_out./sqrt(21);
end
end
end
end
methods (Access=private)
% Cant be seen from outside! So put all your functions here that can/
% shall not be called from outside
end
end

214
Classes/04_DSP/TransmissionPerformance.m Normal file
View File

@@ -0,0 +1,214 @@
classdef TransmissionPerformance
% TransmissionPerformance
%
% This class calculates the best possible net data rate (in bits/s)
% from a given gross rate (bits/s) and a measured channel quality
% parameter (either bit error ratio (BER) or NGMI). The class uses
% lookup tables for three different FEC modes:
%
% 1. SD+HD concatenated (uses overall code rate and an NGMI threshold)
% 2. SD+HD concatenated with punctured LDPC (uses overall code rate and an NGMI threshold)
% 3. HD-only (uses a code rate and a BER threshold)
%
% The basic algorithm is as follows:
% 1. Given a measured quality (NGMI or BER) for each measurement, find in
% the table the first (i.e. “best”) entry that is satisfied by the
% measured value. For NGMI the condition is measured NGMI >= threshold,
% while for BER the condition is measured BER <= threshold.
% 2. Use the corresponding overall code rate to compute:
%
% net rate = gross rate * code rate.
%
% USAGE EXAMPLE (array inputs):
%
% tp = TransmissionPerformance;
% % Suppose we have 3 measurements:
% % Gross rates: 10, 15, and 20 Gbps
% % Measured NGMI: 0.92, 0.88, and 0.85
% % Measured BER: 1e-4, 7e-3, and 2e-2
% netrates = tp.calculatenetrate([10e9 15e9 20e9], ...
% 'NGMI', [0.92 0.88 0.85], ...
% 'BER', [1e-4 7e-3 2e-2]);
%
% The returned structure netrates will contain the net rates (and the
% used code rates and thresholds) for each mode in vector form.
properties (Constant)
%% Lookup Table for SD+HD concatenated (using NGMI)
% Overall code rate and NGMI threshold for each row.
OVERALLCODERATES_SDHD = [0.7519, 0.7602, 0.7684, 0.7766, 0.7850, ...
0.7932, 0.8014, 0.8098, 0.8180, 0.8262, ...
0.8345, 0.8428, 0.8510, 0.8593, 0.8676, ...
0.8733, 0.8790, 0.8848, 0.8905, 0.8962, ...
0.9019, 0.9077, 0.9134, 0.9191, 0.9248, ...
0.9306, 0.9363, 0.9420, 0.9477, 0.9535];
NGMITHRESHOLDS_SDHD = [0.8116, 0.8167, 0.8241, 0.8317, 0.8401, ...
0.8459, 0.8512, 0.8574, 0.8685, 0.8746, ...
0.8829, 0.8892, 0.8958, 0.9022, 0.9090,...
0.9150, 0.9210, 0.9270, 0.9330, 0.9390, ...
0.9450, 0.9510, 0.9570, 0.9630, 0.9690, ...
0.9750, 0.9810, 0.9870, 0.9930, 0.9990];
ISPUNCT_LDPI = [zeros(1,15),ones(1,15)];
%% Lookup Table for HD-only mode (using BER)
% For HD-only, the table gives the code rate and the maximum acceptable BER.
CODE_RATES_HD = [0.7500, 0.8000, 0.8123, 0.8333, 0.8571, ...
0.8750, 0.8889, 0.9000, 0.9091, 0.9167, 0.9412];
BERTHRESHOLDS_HD = [2.12e-2, 1.76e-2, 1.62e-2, 1.44e-2, 1.25e-2, ...
1.03e-2, 9.29e-3, 8.33e-3, 7.54e-3, 7.04e-3, 4.70e-3];
%% LUT for KP4-FEC and Inner Code https://grouper.ieee.org/groups/802/3/dj/public/23_03/patra_3dj_01b_2303.pdf
CODE_RATE_KP4_AND_INNER = [0.885799];
BERTHRESHOLDS_KP4_AND_INNER = 4.85e-3;
end
methods
function netrates = calculateNetRate(obj, grossRate, varargin)
% calculatenetrate Calculate the net data rate(s) from measured data.
%
% netrates = tp.calculatenetrate(grossRate, 'BER', berValue, 'NGMI', ngmiValue)
%
% INPUTS:
% grossRate - (scalar or vector) gross data rate(s) [bits/s]
%
% Optional name/value pairs:
% 'BER' - (scalar or vector) measured bit error ratio(s)
% 'NGMI' - (scalar or vector) measured NGMI value(s)
%
% At least one of 'BER' or 'NGMI' must be provided.
%
% OUTPUT:
% netrates - a structure with the following fields (if available):
% .SDHD - for SD+HD concatenated (using NGMI)
% .SDHD_LDPC - for SD+HD with punctured LDPC (using NGMI)
% .HD - for HD-only (using BER)
%
% Each sub-structure contains fields:
% .CodeRate - the chosen overall code rate from the table (vector)
% .Threshold - the threshold value used from the table (vector)
% .<BER/NGMI> - net rate computed as grossRate * CodeRate (vector)
% Parse inputs.
p = inputParser;
addRequired(p, 'grossRate', @(x) isnumeric(x));
addParameter(p, 'BER', [], @(x) isnumeric(x));
addParameter(p, 'NGMI', [], @(x) isnumeric(x));
parse(p, grossRate, varargin{:});
ber = p.Results.BER;
ngmi = p.Results.NGMI;
if isempty(ber) && isempty(ngmi)
error('At least one of ''BER'' or ''NGMI'' must be provided.');
end
% Determine the number of measurements.
numMeasurements = max([numel(grossRate), numel(ngmi), numel(ber)]);
% Broadcast scalars if needed.
if isscalar(grossRate) && numMeasurements > 1
grossRate = repmat(grossRate, 1, numMeasurements);
elseif numel(grossRate) ~= numMeasurements
error('grossRate must be scalar or have %d elements.', numMeasurements);
end
if ~isempty(ngmi)
if isscalar(ngmi) && numMeasurements > 1
ngmi = repmat(ngmi, 1, numMeasurements);
elseif numel(ngmi) ~= numMeasurements
error('NGMI must be scalar or have %d elements.', numMeasurements);
end
end
if ~isempty(ber)
if isscalar(ber) && numMeasurements > 1
ber = repmat(ber, 1, numMeasurements);
elseif numel(ber) ~= numMeasurements
error('BER must be scalar or have %d elements.', numMeasurements);
end
end
% Initialize the output structure.
netrates = struct;
if ~isempty(ngmi)
netrates.SDHD.GrossRate = NaN(1, numMeasurements);
netrates.SDHD.NetRate = NaN(1, numMeasurements);
netrates.SDHD.punctLDPI = NaN(1, numMeasurements);
netrates.SDHD.CodeRate = NaN(1, numMeasurements);
netrates.SDHD.Threshold = NaN(1, numMeasurements);
end
if ~isempty(ber)
netrates.HD.GrossRate = NaN(1, numMeasurements);
netrates.HD.NetRate = NaN(1, numMeasurements);
netrates.HD.CodeRate = NaN(1, numMeasurements);
netrates.HD.Threshold = NaN(1, numMeasurements);
netrates.KP4_hamming.GrossRate = NaN(1, numMeasurements);
netrates.KP4_hamming.NetRate = NaN(1, numMeasurements);
netrates.KP4_hamming.CodeRate = NaN(1, numMeasurements);
netrates.KP4_hamming.Threshold = NaN(1, numMeasurements);
end
% Process each measurement individually.
for i = 1:numMeasurements
% --- NGMI-based modes ---
if ~isempty(ngmi)
% SD+HD concatenated mode.
idx = find(obj.NGMITHRESHOLDS_SDHD <= ngmi(i), 1, 'last');
if ~isempty(idx)
codeRate = obj.OVERALLCODERATES_SDHD(idx);
netrates.SDHD.NetRate(i) = grossRate(i) * codeRate;
netrates.SDHD.GrossRate(i) = grossRate(i) ;
netrates.SDHD.CodeRate(i) = codeRate;
netrates.SDHD.Threshold(i) = obj.NGMITHRESHOLDS_SDHD(idx);
netrates.SDHD.punctLDPI(i) = obj.ISPUNCT_LDPI(idx);
end
end
% --- BER-based mode (HD-only) ---
if ~isempty(ber)
% Loop through the BER thresholds from the best (highest code rate)
% to the worst until the measured BER is acceptable.
idxBER = [];
for j = length(obj.BERTHRESHOLDS_HD):-1:1
if ber(i) <= obj.BERTHRESHOLDS_HD(j)
idxBER = j;
break;
end
end
if ~isempty(idxBER)
codeRate = obj.CODE_RATES_HD(idxBER);
netrates.HD.NetRate(i) = grossRate(i) * codeRate;
netrates.HD.GrossRate(i) = grossRate(i) ;
netrates.HD.CodeRate(i) = codeRate;
netrates.HD.Threshold(i) = obj.BERTHRESHOLDS_HD(idxBER);
end
idxBER = [];
for j = length(obj.BERTHRESHOLDS_KP4_AND_INNER):-1:1
if ber(i) <= obj.BERTHRESHOLDS_KP4_AND_INNER(j)
idxBER = j;
break;
end
end
if ~isempty(idxBER)
codeRate = obj.CODE_RATE_KP4_AND_INNER(idxBER);
netrates.KP4_hamming.NetRate(i) = grossRate(i) * codeRate;
netrates.KP4_hamming.GrossRate(i) = grossRate(i) ;
netrates.KP4_hamming.CodeRate(i) = codeRate;
netrates.KP4_hamming.Threshold(i) = obj.BERTHRESHOLDS_KP4_AND_INNER(idxBER);
end
end
end
end
end
end

660
Classes/DataBaseHandler/DBHandler.m Normal file
View File

@@ -0,0 +1,660 @@
classdef DBHandler < handle
% DBHANDLER Class to handle database queries
% This class provides methods to interact with an SQLite database, including
% inserting data, retrieving table names, and appending new rows.
properties
conn % Database connection object
pathToDB % Path to the SQLite database
tableNames % Cell array containing names of all tables in the database
tables = struct(); % Structure containing MATLAB tables for each database table
distinctValues
end
methods
function obj = DBHandler(options)
% DBHANDLER Constructor for the DBHandler class
% Initializes the database connection and retrieves table and field names.
%
% Usage:
% obj = DBHandler('pathToDB', 'path/to/database.db');
arguments
options.pathToDB = ""; % Default value for pathToDB if not provided
end
% Assign values to class properties based on input arguments
fn = fieldnames(options);
for n = 1:numel(fn)
try
obj.(fn{n}) = options.(fn{n});
end
end
% Establish a connection to the SQLite database
try
obj.conn = sqlite(obj.pathToDB);
catch e
error('Failed to connect to the database: %s', e.message);
end
if obj.dbIsHealthy
obj.refresh();
else
error('DB seems to be corrupt')
end
end
function obj = refresh(obj)
% Get table names and the first rows of each table to understand the structure
obj.getTableNames();
obj.getTables();
obj.getDistinctValues();
end
function obj = getTableNames(obj)
% Get all table names from the database
try
result = fetch(obj.conn, 'SELECT name FROM sqlite_master WHERE type="table"');
obj.tableNames = result.name;
catch e
error('Failed to retrieve table names: %s', e.message);
end
end
function obj = getTables(obj)
% Get a preview (first row) of each table to understand its structure
for i = 1:numel(obj.tableNames)
try
tableName = obj.tableNames{i};
results = fetch(obj.conn, sprintf('SELECT * FROM %s WHERE 1 = 2', tableName));
% Matlab cant handle if there is a NULL in a returned
% datarow... therefore do not return a row using the
% above condition which is never true
% results = sqlread(obj.conn, tableName, MaxRows=1);
for l = 1:numel(results.Properties.VariableNames)
varName = results.Properties.VariableNames{l};
obj.tables.(tableName).(varName) = []; % Store the preview as a reference
end
catch e
warning('Failed to read the table %s: %s', tableName, e.message);
end
end
end
function obj = getDistinctValues(obj)
% getDistinctValues Retrieves distinct values for each relevant field in all tables
% excluding fields ending with "_id". Stores distinct values in the 'distinctValues'
% property.
% Initialize a structure to store distinct values for each table
obj.distinctValues = struct();
% Iterate over each table in obj.tables
tableNames_ = fieldnames(obj.tables);
for i = 1:numel(tableNames_)
tableName = tableNames_{i};
% Initialize a sub-struct to store distinct values for each field in the table
obj.distinctValues.(tableName) = struct();
% Get all fields of the current table
fieldNames = fieldnames(obj.tables.(tableName));
% Iterate over each field
for j = 1:numel(fieldNames)
fieldName = fieldNames{j};
% % Skip fields ending with '_id' as they don't contain useful distinct values
% if endsWith(fieldName, '_id')
% continue;
% end
% Construct SQL to get distinct values for the current field
query = sprintf('SELECT DISTINCT %s FROM %s', fieldName, tableName);
% Execute query and fetch distinct values
try
result = fetch(obj.conn, query);
% Store the distinct values in the structure
if ~isempty(result)
distinctValues = table2array(result);
else
distinctValues = [];
end
obj.distinctValues.(tableName).(fieldName) = distinctValues;
catch e
% warning('Failed to retrieve distinct values for %s.%s: %s', tableName, fieldName, e.message);
obj.distinctValues.(tableName).(fieldName) = [];
end
end
end
end
function healthyDB = dbIsHealthy(obj)
healthyDB = false;
num_runs = obj.fetch('SELECT COUNT(*) AS total_runs FROM Runs');
num_configs = obj.fetch('SELECT COUNT(*) AS total_configurations FROM Configurations');
num_meas = obj.fetch('SELECT COUNT(*) AS total_measurements FROM Measurements');
assert((num_runs{1,1}==num_configs{1,1})&&(num_configs{1,1}==num_meas{1,1}),'Different num of entries per table');
%Check for any duplicate paths
duplictae_raw = obj.fetch("SELECT COALESCE(Runs.rx_raw_path,'NaN') AS rx_raw_path, COUNT(*) AS occurrences FROM Runs GROUP BY rx_raw_path HAVING COUNT(*) > 1");
duplictae_sync = obj.fetch("SELECT COALESCE(Runs.rx_sync_path,'NaN') AS rx_sync_path, COUNT(*) AS occurrences FROM Runs GROUP BY rx_sync_path HAVING COUNT(*) > 1");
if size(duplictae_raw,2) > 1
for i = 1:size(duplictae_raw,1)
fprintf('Raw Rx Paths: Found %d duplictaes of %s \n',duplictae_raw.occurrences(i),duplictae_raw.rx_raw_path(i));
end
end
healthyDB = true;
end
function lastID = appendToTable(obj, tableName, newRow)
% appendToTable Appends a new row to the specified table
%
% Usage:
% appendToTable(tableName, newRow)
%
% Inputs:
% tableName: The name of the table to append data to.
% newRow: A MATLAB table or struct containing the new row to be appended.
% Check if the table exists in the fetched tables
if ~isfield(obj.tables, tableName)
error('Table %s does not exist in the database or has not been fetched.', tableName);
end
% Convert newRow to a table if it is a struct
if isstruct(newRow)
fields = fieldnames(newRow);
emptyFields = structfun(@isempty,newRow);
if sum(emptyFields)>0
newRow.(fields{emptyFields==1}) = NaN;
disp(['In Table: ',tableName,': ',fields{emptyFields==1},' was empty, is now NaN ',newRow.(fields{emptyFields==1})])
end
newRow = struct2table(newRow);
end
% Ensure the new row matches the structure of the existing table
existingTableStructure = obj.tables.(tableName);
% Perform data type checks and conversions
for colName = newRow.Properties.VariableNames
% Extract the value and its intended column type
value = newRow.(colName{1});
existingValue = existingTableStructure.(colName{1});
% If the value is a class object, convert it to JSON format
if isobject(value) && ~isdatetime(value) && ~isa(value,"string")
newRow.(colName{1}) = string(jsonencode(value));
% If the value is a character array, convert it to a string
elseif ischar(value)
newRow.(colName{1}) = string(value);
end
end
% Append the new row to the database table
try
sqlwrite(obj.conn, tableName, newRow);
% disp(['Successfully appended new row to the table ', tableName]);
catch e
error('Failed to append to the table %s: %s', tableName, e.message);
end
% Retrieve the measurement_id of the newly inserted row for linking other tables
result = fetch(obj.conn, 'SELECT last_insert_rowid()');
lastID = result{1, 1}; % Access the value directly from the table
end
function exists = checkIfRunExists(obj, table2check, column2check, value2check)
% checkIfRunExists Checks if a specific value exists in a specified column of a table
%
% Usage:
% exists = checkIfRunExists(table2check, column2check, value2check)
%
% Inputs:
% table2check: The name of the table to check for duplicates.
% column2check: The name of the column to check within the specified table.
% value2check: The value to check for in the specified column.
%
% Outputs:
% exists: Boolean indicating whether the value already exists in the table.
% Ensure the specified table and column exist in the database
if ~isfield(obj.tables, table2check)
error('Table %s does not exist in the database.', table2check);
end
% Ensure the specified column exists in the table structure
if ~isfield(obj.tables.(table2check), column2check)
error('Column %s does not exist in the table %s.', column2check, table2check);
end
% Construct the query to check for the value in the specified column
query = sprintf('SELECT COUNT(*) FROM %s WHERE %s = "%s"', table2check, column2check, value2check);
% Execute the query and pass the value2check to avoid SQL injection issues
try
result = fetch(obj.conn, query);
count = result{1, 1}; % Extract the count from the result
catch e
error('Failed to execute the duplicate check query: %s', e.message);
end
% If count is greater than 0, then the value exists in the table
exists = count > 0;
if exists
disp(['The value "', value2check, '" already exists in the column "', column2check, '" of the table "', table2check, '".']);
else
% disp(['The value "', value2check, '" does not exist in the column "', column2check, '" of the table "', table2check, '".']);
end
end
function addBEREntry(obj, berValue, occurrence, runID, ffe, dfe, mlse, pf, eqType, ffe_order, dfe_order, len_tr, mu_ffe, mu_dfe, mu_dc, comment)
% addBEREntry Adds a BER entry linked to an existing or new Equalizer entry.
% Usage:
% addBEREntry(runID, eq, ffe, dfe, mlse, pf, eqType, ffe_order, dfe_order, len_tr, mu_ffe, mu_dfe, mu_dc, berValue, comment)
if isempty(pf)
postfilter_taps = [];
else
postfilter_taps = pf.burg_coeff;
end
% Create equalizer data struct for searching and adding if necessary
equalizerData = struct( ...
'ffe', jsonencode(ffe), ...
'dfe', jsonencode(dfe), ...
'mlse', jsonencode(mlse), ...
'pf', jsonencode(pf), ...
'eq_type', string(eqType), ...
'ffe_order', jsonencode(ffe_order), ...
'dfe_order', jsonencode(dfe_order), ...
'postfilter_taps',jsonencode(postfilter_taps),...
'len_tr', len_tr, ...
'mu_ffe', jsonencode(mu_ffe), ...
'mu_dfe', mu_dfe, ...
'mu_dc', mu_dc, ...
'comment', comment ...
);
% Check if exact Equalizer and BER entries already exist in the DB ...
selectedFields = {'Runs.run_id','BERs.ber_id','Equalizer.eq_id','BERs.ber',['BERs.occurrence' ...
'']};
filterParams = obj.tables;
filterParams.Equalizer = equalizerData;
[dataTable,sql_query] = obj.queryDB(filterParams, selectedFields);
% get or insert Equalizer
if ~isempty(dataTable)
% Equalizer entry already exists, use the existing eq_id
cur_eq_id = dataTable.eq_id;
else
% Insert the new Equalizer entry
cur_eq_id = obj.appendToTable('Equalizer', equalizerData);
end
% skip if already here or insert BER entry
if ~isempty(dataTable)
% A BER entry with the same eq_id, run_id, and occurrence already exists
existingBERValue = dataTable.ber;
% Compare the existing BER value with the new BER value
if existingBERValue == berValue
fprintf('The BER entry %.2e || -- eq_id: %d -- run_id: %d -- occurrence: %d already exists. \n',berValue, cur_eq_id, runID, occurrence);
else
fprintf('Already found BER for EQ: %.2e ~= %.2e || -- eq_id: %d -- run_id: %d -- occurrence: %d already exists.\n', berValue, existingBERValue, cur_eq_id, runID, occurrence);
end
else
% No such BER entry exists, insert the new BER entry
berData = struct( ...
'run_id', runID, ...
'eq_id', cur_eq_id, ...
'ber', berValue, ...
'occurrence', occurrence ...
);
obj.appendToTable('BERs', berData);
end
end
function executeSQL(obj, query)
% This method executes an SQL statement using MATLAB's execute function.
execute(obj.conn, query);
end
function answer = fetch(obj,query)
answer = fetch(obj.conn,query);
end
function [result,query] = queryDB(obj, filterParams, selectedFields)
% getPathsWithFlexibleFilter Retrieves values from Runs table with flexible filtering
% and lets the user select which fields to include in the SELECT statement.
%
% Usage:
% [rxRawPaths, filteredValues] = getPathsWithFlexibleFilter(filterParams)
%
% Inputs:
% filterParams: A structure containing the parameters with their values.
% If left empty, two popup windows will prompt the user for input.
%
% Outputs:
% result: table with sql return
% query: this was send to SQL DB
arguments
obj
filterParams = [];
selectedFields = [];
end
% Step 1: Prompt the user to input filter parameters if not provided
if isempty(filterParams)
filterParams = obj.promptFilterParameters();
end
% Step 2: Prompt the user to select fields to include in the SELECT statement
if isempty(selectedFields)
selectedFields = obj.promptSelectFields();
else
if iscell(selectedFields)
elseif isstruct(selectedFields)
end
end
% Step 3: Construct the SQL query based on the inputs
query = obj.constructSQLQuery(filterParams, selectedFields);
% Step 4: Execute the query and handle results
result = obj.fetch(query);
end
function query = constructSQLQuery(obj, filterParams, selectedFields)
% constructSQLQuery Constructs the SQL query based on filter parameters and selected fields.
% Construct the SELECT clause dynamically based on user selection
selectClause = 'SELECT DISTINCT ';
for i = 1:numel(selectedFields)
fieldParts = strsplit(selectedFields{i}, '.');
tableName = fieldParts{1};
fieldName = fieldParts{2};
if isnumeric(obj.tables.(tableName).(fieldName))
selectClause = [selectClause, 'COALESCE(', selectedFields{i}, ', ''NaN'') AS ', fieldName];
else
selectClause = [selectClause, 'COALESCE(', selectedFields{i}, ', '''') AS ', fieldName];
end
if i < numel(selectedFields)
selectClause = [selectClause, ', '];
else
selectClause = [selectClause, ' '];
end
end
% Construct the FROM and WHERE clause
baseQuery = [selectClause, 'FROM Runs ' ...
'LEFT JOIN Configurations ON Runs.run_id = Configurations.run_id ' ...
'LEFT JOIN Measurements ON Runs.run_id = Measurements.run_id ' ...
'WHERE '];
% 'LEFT JOIN BERs ON Runs.run_id = BERs.run_id ' ...
% 'LEFT JOIN Equalizer ON BERs.eq_id = Equalizer.eq_id ' ...
% Loop through each table in filterParams
filterClauses = [];
tableNames_ = fieldnames(filterParams);
for t = 1:numel(tableNames_)
tableName = tableNames_{t};
tableParams = filterParams.(tableName);
% Loop through each parameter in the table
fieldNames = fieldnames(tableParams);
for i = 1:numel(fieldNames)
fieldName = fieldNames{i};
value = tableParams.(fieldName);
% Construct the full column name in the format "tableName.fieldName"
fullName = sprintf('%s.%s', tableName, fieldName);
% Handle different types of values for SQL query construction
if isempty(value)
% Skip this parameter if it is empty (include all values)
continue;
elseif isnumeric(value) && isnan(value)
% If value is NaN, use IS NULL in SQL
filterClause = sprintf('%s IS NULL', fullName);
elseif isnumeric(value) && ~isEnumeration(value)
filterClause = sprintf('%s = %f', fullName, value);
elseif islogical(value) || (isnumeric(value) && ismember(value, [0, 1])) && ~isEnumeration(value)
filterClause = sprintf('%s = %d', fullName, value);
elseif ischar(value) || isstring(value)
filterClause = sprintf('%s = ''%s''', fullName, char(value)); %nicht nach string suchen sondern nach chararray -> 'bla' statt "bla"
elseif isEnumeration(value)
filterClause = sprintf('%s = ''%s''', fullName, value);
else
error('Unsupported data type for field "%s".', fullName);
end
% Add the constructed filter clause to the list
filterClauses = [filterClauses, filterClause, ' AND '];
end
end
% Remove trailing ' AND ' from the filter clauses if any filters were added
if ~isempty(filterClauses)
filterClauses = filterClauses(1:end-5);
end
% Construct the final SQL query
if isempty(filterClauses)
query = [selectClause, 'FROM Runs ' ...
'LEFT JOIN Configurations ON Runs.run_id = Configurations.run_id ' ...
'LEFT JOIN Measurements ON Runs.run_id = Measurements.run_id ' ...
'LEFT JOIN BERs ON Runs.run_id = BERs.run_id'];
else
query = [baseQuery, filterClauses];
end
end
function selectedFields = promptSelectFields(obj)
% promptSelectFields Prompts the user to select fields from multiple tables to include in the SELECT statement using settingsdlg.
% Get all possible fields from all tables (excluding sqlite_sequence)
tableNames = fieldnames(obj.tables);
tableNames = setdiff(tableNames, {'sqlite_sequence'}); % Remove sqlite_sequence
% Prepare the inputs for settingsdlg
promptSettings = {};
allFieldsFullName = {};
convertedFieldNames = {};
for i = 1:numel(tableNames)
tableFields = fieldnames(obj.tables.(tableNames{i}));
for j = 1:numel(tableFields)
fieldName = tableFields{j};
fullName = sprintf('%s.%s', tableNames{i}, fieldName);
convertedName = strrep(fullName, '.', '_'); % Replace '.' with '_'
allFieldsFullName{end + 1} = fullName; % Add full name to the list
convertedFieldNames{end + 1} = convertedName; % Store the converted name
% Add the field name and checkbox setting to the prompt
promptSettings{end + 1} = {sprintf('Include %s', fullName), convertedName};
promptSettings{end + 1} = false; % Default: not selected
end
end
% Create the settings dialog
[settings, button] = settingsdlg(...
'title', 'Select Fields for the SQL Query', ...
'description', 'Check the boxes for the fields you want to include in the SELECT statement.', ...
promptSettings{:} ...
);
% If the user cancels, default to selecting all fields
if strcmp(button, 'cancel')
selectedFields = allFieldsFullName;
return;
end
% Parse user input into selectedFields
selectedFields = {};
for i = 1:numel(allFieldsFullName)
convertedName = convertedFieldNames{i};
if isfield(settings, convertedName) && settings.(convertedName) % Add to selectedFields if the checkbox was selected
selectedFields{end + 1} = allFieldsFullName{i}; %#ok<AGROW>
end
end
% If no fields are selected, default to selecting all fields
if isempty(selectedFields)
selectedFields = allFieldsFullName;
end
end
function filterParams = promptFilterParameters(obj)
% promptFilterParameters Prompts the user to enter filter parameters using the settingsdlg framework.
% Get all possible parameters from all tables (excluding sqlite_sequence)
tableNames_ = fieldnames(obj.tables);
tableNames_ = setdiff(tableNames_, {'sqlite_sequence'}); % Remove sqlite_sequence
% Prepare the inputs for settingsdlg with sections and separators
promptSettings = {};
allFieldsFullName = {};
convertedFieldNames = {};
for i = 1:numel(tableNames_)
% Add a separator for each table section
promptSettings{end + 1} = 'separator';
promptSettings{end + 1} = tableNames_{i};
% Get all fields from the current table
tableFields = fieldnames(obj.tables.(tableNames_{i}));
% Prepare each field to be added to the dialog
for j = 1:numel(tableFields)
fieldName = tableFields{j};
fullName = sprintf('%s.%s', tableNames_{i}, fieldName);
convertedName = strrep(fullName, '.', '_'); % Replace '.' with '_'
% Skip fields that do not have distinct values stored
if ~isfield(obj.distinctValues.(tableNames_{i}), fieldName)
continue;
end
% Get the distinct values for the field
distinctValues_ = obj.distinctValues.(tableNames_{i}).(fieldName);
% Prepare distinct values for dropdown
if isempty(distinctValues_)
% If there are no distinct values, use only an "All" entry
distinctValues_ = {'All'};
else
% Ensure distinctValues is a cell array of strings
if isnumeric(distinctValues_)
distinctValues_ = arrayfun(@(x) num2str(x), distinctValues_, 'UniformOutput', false);
elseif isstring(distinctValues_)
distinctValues_ = cellstr(distinctValues_);
elseif iscell(distinctValues_) && ~iscellstr(distinctValues_)
distinctValues_ = cellfun(@num2str, distinctValues_, 'UniformOutput', false);
end
% Add an "All" option at the beginning of the distinct values list
distinctValues_ = [{'All'}; distinctValues_];
end
allFieldsFullName{end + 1} = fullName; % Add full name to the list
convertedFieldNames{end + 1} = convertedName; % Store the converted name
% Add the field name and value setting to the prompt
promptSettings{end + 1} = {sprintf('%s', fullName), convertedName};
promptSettings{end + 1} = distinctValues_; % Add distinct values as dropdown options
end
end
% Create the settings dialog
[settings, button] = settingsdlg(...
'title', 'Input Parameters for Filtering', ...
'description', 'Enter the values for each field to filter. Select "All" to include all values.', ...
promptSettings{:} ...
);
% If the user cancels, return an empty struct
if strcmp(button, 'cancel')
filterParams = struct();
return;
end
% Parse user input into filterParams structure
filterParams = struct();
for i = 1:numel(allFieldsFullName)
value = settings.(convertedFieldNames{i});
% Split full name to get table and field names
fieldParts = strsplit(allFieldsFullName{i}, '.');
tableName = fieldParts{1};
fieldName = fieldParts{2};
% If the table does not exist in the filterParams struct, create it
if ~isfield(filterParams, tableName)
filterParams.(tableName) = struct();
end
% Assign values to the respective fields under each table
if strcmp(value, 'All')
filterParams.(tableName).(fieldName) = []; % Set to empty to include all values
elseif isnumeric(value) && isnan(value)
filterParams.(tableName).(fieldName) = NaN; % Use NaN to handle as NULL
else
filterParams.(tableName).(fieldName) = value; % Use the entered value
end
end
end
end
end

76
Classes/DataBaseHandler/copyStylingFrom.m Normal file
View File

@@ -0,0 +1,76 @@
function copyStylingFrom(figNumSource, figNumTgt)
% Get handles to the source and target figures
sourceFig = figure(figNumSource);
targetFig = figure(figNumTgt);
% Get axes of source and target figures
sourceAxes = findall(sourceFig, 'type', 'axes');
targetAxes = findall(targetFig, 'type', 'axes');
% Ensure the number of axes match
if length(sourceAxes) ~= length(targetAxes)
error('Number of axes in source and target figures must be the same.');
end
% Loop through each pair of axes and copy styling properties
for i = 1:length(sourceAxes)
copyAxesProperties(sourceAxes(i), targetAxes(i));
end
% Apply general figure properties if desired
targetFig.Color = sourceFig.Color; % Background color
end
function copyAxesProperties(sourceAx, targetAx)
% List of properties to copy from source to target axes
propsToCopy = {'XColor', 'YColor', 'ZColor', 'FontSize', 'FontName', ...
'GridColor', 'GridLineStyle', 'MinorGridColor', 'Box', ...
'XGrid', 'YGrid', 'ZGrid', 'XMinorGrid', 'YMinorGrid', 'ZMinorGrid', ...
'LineWidth', 'TitleFontSizeMultiplier', 'LabelFontSizeMultiplier'};
% Copy properties from source to target
for i = 1:length(propsToCopy)
try
targetAx.(propsToCopy{i}) = sourceAx.(propsToCopy{i});
catch
% Skip property if it doesn't exist or can't be copied
end
end
% Copy axis labels and titles
targetAx.Title.String = sourceAx.Title.String;
targetAx.XLabel.String = sourceAx.XLabel.String;
targetAx.YLabel.String = sourceAx.YLabel.String;
targetAx.ZLabel.String = sourceAx.ZLabel.String;
% Copy children elements like lines, patches, etc.
sourceChildren = allchild(sourceAx);
targetChildren = allchild(targetAx);
% Ensure the number of children elements match
if length(sourceChildren) ~= length(targetChildren)
warning('Number of elements in source and target axes differ. Styling may not be applied completely.');
end
% Copy properties of children (like lines, patches, etc.), except colors and legends
for i = 1:min(length(sourceChildren), length(targetChildren))
copyObjectProperties(sourceChildren(i), targetChildren(i));
end
end
function copyObjectProperties(sourceObj, targetObj)
% List of common properties to copy for plot elements (lines, patches, etc.)
propsToCopy = {'LineStyle', 'LineWidth', 'Marker', 'MarkerSize', ...
'MarkerEdgeColor', 'MarkerFaceColor', 'DisplayName'};
% Copy properties from source to target, excluding colors
for i = 1:length(propsToCopy)
try
if ~contains(propsToCopy{i}, 'Color') % Skip color properties
targetObj.(propsToCopy{i}) = sourceObj.(propsToCopy{i});
end
catch
% Skip property if it doesn't exist or can't be copied
end
end
end

78
Classes/Warehouse_class/classes/DataStorage.m
View File

@@ -114,6 +114,12 @@ classdef DataStorage < handle
end end
function addValueToStorageByLinIdx(obj, valueToStore ,storageVarName, lin_idx)
obj.sto.(storageVarName){lin_idx} = valueToStore;
end
% Access Value(s) % Access Value(s)
function value = getStoValue(obj,storageVarName, varargin) function value = getStoValue(obj,storageVarName, varargin)
@@ -129,7 +135,9 @@ classdef DataStorage < handle
tmp = obj.sto.(storageVarName){lin_idx(i)}; tmp = obj.sto.(storageVarName){lin_idx(i)};
if ~isempty(tmp) if ~isempty(tmp)
if isa(tmp,'Signal') || isa(tmp,'struct') || isa(tmp,'Exfo_laser') || isa(tmp,'DC_supply') if isa(tmp,'double')
value(i) = tmp ;
elseif isa(tmp,'Signal') || isa(tmp,'struct') || isa(tmp,'Exfo_laser') || isa(tmp,'DC_supply')
if i == 1 if i == 1
value = {}; value = {};
end end
@@ -148,7 +156,10 @@ classdef DataStorage < handle
else else
try try
value(i,:) = tmp ; if i == 1
value = {};
end
value{i} = tmp ;
catch catch
% value(i,:) = tmp(1:size(value,2)) ; % value(i,:) = tmp(1:size(value,2)) ;
@@ -265,6 +276,69 @@ classdef DataStorage < handle
function [phys_indices,param_name] = getPhysIndicesByLinIndex(obj, lin_idx)
% Converts a linear index into the corresponding physical parameter values
% Inputs:
% - lin_idx: The linear index within the storage array
% Output:
% - phys_indices: A cell array containing the physical parameter values for each dimension
% Initialize output cell array
phys_indices = cell(1, numel(obj.fn));
% Convert linear index to subscript indices
[subscripts{1:numel(obj.dim)}] = ind2sub(obj.dim, lin_idx);
% Map subscripts to physical values for each parameter
for i = 1:numel(obj.fn)
param_name{i} = obj.fn(i);
phys_indices{i} = obj.parameter.(param_name{i}).getPhysForIndex(subscripts{i});
end
end
function [physStruct, stored_value] = getPhysAndValueByLinIndex(obj, storageVarName, lin_idx)
% Retrieves a structure with physical parameter values as fieldnames,
% their corresponding parameter names as values, and the stored value
% for a given linear index.
% Inputs:
% - storageVarName: Name of the storage variable in obj.sto
% - lin_idx: The linear index within the storage array
% Outputs:
% - physStruct: A structure with physical parameter values as fieldnames
% and parameter names as values
% - stored_value: The value stored at the given linear index in the
% specified storage variable
% Initialize an empty structure
physStruct = struct();
% Convert linear index to subscript indices
[subscripts{1:numel(obj.dim)}] = ind2sub(obj.dim, lin_idx);
% Map subscripts to physical values and parameter names for each dimension
for i = 1:numel(obj.fn)
param_name = obj.fn(i);
phys_value = obj.parameter.(param_name).getPhysForIndex(subscripts{i});
% Add to the structure with phys_value as the fieldname and param_name as the value
physStruct.(param_name) = phys_value;
end
% Retrieve the stored value at the given linear index
stored_value = obj.sto.(storageVarName){lin_idx};
end
function num_elements = getLastLinIndice(obj)
% Returns all possible linear indices for the data structure
% Output:
% - lin_indices: A column vector containing all linear indices for the storage array
% Calculate the total number of elements in the storage array
num_elements = prod(obj.dim);
end
end end
end end

114
Classes/Warehouse_class/classes/DataStorage2.m
View File

@@ -1,114 +0,0 @@
classdef DataStorage2 < handle
% DATASTORAGE: Stores data with physical parameter mappings
properties
inputParams = struct;
parameter = struct;
fn = [];
dim = [];
sto = struct;
end
methods
function obj = DataStorage2(inputParams)
% Constructor to initialize the DataStorage object
if nargin > 0
obj.inputParams = inputParams;
obj.fn = string(fieldnames(inputParams));
obj = obj.buildParameter();
obj.dim = obj.getDimension();
obj.sto = struct;
else
error('Input parameters are required.');
end
end
function showInfo(obj)
% Displays information about the storage and its dimensions
disp("Data Structure with fields:");
fprintf('%-12s | %-8s | %-12s\n', 'Name', 'Dimension', 'Physical Values');
disp('-------------------------------------------------------');
for i = 1:numel(obj.fn)
fprintf('%-12s | %-8d | %-12s\n', ...
char(obj.fn(i)), obj.dim(i), ...
strjoin(string(obj.parameter.(obj.fn(i)).values), ', '));
end
disp('-------------------------------------------------------');
end
function dim = getDimension(obj)
% Get the dimensions based on the length of parameters
dim = zeros(1, numel(obj.fn));
for p = 1:numel(obj.fn)
dim(p) = obj.parameter.(obj.fn(p)).length;
end
end
function obj = buildParameter(obj)
% Build the Parameter objects for each input parameter
for p = 1:numel(obj.fn)
name = obj.fn(p);
values = obj.inputParams.(name);
obj.parameter.(name) = Parameter2(name, values);
end
end
function addStorage(obj, varName)
% Create an empty storage for a specific variable name
obj.sto.(string(varName)) = cell(obj.dim);
end
function addValueToStorage(obj, valueToStore, storageVarName, varargin)
% Add a value to the storage at the specified indices
if nargin - 3 == numel(obj.fn)
lin_idx = obj.getIndicesByPhys(varargin);
obj.sto.(storageVarName){lin_idx} = valueToStore;
else
error('Please provide all indices for the storage.');
end
end
function value = getStoValue(obj, storageVarName, varargin)
% Retrieve a value from storage based on physical parameters
if nargin - 2 == numel(obj.fn)
lin_idx = obj.getIndicesByPhys(varargin);
value = cell(1, numel(lin_idx));
for i = 1:numel(lin_idx)
value{i} = obj.sto.(storageVarName){lin_idx(i)};
end
value = value(~cellfun('isempty', value)); % Remove empty entries
else
error('Please provide all physical parameters.');
end
end
function lin_idx = getIndicesByPhys(obj, varargin)
% Unpack nested cell array if needed
if numel(varargin) == 1 && iscell(varargin{1})
varargin = varargin{1}; % Unpack if single cell array is passed
end
indices = cell(1, numel(obj.fn));
% Loop through each parameter (e.g., L, D)
for p = 1:numel(obj.fn)
% Unwrap if it's a cell
if iscell(varargin{p})
physVal = varargin{p}{1}; % Extract scalar from cell
else
physVal = varargin{p}; % It's already a scalar
end
paramName = obj.fn(p); % Get the parameter name (e.g., 'L' or 'D')
% Call getIndexByPhys on the corresponding Parameter2 object
indices{p} = obj.parameter.(paramName).getIndexByPhys(physVal);
end
% Convert subscript indices to a linear index
lin_idx = sub2ind(obj.dim, indices{:});
end
end
end

51
Classes/Warehouse_class/classes/Parameter2.m
View File

@@ -1,51 +0,0 @@
classdef Parameter2 < handle
% PARAMETER2: Represents a physical parameter with mappings between values and indices
properties
name
values
length
physToIndexMap % Rename this from 'getPhysForIndex'
indexToPhysMap % Rename this from 'getIndexForPhys'
end
methods
function obj = Parameter2(name, values)
% Constructor to initialize the Parameter2 object
obj.name = name;
obj.values = values;
obj.length = numel(values);
% Initialize the mappings
obj.physToIndexMap = containers.Map('KeyType', 'double', 'ValueType', 'any');
obj.indexToPhysMap = containers.Map('KeyType', 'double', 'ValueType', 'any');
obj = obj.buildMappings();
end
function obj = buildMappings(obj)
% Build mappings between physical values and indices
for idx = 1:obj.length
obj.indexToPhysMap(idx) = obj.values(idx);
obj.physToIndexMap(obj.values(idx)) = idx;
end
end
function physVal = getPhysForIndex(obj, idx)
% Return the physical value corresponding to the index
if isKey(obj.indexToPhysMap, idx)
physVal = obj.indexToPhysMap(idx);
else
error('Index out of range for parameter %s', obj.name);
end
end
function idx = getIndexByPhys(obj, physVal)
% Return the index corresponding to the physical value
if isKey(obj.physToIndexMap, physVal)
idx = obj.physToIndexMap(physVal);
else
error('Physical value %g not found in parameter %s', physVal, obj.name);
end
end
end
end

45
Classes/Warehouse_class/classes/minimalExample_gen2.m
View File

@@ -1,45 +0,0 @@
% Define input parameters for the DataStorage2
inputParams.L = [1, 2, 10, 80]; % Length in kilometers
inputParams.D = [16, 17, 18]; % Diameter in millimeters
% Create a DataStorage2 instance with the input parameters
dataStorage = DataStorage2(inputParams); % Using DataStorage2 class
% Display the current information about the data storage structure
dataStorage.showInfo();
% Add a storage variable named 'testStorage'
dataStorage.addStorage('testStorage');
% Add a value (e.g., 100) to the storage at specific physical parameter values
% For example, we store the value 100 at L = 10 km and D = 17 mm
dataStorage.addValueToStorage(100, 'testStorage', 10, 16);
dataStorage.addValueToStorage(100, 'testStorage', 10, 17);
dataStorage.addValueToStorage(100, 'testStorage', 10, 18);
% Retrieve the value from the storage at the same physical parameter values
storedValue = dataStorage.getStoValue('testStorage', 10, 16:18);
disp('Retrieved value from storage:');
disp(storedValue);
% Retrieve another value at a non-existent location (L = 2 km, D = 8 mm)
% This will show how the function handles empty storage entries
nonExistentValue = dataStorage.getStoValue('testStorage', 2, 8);
disp('Retrieved value from empty location:');
disp(nonExistentValue);
% Use the internal mappings to check how physical values map to indices
% Get the linear index for physical values L = 10 km and D = 17 mm
lin_idx = dataStorage.getIndicesByPhys(10, 17);
disp('Linear index for L=10 km and D=17 mm:');
disp(lin_idx);
% Check the reverse mapping: physical value for index 2 of parameter L
physValForIndex = dataStorage.parameter.L.getPhysForIndex(2);
disp('Physical value for index 2 of parameter L:');
disp(physValForIndex);
% Check the mapping: index for physical value D = 21 mm
indexForPhys = dataStorage.parameter.D.getIndexByPhys(21);
disp('Index for physical value D=21 mm:');
disp(indexForPhys);

BIN
Classes/matlab.mat
View File

Binary file not shown.

11
Datatypes/db_mode.m Normal file
View File

@@ -0,0 +1,11 @@
classdef db_mode < int32
enumeration
no_db (0)
db_precoded (1) %sequence is precoded
db_encoded (2) %sequence is precoded and encoded
db_emulate (3) %for eq'ing: emulate precode on sequence that was not precoded at tx
db_discard (4) %for eq'ing: discard precode on sequence that was precoded at tx
end
end

11
Datatypes/equalizer_structure.m Normal file
View File

@@ -0,0 +1,11 @@
classdef equalizer_structure < int32
enumeration
ffe (0)
vnle (1)
vnle_pf_mlse (2)
db_precoded (3)
db_encoded (4)
end
end

11
Datatypes/isEnumeration.m Normal file
View File

@@ -0,0 +1,11 @@
function result = isEnumeration(value)
% isEnumeration Checks if a given value is an instance of an enumeration class
% Usage:
% result = isEnumeration(value);
% Get meta information about the class of the value
metaInfo = metaclass(value);
% Check if the meta information indicates an enumeration class
result = metaInfo.Enumeration;
end

18
Functions/EQ_structures/duobinary_signaling.m Normal file
View File

@@ -0,0 +1,18 @@
function [eq_package] = duobinary_signaling(eq_, mlse_,M ,rx_signal, tx_symbols, tx_bits)
%Duobinary Signaling
[eq_signal, eq_noise] = eq_.process(rx_signal,tx_symbols);
eq_signal = mlse_.process(eq_signal);
eq_signal = Duobinary().encode(eq_signal);
eq_signal = Duobinary().decode(eq_signal);
% M = numel(unique(eq_signal.signal));
rx_bits = PAMmapper(M,0).demap(eq_signal);
[~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
eq_package.ber = ber;
end

77
Functions/EQ_structures/duobinary_target.m Normal file
View File

@@ -0,0 +1,77 @@
function [eq_package] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx_bits, options)
arguments
eq_
mlse_
M
rx_signal
tx_symbols
tx_bits
options.precode_mode db_mode
options.showAnalysis = 0;
end
%Duobinary Targeting
[eq_signal, eq_noise] = eq_.process(rx_signal,Duobinary().encode(tx_symbols));
% dir = [1,1];
mlse_sig_sd = mlse_.process(eq_signal);
mlse_sig_hd = PAMmapper(M,0).quantize(mlse_sig_sd);
% precoding to mitigate error propagation, most prominently used in
% combination with duobinary signaling to avoid catastrophic error
% behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling)
% takes:
% -> eq_signal_hd: hard decision signal after eq
% -> tx_symbols: that where used as reference for eq
switch options.precode_mode
case db_mode.db_emulate
mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M);
mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M);
tx_symbols_precoded = Duobinary().encode(tx_symbols);
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
tx_bits = PAMmapper(M,0).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);
case db_mode.db_encoded
% normal DB encoded data (only for 10KM)
case db_mode.db_precoded
mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M);
mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M);
end
% M = numel(unique(tx_symbols.signal));
rx_bits = PAMmapper(M,0).demap(mlse_sig_hd);
[~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
eq_package.ber = ber;
if options.showAnalysis
eq_noise = eq_noise - mean(eq_noise.signal);
rx_signal.spectrum("normalizeTo0dB",1,"fignum",250);
showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise');
Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",250);
end
end

86
Functions/EQ_structures/vnle.m Normal file
View File

@@ -0,0 +1,86 @@
function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options)
%VNLE Apply an equalization algorithm to the received signal and calculate BER
% This function takes an equalizer object, a received signal, and the
% transmitted symbols to apply equalization, map the received signal back to bits,
% and compute the bit error rate (BER).
%
% Inputs:
% EQ - Equalizer object that provides the equalization method
% rx_signal - Received signal that needs to be equalized
% tx_symbols - Transmitted symbols used as a reference for BER calculation
%
% Outputs:
% eq_signal - Equalized version of the received signal
% ber - Bit error rate after equalization
% numErrors - Number of bit errors detected
arguments
eq_
M
rx_signal
tx_symbols
tx_bits
options.precode_mode db_mode
options.showAnalysis = 0
end
%FFE or VNLE
[eq_signal_sd,eq_noise] = eq_.process(rx_signal,tx_symbols);
eq_signal_hd = PAMmapper(M,0).quantize(eq_signal_sd);
% precoding to mitigate error propagation, most prominently used in
% combination with duobinary signaling to avoid catastrophic error
% behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling)
% takes:
% -> eq_signal_hd: hard decision signal after eq
% -> tx_symbols: that where used as reference for eq
switch options.precode_mode
case db_mode.db_emulate
% re
eq_signal_hd = Duobinary().encode(eq_signal_hd,"M",M);
eq_signal_hd = Duobinary().decode(eq_signal_hd,"M",M);
tx_symbols_precoded = Duobinary().encode(tx_symbols);
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
tx_bits = PAMmapper(M,0).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);
case db_mode.db_encoded
% normal DB encoded data (only for 10KM)
case db_mode.db_precoded
eq_signal_hd = Duobinary().encode(eq_signal_hd,"M",M);
eq_signal_hd = Duobinary().decode(eq_signal_hd,"M",M);
end
rx_bits = PAMmapper(M,0).demap(eq_signal_hd);
[~,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);
eq_package.ber_vnle = ber;
eq_package.evm_total = evm_total;
eq_package.evm_lvl = evm_lvl;
eq_package.inf_rate_vnle = inf_rate;
if options.showAnalysis
fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl);
fprintf('VNLE BER: %.2e \n',ber);
end
end

125
Functions/EQ_structures/vnle_postfilter_mlse.m Normal file
View File

@@ -0,0 +1,125 @@
function [eq_package] = vnle_postfilter_mlse(eq_,pf_,mlse_,M,rx_signal,tx_symbols,tx_bits,options)
arguments
eq_
pf_
mlse_
M
rx_signal
tx_symbols
tx_bits
options.precode_mode db_mode
options.showAnalysis = 0;
end
%FFE or VNLE
[eq_signal_sd,eq_noise] = eq_.process(rx_signal,tx_symbols);
eq_signal_hd = PAMmapper(M,0).quantize(eq_signal_sd);
mlse_sig_sd = pf_.process(eq_signal_sd,eq_noise);
mlse_.DIR = pf_.coefficients;
% [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);
% precoding to mitigate error propagation, most prominently used in
% combination with duobinary signaling to avoid catastrophic error
% behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling)
% takes:
% -> M
% -> eq_signal_hd: hard decision signal after eq
% -> tx_symbols: that where used as reference for eq
switch options.precode_mode
case db_mode.db_emulate
% re
eq_signal_hd = Duobinary().encode(eq_signal_hd,"M",M);
eq_signal_hd = Duobinary().decode(eq_signal_hd,"M",M);
mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M);
mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M);
tx_symbols_precoded = Duobinary().encode(tx_symbols);
tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
tx_bits = PAMmapper(M,0).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);
case db_mode.db_encoded
% normal DB encoded data (only for 10KM)
case db_mode.db_precoded
eq_signal_hd = Duobinary().encode(eq_signal_hd,"M",M);
eq_signal_hd = Duobinary().decode(eq_signal_hd,"M",M);
mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M);
mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M);
end
% METRICS OF VNLE %
rx_bits_vnle = PAMmapper(M,0).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
[inf_rate_vnle] = calc_air(eq_signal_sd,tx_symbols,"skip_front",10000,"skip_end",10000);
[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);
[~,~,ber_mlse,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
eq_package.ber_mlse = ber_mlse;
eq_package.ber_vnle = ber_vnle;
eq_package.evm_vnle_total = evm_vnle_total;
eq_package.evm_vnle_lvl = evm_vnle_lvl;
eq_package.air = inf_rate_vnle;
eq_package.eq = eq_;
eq_package.pf = pf_;
eq_package.mlse = mlse_;
if options.showAnalysis
% fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl);
fprintf('VNLE BER: %.2e \n',ber_vnle);
fprintf('MLSE BER: %.2e \n',ber_mlse);
showEQNoisePSD(eq_noise,"fignum",336,"displayname",'VNLE+DFE','postfilter_taps',pf_.coefficients);
rx_signal.spectrum("normalizeTo0dB",1,"fignum",337,"displayname",'Rx Signal');
tx_symbols.spectrum("normalizeTo0dB",1,"fignum",337,'displayname','Tx Signal');
showLevelHistogram(eq_signal_sd,tx_symbols)
% showLevelHistogram(mlse_sig_sd,tx_symbols)
showEQcoefficients(eq_.e,eq_.e2,eq_.e3,"displayname",'Coefficients');
showEQNoiseSNR(tx_symbols,eq_noise,"displayname",'vnle snr','fignum',101);
%%% EQ SNR Spectrum %230
%snr
snr_vnle = snr(tx_symbols.signal,eq_noise.signal);
% showErrorBurstCount(eq_signal_sd,tx_symbols)
end
end

150
Functions/EQ_visuals/analyzeEQperformance.m Normal file
View File

@@ -0,0 +1,150 @@
function analyzeEQperformance(ref_bits,ref_symbols,rx_signal,eq_signal,eq_decisions,fsym,M,options)
arguments
ref_bits
ref_symbols
rx_signal
eq_signal
eq_decisions
fsym
M
options.postfilterclass
options.eqclass
options.mlseclass
options.db_precoded
options.displayname
end
% toolkit to visualize stuff related to DSP of IM/DD
%%% Preps
if isempty(eq_decisions)
eq_decisions = PAMmapper(M,0).quantize(eq_signal);
end
%%% Demap eqlzd signal to determine BER
if options.db_precoded
eq_hd = PAMmapper(M,0).quantize(eq_signal);
eq_hd = Duobinary().encode(eq_hd);
eq_hd = Duobinary().decode(eq_hd,"M",M);
rx_bits = PAMmapper(M,0).demap(eq_hd);
else
rx_bits = PAMmapper(M,0).demap(eq_signal);
end
[~,numerr,ber_sd,errpos] = calc_ber(rx_bits.signal,ref_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% fprintf('SD BER: %.2e \n',ber_sd);
%%% Demap provided decisions to determine BER
rx_bits = PAMmapper(M,0).demap(eq_decisions);
[~,numerr,ber_mlse,errpos] = calc_ber(rx_bits.signal,ref_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% fprintf('MLSE BER: %.2e \n',ber_mlse);
%%% Noise prior to DSP (is thius accurate with resampling?)
rx_resampled = rx_signal.normalize("mode","rms").resample("fs_out",fsym);
rx_noise = rx_resampled-ref_symbols;
%%% Noise after to soft-decision DSP
eq_noise = eq_signal-ref_symbols;
col = cbrewer2('paired',12);
lines = numel(findall(figure(200), 'Type', 'Line'))+1;
darkcoloridx = max(mod(2*lines,12),2);
lightcoloridx = max(mod(2*lines-1,12),1);
%%% Separate Classes
constellation = unique(ref_symbols.signal);
received_sd = NaN(numel(constellation),length(ref_symbols));
received_hd = NaN(numel(constellation),length(ref_symbols));
lvlcol = cbrewer2('Set1',numel(constellation));
for lvl = 1:numel(constellation)
%Separate the equalized signal into the
%respective levels based on the actually
%transmitted level!
received_sd(lvl,ref_symbols.signal==constellation(lvl)) = eq_signal.signal(ref_symbols.signal==constellation(lvl));
received_hd(lvl,ref_symbols.signal==constellation(lvl)) = eq_decisions.signal(ref_symbols.signal==constellation(lvl));
end
%%% bursts
% Find differences between consecutive elements
diff_indices = diff(errpos);
% Identify the start of new sequences (when the difference is not 1)
sequence_starts = [1, find(diff_indices ~= 1) + 1]; % Include the first index
sequence_ends = [sequence_starts(2:end) - 1, length(errpos)]; % Calculate end indices
% Initialize burst count and print bursts longer than 10
burst_len = 1:10;
burst_count = zeros(length(burst_len),1);
for t = 1:numel(burst_len)
for i = 1:length(sequence_starts)
% Extract current sequence
current_burst = errpos(sequence_starts(i):sequence_ends(i));
% Check if the sequence length matches criterion
if length(current_burst) == burst_len(t)
burst_count(t) = burst_count(t) + 1;
end
end
end
burst_symbols = burst_count .* burst_len';
burst_rate = burst_symbols ;%./ length(rx_signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Rx Spectrum %200
rx_signal.spectrum("displayname",sprintf('Rx %d GBd PAM%d',fsym.*1e-9,M),'fignum',200,'normalizeTo0dB',1,'color',col(darkcoloridx,:));
xline([-fsym/2,fsym/2].*1e-9,'Color',col(mod(lines,12)+2,:),'HandleVisibility','off');
ylim([-20,3]);
%%% EQ Spectrum
%
%%% EQ Time Series %210
showEQTimeSignal(eq_signal,ref_symbols)
%%% EQ Noise Spectrum + inverted Postfilter %220
showEQNoisePSD(eq_noise,options.postfilterclass.burg_coeff,"fignum",220,"color",col(darkcoloridx,:));
%%% EQ SNR Spectrum %230
fft_length = 2^13;
[s_lin,w] = pwelch(eq_signal.signal,hanning(fft_length),fft_length/2,fft_length,eq_signal.fs,"centered","psd","mean");
[n_lin,w] = pwelch(eq_noise.signal,hanning(fft_length),fft_length/2,fft_length,eq_noise.fs,"centered","psd","mean");
w = w.*1e-9;
snr_dbm = 10*log10(s_lin./n_lin);
figure(230)
hold on
plot(w,snr_dbm,'DisplayName','SNR after EQ','LineWidth',1,'Color',col(darkcoloridx,:));
xlabel("Frequency in GHz");
edgetick = 2^(nextpow2(eq_signal.fs*1e-9));
xticks(-edgetick:16:edgetick);
xlim([-128 128]);
ylim([-20,35]);
grid minor
yticks(-200:10:100);
grid on; grid minor;
legend('Interpreter','none');
title('Noise of soft decision signal (not MLSE)');
%%% FFE histogram %240
showLevelHistogram(eq_signal,ref_symbols)
%%% Confusion Matrix
showLevelConfusionMatrix(eq_decisions,ref_symbols,"fignum",250)
%%% Burst Count
figure(260)
hold on
plot(burst_len,burst_rate,'Marker','x','Color',col(darkcoloridx,:),"displayname",sprintf('Rx %d GBd PAM%d; %s',fsym.*1e-9,M,options.displayname));
xlabel('length of error burst');
ylabel('occurences')
grid on
set(gca, 'YScale', 'log');
autoArrangeFigures
end

44
Functions/EQ_visuals/showEQNoisePSD.m Normal file
View File

@@ -0,0 +1,44 @@
function showEQNoisePSD(eq_noise, options)
arguments
eq_noise
options.postfilter_taps = NaN
options.fignum (1,1) double = NaN % Default to NaN if not provided
options.displayname (1,:) char = '' % Default to an empty string if not provided
options.color = [0.2157 0.4941 0.7216];
end
% Determine the figure number to use or create a new figure
if isnan(options.fignum)
fig = figure; % Create a new figure and get its handle
else
fig = figure(options.fignum); % Use the specified figure number
end
hold on
ax = gca;
% N = numel(ax.Children);
N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children));
cmap = linspecer(8);
options.color = cmap(mod(N, size(cmap, 1)) + 1, :);
% Ensure the figure is ready before calling spectrum
eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 1,"color",options.color);
title('Noise of soft decision signal (not MLSE)')
if ~isnan(options.postfilter_taps)
% Hold on to the figure for further plotting
hold on;
% Compute the frequency response of the postfilter
[h, w] = freqz(1, options.postfilter_taps, length(eq_noise), "whole", eq_noise.fs);
h = h / max(abs(h)); % Normalize the filter response
% Adjust frequency axis to center at 0
w_ = (w - eq_noise.fs / 2);
% Plot the inverted postfilter response
plot(w_ * 1e-9, 20 * log10(fftshift(abs(h))), 'DisplayName', ['Burg Coeffs: ', num2str(round(options.postfilter_taps, 2)), ' '], 'LineWidth', 1,'Color',options.color,'LineStyle','--');
% Ensure a legend is displayed
legend('show');
end
end

70
Functions/EQ_visuals/showEQNoiseSNR.m Normal file
View File

@@ -0,0 +1,70 @@
function showEQNoiseSNR(tx_signal, rx_signal, options)
arguments
tx_signal
rx_signal
options.fs_tx
options.fs_rx
options.fignum (1,1) double = NaN % Default to NaN if not provided
options.displayname (1,:) char = '' % Default to an empty string if not provided
options.color = [0.2157 0.4941 0.7216];
end
% Determine the figure number to use or create a new figure
if isnan(options.fignum)
fig = figure; % Create a new figure and get its handle
else
fig = figure(options.fignum); % Use the specified figure number
end
if isa(tx_signal,'Signal')
options.fs_tx = tx_signal.fs;
tx_signal = tx_signal.signal;
end
if isa(rx_signal,'Signal')
options.fs_rx = rx_signal.fs;
rx_signal = rx_signal.signal;
end
hold on
ax = gca;
% N = numel(ax.Children);
N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children));
cmap = linspecer(8);
options.color = cmap(mod(N, size(cmap, 1)) + 1, :);
% Ensure the figure is ready before calling spectrum
title('SNR of received Signal')
fft_length = 2^(nextpow2(length(tx_signal))-7);
[s_lin,w] = pwelch(tx_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_tx,"centered","psd","mean");
[n_lin,w] = pwelch(rx_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_rx,"centered","psd","mean");
w = w.*1e-9;
snr_dbm = 10*log10(s_lin./n_lin);
% figure(231)
hold on
plot(w,snr_dbm,'DisplayName','SNR','LineWidth',0.5,'Color',options.color);
xlabel("Frequency in GHz");
edgetick = 2^(nextpow2(options.fs_tx*1e-9));
ticks = -edgetick:16:edgetick;
xticks(ticks);
[~,b]=min(abs((-edgetick:16:edgetick)-max(w)));
xlim([-ticks(b+1) ticks(b+1)]);
max_snr = ceil(max(snr_dbm)/10)*10;
min_snr = floor(min(snr_dbm)/10)*10;
ylim([min_snr,max_snr]);
yticks(-200:10:100);
grid on; grid minor;
legend('Interpreter','none');
title('Noise of soft decision signal (not MLSE)');
end

54
Functions/EQ_visuals/showEQTimeSignal.m Normal file
View File

@@ -0,0 +1,54 @@
function showEQTimeSignal(eq_signal,ref_symbols,options)
arguments
eq_signal
ref_symbols
options.fignum (1,1) double = NaN % Default to NaN if not provided
options.displayname (1,:) char = '' % Default to an empty string if not provided
options.color = [0.2157 0.4941 0.7216];
end
% Determine the figure number to use or create a new figure
if isnan(options.fignum)
fig = figure; % Create a new figure and get its handle
else
fig = figure(options.fignum); % Use the specified figure number
end
M = numel(unique(ref_symbols.signal));
lvlcol = cbrewer2('Set1',M);
col = cbrewer2('paired',2);
eq_decisions = PAMmapper(M,0).quantize(eq_signal);
rx_bits = PAMmapper(M,0).demap(eq_signal);
ref_bits = PAMmapper(M,0).demap(ref_symbols);
%%% Separate Classes
constellation = unique(ref_symbols.signal);
received_sd = NaN(numel(constellation),length(ref_symbols));
received_hd = NaN(numel(constellation),length(ref_symbols));
lvlcol = cbrewer2('Set1',numel(constellation));
for lvl = 1:numel(constellation)
%Separate the equalized signal into the
%respective levels based on the actually
%transmitted level!
received_sd(lvl,ref_symbols.signal==constellation(lvl)) = eq_signal.signal(ref_symbols.signal==constellation(lvl));
received_hd(lvl,ref_symbols.signal==constellation(lvl)) = eq_decisions.signal(ref_symbols.signal==constellation(lvl));
end
[~,numerr,ber_sd,errpos] = calc_ber(rx_bits.signal,ref_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
%%% EQ Time Series %210
eq_signal.plot("fignum",fig.Number,"displayname",'Equalized Signal','color',col(1,:),'clear',1);
hold on
try
yline(PAMmapper(M,0).get_demodulation_thresholds,'HandleVisibility','off','LineStyle','--');
for c = 1:M
scatter(errpos./eq_signal.fs,received_sd(c,errpos),1,'x','MarkerEdgeColor',lvlcol(c,:),'LineWidth',2,'DisplayName',sprintf('Tx Lvl: %d',c));
end
end
end

58
Functions/EQ_visuals/showEQcoefficients.m Normal file
View File

@@ -0,0 +1,58 @@
function showEQcoefficients(n1, n2, n3, options)
% Show filter coefficients as stem plot
% n1, n2, and n3 in different subplots
% Scale all y-axis to -1 and 1
arguments
n1
n2
n3
options.fignum (1,1) double = NaN % Default to NaN if not provided
options.displayname (1,:) char = '' % Default to an empty string if not provided
options.color = [0.2157, 0.4941, 0.7216];
options.clf = 0; % Clear figure before plotting new
end
% Determine the figure number to use or create a new figure
if isnan(options.fignum)
fig = figure; % Create a new figure and get its handle
else
fig = figure(options.fignum); % Use the specified figure number
end
if options.clf
clf(fig); % Clear the figure if requested
end
hold on
ax = gca;
N = numel(ax.Children);
% Set up a colormap for consistent coloring
cmap = linspecer(8);
options.color = cmap(mod(N, size(cmap, 1)) + 1, :);
% Create subplots for n1, n2, n3
for i = 1:3
subplot(3, 1, i);
switch i
case 1
stem(n1, 'Color', options.color, 'LineWidth', 1,'Marker','.','MarkerSize',10);
title(sprintf('1st order Filter Coefficients: %d',numel(n1)));
case 2
stem(n2, 'Color', options.color, 'LineWidth', 1,'Marker','.','MarkerSize',10);
title(sprintf('2nd order Filter Coefficients: %d',numel(n2)));
case 3
stem(n3, 'Color', options.color, 'LineWidth', 1,'Marker','.','MarkerSize',10);
title(sprintf('3rd order Filter Coefficients: %d',numel(n3)));
end
ylim([-1, 1]); % Scale y-axis to -1 and 1
grid on;
grid minor
xlabel('Coefficient Index');
ylabel('Amplitude');
end
% Ensure the layout is tight for better visibility
sgtitle('Filter Coefficients'); % Overall title
end

47
Functions/EQ_visuals/showErrorBurstCount.m Normal file
View File

@@ -0,0 +1,47 @@
function showErrorBurstCount(eq_signal,ref_symbols,options)
arguments
eq_signal
ref_symbols
options.fignum (1,1) double = NaN % Default to NaN if not provided
options.displayname (1,:) char = '' % Default to an empty string if not provided
end
% Determine the figure number to use or create a new figure
if isnan(options.fignum)
fig = figure; % Create a new figure and get its handle
else
fig = figure(options.fignum); % Use the specified figure number
end
if numel(unique(eq_signal.signal)) > 20
M = numel(unique(ref_symbols.signal));
eq_signal = PAMmapper(M,0).quantize(eq_signal);
end
diff_indices = eq_signal.signal == ref_symbols.signal;
% Identify the start of new sequences (when the difference is not 1)
sequence_starts = [1, find(diff_indices ~= 1) + 1]; % Include the first index
sequence_ends = [sequence_starts(2:end) - 1, length(errpos)]; % Calculate end indices
% Initialize burst count and print bursts longer than 10
burst_len = 1:10;
burst_count = zeros(length(burst_len),1);
for t = 1:numel(burst_len)
for i = 1:length(sequence_starts)
% Extract current sequence
current_burst = errpos(sequence_starts(i):sequence_ends(i));
% Check if the sequence length matches criterion
if length(current_burst) == burst_len(t)
burst_count(t) = burst_count(t) + 1;
end
end
end
burst_symbols = burst_count .* burst_len';
burst_rate = burst_symbols ;%./ length(rx_signal);
end

27
Functions/EQ_visuals/showLevelConfusionMatrix.m Normal file
View File

@@ -0,0 +1,27 @@
function showLevelConfusionMatrix(decided_symbols,ref_symbols,options)
arguments
decided_symbols
ref_symbols
options.M = NaN
options.fignum (1,1) double = NaN % Default to NaN if not provided
options.displayname (1,:) char = '' % Default to an empty string if not provided
end
% Determine the figure number to use or create a new figure
if isnan(options.fignum)
fig = figure; % Create a new figure and get its handle
else
fig = figure(options.fignum); % Use the specified figure number
end
if length(unique(decided_symbols.signal))>20
assert(~isnan(options.M),'Provide either decided symbol sequence or modulation order M (PAM-4 -> M=4)');
decided_symbols = PAMmapper(options.M,0).quantize(decided_symbols);
end
%%% Confusion Matrix
cm = confusionchart(ref_symbols.signal,decided_symbols.signal,'RowSummary','row-normalized','ColumnSummary','column-normalized','Title','Confusion Matrix','XLabel','Decisions','YLabel','Transmitted');
end

50
Functions/EQ_visuals/showLevelHistogram.m Normal file
View File

@@ -0,0 +1,50 @@
function showLevelHistogram(eq_signal,ref_symbols,options)
arguments
eq_signal
ref_symbols
options.fignum (1,1) double = NaN % Default to NaN if not provided
options.displayname (1,:) char = '' % Default to an empty string if not provided
end
if isa(eq_signal,'Signal')
eq_signal = eq_signal.signal;
end
if isa(ref_symbols,'Signal')
ref_symbols = ref_symbols.signal;
end
% Determine the figure number to use or create a new figure
if isnan(options.fignum)
fig = figure; % Create a new figure and get its handle
else
fig = figure(options.fignum); % Use the specified figure number
end
%%% Separate Classes
constellation = unique(ref_symbols);
received_sd = NaN(numel(constellation),length(ref_symbols));
lvlcol = cbrewer2('Set1',numel(constellation));
for lvl = 1:numel(constellation)
%Separate the equalized signal into the
%respective levels based on the actually
%transmitted level!
received_sd(lvl,ref_symbols==constellation(lvl)) = eq_signal(ref_symbols==constellation(lvl));
end
%%% FFE histogram
clf
for lvl = 1:numel(constellation)
intermediate = received_sd(lvl,:);
cnt(lvl) = round(numel(intermediate(~isnan(intermediate)))./length(eq_signal),3).*100;
hold on
histogram(received_sd(lvl,:),1000,"EdgeAlpha",0,'DisplayName',['Lvl ',num2str(lvl),' | ',num2str(cnt(lvl)),' %'],'FaceColor',lvlcol(lvl,:),'Normalization','pdf');
end
legend
grid on
end

0
Functions/holdAndShowValue.m → Functions/Lab_helper/holdAndShowValue.m
View File

5
Functions/showCurrentMeasurement.m → Functions/Lab_helper/showCurrentMeasurement.m
View File

@@ -1,5 +1,8 @@
function showCurrentMeasurement(varargin) function showCurrentMeasurement(varargin)
% showCurrentMeasurement displays measurement data in a figure with variable names
% USEFUL FOR LAB!
% showCurrentMeasurement displays lab measurement data in a table figure with variable names
% as column headers and values listed below. Calling the function multiple times % as column headers and values listed below. Calling the function multiple times
% with the same variable names but different values adds more data points. % with the same variable names but different values adds more data points.
% %

0
Functions/waitUntilClick.m → Functions/Lab_helper/waitUntilClick.m
View File

156
Functions/Metrics/air_garcia_implementation.m Normal file
View File

@@ -0,0 +1,156 @@
function air = air_garcia_implementation(x,r,idx_tx,Px,M_training)
MAX_MEMORY = 200e6; % maximum allowed size for a matrix
if nargin == 3
Px = [];
M_training = [];
end
if nargin == 4
M_training = [];
end
% if input is complex, separate into real and imaginary parts
if any(imag(x(:))~=0) || any(imag(r(:))~=0)
x = [real(x); imag(x)];
r = [real(r); imag(r)];
end
D = size(x, 1); % D = 2 if complex x
N = size(x, 2); % number of constellation points
M = size(r, 2); % number of samples
% set default training set size
if isempty(M_training)
M_training = ceil(0.3*M);
end
M_testing = M - M_training;
% Training: estimate parameters of the conditionally Gaussian model
% sort according to transmit index
[idx_tx_training, idx_sort] = sort(idx_tx(1:M_training));
r_training = r(:, idx_sort);
i_bounds = zeros(1, N+1);
% compute conditional means and covariance matrices
C_n = zeros(D, D, N);
det_n = zeros(1, N);
for n=1:N
% find how many times x(:, n) was transmitted and update i_bounds
N_current_x = find(idx_tx_training((i_bounds(n)+1):end)==n, 1, 'last');
if isempty(N_current_x), N_current_x=0; end
i_bounds(n+1) = i_bounds(n) + N_current_x;
if N_current_x > 0
% Compute mu_n=E[Y|X=x_n] according to Eq. (14) and store it in
% x(:, n) to save space
x(:, n) = sum(r_training(:, (i_bounds(n)+1):i_bounds(n+1)), 2)/(i_bounds(n+1)-i_bounds(n));
% compute C_n=cov[Y|X=x_n] according to Eq. (15)
r_meanfree = r_training(:, (i_bounds(n)+1):i_bounds(n+1)) - x(:, n);
C_n(:, :, n) = (r_meanfree*r_meanfree')/(i_bounds(n+1)-i_bounds(n));
% store also the determinant of C(:, :, n)
det_n(n) = det(C_n(:, :, n));
% if the determinant is 0, or if the matrix is badly conditioned,
% regularize by adding a small identity matrix. Note that we do
% need the check for 0 determinant, in case a cloud has exactly 0
% variance according to the training set
if det_n(n)==0 || cond(C_n(:, :, n))>1e16
C_n(:, :, n) = C_n(:, :, n) + 5 * eps * eye(D);
det_n(n) = (5*eps)^D;
end
end
end
% uniform input pmf Px if not provided
if isempty(Px)
Px = repmat(1/N, [1, N]);
end
% extract testing set and sort it according to transmit index
[idx_tx_testing, idx_sort] = sort(idx_tx((M_training+1):M));
r_testing = r(:, M_training+idx_sort);
% computation of h(Y|X)
h_Y_X = 0;
i_bounds_testing = zeros(1, N+1);
% loop over constellation points to compute h(Y|X)
for n = 1:N
% find how many times x(:, n) was transmitted and update
% i_bounds_testing
N_current_x = find(idx_tx_testing((i_bounds_testing(n)+1):end)==n, 1, 'last');
if isempty(N_current_x), N_current_x=0; end
i_bounds_testing(n+1) = i_bounds_testing(n) + N_current_x;
% add the corresponding contribution to the mutual information (two
% first lines of Eq. (17)). This, together with
% D/2*log2(2*pi) after the end of the loop, gives h(Y|X)
h_Y_X = h_Y_X + N_current_x * log2(det_n(n))/2+...
sum(sum(conj(r_testing(:, (i_bounds_testing(n)+1):i_bounds_testing(n+1))-x(:, n)).*(C_n(:, :, n)\(r_testing(:, (i_bounds_testing(n)+1):i_bounds_testing(n+1))-x(:, n)))))/2/log(2);
end
h_Y_X = D/2*log2(2*pi) + h_Y_X/M_testing;
% When computing log(py), we might run out of memory. If necessary, we
% doe the computation in blocks
logpy = zeros(1, M_testing);
BLOCK_SIZE = floor(MAX_MEMORY/N);
N_blocks = ceil(M_testing/BLOCK_SIZE);
% loop over blocks of symbols. This loop can be replaced by parfor to allow
% parallel computation
for i_block = 1:N_blocks
logpy_cur = zeros(1, M_testing);
% beginning of block
i_start = (i_block-1) * BLOCK_SIZE + 1;
% end of block
i_end = min(M_testing, i_block*BLOCK_SIZE);
% block size
current_block_size = i_end-i_start+1;
% compute exponents of third line of (17)
exponents = zeros(N, current_block_size);
for n = 1:N
exponents(n, :) = -log(det_n(n))/2-real(sum(conj(r_testing(:, i_start:i_end)-x(:, n)).*(C_n(:, :, n)\(r_testing(:, i_start:i_end)-x(:, n))), 1))/2;
%sum über 2 einträge von r
end
% compute third line of Eq. (17). Use a custom function
% that computes log(sum(exp(x))) avoiding overflow errors
logpy_cur(i_start:i_end) = math_logsumexp(log(Px(:))+exponents, 1);
logpy = logpy + logpy_cur;
end
% output entropy h(Y)
h_Y = D/2*log2(2*pi) - mean(logpy)/log(2);%log basis change
% compute mutual information
air = h_Y - h_Y_X;
end
function [y] = math_logsumexp(x, dim)
%[y] = math_logsumexp(x, dim)
% Computes log(sum(exp(x), dim)), avoiding overflow errors when one of the
% x is large.
if nargin<2 || isempty(dim)
m = max(x);
y = m + log(sum(exp(x-m)));
else
m = max(x, [], dim);
y = m + log(sum(exp(x-m), dim));
end
end

48
Functions/Metrics/calc_air.m Normal file
View File

@@ -0,0 +1,48 @@
function [ach_inf_rate] = calc_air(test_signal,reference_signal,options)
% Calculation of AIR acc. to J. Kozesnik, „Numerically Computing Achievable Rates of Memoryless Channels“, Francisco Javier Garcıa-Gomez, doi: 10.1007/978-94-009-9857-5.
% Implementation is not accessible, I mailed TUM to get the code...
arguments(Input)
test_signal;
reference_signal;
options.skip_front = 0;
options.skip_end = 0;
options.returnErrorLocation = 0;
end
options.skip_end = abs(options.skip_end);
options.skip_front = abs(options.skip_front);
assert((options.skip_end+options.skip_front)<length(test_signal),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
if isa(reference_signal,'Signal')
reference_signal = reference_signal.signal;
end
if isa(test_signal,'Signal')
test_signal = test_signal.signal;
end
% TRIM
[test_signal,reference_signal]=trimseq(test_signal,reference_signal,options.skip_front,options.skip_end);
% CALC EVM
%%% new implementation of AIR
constellation = unique(reference_signal);
reference_idx = arrayfun(@(x) find(constellation == x, 1), reference_signal);
ach_inf_rate = air_garcia_implementation(constellation',test_signal',reference_idx');
function [data_,reference_]=trimseq(data,reference,skipstart,skip_end)
data_ = data(skipstart+1:end-skip_end,:);
delta_bits = length(reference) - length(data);
skip_end = delta_bits + skip_end;
reference_ = reference(skipstart+1:end-skip_end,:);
end
end

14
Functions/calc_ber.m → Functions/Metrics/calc_ber.m
View File

@@ -12,9 +12,6 @@ options.skip_front = abs(options.skip_front);
assert((options.skip_end+options.skip_front)<length(data_in),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in"); assert((options.skip_end+options.skip_front)<length(data_in),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
bits = 0;
errors= 0;
ber= 0;
errorIndice= []; errorIndice= [];
% trim sequence to given start; stop. trim reference if signal is shorter % trim sequence to given start; stop. trim reference if signal is shorter
@@ -22,21 +19,15 @@ errorIndice= [];
if length(data_ref) == length(data_in) if length(data_ref) == length(data_in)
bits = 0; bits = numel(data_in(:,options.skip_front+1:end));
bits = bits+numel(data_in(:,options.skip_front+1:end));
try
if options.returnErrorLocation == 0 if options.returnErrorLocation == 0
errors = sum( data_in ~= data_ref,"all" ); errors = sum( data_in ~= data_ref,"all" );
else else
errorIndice = sum(data_in ~= data_ref,1); errorIndice = sum(data_in ~= data_ref,1);
errors = sum(errorIndice ,"all" ); errors = sum(errorIndice ,"all" );
[~,errorIndice] = find(errorIndice~=0); [~,errorIndice] = find(errorIndice~=0);
end errorIndice = errorIndice+options.skip_front;
catch
warning('BER calculation not optimal: Arrays have incompatible sizes for this operation.')
errors = NaN;
end end
% Determine BER % Determine BER
@@ -52,7 +43,6 @@ end
delta_bits = length(reference) - length(data); delta_bits = length(reference) - length(data);
% skip_end = max(skip_end,delta_bits);
skip_end = delta_bits + skip_end; skip_end = delta_bits + skip_end;
reference_ = logical(reference(skipstart+1:end-skip_end,:))'; reference_ = logical(reference(skipstart+1:end-skip_end,:))';

66
Functions/Metrics/calc_evm.m Normal file
View File

@@ -0,0 +1,66 @@
function [evm_total,evm_lvl] = calc_evm(test_signal,reference_signal,options)
arguments(Input)
test_signal;
reference_signal;
options.skip_front = 0;
options.skip_end = 0;
options.returnErrorLocation = 0;
end
options.skip_end = abs(options.skip_end);
options.skip_front = abs(options.skip_front);
assert((options.skip_end+options.skip_front)<length(test_signal),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
if isa(reference_signal,'Signal')
reference_signal = reference_signal.signal;
end
if isa(test_signal,'Signal')
test_signal = test_signal.signal;
end
% TRIM
[test_signal,reference_signal]=trimseq(test_signal,reference_signal,options.skip_front,options.skip_end);
% CALC EVM
[evm_total,evm_lvl] = calc_evm_(test_signal,reference_signal);
function [evm_total,evm_lvl] = calc_evm_(test_signal,reference_signal)
assert(length(test_signal) == length(reference_signal),"Sequence length does not match");
error_vector = (test_signal-reference_signal);
%%% Overall EVM
evm_total = rms(error_vector);
try
%%% Per Level EVM
k = unique(reference_signal);
for lvl = 1:length(k)
lvl_errors = error_vector(reference_signal==k(lvl));
evm_lvl(lvl) = rms(lvl_errors);
end
catch
evm_lvl = NaN;
warning('No EVM per level calculated')
end
end
function [data_,reference_]=trimseq(data,reference,skipstart,skip_end)
data_ = data(skipstart+1:end-skip_end,:);
delta_bits = length(reference) - length(data);
skip_end = delta_bits + skip_end;
reference_ = reference(skipstart+1:end-skip_end,:);
end
end

123
Functions/Metrics/calc_ngmi.m Normal file
View File

@@ -0,0 +1,123 @@
function [GMI,NGMI] = calc_ngmi(test_signal,reference_signal,options)
% Silas implementation of (N)GMI calculation according to: J. Cho, L. Schmalen, und P. J. Winzer,
% „Normalized Generalized Mutual Information as a Forward Error Correction Threshold for Probabilistically Shaped QAM“,
% in 2017 European Conference on Optical Communication (ECOC), Sep. 2017, doi: 10.1109/ECOC.2017.8345872.
% This implementation assumes the same normal distributed noise for each
% channel (sigma2 is calculated once for the whole signal)
arguments(Input)
test_signal;
reference_signal;
options.skip_front = 0;
options.skip_end = 0;
options.returnErrorLocation = 0;
end
options.skip_end = abs(options.skip_end);
options.skip_front = abs(options.skip_front);
assert((options.skip_end+options.skip_front)<length(test_signal),"You can not skip more bits than overall length of data! Set skip_front or skip_end to lower value or check data_in");
if isa(reference_signal,'Signal')
reference_signal = reference_signal.signal;
end
if isa(test_signal,'Signal')
test_signal = test_signal.signal;
end
% TRIM
[test_signal,reference_signal]=trimseq(test_signal,reference_signal,options.skip_front,options.skip_end);
% CALC EVM
[GMI,NGMI] = calc_ngmi_(test_signal,reference_signal);
function [GMI,NGMI] = calc_ngmi_(test_signal,reference_signal)
assert(length(test_signal) == length(reference_signal),"Sequence length does not match");
%%% implemented according to [1] J. Cho, L. Schmalen, und P. J. Winzer,
% „Normalized Generalized Mutual Information as a Forward Error Correction Threshold for Probabilistically Shaped QAM“,
% in 2017 European Conference on Optical Communication (ECOC), Sep. 2017, S. 1–3. doi: 10.1109/ECOC.2017.8345872.
error_vector = (test_signal-reference_signal);
sigma2 = var(error_vector); %noise variance
%%% Separate Classes
constellation = unique(reference_signal);
received_sd = NaN(numel(constellation),length(reference_signal));
lvlcol = cbrewer2('Set1',numel(constellation));
for lvl = 1:numel(constellation)
%Separate the equalized signal into the
%respective levels based on the actually
%transmitted level!
received_sd(lvl,reference_signal==constellation(lvl)) = test_signal(reference_signal==constellation(lvl));
end
N = length(test_signal); % Number of received samples
M = length(constellation); %P
m = log2(M); %bits per symbol
entries = sum(~isnan(received_sd),2)';
P_X = entries./N;
% Parameters
symbols = constellation'; % PAM-4 symbols
gray_bits = PAMmapper(M,0).demap_(constellation); % Gray coding (bits per symbol)
% Conditional probability function for AWGN
q_Y_given_X = @(y, x) (1 / sqrt(2 * pi * sigma2)) * exp(-(y - x).^2 / (2 * sigma2));
% Entropy term
H_X = -sum(P_X .* log2(P_X)); % Entropy of input distribution
% GMI computation
noise_impact_term = 0;
for k = 1:N
y_k = test_signal(k); % Current received sample
[~, closest_symbol_idx] = min(abs(symbols - y_k)); % Closest symbol index
closest_symbol = symbols(closest_symbol_idx); % Closest symbol
for i = 1:m
% Extract i-th bit for each symbol
bit_mask = gray_bits(:, i); % Binary column for i-th bit of all symbols
matching_symbols = symbols(bit_mask == gray_bits(closest_symbol_idx, i));
% Numerator: Sum over x in x_{b_{k, i}}
numerator = sum(q_Y_given_X(y_k, matching_symbols) .* P_X(ismember(symbols, matching_symbols)));
% Denominator: Sum over all x
denominator = sum(q_Y_given_X(y_k, symbols) .* P_X);
% Logarithmic contribution
noise_impact_term = noise_impact_term + log2(numerator / denominator);
end
end
% Normalize the noise impact term by N
noise_impact_term = noise_impact_term / N;
% GMI
GMI = H_X + noise_impact_term;
NGMI = GMI / m;
end
function [data_,reference_]=trimseq(data,reference,skipstart,skip_end)
data_ = data(skipstart+1:end-skip_end,:);
delta_bits = length(reference) - length(data);
skip_end = delta_bits + skip_end;
reference_ = reference(skipstart+1:end-skip_end,:);
end
end

35
Functions/beautifyBERplot.m Normal file
View File

@@ -0,0 +1,35 @@
function beautifyBERplot()
% BEAUTIFYBERPLOT Enhances a BER plot for publication-quality figures.
% Set line properties for all current plot lines
lines = findall(gca, 'Type', 'Line');
markers = {'o', 's', 'd', '^', 'v', '>', '<', 'p', 'h'}; % Define marker styles
num_markers = length(markers);
for i = 1:length(lines)
lines(i).LineWidth = 1.3; % Thicker line width
%lines(i).LineStyle = '-'; % Solid lines for simplicity
lines(i).Marker = markers{mod(i-1, num_markers) + 1}; % Assign markers cyclically
lines(i).MarkerSize = 4; % Marker size
lines(i).MarkerFaceColor = 'auto'; % Use line color for marker face
end
% Change all text interpreters to LaTeX
set(findall(gca, '-property', 'Interpreter'), 'Interpreter', 'latex');
% Set figure background to white
set(gcf, 'Color', 'w');
% Set logarithmic scale for y-axis, but only if it makes sense.
% If this is not always desired, you could condition this on the presence of lines or data.
% set(gca, 'YScale', 'log');
% Customize grid and box appearance
set(gca, 'Box', 'on', 'LineWidth', 0.8); % Thicker border
grid on;
% grid minor;
% Adjust font size and style for better readability
set(gca, 'FontSize', 10, 'FontName', 'Times New Roman');
end

19
Functions/calc_evm.m
View File

@@ -1,19 +0,0 @@
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

5
Functions/channel_structures/awgn_channel.m Normal file
View File

@@ -0,0 +1,5 @@
function signal_out = awgn_channel(signal_in)
signal_out = signal_in;
end

76
Functions/copyStylingFrom.m Normal file
View File

@@ -0,0 +1,76 @@
function copyStylingFrom(figNumSource, figNumTgt)
% Get handles to the source and target figures
sourceFig = figure(figNumSource);
targetFig = figure(figNumTgt);
% Get axes of source and target figures
sourceAxes = findall(sourceFig, 'type', 'axes');
targetAxes = findall(targetFig, 'type', 'axes');
% Ensure the number of axes match
if length(sourceAxes) ~= length(targetAxes)
error('Number of axes in source and target figures must be the same.');
end
% Loop through each pair of axes and copy styling properties
for i = 1:length(sourceAxes)
copyAxesProperties(sourceAxes(i), targetAxes(i));
end
% Apply general figure properties if desired
targetFig.Color = sourceFig.Color; % Background color
end
function copyAxesProperties(sourceAx, targetAx)
% List of properties to copy from source to target axes
propsToCopy = {'XColor', 'YColor', 'ZColor', 'FontSize', 'FontName', ...
'GridColor', 'GridLineStyle', 'MinorGridColor', 'Box', ...
'XGrid', 'YGrid', 'ZGrid', 'XMinorGrid', 'YMinorGrid', 'ZMinorGrid', ...
'LineWidth', 'TitleFontSizeMultiplier', 'LabelFontSizeMultiplier'};
% Copy properties from source to target
for i = 1:length(propsToCopy)
try
targetAx.(propsToCopy{i}) = sourceAx.(propsToCopy{i});
catch
% Skip property if it doesn't exist or can't be copied
end
end
% Copy axis labels and titles
targetAx.Title.String = sourceAx.Title.String;
targetAx.XLabel.String = sourceAx.XLabel.String;
targetAx.YLabel.String = sourceAx.YLabel.String;
targetAx.ZLabel.String = sourceAx.ZLabel.String;
% Copy children elements like lines, patches, etc.
sourceChildren = allchild(sourceAx);
targetChildren = allchild(targetAx);
% Ensure the number of children elements match
if length(sourceChildren) ~= length(targetChildren)
warning('Number of elements in source and target axes differ. Styling may not be applied completely.');
end
% Copy properties of children (like lines, patches, etc.), except colors and legends
for i = 1:min(length(sourceChildren), length(targetChildren))
copyObjectProperties(sourceChildren(i), targetChildren(i));
end
end
function copyObjectProperties(sourceObj, targetObj)
% List of common properties to copy for plot elements (lines, patches, etc.)
propsToCopy = {'LineStyle', 'LineWidth', 'Marker', 'MarkerSize', ...
'MarkerEdgeColor', 'MarkerFaceColor', 'DisplayName'};
% Copy properties from source to target, excluding colors
for i = 1:length(propsToCopy)
try
if ~contains(propsToCopy{i}, 'Color') % Skip color properties
targetObj.(propsToCopy{i}) = sourceObj.(propsToCopy{i});
end
catch
% Skip property if it doesn't exist or can't be copied
end
end
end

260
Functions/helper_functions_community/cbrewer2/cbrewer2.m
View File

@@ -1,260 +0,0 @@
%CBREWER2 Interpolated versions of Cynthia Brewer's ColorBrewer colormaps
% CBREWER2(CNAME, NCOL) returns the colour scheme CNAME with the number
% of colours equal to NCOL. If there is a ColorBrewer scheme with exactly
% this number of colours, the color scheme is returned as-is. If NCOL
% larger (or smaller) than the designed colormaps for this scheme, the
% largest (smallest) one is interpolated to provide enough colours,
% unless the requested colour scheme CNAME is a qualitative palette. For
% a qualitative scheme, the colours are repeated, cycling from the
% beginning again, to output the requested NCOL colours.
%
% CBREWER2(CNAME) without an NCOL input will use the same number of
% colours as the current colormap.
%
% CBREWER2(CNAME, NCOL, INTERP_METHOD) allows you to change the method
% used for the interpolation. The default is 'cubic'.
%
% CBREWER2(CNAME, NCOL, INTERP_METHOD, INTERP_SPACE) allows you to
% change the colorspace used for the interpolation. By default, this is
% in the CIELAB colorspace, which is approximately perceptually uniform.
% Options for INTERP_SPACE are
% 'rgb' : interpolation in sRGB (as used in original CBREWER)
% 'lab' : interpolation in CIELAB (default)
% 'lch' : interpolation in CIELCH_ab (not recommended due to the
% discontinuities at C=0 and H=0)
% Anything else supported by COLORSPACE will also function.
%
% The input format CBREWER2(TYPE, ...) can also be used, where TYPE is
% one of 'seq', 'div', 'qual'. This input is redandant and will be
% ignored. This input format is provided for backwards compatibility with
% the original CBREWER.
%
% Example 1 (sequential heatmap):
% C = [0 2 4 6; 8 10 12 14; 16 18 20 22];
% imagesc(C);
% colorbar;
% colormap(cbrewer('YlOrRd', 256);
%
% Example 2 (line plot):
% x = 0:0.01:2;
% sc = [0.5; 1; 2];
% t0 = [0; 0.2; 0.4];
% t = bsxfun(@rdivide, bsxfun(@plus, x, t0), sc);
% y = sin(t * 2 * pi);
% cmap = cbrewer2('Set1', numel(sc));
% axes('ColorOrder', cmap, 'NextPlot', 'ReplaceChildren');
% plot(x, y);
%
% Example 3 (divergent heatmap):
% [X,Y,Z] = peaks(30);
% surfc(X,Y,Z);
% colormap(cbrewer2('RdBu'));
%
% This product includes color specifications and designs developed by
% Cynthia Brewer (http://colorbrewer.org/). For more information on
% ColorBrewer, please visit http://colorbrewer.org/.
%
% CBREWER2 uses a cached copy of the Cynthia Brewer color schemes which
% was converted to .mat format by Charles Robert for use with CBREWER.
% CBREWER is available from the MATLAB FileExchange under the MIT license.
%
% See also CBREWER, BREWERMAP, COLORSPACE, INTERP1.
% Copyright (c) 2016 Scott Lowe
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
function colormap = cbrewer2(...
cname, ncol, interp_method, interp_space, varargin)
% Definitions -------------------------------------------------------------
% List of all of Cynthia Brewer's colormaps and their types
% seq: sequential
% div: divergent
% qual: qualitative
cbdict = {...
'Blues', 'seq'; ...
'BuGn', 'seq'; ...
'BuPu', 'seq'; ...
'GnBu', 'seq'; ...
'Greens', 'seq'; ...
'Greys', 'seq'; ...
'Oranges', 'seq'; ...
'OrRd', 'seq'; ...
'PuBu', 'seq'; ...
'PuBuGn', 'seq'; ...
'PuRd', 'seq'; ...
'Purples', 'seq'; ...
'RdPu', 'seq'; ...
'Reds', 'seq'; ...
'YlGn', 'seq'; ...
'YlGnBu', 'seq'; ...
'YlOrBr', 'seq'; ...
'YlOrRd', 'seq'; ...
'BrBG', 'div'; ...
'PiYG', 'div'; ...
'PRGn', 'div'; ...
'PuOr', 'div'; ...
'RdBu', 'div'; ...
'RdGy', 'div'; ...
'RdYlBu', 'div'; ...
'RdYlGn', 'div'; ...
'Spectral', 'div'; ...
'Accent', 'qual'; ...
'Dark2', 'qual'; ...
'Paired', 'qual'; ...
'Pastel1', 'qual'; ...
'Pastel2', 'qual'; ...
'Set1', 'qual'; ...
'Set2', 'qual'; ...
'Set3', 'qual'; ...
};
% Input handling ----------------------------------------------------------
narginchk(1, 5);
% Initialise variables if not supplied
if nargin<2
ncol = [];
end
if nargin<3
interp_method = [];
end
if nargin<4
interp_space = [];
end
if nargin<5
varargin = {[]};
end
% Check if the colormap type was unnecessarily input
types = unique(cbdict(:, 2));
if nargin > 1 && ischar(cname) && ischar(ncol)
LI = ismember({cname ncol}, types);
if ~any(LI); error('Number of colors cannot be a string'); end;
if all(LI); error('Incorrect colormap name'); end;
if LI(1)
vgn = {cname; ncol; interp_method; interp_space};
cname = vgn{2};
ncol = vgn{3};
interp_method = vgn{4};
interp_space = varargin{1};
ctype_input = vgn{1};
elseif LI(2)
vgn = {cname; ncol; interp_method; interp_space};
cname = vgn{1};
ncol = vgn{3};
interp_method = vgn{4};
interp_space = varargin{1};
ctype_input = vgn{2};
end
else
ctype_input = '';
end
% Default values
if isempty(ncol)
% Number of colours in the colormap
ncol = size(get(gcf,'colormap'), 1);
end
if isempty(interp_method)
interp_method = 'pchip';
end
if isempty(interp_space)
interp_space = 'lab';
end
% Load colorbrewer data ---------------------------------------------------
Tmp = load('colorbrewer.mat');
colorbrewer = Tmp.colorbrewer;
[TF, idict] = ismember(lower(cname), lower(cbdict(:, 1)));
if ~TF
error('%s is not a recognised Brewer colormap',cname);
end
cname = cbdict{idict, 1};
ctype = cbdict{idict, 2};
if (~isfield(colorbrewer.(ctype), cname))
error('Colormap %s is not present in loaded data',cname);
end
% Main script -------------------------------------------------------------
if ncol > length(colorbrewer.(ctype).(cname))
% If we specified too many colours, we take the maximum and interpolate
colormap = colorbrewer.(ctype).(cname){length(colorbrewer.(ctype).(cname))};
colormap = colormap ./ 255;
elseif isempty(colorbrewer.(ctype).(cname){ncol})
% If we specified too few colours, we take the minimum and interpolate
nmin = find(~cellfun(@isempty, colorbrewer.(ctype).(cname)), 1);
colormap = colorbrewer.(ctype).(cname){nmin};
colormap = colormap./255;
else
% If we specified a number of colours in the pre-designed range, no
% need to interpolate
colormap = (colorbrewer.(ctype).(cname){ncol}) ./ 255;
return;
end
% Don't interpolate if qualitative type
if strcmp(ctype,'qual')
if size(colormap, 1) >= ncol
colormap = colormap(1:ncol, :);
return;
end
warning('CBREWER2:QualTooManyColors', ...
['Too many colors requested: cannot interpolate a qualitative' ...
' colorscheme']);
% Cycle the colours from the beginning again, so we have enough to
% return
colormap = repmat(colormap, ceil(ncol / size(colormap, 1)), 1);
colormap = colormap(1:ncol, :);
return;
end
% Make sure we have colorspace downloaded from the FEX
if ~strcmpi(interp_space, 'rgb') && ~exist('colorspace.m', 'file')
P = requireFEXpackage(28790);
if isempty(P);
error(...
['You need to download COLORSPACE from the MATLAB FEX and' ...
' add it to the MATLAB path.']);
end;
end
% Move to perceptually uniform space
if ~strcmpi(interp_space,'rgb')
colormap = colorspace(['rgb->' interp_space], colormap);
end
% Linearly interpolate
X = linspace(0, 1, size(colormap, 1));
XI = linspace(0, 1, ncol);
colormap = interp1(X, colormap, XI, interp_method);
% Move from perceptually uniform space back to sRGB
if ~strcmpi(interp_space,'rgb')
colormap = colorspace(['rgb<-' interp_space], colormap);
end
end

BIN
Functions/helper_functions_community/cbrewer2/colorbrewer.mat
View File

Binary file not shown.

194
Functions/helper_functions_community/cbrewer2/requireFEXpackage.m
View File

@@ -1,194 +0,0 @@
function AddedPath = requireFEXpackage(FEXSubmissionID)
%Function requireFEXpackage -
%installs Matlab Central File Exchange (FEX) submission
%with given ID into the directory chosen by the user.
%A new FEX submissions may use previous FEX submissions as its part.
%The function 'requireFEXpackage' helps in adding those previous
%submissions to the user's MATLAB installation.
%
%This function is a part of File Exchange submission 31069.
%Download the entire submission:
%http://www.mathworks.com/matlabcentral/fileexchange/31069
%
% SYNTAX:
% AddedPath = requireFEXpackage(FEXSubmissionID)
%
% INPUT:
% ID of the required submission to File Exchange
%
% OUTPUT:
% the path to that submission added to the user's MATLAB path.
%
% HOW TO CALL:
% The command
% P = requireFEXpackage(8277)
% will download and install the package with ID 8277
% (namely, nice 'fminsearchbnd' by John D'Errico)
%
% EXAMPLES -- HOW TO USE:
%
% EXAMPLE 1 (using 'exist' command):
%
% % first, somewhere in the very beginning of your code,
% % check if the function 'fminsearchbnd' from the FEX package 8277
% % is on your MATLAB path, and if it is not there,
% % require the FEX package 8277:
% if ~(exist('fminsearchbnd', 'file') == 2)
% P = requireFEXpackage(8277); % fminsearchbnd is part of 8277
% end
%
% % Then just use 'fminsearchbnd' where you need it:
% syms x
% RosenbrockBananaFunction = @(x) (1-x(1)).^2 + 100*(x(2)-x(1).^2).^2;
% x = fminsearchbnd(RosenbrockBananaFunction,[3 3])
% EXAMPLE 2 (using 'try-catch' command):
%
% syms x
% RosenbrockBananaFunction = @(x) (1-x(1)).^2+100*(x(2)-x(1).^2).^2;
% try
% % if function 'fminsearchbnd' already exists in your MATLAB
% % installation, just use it:
% x = fminsearchbnd(RosenbrockBananaFunction,[3 3])
% catch
% % if function 'fminsearchbnd' is not present in your MATLAB
% % installation, first get the package 8277 (to which it belongs)
% % from the MATLAB Central File Exchange (FEX)
% P = requireFEXpackage(8277); % fminsearchbnd is part of 8277
% % and then use that function:
% x = fminsearchbnd(RosenbrockBananaFunction,[3 3])
% end
%
%
% NOTE: on Mac platform, the title of the dialog box for
% choosing the directory for installing the required FEX package
% is not shown; this is not a bug, this is how UIGETDIR works on Macs --
% see the documentation for UIGETDIR
% http://www.mathworks.com/help/techdoc/ref/uigetdir.html
%
% (C) Igor Podlubny, 2011
% Copyright (c) 2011, Igor Podlubny
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
ID = num2str(FEXSubmissionID);
% Ask user for the confirmation of the installation
% of the required FEX package
yes = ['YES, Install package ' ID];
no = 'NO, do not install';
userchoice = questdlg(['The Matlab function/toolbox, which you are running, ' ...
'requires the presence of the package ' ID ...
' from Matlab Central File Exchange.' ...
sprintf('\n\n') ...
'Would you like to install the FEX package ' ID ' now?'] , ...
['Required package ' ID], ...
yes, no, yes);
% Handle response
switch userchoice
case yes,
install = 1;
case no,
install = 0;
otherwise,
install = 0;
end
if install == 1
baseURL = 'http://www.mathworks.com/matlabcentral/fileexchange/';
query = '?download=true';
location = uigetdir(pwd, ['Select the directory for installing the required FEX package' ID ]);
if location ~= 0
% download package 'ID' from Matlab Central File Exchange
filetosave = [location filesep ID '.zip'];
FEXpackage = [baseURL ID query];
[f,status] = urlwrite(FEXpackage,filetosave);
if status==0
warndlg(['No connection to Matlab Central File Exchange,' ' or package ' ID ' does not exist.' ...
' Package ' ID ' has not been installed. ' ...
' Check you internet settings and the ID of the required package, and try again. '] , ...
['No connection to Matlab Central File Exchange' ' or package ' ID ' does not exist'], ...
'modal');
AddedPath = '';
return
end
% unzip the downloaded file to the subdirectory 'ID'
todir = [location filesep ID];
% if the directory 'ID' doesn't exist at given location, create it
if ~(exist([location filesep ID], 'dir') == 7)
mkdir(location, ID);
end
try
unzip(filetosave, todir);
% after unzipping, delete the downloaded ZIP file
delete(filetosave);
% prepend the paths to the downloaded package to the MATLAB path
P = genpath([location filesep ID]);
path(P,path);
catch
% if the FEX package is not ZIP, then it is a single m-file
% just move the file to the ID directory
[pathstr, name, ext] = fileparts(filetosave);
movefile(filetosave, [todir filesep name '.m']);
P = genpath([location filesep ID]);
path(P,path);
end
else
P = '';
end
AddedPath = P;
else
AddedPath = '';
end
if install == 1,
% Ask user about reviewing and saving the modified MATLAB path,
% and take him to PATHTOOL, if the user wants to save the modified path
yes = 'YES, I want to review and save the MATLAB path';
no = 'NO, I don''t want to save the path permanently';
userchoice = questdlg(['After adding the package ' ID ...
' from Matlab Central File Exchange to your MATLAB installation,' ...
' the MATLAB path has been modified accordingly. ', ...
'Would you like to review and save the modified MATLAB path?'] , ...
'Review and save the modified MATLAB path for future use?', ...
yes, no, yes);
% Handle response
switch userchoice
case yes,
pathtool;
case no,
otherwise,
end
end

261
Functions/helper_functions_community/linspecer.m
View File

@@ -1,261 +0,0 @@
% function lineStyles = linspecer(N)
% This function creates an Nx3 array of N [R B G] colors
% These can be used to plot lots of lines with distinguishable and nice
% looking colors.
%
% lineStyles = linspecer(N); makes N colors for you to use: lineStyles(ii,:)
%
% colormap(linspecer); set your colormap to have easily distinguishable
% colors and a pleasing aesthetic
%
% lineStyles = linspecer(N,'qualitative'); forces the colors to all be distinguishable (up to 12)
% lineStyles = linspecer(N,'sequential'); forces the colors to vary along a spectrum
%
% % Examples demonstrating the colors.
%
% LINE COLORS
% N=6;
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sin(x+2*n*pi/N), X.', 1:N);
% C = linspecer(N);
% axes('NextPlot','replacechildren', 'ColorOrder',C);
% plot(X,Y,'linewidth',5)
% ylim([-1.1 1.1]);
%
% SIMPLER LINE COLOR EXAMPLE
% N = 6; X = linspace(0,pi*3,1000);
% C = linspecer(N)
% hold off;
% for ii=1:N
% Y = sin(X+2*ii*pi/N);
% plot(X,Y,'color',C(ii,:),'linewidth',3);
% hold on;
% end
%
% COLORMAP EXAMPLE
% A = rand(15);
% figure; imagesc(A); % default colormap
% figure; imagesc(A); colormap(linspecer); % linspecer colormap
%
% See also NDHIST, NHIST, PLOT, COLORMAP, 43700-cubehelix-colormaps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% by Jonathan Lansey, March 2009-2013 – Lansey at gmail.com %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% credits and where the function came from
% The colors are largely taken from:
% http://colorbrewer2.org and Cynthia Brewer, Mark Harrower and The Pennsylvania State University
%
%
% She studied this from a phsychometric perspective and crafted the colors
% beautifully.
%
% I made choices from the many there to decide the nicest once for plotting
% lines in Matlab. I also made a small change to one of the colors I
% thought was a bit too bright. In addition some interpolation is going on
% for the sequential line styles.
%
%
%%
function lineStyles=linspecer(N,varargin)
if nargin==0 % return a colormap
lineStyles = linspecer(128);
return;
end
if ischar(N)
lineStyles = linspecer(128,N);
return;
end
if N<=0 % its empty, nothing else to do here
lineStyles=[];
return;
end
% interperet varagin
qualFlag = 0;
colorblindFlag = 0;
if ~isempty(varargin)>0 % you set a parameter?
switch lower(varargin{1})
case {'qualitative','qua'}
if N>12 % go home, you just can't get this.
warning('qualitiative is not possible for greater than 12 items, please reconsider');
else
if N>9
warning(['Default may be nicer for ' num2str(N) ' for clearer colors use: whitebg(''black''); ']);
end
end
qualFlag = 1;
case {'sequential','seq'}
lineStyles = colorm(N);
return;
case {'white','whitefade'}
lineStyles = whiteFade(N);return;
case 'red'
lineStyles = whiteFade(N,'red');return;
case 'blue'
lineStyles = whiteFade(N,'blue');return;
case 'green'
lineStyles = whiteFade(N,'green');return;
case {'gray','grey'}
lineStyles = whiteFade(N,'gray');return;
case {'colorblind'}
colorblindFlag = 1;
otherwise
warning(['parameter ''' varargin{1} ''' not recognized']);
end
end
% *.95
% predefine some colormaps
set3 = colorBrew2mat({[141, 211, 199];[ 255, 237, 111];[ 190, 186, 218];[ 251, 128, 114];[ 128, 177, 211];[ 253, 180, 98];[ 179, 222, 105];[ 188, 128, 189];[ 217, 217, 217];[ 204, 235, 197];[ 252, 205, 229];[ 255, 255, 179]}');
set1JL = brighten(colorBrew2mat({[228, 26, 28];[ 55, 126, 184]; [ 77, 175, 74];[ 255, 127, 0];[ 255, 237, 111]*.85;[ 166, 86, 40];[ 247, 129, 191];[ 153, 153, 153];[ 152, 78, 163]}'));
set1 = brighten(colorBrew2mat({[ 55, 126, 184]*.85;[228, 26, 28];[ 77, 175, 74];[ 255, 127, 0];[ 152, 78, 163]}),.8);
% colorblindSet = {[215,25,28];[253,174,97];[171,217,233];[44,123,182]};
colorblindSet = {[215,25,28];[253,174,97];[171,217,233]*.8;[44,123,182]*.8};
set3 = dim(set3,.93);
if colorblindFlag
switch N
% sorry about this line folks. kind of legacy here because I used to
% use individual 1x3 cells instead of nx3 arrays
case 4
lineStyles = colorBrew2mat(colorblindSet);
otherwise
colorblindFlag = false;
warning('sorry unsupported colorblind set for this number, using regular types');
end
end
if ~colorblindFlag
switch N
case 1
lineStyles = { [ 55, 126, 184]/255};
case {2, 3, 4, 5 }
lineStyles = set1(1:N);
case {6 , 7, 8, 9}
lineStyles = set1JL(1:N)';
case {10, 11, 12}
if qualFlag % force qualitative graphs
lineStyles = set3(1:N)';
else % 10 is a good number to start with the sequential ones.
lineStyles = cmap2linspecer(colorm(N));
end
otherwise % any old case where I need a quick job done.
lineStyles = cmap2linspecer(colorm(N));
end
end
lineStyles = cell2mat(lineStyles);
end
% extra functions
function varIn = colorBrew2mat(varIn)
for ii=1:length(varIn) % just divide by 255
varIn{ii}=varIn{ii}/255;
end
end
function varIn = brighten(varIn,varargin) % increase the brightness
if isempty(varargin),
frac = .9;
else
frac = varargin{1};
end
for ii=1:length(varIn)
varIn{ii}=varIn{ii}*frac+(1-frac);
end
end
function varIn = dim(varIn,f)
for ii=1:length(varIn)
varIn{ii} = f*varIn{ii};
end
end
function vOut = cmap2linspecer(vIn) % changes the format from a double array to a cell array with the right format
vOut = cell(size(vIn,1),1);
for ii=1:size(vIn,1)
vOut{ii} = vIn(ii,:);
end
end
%%
% colorm returns a colormap which is really good for creating informative
% heatmap style figures.
% No particular color stands out and it doesn't do too badly for colorblind people either.
% It works by interpolating the data from the
% 'spectral' setting on http://colorbrewer2.org/ set to 11 colors
% It is modified a little to make the brightest yellow a little less bright.
function cmap = colorm(varargin)
n = 100;
if ~isempty(varargin)
n = varargin{1};
end
if n==1
cmap = [0.2005 0.5593 0.7380];
return;
end
if n==2
cmap = [0.2005 0.5593 0.7380;
0.9684 0.4799 0.2723];
return;
end
frac=.95; % Slight modification from colorbrewer here to make the yellows in the center just a bit darker
cmapp = [158, 1, 66; 213, 62, 79; 244, 109, 67; 253, 174, 97; 254, 224, 139; 255*frac, 255*frac, 191*frac; 230, 245, 152; 171, 221, 164; 102, 194, 165; 50, 136, 189; 94, 79, 162];
x = linspace(1,n,size(cmapp,1));
xi = 1:n;
cmap = zeros(n,3);
for ii=1:3
cmap(:,ii) = pchip(x,cmapp(:,ii),xi);
end
cmap = flipud(cmap/255);
end
function cmap = whiteFade(varargin)
n = 100;
if nargin>0
n = varargin{1};
end
thisColor = 'blue';
if nargin>1
thisColor = varargin{2};
end
switch thisColor
case {'gray','grey'}
cmapp = [255,255,255;240,240,240;217,217,217;189,189,189;150,150,150;115,115,115;82,82,82;37,37,37;0,0,0];
case 'green'
cmapp = [247,252,245;229,245,224;199,233,192;161,217,155;116,196,118;65,171,93;35,139,69;0,109,44;0,68,27];
case 'blue'
cmapp = [247,251,255;222,235,247;198,219,239;158,202,225;107,174,214;66,146,198;33,113,181;8,81,156;8,48,107];
case 'red'
cmapp = [255,245,240;254,224,210;252,187,161;252,146,114;251,106,74;239,59,44;203,24,29;165,15,21;103,0,13];
otherwise
warning(['sorry your color argument ' thisColor ' was not recognized']);
end
cmap = interpomap(n,cmapp);
end
% Eat a approximate colormap, then interpolate the rest of it up.
function cmap = interpomap(n,cmapp)
x = linspace(1,n,size(cmapp,1));
xi = 1:n;
cmap = zeros(n,3);
for ii=1:3
cmap(:,ii) = pchip(x,cmapp(:,ii),xi);
end
cmap = (cmap/255); % flipud??
end

47
Functions/removeDotsFromFileNames.m Normal file
View File

@@ -0,0 +1,47 @@
folderPath = "/Volumes/NT-Labor/2024/sioe/High Speed Messungen Oktober/";
% Get list of all files in the folder and subfolders
fileList = dir(fullfile(folderPath, '**', '*'));
% Loop through each file
for i = 1:length(fileList)
% Get the file name and path
oldFileName = fileList(i).name;
oldFilePath = fullfile(fileList(i).folder, oldFileName);
% Skip if it’s a directory or hidden file
if fileList(i).isdir || startsWith(oldFileName, '.')
continue;
end
ismat = strfind(oldFileName, '.mat');
if ismat
continue;
else
% Find position of last dot (for file extension)
dotIndex = strfind(oldFileName, '.');
% Only proceed if there is more than one dot in the file name
if ~isempty(dotIndex)
% Replace all dots in the "base name" with underscores
sanitizedBaseName = strrep(oldFileName, '.', '_');
% Create the new file name with the corrected extension
newFileName = [sanitizedBaseName, '.mat'];
% Full path of the new file
newFilePath = fullfile(fileList(i).folder, newFileName);
% Rename the file
movefile(oldFilePath, newFilePath);
fprintf('Renamed: %s -> %s\n', oldFilePath, newFilePath);
end
end
end

0
Functions/helper_functions_community/Violinplot-Matlab-master/LICENSE → Libs/Violinplot-Matlab-master/LICENSE
View File

0
Functions/helper_functions_community/Violinplot-Matlab-master/README.md → Libs/Violinplot-Matlab-master/README.md
View File

0
Functions/helper_functions_community/Violinplot-Matlab-master/Violin.m → Libs/Violinplot-Matlab-master/Violin.m
View File

0
Functions/helper_functions_community/Violinplot-Matlab-master/example.png → Libs/Violinplot-Matlab-master/example.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 42 KiB

After

Width:  |  Height:  |  Size: 42 KiB

0
Functions/helper_functions_community/Violinplot-Matlab-master/example2.png → Libs/Violinplot-Matlab-master/example2.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 34 KiB

After

Width:  |  Height:  |  Size: 34 KiB

0
Functions/helper_functions_community/Violinplot-Matlab-master/test_cases/testviolinplot.m → Libs/Violinplot-Matlab-master/test_cases/testviolinplot.m
View File

0
Functions/helper_functions_community/Violinplot-Matlab-master/violinplot.m → Libs/Violinplot-Matlab-master/violinplot.m
View File

0
Functions/helper_functions_community/autoArrangeFigures.m → Libs/autoArrangeFigures.m
View File

0
Functions/helper_functions_community/boundedlines/.gitignore → Libs/boundedlines/.gitignore vendored
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/.gitignore → Libs/boundedlines/Inpaint_nans/.gitignore vendored
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo.html → Libs/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo.html
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo.png → Libs/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 3.1 KiB

After

Width:  |  Height:  |  Size: 3.1 KiB

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_01.png → Libs/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_01.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 16 KiB

After

Width:  |  Height:  |  Size: 16 KiB

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_02.png → Libs/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_02.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 17 KiB

After

Width:  |  Height:  |  Size: 17 KiB

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_03.png → Libs/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_03.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 16 KiB

After

Width:  |  Height:  |  Size: 16 KiB

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_04.png → Libs/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_04.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 16 KiB

After

Width:  |  Height:  |  Size: 16 KiB

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_05.png → Libs/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_05.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 13 KiB

After

Width:  |  Height:  |  Size: 13 KiB

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_06.png → Libs/boundedlines/Inpaint_nans/demo/html/inpaint_nans_demo_06.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 9.2 KiB

After

Width:  |  Height:  |  Size: 9.2 KiB

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/demo/inpaint_nans_demo_old.m → Libs/boundedlines/Inpaint_nans/demo/inpaint_nans_demo_old.m
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/doc/Nomination comments.rtf → Libs/boundedlines/Inpaint_nans/doc/Nomination comments.rtf
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/doc/methods_of_inpaint_nans.m → Libs/boundedlines/Inpaint_nans/doc/methods_of_inpaint_nans.m
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/garden50.jpg → Libs/boundedlines/Inpaint_nans/garden50.jpg
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 22 KiB

After

Width:  |  Height:  |  Size: 22 KiB

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/inpaint_nans.m → Libs/boundedlines/Inpaint_nans/inpaint_nans.m
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/inpaint_nans_bc.m → Libs/boundedlines/Inpaint_nans/inpaint_nans_bc.m
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/inpaint_nans_demo.m → Libs/boundedlines/Inpaint_nans/inpaint_nans_demo.m
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/license.txt → Libs/boundedlines/Inpaint_nans/license.txt
View File

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/monet_adresse.jpg → Libs/boundedlines/Inpaint_nans/monet_adresse.jpg
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 145 KiB

After

Width:  |  Height:  |  Size: 145 KiB

0
Functions/helper_functions_community/boundedlines/Inpaint_nans/test/test_main.m → Libs/boundedlines/Inpaint_nans/test/test_main.m
View File

0
Functions/helper_functions_community/boundedlines/LICENSE.txt → Libs/boundedlines/LICENSE.txt
View File

0
Functions/helper_functions_community/boundedlines/README.html → Libs/boundedlines/README.html
View File

0
Functions/helper_functions_community/boundedlines/README.m → Libs/boundedlines/README.m
View File

0
Functions/helper_functions_community/boundedlines/README.md → Libs/boundedlines/README.md
View File

0
Functions/helper_functions_community/boundedlines/boundedline/.gitignore → Libs/boundedlines/boundedline/.gitignore vendored
View File

0
Functions/helper_functions_community/boundedlines/boundedline/boundedline.m → Libs/boundedlines/boundedline/boundedline.m
View File

0
Functions/helper_functions_community/boundedlines/boundedline/outlinebounds.m → Libs/boundedlines/boundedline/outlinebounds.m
View File

0
Functions/helper_functions_community/boundedlines/catuneven/.gitignore → Libs/boundedlines/catuneven/.gitignore vendored
View File

0
Functions/helper_functions_community/boundedlines/catuneven/catuneven.m → Libs/boundedlines/catuneven/catuneven.m
View File

0
Functions/helper_functions_community/boundedlines/readmeExtras/README_01.png → Libs/boundedlines/readmeExtras/README_01.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 17 KiB

After

Width:  |  Height:  |  Size: 17 KiB

0
Functions/helper_functions_community/boundedlines/readmeExtras/README_02.png → Libs/boundedlines/readmeExtras/README_02.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 26 KiB

After

Width:  |  Height:  |  Size: 26 KiB

0
Functions/helper_functions_community/boundedlines/readmeExtras/README_03.png → Libs/boundedlines/readmeExtras/README_03.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 34 KiB

After

Width:  |  Height:  |  Size: 34 KiB

0
Functions/helper_functions_community/boundedlines/readmeExtras/README_04.png → Libs/boundedlines/readmeExtras/README_04.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 44 KiB

After

Width:  |  Height:  |  Size: 44 KiB

0
Functions/helper_functions_community/boundedlines/readmeExtras/README_05.png → Libs/boundedlines/readmeExtras/README_05.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 16 KiB

After

Width:  |  Height:  |  Size: 16 KiB

0
Functions/helper_functions_community/boundedlines/readmeExtras/README_06.png → Libs/boundedlines/readmeExtras/README_06.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 24 KiB

After

Width:  |  Height:  |  Size: 24 KiB

0
Functions/helper_functions_community/boundedlines/readmeExtras/README_07.png → Libs/boundedlines/readmeExtras/README_07.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 25 KiB

After

Width:  |  Height:  |  Size: 25 KiB

0
Functions/helper_functions_community/boundedlines/readmeExtras/README_08.png → Libs/boundedlines/readmeExtras/README_08.png
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 24 KiB

After

Width:  |  Height:  |  Size: 24 KiB

0
Functions/helper_functions_community/boundedlines/singlepatch/.gitignore → Libs/boundedlines/singlepatch/.gitignore vendored
View File

0
Functions/helper_functions_community/boundedlines/singlepatch/singlepatch.m → Libs/boundedlines/singlepatch/singlepatch.m
View File

Some files were not shown because too many files have changed in this diff Show More