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

- bring FWM plots back to life!

Browse Source 
- some dispersion plots with the great help of chatGPT :-D
This commit is contained in:
Silas Oettinghaus
2025-10-17 11:50:27 +02:00
parent 99898da519
commit 3ea6439947
13 changed files with 720 additions and 62 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
Classes/DataBaseHandler/DBHandler.m
View File

@@ -588,14 +588,13 @@ classdef DBHandler < handle
return
catch ME
if attempt < maxFast
pause(1)
pause(0.1)
else
pause(10)
pause(1)
end
lastErr = ME;
end
end
pause(60)
error('Database fetch failed after %d attempts:\n%s', maxSlow, lastErr.getReport())
end

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

@@ -137,7 +137,7 @@ classdef DataStorage < handle
if ~isempty(tmp)
if isa(tmp,'double')
value(i) = tmp ;
value(i,:) = tmp ;
elseif isa(tmp,'Signal') || isa(tmp,'struct') || isa(tmp,'Exfo_laser') || isa(tmp,'DC_supply')
if i == 1
value = {};

7
Classes/Warehouse_class/classes/Parameter.m Normal file
View File

@@ -0,0 +1,7 @@
classdef Parameter < StorageParameter
methods
function obj = Parameter(varargin)
obj@StorageParameter(varargin{:});
end
end
end

69
Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m
View File

@@ -5,6 +5,16 @@
wh = load([path filesep file]);
wh = wh.wh;
% fields = fieldnames(wh.parameter);
% for k = 1:numel(fields)
% oldParam = wh.parameter.(fields{k});
% % copy over the properties to your new class
% wh.parameter.(fields{k}) = StorageParameter(...
% oldParam.Name, oldParam.values);
% end
%
plotJob = struct();
width = 350;
height = 200;
@@ -17,15 +27,15 @@ plotJob.d = 0;
plotJob.sgm = 0;
plotJob.pol = "copolarized";
plotJob.p_in = 3;
plotJob.gamma = 0;
plotJob.pmd = 0;
plotJob.gamma = 0.0023;
plotJob.pmd = 0.1;
plotJob.channelspacing = 400e9;
plotJob.randzdw = 0;
plotJob.plot_ber_curve = 1;
plotJob.plot_ber_curve = 0;
plotJob.plot_3dber_curve = 0;
plotJob.plot_violin = 0;
plotJob.plot_violin = 1;
plotJob.plot_wavelength_sweep = 0;
plotJob.plot_wavelength_sweep_failure_rate = 0;
@@ -34,13 +44,20 @@ plotJob.dataStatArg = 'Lineplot with quartiles';
plotJob.plotTypeArg = 'Lines';
plotJob.displayname = 'a';
plotJob.title = 'title';
plotJob.figName = '16 Chann__';
plotJob.figName = '16 Chann';
plotJob.xAxisLabel = 'ROP per Channel in dBm';
plotJob.yAxisLabel = 'BER';
% createbercurves(wh,plotJob)
createviolinplots(wh,plotJob);
% createsweepplots(wh,plotJob);
%%
% createbercurves(wh,plotJob)
P = [3,6];
for i = 1:2
plotJob.p_in = P(i);
createviolinplots(wh,plotJob);
end
% createsweepplots(wh,plotJob);
%% 1
@@ -64,7 +81,7 @@ D = [0,0,0,3,0,0,0,3];
Sgm = [0,0,0,1,0,0,0,1];
colidx = [4,8,6];
P_launch = [0,3,6];
P_launch = [3,6];
fig = figure('Name',plotJob.figName);
fig.Position = plotJob.Position;
@@ -81,7 +98,7 @@ for idx = 1:(numRows * numCols)
plotJob.sgm = Sgm(idx);
plotJob.randzdw = 1;
for i = 1:3
for i = 1:length(P_launch)
plotJob.p_in = P_launch(i);
plotJob.color = cols(colidx(i),:);
@@ -199,6 +216,7 @@ end
%% 2
function createviolinplots(wh,plotJob)
width = 350;
height = 200;
s = 100;
@@ -208,30 +226,34 @@ cols = cbrewer2("paired",12);
numRows = 1;
numCols = 4;
plotJob.ch = 16;
plotJob.p_in = 3;
plotJob.randzdw = 0;
plotJob.randzdw = 1;
Pol = ["copolarized","copolarized","alternated","paired",];
Title = ["Co Pol.","Link Segmentation","Paired Pol. Interl.","Alternating Pol. Interl."];
D = [0,3,0,0];
Sgm = [0,1,0,0];
colidx = [2];
Len = [2];
colidx = [3];
Len = [10];
plotJob.figName = ['_Violin',num2str(plotJob.ch),' Channels; ',num2str(plotJob.channelspacing*1e-9),' GHz; ',num2str(plotJob.p_in),' dBm; randomized: ', num2str(plotJob.randzdw)];
fig = figure('Name',plotJob.figName);
fig.Position = plotJob.Position;
fig.Units = "centimeters";
fig.Position = [0 0 18 7];
fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
if isvalid(fig)
figure(fig)
% fig = get(fig);
AxesMain = fig.CurrentAxes;
hold on
% t = tiledlayout(1,4,'TileSpacing','compact','Padding','compact');
else
fig = figure('name',char(plotJob.figName));
AxesMain = gca;
hold on; grid on;
% t = tiledlayout(1,4,'TileSpacing','compact','Padding','compact');
end
t = tiledlayout(numRows,numCols,'TileSpacing','compact','Padding','compact');
for idx = 1:(numRows * numCols)
% Create subplot
%subplot(numRows, numCols, idx);
nexttile;
subplot(numRows, numCols, idx);
plotJob.pol = Pol(idx);
plotJob.d = D(idx);
@@ -295,7 +317,6 @@ annotation(fig,'textbox',...
'FontSize',8,...
'FitBoxToText','off');
copygraphics(t,'BackgroundColor','none');
% lgd = legend('$P_{\mathrm{in}}=0$ dBm','$P_{\mathrm{in}}=3$ dBm','$P_{\mathrm{in}}=6$ dBm','Interpreter','latex');
% lgd.NumColumns = 3;
% lgd.Layout.Tile = 'south';

137
Classes/Warehouse_class/functions/fwm_plots/automate_PTL_plot_new.m Normal file
View File

@@ -0,0 +1,137 @@
% Select dataset
[file, path] = uigetfile("C:\Users\Silas\Documents\MATLAB\Datensätze\Raw_Cluster_Simulations\session_februar_24\wh_mi_nacht.mat");
wh = load(fullfile(path, file));
wh = wh.wh;
%% --- Plot Settings ---
cols = cbrewer2("Paired", 12);
plotJob = struct();
plotJob.Position = [100 100 600 400];
plotJob.channelspacing = 400e9;
plotJob.ch = 16;
plotJob.d = 0;
plotJob.sgm = 0;
plotJob.gamma = 0.0023;
plotJob.pmd = 0.1;
plotJob.randzdw = 1;
plotJob.plot_ber_curve = 1;
plotJob.xAxisLabel = 'ROP per $\lambda$ [dBm]';
plotJob.yAxisLabel = 'BER';
plotJob.figName = 'avg BER_vs_Plaunch_combined';
plotJob.dataStatArg = 'Lineplot with quartiles';%'All Channels; mean(PMD Realizations)';Lineplot with quartiles
plotJob.plotTypeArg = 'Lines';
plotJob.displayname = 'bla';
% --- Parameter combinations ---
Len = [2, 10];
Pol = ["copolarized", "copolarized", "alternated", "paired"];
Title = ["CoPol","LS", "API", "PPI"];
D = [0, 3, 0, 0];
Sgm = [0, 1, 0, 0];
% Pol = ["copolarized", "copolarized"];
% Title = ["CoPol","LS"];
% D = [0, 3];
% Sgm = [0, 1];
colidx = [6,4,2,2]; % color indices for different schemes
P_launch = [0,3,6]; % input power sweep
%% --- Create Figure ---
fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
if isvalid(fig)
figure(fig)
% fig = get(fig);
AxesMain = fig.CurrentAxes;
hold on
% t = tiledlayout(1,4,'TileSpacing','compact','Padding','compact');
else
fig = figure('name',char(plotJob.figName));
AxesMain = gca;
hold on; grid on;
% t = tiledlayout(1,4,'TileSpacing','compact','Padding','compact');
end
dsa = ["AVG"];
for d = 1
plotJob.dataStatArg = dsa(d);
cnt = 1;
for l = 1:numel(Len)
subplot(1,2,cnt);
cnt = cnt+1;
for s = 1:length(P_launch)
for p = 1:numel(Title)
plotJob.l = Len(l);
plotJob.pol = Pol(p);
plotJob.d = D(p);
plotJob.sgm = Sgm(p);
% color + style per length
baseColor = cols(colidx(p), :);
if s == 1
plotJob.linestyle = '-';
elseif s == 2
plotJob.linestyle = '--';
else
plotJob.linestyle = ':';
end
if p == 1
% plotJob.linestyle = '-';
plotJob.markerstyle = 'o';
plotJob.markersize = 2;
elseif p == 2
% plotJob.linestyle = ':';
plotJob.markerstyle = 'square';
plotJob.markersize = 2;
elseif p == 3
% plotJob.linestyle = '-';
plotJob.markerstyle = 'x';
plotJob.markersize = 6;
else
% plotJob.linestyle = '-';
plotJob.markerstyle = 'diamond';
plotJob.markersize = 2;
end
% if d == 1
% plotJob.linestyle = '-';
% else
% plotJob.linestyle = ':';
% plotJob.markerstyle = 'none';
% end
plotJob.p_in = P_launch(s);
plotJob.displayname = sprintf('%s',Title(p));
plotJob.color = baseColor;% * (1 - 0.15*(p-1)); % slight shade for powers
plotCurve(wh, plotJob);
% h = findobj(gca,'Type','Line','-not','Tag','FEC');
% set(h(p),'DisplayName',sprintf('%s (%.0f km, %.0f dBm)',Title(p),Len(l),P_launch(s)));
% title(sprintf('%d km; %d Channels, \Delta f = %.0f GHz', ...
% plotJob.l, plotJob.ch, plotJob.channelspacing*1e-9));
end
end
end
end
set(gca, 'YScale', 'log');
xlabel(plotJob.xAxisLabel);
ylabel(plotJob.yAxisLabel);
% legend('Interpreter','latex','NumColumns',2,'Location','southoutside');
grid on; box on;
copygraphics(fig, 'BackgroundColor','none');

62
Classes/Warehouse_class/functions/fwm_plots/plotCurve.m
View File

@@ -6,7 +6,7 @@ fig = findall(groot, 'Type', 'figure', 'Name', plotJob.figName);
if isvalid(fig)
figure(fig)
fig = get(fig);
% fig = get(fig);
AxesMain = fig.CurrentAxes;
hold on
else
@@ -28,32 +28,39 @@ realization = wh.parameter.realization.values(1:end);
% get all xAxis values
xAxis = wh.parameter.p_out.values;
markerstyle = 'o';
linestyle = '-';
% Fetch Data from Warehouse
for xl = 1:numel(xAxis)
p_out = xAxis(xl);
if string(plotJob.dataStatArg) == "Worst"
temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
temp = removeZeros(temp);
ber(xl) = max(temp,[],'all');
ber(xl) = quantile(temp,0.9,"all");
% ber(xl) = max(temp,[],'all');
linew = 1.0;
markersz = 3;
linestyle = '-';
markersz = plotJob.markersize;
markerstyle = plotJob.markerstyle;
linestyle = plotJob.linestyle;
elseif string(plotJob.dataStatArg) == "AVG"
temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
temp = removeZeros(temp);
ber(xl) = mean(temp,'all');
linew = 1.0;
markersz = 3;
linestyle = '--';
markersz = plotJob.markersize;
markerstyle = plotJob.markerstyle;
linestyle = plotJob.linestyle;
elseif string(plotJob.dataStatArg) == "Best"
temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
temp = removeZeros(temp);
ber(xl) = min(temp,[],'all');
linew = 1.0;
markersz = 3;
linestyle = '-';
markersz = plotJob.markersize;
markerstyle = plotJob.markerstyle;
linestyle = plotJob.linestyle;
elseif string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)"
temp = wh.getStoValue('ber',plotJob.l,plotJob.d,plotJob.sgm,string(plotJob.pol),plotJob.p_in,p_out,plotJob.pmd,plotJob.gamma,realization,plotJob.ch,1310,plotJob.channelspacing,plotJob.randzdw);
@@ -63,7 +70,8 @@ for xl = 1:numel(xAxis)
linew = 1;
markersz = 1;
linestyle = '-';
markerstyle = plotJob.markerstyle;
linestyle = plotJob.linestyle;
elseif string(plotJob.dataStatArg) == "Lineplot with quartiles"
@@ -75,14 +83,16 @@ for xl = 1:numel(xAxis)
upperq(xl) = quantile(temp,0.99,"all");
lowerq(xl) = quantile(temp,0.04,"all");
if lowerq(xl) == 0
lowerq(xl) = lowerq(xl-1);
end
% upperq(xl) = 0.5*std(tmp,1,'all','omitnan');
% lowerq(xl) = 0.5*std(tmp,1,'all','omitnan');
% upperq(xl) = max(dataNoNans);
% lowerq(xl) = min(dataNoNans);
upperq(xl) = max(temp(:));
lowerq(xl) = min(temp(:));
if lowerq(xl) == 0
lowerq(xl) = 1e-8;
end
% upperq(xl) = mean(tmp,"all","omitnan") + 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp)));
% lowerq(xl) = mean(tmp,"all","omitnan") - 1.96 * (std(tmp,1,'all','omitnan')/sqrt(numel(tmp)));
@@ -149,6 +159,7 @@ elseif string(plotJob.plotTypeArg) == "Lines"
% [xAxis,ber] = interpCurve(xAxis, ber);
% end
cols = cbrewer2('RdBu',size(ber,1));
for rlz = 1:size(ber,1)
%
if (string(plotJob.dataStatArg) == "All Channels ;All PMD Realizations")||(string(plotJob.dataStatArg) == "All Channels; mean(PMD Realizations)")
@@ -160,6 +171,7 @@ elseif string(plotJob.plotTypeArg) == "Lines"
if rlz < size(ber,1)
col = cols(rlz,:);
s = plot(xAxis,ber(rlz,:),linestyle,'Marker',"none",'MarkerSize',markersz,'LineWidth',linew,'Color',col,'Parent', AxesMain,'HandleVisibility','off');
@@ -169,22 +181,25 @@ elseif string(plotJob.plotTypeArg) == "Lines"
s.DataTipTemplate.DataTipRows(1);
s.DataTipTemplate.DataTipRows(2) = [];
else
if string(plotJob.dataStatArg) == "Lineplot with quartiles"
[hl,hp] = boundedline(xAxis,ber(rlz,:),([(ber(rlz,:)-lowerq(rlz,:));upperq(rlz,:)-ber(rlz,:)]'),'-o','alpha','Color',col,'transparency', 0.1,'linewidth',0.7);
[hl,hp] = boundedline(xAxis,ber(rlz,:),([(ber(rlz,:)-lowerq(rlz,:));upperq(rlz,:)-ber(rlz,:)]'),'-o','alpha','Color',col,'transparency', 0.06,'linewidth',0.7);
hl.MarkerFaceColor = col;
hl.MarkerSize = 3;
hl.MarkerSize = 2;
set(hp,'HandleVisibility','off');
%hp.LineWidth = 1.2;
ho = outlinebounds(hl,hp);
set(ho, 'linestyle', ':', 'color', col,'Linewidth',0.5);
set(ho, 'linestyle', ':', 'color', col,'Linewidth',0.6);
set(ho,'HandleVisibility','off');
% errorbar(xAxis,ber(rlz,:),ber(rlz,:)-lowerq(rlz,:),upperq(rlz,:)-ber(rlz,:),'-o','Color',col,'linewidth',0.7);
else
s = plot(xAxis,ber(rlz,:),linestyle,'Marker',"o",'MarkerFaceColor',col,'MarkerSize',markersz,'LineWidth',linew,'Color',col,'Parent', AxesMain,'DisplayName',[plotJob.displayname]);
s = plot(xAxis,ber(rlz,:),linestyle,'Marker',markerstyle,'MarkerFaceColor',col,'MarkerSize',markersz,'LineWidth',linew,'Color',col,'Parent', AxesMain,'DisplayName',[plotJob.displayname]);
s.DataTipTemplate.Interpreter = "latex";
s.DataTipTemplate.DataTipRows(1).Label = "Ch: ";
@@ -245,9 +260,10 @@ grid minor
fontsize(AxesMain,8,"points")
% fontname(AxesMain,"Arial")
fig.Position = plotJob.Position;
fig.Units = "centimeters";
fig.Position = [2 2 8.5 7];
% fig.Position = plotJob.Position;
% fig.Units = "centimeters";
% fig.Position = [2 2 8.5 7];
set(AxesMain,'TickLabelInterpreter','latex')
@@ -257,7 +273,7 @@ set(AxesMain.Legend,'Interpreter','latex')
ylim([1e-5,0.3]);
xlim([min(xAxis),-3]);
xlim([-10,-4]);
% annotation('textbox', [0.125, 0.32, 0.1, 0.1], 'String', "FEC 3.8e-3","LineStyle","none","FontSize",8,"FontUnits","points")

2
Classes/Warehouse_class/functions/fwm_plots/plotViolin.m
View File

@@ -155,7 +155,7 @@ ylabel('Penalty in dB');
xlabel('Channel Number');
ylim([-9.3,-3]);
xlim([0,plotJob.ch+1]);
grid minor;
% grid minor;
set(gca, 'color', 'none');
legend = [];

5
Functions/Theory/dispersion_contour.m
View File

@@ -1,14 +1,15 @@
% Gitter für lambda0 und S0
lambda0_vec = linspace(1300,1320,200);
lambda0_vec = linspace(1260,1360,200);
S0_vec = linspace(0.06,0.1,200);
[Lambda0, S0] = meshgrid(lambda0_vec, S0_vec);
% Festen Betriebsparameter
lambda = 1293; % nm
L = 10; % km
L = 1; % km
% Dispersion berechnen (lineare Näherung)
D = S0 .* ( lambda - Lambda0 ) * L;
% D = (S0./4) .* ( lambda - (Lambda0.^4)./(lambda^3) ) * L;
%% 2D-Konturplot nur mit Linien und Text
figure('Color','w');

146
Functions/Theory/dispersion_contour_bandwidth_lambda.m Normal file
View File

@@ -0,0 +1,146 @@
%% ------------------------------------------------------------
% Contour plot: λ_null as function of bandwidth (f_target) and reach (L)
% ------------------------------------------------------------
% Parameters
lambda0 = 1310e-9; % [m]
S0 = 0.08; % [ps/(nm²·km)]
c = physconst('lightspeed');
% Sweep dimensions
f_targets = linspace(50e9, 120e9, 100); % [Hz] (x-axis)
L_values = linspace(0.5e3, 10e3, 100); % [m] (y-axis)
lambda_surface = zeros(numel(L_values), numel(f_targets));
Dacc_surface = zeros(numel(L_values), numel(f_targets));
% Outer loop over fiber length (since L must be scalar)
for iL = 1:numel(L_values)
L = L_values(iL);
[lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_targets, L, lambda0, S0);
lambda_vec = 2*abs(lambda0 - lambda_vec);
if 0
fprintf('\n- %d km ------------------------------------\n',L);
fprintf(' f_null [GHz] lambda [nm] Dacc [ps/nm]\n');
fprintf('----------------------------------------------\n');
fprintf('%10.1f %8.2f %+8.3f\n',[f_targets(:)/1e9, lambda_vec(:)*1e9, Dacc_vec(:)].');
fprintf('----------------------------------------------\n\n');
end
lambda_surface(iL, :) = lambda_vec; % λ for each f_target
Dacc_surface(iL, :) = Dacc_vec; % corresponding accumulated dispersion
end
% Convert for plotting
lambda_surface_nm = lambda_surface * 1e9; % [nm]
L_km = L_values / 1000; % [km]
f_GHz = f_targets / 1e9; % [GHz]
%% Contour plot
figure('Color','w');
% Define wavelength contour levels [nm]
lambda_levels = [1260:10:1290, 1290:5:1300, 1300:2.5:1310];
lambda_levels = [100:-20:50, 50:-10:30,30:-5:0];
% Contour plot
contour(f_GHz, L_km, lambda_surface_nm, lambda_levels, ...
'LineWidth', 1.5, ...
'ShowText', 'on', ...
'LabelFormat', '%.0f nm');
% Colormap and colorbar
colormap((cbrewer2('RdYlGn',100)));
colorbar;
clim([0 100]);
% Axis formatting
xlabel('Signal Bandwidth [GHz]');
ylabel('Fiber length [km]');
legend('$\Delta \lambda$')
% X-axis ticks at 56 : 16 : 150 GHz
xticks(56:8:150);
grid on; box on;
%% Optional: overlay accumulated-dispersion contours
if 0
hold on;
[CS, h] = contour(f_GHz, L_km, Dacc_surface, 10, 'k--', 'LineWidth', 0.8);
clabel(CS, h, 'Color','k', 'FontSize',8);
end
function [lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_target, L, lambda0, S0)
% lambda_for_first_null_full (stable, single-branch + validity checks)
% --------------------------------------------------------------------
% Computes the wavelength(s) at which the first IM/DD fading null
% occurs at frequency/ies f_target using the full dispersion model:
%
% D(lambda) = (S0/4)*(lambda - lambda0^4 / lambda^3)
%
% Restricted to the NORMAL-dispersion branch (λ < λ0),
% and valid only in the O-band (1260–1360 nm).
%
% Inputs:
% f_target - scalar or vector of target null frequencies [Hz]
% L - fiber length [m]
% lambda0 - zero-dispersion wavelength (ZDW) [m]
% S0 - dispersion slope at ZDW [ps/(nm²·km)]
%
% Outputs:
% lambda_vec - wavelength(s) [m] where first null occurs (clamped to O-band)
% Dacc_vec - accumulated dispersion(s) [ps/nm] (NaN if out of valid range)
% --------------------------------------------------------------------
c = physconst('lightspeed');
S0_si = S0 * 1e3; % ps/(nm²·km) -> s/(m³)
% Define O-band boundaries (in meters)
lambda_min = 1255e-9;
lambda_max = 1361e-9;
% Force column vector
f_target = f_target(:);
N = numel(f_target);
lambda_vec = NaN(N,1);
Dacc_vec = NaN(N,1);
for k = 1:N
RHS = c * 0.5 / (f_target(k)^2 * L);
% Normal-dispersion branch (λ < λ0)
fun = @(lambda) -(S0_si/4).*(lambda - (lambda0^4)./(lambda.^3)).*lambda.^2 - RHS;
% Limit the search to [λ_min, λ0)
try
lambda_sol = fzero(fun, [lambda_min, lambda0 * 0.999]);
catch
% If the zero is not within bounds, skip this point
lambda_sol = NaN;
end
% Validate solution
if isnan(lambda_sol) || lambda_sol < lambda_min || lambda_sol > lambda_max
lambda_vec(k) = NaN;
Dacc_vec(k) = NaN;
continue
end
% Compute D(lambda) and accumulated dispersion
D_lambda = (S0_si/4) * (lambda_sol - (lambda0^4)/(lambda_sol^3)) / 1e-6; % ps/(nm·km)
Dacc_val = D_lambda * (L/1000); % ps/nm
% Sanity bound on dispersion (avoid unphysical > ±100 ps/nm)
if abs(Dacc_val) > 100
lambda_vec(k) = NaN;
Dacc_vec(k) = NaN;
else
lambda_vec(k) = lambda_sol;
Dacc_vec(k) = Dacc_val;
end
end
end

38
Functions/Theory/dispersion_first_notch_10km.m Normal file
View File

@@ -0,0 +1,38 @@
%% ------------------------------------------------------------
% Plot: Maximum usable IM/DD bandwidth vs wavelength
% ------------------------------------------------------------
% Fiber and dispersion parameters
lambda0 = 1310e-9; % [m]
S0 = 0.08; % [ps/(nm²·km)]
L = 10000; % [m]
c = physconst('lightspeed');
% Wavelength range around ZDW
lambda_vec = linspace(1250e-9, 1350e-9, 200); % [m]
% Compute D(lambda) using full model
lambda_nm = lambda_vec * 1e9;
lambda0_nm = lambda0 * 1e9;
D_lambda = (S0/4) .* (lambda_nm - (lambda0_nm.^4) ./ (lambda_nm.^3)); % [ps/(nm·km)]
% Convert D to [s/m²]
D_si = D_lambda * 1e-6;
% Compute first null frequency (f₀) for each wavelength
f_null = sqrt(c*(0.5) ./ (abs(D_si).*lambda_vec.^2*L)); % [Hz]
% Plot
figure('Color','w');
plot(lambda_vec*1e9, f_null/1e9, 'LineWidth', 1.6);
grid on; box on;
xlabel('Wavelength [nm]');
ylabel('First Fading Null Frequency [GHz]');
title(sprintf('IM/DD Bandwidth Limit vs. Wavelength (L = %.1f km)', L/1000));
% Highlight useful bandwidth thresholds
yline(25, '--', '25 GHz','Color',[0.4 0.4 0.4],'LabelHorizontalAlignment','left');
yline(50, '--', '50 GHz','Color',[0.2 0.6 0.2],'LabelHorizontalAlignment','left');
yline(100,'--', '100 GHz','Color',[0.6 0.2 0.2],'LabelHorizontalAlignment','left');
legend('First fading notch (f_{null})','Location','best');

165
Functions/Theory/dispersion_power_fading.m Normal file
View File

@@ -0,0 +1,165 @@
%% Chromatic Dispersion Power Fading Demonstration
% ------------------------------------------------------------
% This script computes and visualizes power fading after
% photodiode detection caused by chromatic dispersion in IM/DD links.
%
% It also determines the wavelength λ that produces the first
% fading null at a specified RF frequency f_target using the
% full physical dispersion model:
%
% D(λ) = (S0/4) * (λ - λ0^4 / λ^3)
%
% and compares the analytic null frequency with simulation.
% ------------------------------------------------------------
% clear; close all; clc;
%% Fiber and wavelength parameters
lambda0 = 1310e-9; % Zero-dispersion wavelength (ZDW) [m]
S0 = 0.08; % Dispersion slope at ZDW [ps/(nm^2·km)]
L = 10000; % Fiber length [m]
alpha_dB = 0; % Attenuation [dB/m] (ignored here)
%% Target null frequency
f_targets = linspace(55e9,58e9,10);
f_targets = 56e9;
% f_targets = 80e9;
% Compute wavelength that gives the first null at f_target
[lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_targets, L, lambda0, S0);
% lambda_vec = 1293e-9;
fprintf('\n----------------------------------------------\n');
fprintf(' f_null [GHz] lambda [nm] Dacc [ps/nm]\n');
fprintf('----------------------------------------------\n');
fprintf('%10.1f %8.2f %+8.3f\n',[f_targets(:)/1e9, lambda_vec(:)*1e9, Dacc_vec(:)].');
fprintf('----------------------------------------------\n\n');
%% Frequency grid
f_simu = 500e9; % Simulation bandwidth [Hz]
N_freq = 500000;
faxis = linspace(-f_simu/2, f_simu/2, N_freq);
%% Derived fiber parameters
c = physconst('lightspeed');
S0_si = S0 * 1e3; % ps/(nm²·km) -> s/m³
% Convert wavelengths to nm for the D(lambda) model
lambda_nm = lambda_vec(end) * 1e9;
lambda0_nm = lambda0 * 1e9;
% Dispersion parameter [ps/(nm·km)]
D_lambda = (S0/4) * (lambda_nm - (lambda0_nm^4)/(lambda_nm^3));
% Convert to [s/m²]
D_si = D_lambda * 1e-6;
% β2 in [s²/m]
b2 = -D_si * lambda_vec(end)^2 / (2*pi*c);
%% IM/DD intensity response (simulation)
phi = 2*pi^2*b2*faxis.^2*L;
H_field_pos = exp(-1j*phi); % +f sideband
H_field_neg = exp(+1j*phi); % -f sideband
H_intensity = 0.5 * (H_field_pos + H_field_neg); % PD beating term
H_sim = abs(H_intensity);
%% Theoretical analytical IM/DD response
phi = 2*pi^2 * abs(b2) * faxis.^2 * L;
H_theoretical = abs(cos(phi));
%% Analytic first null (for verification)
f_null_analytic = sqrt(c*(0.5)/(abs(D_si)*lambda_vec(end)^2*L));
fprintf('Analytic first null from D,λ,L: %.2f GHz\n\n', f_null_analytic/1e9);
%% Plot
cols = linspecer(5);
figure('Color','w'); hold on; grid on; box on;
plot(faxis*1e-9, 10*log10(H_sim), 'DisplayName','$|H_{sim}|$ (IM/DD simulation)','Color',cols(1,:));
plot(faxis*1e-9, 10*log10(H_theoretical), 'DisplayName','|cos($\phi$)| (theory)','Color',cols(2,:),'LineStyle','--');
xline(f_targets(end)/1e9,'k:','LineWidth',1.2,'DisplayName','Target null (56 GHz)');
xline(f_null_analytic/1e9,'Color',[0.2 0.6 0.2],'LineStyle','-.','LineWidth',1.2,'DisplayName','Analytic null');
xlabel('Frequency [GHz]');
ylabel('Magnitude [dB]');
title(sprintf('Power Fading for %.2f nm, L = %.1f km',lambda_nm,L/1000));
legend('Location','best'); ylim([-30 0]);
%% Plot Bandwidth vs Lambda max
figure();
hold on;
plot(lambda_vec.*1e6,f_targets.*1e-9)
xlabel('wavelength');
ylabel('max. Bandwidth')
function [lambda_vec, Dacc_vec] = lambda_for_first_null_full(f_target, L, lambda0, S0)
% lambda_for_first_null_full (stable, single-branch + validity checks)
% --------------------------------------------------------------------
% Computes the wavelength(s) at which the first IM/DD fading null
% occurs at frequency/ies f_target using the full dispersion model:
%
% D(lambda) = (S0/4)*(lambda - lambda0^4 / lambda^3)
%
% Restricted to the NORMAL-dispersion branch (λ < λ0),
% and valid only in the O-band (1260–1360 nm).
%
% Inputs:
% f_target - scalar or vector of target null frequencies [Hz]
% L - fiber length [m]
% lambda0 - zero-dispersion wavelength (ZDW) [m]
% S0 - dispersion slope at ZDW [ps/(nm²·km)]
%
% Outputs:
% lambda_vec - wavelength(s) [m] where first null occurs (clamped to O-band)
% Dacc_vec - accumulated dispersion(s) [ps/nm] (NaN if out of valid range)
% --------------------------------------------------------------------
c = physconst('lightspeed');
S0_si = S0 * 1e3; % ps/(nm²·km) -> s/(m³)
% Define O-band boundaries (in meters)
lambda_min = 1255e-9;
lambda_max = 1361e-9;
% Force column vector
f_target = f_target(:);
N = numel(f_target);
lambda_vec = NaN(N,1);
Dacc_vec = NaN(N,1);
for k = 1:N
RHS = c * 0.5 / (f_target(k)^2 * L);
% Normal-dispersion branch (λ < λ0)
fun = @(lambda) -(S0_si/4).*(lambda - (lambda0^4)./(lambda.^3)).*lambda.^2 - RHS;
% Limit the search to [λ_min, λ0)
try
lambda_sol = fzero(fun, [lambda_min, lambda0 * 0.999]);
catch
% If the zero is not within bounds, skip this point
lambda_sol = NaN;
end
% Validate solution
if isnan(lambda_sol) || lambda_sol < lambda_min || lambda_sol > lambda_max
lambda_vec(k) = NaN;
Dacc_vec(k) = NaN;
continue
end
% Compute D(lambda) and accumulated dispersion
D_lambda = (S0_si/4) * (lambda_sol - (lambda0^4)/(lambda_sol^3)) / 1e-6; % ps/(nm·km)
Dacc_val = D_lambda * (L/1000); % ps/nm
% Sanity bound on dispersion (avoid unphysical > ±100 ps/nm)
if abs(Dacc_val) > 100
lambda_vec(k) = NaN;
Dacc_vec(k) = NaN;
else
lambda_vec(k) = lambda_sol;
Dacc_vec(k) = Dacc_val;
end
end
end

120
Functions/Theory/matched_filter_rrc.m Normal file
View File

@@ -0,0 +1,120 @@
%% Matched Filter SNR Demonstration (Correct Timing)
% clear; close all; clc;
%% Parameters
M = 4; % QPSK
numSymbols = 1e6;
sps = 25; % samples per symbol
rolloff = 0.5;
EbNo_dB = 10;
%% Generate random data
data = randi([0 M-1], numSymbols, 1);
txSym = qammod(data, M, 'UnitAveragePower', true);
%% Root Raised Cosine filters
span = 64; % filter span in symbols
rrcTx = rcosdesign(rolloff, span, sps, 'sqrt');
rrcRx = rrcTx; % matched filter
txSignal2 = ifft(fft(rrcTx).*fft(txSym));
%% Transmit filtering (includes upsampling)
txSignal = upfirdn(txSym, rrcTx, sps, 1);
%% AWGN channel
rxSignal = awgn(txSignal, EbNo_dB + 10*log10(sps), 'measured');
%% Receiver matched filter
rxFilt = conv(rxSignal, rrcRx, 'same');
%% Symbol timing (group delay compensation)
delay = span * sps / 2; % total delay per filter is span*sps/2
rxAligned = rxFilt(delay+1 : end-delay);
%% Downsample to symbol rate
rxSampled = rxAligned(1:sps:end);
%% Align lengths
L = min(length(rxSampled), length(txSym));
rxSampled = rxSampled(1:L);
txSym = txSym(1:L);
%% Decision and BER
rxSym = qamdemod(rxSampled, M, 'UnitAveragePower', true);
[~, ber] = biterr(data(1:L), rxSym);
%% Compute effective SNR
snr_meas = 10*log10(mean(abs(txSym).^2) / mean(abs(txSym - rxSampled).^2));
fprintf('Measured BER: %.3e | Effective SNR: %.2f dB\n', ber, snr_meas);
%% Eye diagrams
eyediagram(rxSignal(1:4000), 2*sps);
title('Received Signal (Before Matched Filter)');
eyediagram(rxFilt(1:4000), 2*sps);
title('After Matched Filter (RRC)');
%% --------------------------------------------------------------
%% Spectrum analysis of shaped and filtered signals
%% --------------------------------------------------------------
Fs = sps; % normalized sample rate (symbol rate = 1)
Nfft = 2^16; % FFT size for high resolution
f = (-Nfft/2:Nfft/2-1)/Nfft * Fs; % normalized frequency axis (symbol-rate units)
% Spectra
S_tx = 20*log10(abs(fftshift(fft(txSignal, Nfft)))/max(abs(fft(txSignal, Nfft))));
S_rx = 20*log10(abs(fftshift(fft(rxFilt, Nfft)))/max(abs(fft(rxFilt, Nfft))));
% Unshaped (rectangular pulse) for comparison
txRect_unf = upfirdn(txSym, ones(1, sps), sps, 1);
S_rect = 20*log10(abs(fftshift(fft(txRect_unf, Nfft)))/max(abs(fft(txRect_unf, Nfft))));
% Plot
figure('Name','Spectrum after Pulse Shaping');
plot(f, S_rect, '--', 'DisplayName','Rectangular pulse');
hold on;
plot(f, S_tx, 'LineWidth',1.4, 'DisplayName','RRC (TX)');
plot(f, S_rx, 'LineWidth',1.4, 'DisplayName','After Matched Filter');
grid on;
xlabel('Normalized frequency (× symbol rate)');
ylabel('Magnitude [dB]');
title('Spectra Before and After RRC Pulse Shaping');
legend('Location','best');
xlim([-1.5 1.5]);
ylim([-60 0]);
%% --------------------------------------------------------------
%% Visualization: RRC and Raised-Cosine Frequency Responses
%% --------------------------------------------------------------
% Frequency axis for plotting (normalized to symbol rate)
Nfft = 4096;
H_rrc = fftshift(fft(rrcTx, Nfft));
H_rc = H_rrc .* H_rrc; % cascade of TX and RX RRC = full RC
f = linspace(-0.5, 0.5, Nfft); % normalized frequency (symbol-rate units)
figure('Name','Raised Cosine Filter Characteristics');
subplot(2,1,1);
plot(f, 20*log10(abs(H_rrc)/max(abs(H_rrc))), 'LineWidth', 1.5);
hold on;
plot(f, 20*log10(abs(H_rc)/max(abs(H_rc))), '--', 'LineWidth', 1.5);
grid on;
xlabel('Normalized frequency (× symbol rate)');
ylabel('Magnitude [dB]');
title(sprintf('RRC (rolloff = %.2f) and Full RC Spectrum', rolloff));
legend('Root Raised Cosine','Raised Cosine (TX×RX)','Location','best');
ylim([-60 5]);
subplot(2,1,2);
t = (-span*sps/2 : span*sps/2) / sps; % time axis in symbol durations
plot(t, rrcTx, 'LineWidth', 1.5);
grid on;
xlabel('Time [symbols]');
ylabel('Amplitude');
title('RRC Impulse Response');

24
projects/ECOC_2025/plots_from_database/plot_mpi_trial.m
View File

@@ -43,12 +43,13 @@
db = DBHandler("type","mysql","dataBase",'labor');
fp = QueryFilter();
% fp.where('Runs', 'loop_id','EQUALS', 209);
% fp.where('Runs', 'sir','EQUALS', 21);
fp.where('Runs', 'pam_level','EQUALS', 4);
fn = [db.getTableFieldNames('Runs');db.getTableFieldNames('Results')];
% fp.where('mpi_superview', 'loop_id','EQUALS', 209);
fp.where('mpi_superview', 'symbolrate','EQUALS', 112e9);
fp.where('mpi_superview', 'pam_level','EQUALS', 4);
fn = [db.getTableFieldNames('mpi_superview')];
[dataTable,sql_query] = db.queryDB(fp,fn);
%%
dataTable_clean = dataTable;
@@ -74,7 +75,7 @@ for int_len = [0,50,300,1000]
mode = 4;
for mode = [1,4]
for mode = [1,2]
hold on;
dataTable = dataTable_clean;
@@ -132,7 +133,7 @@ for int_len = [0,50,300,1000]
method = 'ideal dc tracking';
end
dataTable(dataTable.eq_id==0,:) = [];
% dataTable(dataTable.eq_id==0,:) = [];
dataTable(dataTable.equalizer_structure~=1,:) = [];
% Modify values in 'interference_path_length' where the condition is met
@@ -190,6 +191,10 @@ for int_len = [0,50,300,1000]
dataTableGrpd_min = groupIt(fixedVars, dataTable, @min);
dataTableGrpd_max = groupIt(fixedVars, dataTable, @max);
% dataTableGrpd_mean(dataTableGrpd_mean.nRows<50,:) = [];
% dataTableGrpd_min(dataTableGrpd_min.nRows<50,:) = [];
% dataTableGrpd_max(dataTableGrpd_max.nRows<50,:) = [];
% Create a new figure
hold on
@@ -244,6 +249,8 @@ for int_len = [0,50,300,1000]
% Hide patch (shaded area) from legend
set(hp, 'HandleVisibility', 'off','LineStyle',':','LineWidth',0.5,'Marker','none');
% Fit a 4th-order polynomial to log10(BER)
p = polyfit(x_values, log10(y_mean), 3); % 4 is fitting order, adjust as needed
@@ -294,9 +301,10 @@ for int_len = [0,50,300,1000]
'LineWidth', 0.5, 'HandleVisibility', 'off', 'DisplayName', string(dispname));
pair_one = {'Run ID', dataTable.run_id(loopFiltSingle, :)};
pair_two = {'Rate', dataTable.bitrate(loopFiltSingle, :) * 1e-9};
pair_two = {'Baud', dataTable.symbolrate(loopFiltSingle, :) * 1e-9};
pair_three = {'PD in', round(dataTable.power_pd_in(loopFiltSingle, :), 2)};
addDatatips(sc, pair_one, pair_two, pair_three);
pair_four = {'#bits', round(dataTable.numBits(loopFiltSingle, :), 2)};
addDatatips(sc, pair_one, pair_two, pair_three,pair_four);
end
end