From 3ea64399478515347548640add7d75993f84e10c Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Fri, 17 Oct 2025 11:50:27 +0200 Subject: [PATCH] - bring FWM plots back to life! - some dispersion plots with the great help of chatGPT :-D --- Classes/DataBaseHandler/DBHandler.m | 5 +- Classes/Warehouse_class/classes/DataStorage.m | 2 +- Classes/Warehouse_class/classes/Parameter.m | 7 + .../functions/fwm_plots/automate_JLT_plots.m | 69 +++++--- .../fwm_plots/automate_PTL_plot_new.m | 137 +++++++++++++++ .../functions/fwm_plots/plotCurve.m | 62 ++++--- .../functions/fwm_plots/plotViolin.m | 2 +- Functions/Theory/dispersion_contour.m | 5 +- .../dispersion_contour_bandwidth_lambda.m | 146 ++++++++++++++++ .../Theory/dispersion_first_notch_10km.m | 38 ++++ Functions/Theory/dispersion_power_fading.m | 165 ++++++++++++++++++ Functions/Theory/matched_filter_rrc.m | 120 +++++++++++++ .../plots_from_database/plot_mpi_trial.m | 24 ++- 13 files changed, 720 insertions(+), 62 deletions(-) create mode 100644 Classes/Warehouse_class/classes/Parameter.m create mode 100644 Classes/Warehouse_class/functions/fwm_plots/automate_PTL_plot_new.m create mode 100644 Functions/Theory/dispersion_contour_bandwidth_lambda.m create mode 100644 Functions/Theory/dispersion_first_notch_10km.m create mode 100644 Functions/Theory/dispersion_power_fading.m create mode 100644 Functions/Theory/matched_filter_rrc.m diff --git a/Classes/DataBaseHandler/DBHandler.m b/Classes/DataBaseHandler/DBHandler.m index f1ca74e..5569746 100644 --- a/Classes/DataBaseHandler/DBHandler.m +++ b/Classes/DataBaseHandler/DBHandler.m @@ -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 diff --git a/Classes/Warehouse_class/classes/DataStorage.m b/Classes/Warehouse_class/classes/DataStorage.m index 5fd288d..93a7cf4 100644 --- a/Classes/Warehouse_class/classes/DataStorage.m +++ b/Classes/Warehouse_class/classes/DataStorage.m @@ -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 = {}; diff --git a/Classes/Warehouse_class/classes/Parameter.m b/Classes/Warehouse_class/classes/Parameter.m new file mode 100644 index 0000000..1673241 --- /dev/null +++ b/Classes/Warehouse_class/classes/Parameter.m @@ -0,0 +1,7 @@ +classdef Parameter < StorageParameter + methods + function obj = Parameter(varargin) + obj@StorageParameter(varargin{:}); + end + end +end diff --git a/Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m b/Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m index 68b77b6..cc2c8d6 100644 --- a/Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m +++ b/Classes/Warehouse_class/functions/fwm_plots/automate_JLT_plots.m @@ -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'; diff --git a/Classes/Warehouse_class/functions/fwm_plots/automate_PTL_plot_new.m b/Classes/Warehouse_class/functions/fwm_plots/automate_PTL_plot_new.m new file mode 100644 index 0000000..f92f86f --- /dev/null +++ b/Classes/Warehouse_class/functions/fwm_plots/automate_PTL_plot_new.m @@ -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'); diff --git a/Classes/Warehouse_class/functions/fwm_plots/plotCurve.m b/Classes/Warehouse_class/functions/fwm_plots/plotCurve.m index a7134c9..e06e7e5 100644 --- a/Classes/Warehouse_class/functions/fwm_plots/plotCurve.m +++ b/Classes/Warehouse_class/functions/fwm_plots/plotCurve.m @@ -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") diff --git a/Classes/Warehouse_class/functions/fwm_plots/plotViolin.m b/Classes/Warehouse_class/functions/fwm_plots/plotViolin.m index 2453744..93cfc10 100644 --- a/Classes/Warehouse_class/functions/fwm_plots/plotViolin.m +++ b/Classes/Warehouse_class/functions/fwm_plots/plotViolin.m @@ -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 = []; diff --git a/Functions/Theory/dispersion_contour.m b/Functions/Theory/dispersion_contour.m index 3ed54aa..324dc2d 100644 --- a/Functions/Theory/dispersion_contour.m +++ b/Functions/Theory/dispersion_contour.m @@ -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'); diff --git a/Functions/Theory/dispersion_contour_bandwidth_lambda.m b/Functions/Theory/dispersion_contour_bandwidth_lambda.m new file mode 100644 index 0000000..1f46f0e --- /dev/null +++ b/Functions/Theory/dispersion_contour_bandwidth_lambda.m @@ -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 diff --git a/Functions/Theory/dispersion_first_notch_10km.m b/Functions/Theory/dispersion_first_notch_10km.m new file mode 100644 index 0000000..96a2149 --- /dev/null +++ b/Functions/Theory/dispersion_first_notch_10km.m @@ -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'); diff --git a/Functions/Theory/dispersion_power_fading.m b/Functions/Theory/dispersion_power_fading.m new file mode 100644 index 0000000..6ed8fab --- /dev/null +++ b/Functions/Theory/dispersion_power_fading.m @@ -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 \ No newline at end of file diff --git a/Functions/Theory/matched_filter_rrc.m b/Functions/Theory/matched_filter_rrc.m new file mode 100644 index 0000000..577ec10 --- /dev/null +++ b/Functions/Theory/matched_filter_rrc.m @@ -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'); diff --git a/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m index 68ac774..276c7d2 100644 --- a/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m +++ b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m @@ -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