Compare commits
1 Commits
cae81c0dae
...
revert-798
| Author | SHA1 | Date | |
|---|---|---|---|
|
9093fb2452 |
@@ -172,9 +172,11 @@ classdef Signal
|
||||
|
||||
hold on;
|
||||
if isempty(options.color)
|
||||
plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
|
||||
% plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
|
||||
plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1);
|
||||
else
|
||||
plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color);
|
||||
% plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color);
|
||||
plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Color',options.color);
|
||||
end
|
||||
% 2 c)
|
||||
% - xlabel if not already here: time in readable format (1 ms and not 1e-3 s)
|
||||
|
||||
@@ -24,8 +24,6 @@ classdef DP_Fiber
|
||||
SS_dzmax % [m] max dz (adaptive SSFM)
|
||||
SS_dzmin % [m] min dz (adaptive SSFM)
|
||||
n_waveplates % number of PMD waveplates
|
||||
useGPU % GPU acceleration: true, false, or 'auto' (default)
|
||||
useSingle % Use single precision on GPU (default: false)
|
||||
|
||||
% ---- Internal state (persistent between calls) ----
|
||||
state % struct mirroring legacy 'state'
|
||||
@@ -58,8 +56,6 @@ classdef DP_Fiber
|
||||
options.SS_dzmax = 2e4 % m
|
||||
options.SS_dzmin = 100 % m
|
||||
options.n_waveplates = 100
|
||||
options.useGPU = 'auto' % 'auto', true, or false
|
||||
options.useSingle = false % single precision GPU
|
||||
end
|
||||
|
||||
% Copy provided options into properties
|
||||
@@ -212,7 +208,7 @@ classdef DP_Fiber
|
||||
% Frequency-dependent PMD phase term (legacy form)
|
||||
st.brf.db0 = (R.rand(st.wave_plates,1)*2*pi - pi) * brf_multiplier;
|
||||
st.brf.db1 = sqrt(3*pi/8)*(st.dgd/obj.fa)/st.wave_plates .* st.omega;
|
||||
st.brf.simdgd = 0;
|
||||
st.brf.simdgd = 0;
|
||||
% cumsum used in legacy only for debug; keep compatibility variable:
|
||||
~cumsum(st.brf.db0); % no-op to mirror legacy path
|
||||
|
||||
@@ -232,18 +228,7 @@ classdef DP_Fiber
|
||||
x_in = signal_in(:,1).';
|
||||
y_in = signal_in(:,2).';
|
||||
|
||||
% Determine GPU usage
|
||||
if ischar(obj.useGPU) || isstring(obj.useGPU)
|
||||
if strcmpi(obj.useGPU, 'auto')
|
||||
gpuFlag = []; % Let CNLSE_plain auto-detect
|
||||
else
|
||||
error('DP_Fiber:InvalidGPU', 'useGPU must be true, false, or ''auto''');
|
||||
end
|
||||
else
|
||||
gpuFlag = logical(obj.useGPU);
|
||||
end
|
||||
|
||||
[x_out, y_out, obj.state] = CNLSE_plain(x_in, y_in, obj.state, gpuFlag, obj.useSingle);
|
||||
[x_out, y_out, obj.state] = CNLSE_plain(x_in, y_in, obj.state);
|
||||
|
||||
obj.state.propagated_length = obj.state.propagated_length + obj.state.L;
|
||||
|
||||
|
||||
@@ -42,7 +42,7 @@ classdef Optical_Demultiplex < handle
|
||||
|
||||
function signalclasses_out = process(obj, signalclass_in)
|
||||
|
||||
% ---- Infer wavelength: either given or from input total signal
|
||||
% ---- Infer wavelength: either given or from input total signal
|
||||
if isempty(obj.wavelengthplan)
|
||||
obj.wavelengthplan = signalclass_in.lambda; %meter
|
||||
else
|
||||
@@ -81,7 +81,7 @@ classdef Optical_Demultiplex < handle
|
||||
obj
|
||||
signal_in
|
||||
end
|
||||
|
||||
|
||||
w = obj.fs_out ./ obj.fs_in ;
|
||||
blocklen_in = length(signal_in);
|
||||
blocklen_out = w*blocklen_in;
|
||||
@@ -119,31 +119,30 @@ classdef Optical_Demultiplex < handle
|
||||
N = size(lo,1);
|
||||
C = size(lo,2);
|
||||
|
||||
% ---- VECTORIZED: Process all channels in parallel ----
|
||||
% Batched FFT operates on each column simultaneously on GPU
|
||||
x_envelopes = zeros(N, C, 'like', signal_in);
|
||||
y_envelopes = zeros(N, C, 'like', signal_in);
|
||||
|
||||
% Extract polarization signals
|
||||
s1 = signal_in(:,1); % X polarization [N×1]
|
||||
s2 = signal_in(:,2); % Y polarization [N×1]
|
||||
s1 = signal_in(:,1);
|
||||
s2 = signal_in(:,2);
|
||||
|
||||
% Broadcast signal to all channels and multiply with LO
|
||||
% s1, s2 are [N×1], lo is [N×C] → result is [N×C]
|
||||
x_mixed = att .* s1 .* lo; % [N×C]
|
||||
y_mixed = att .* s2 .* lo; % [N×C]
|
||||
|
||||
% Batched FFT: each column computed in parallel
|
||||
x_freq = fft(x_mixed); % [N×C]
|
||||
y_freq = fft(y_mixed); % [N×C]
|
||||
|
||||
% Apply filter (H is [N×1], broadcasts across columns)
|
||||
x_filtered = x_freq .* H; % [N×C]
|
||||
y_filtered = y_freq .* H; % [N×C]
|
||||
|
||||
% Batched IFFT
|
||||
x_envelopes = ifft(x_filtered); % [N×C]
|
||||
y_envelopes = ifft(y_filtered); % [N×C]
|
||||
% Reusable work buffers (avoid reallocations)
|
||||
wrk_time = zeros(N,1, 'like', signal_in);
|
||||
wrk_freq = zeros(N,1, 'like', signal_in);
|
||||
|
||||
for c = 1:C
|
||||
% ---- X branch ----
|
||||
wrk_time(:) = att .* s1 .* lo(:,c); % N×1
|
||||
wrk_freq(:) = fft(wrk_time); % N×1
|
||||
wrk_freq(:) = wrk_freq .* H; % N×1
|
||||
x_envelopes(:,c) = ifft(wrk_freq); % N×1
|
||||
|
||||
% ---- Y branch ----
|
||||
wrk_time(:) = att .* s2 .* lo(:,c);
|
||||
wrk_freq(:) = fft(wrk_time);
|
||||
wrk_freq(:) = wrk_freq .* H;
|
||||
y_envelopes(:,c) = ifft(wrk_freq);
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
@@ -1,10 +1,10 @@
|
||||
classdef Optical_Multiplex < handle
|
||||
% Takes a cell array of signals
|
||||
% returns a total field signal
|
||||
% WDM spacing is given in wavelength plan OR via delta_F
|
||||
|
||||
% returns a total field signal
|
||||
% WDM spacing is given in wavelength plan OR via delta_F
|
||||
|
||||
% The grid is stored in the output signal -> the demux will ideally
|
||||
% look this up and use this as the demux frequencies...
|
||||
% look this up and use this as the demux frequencies...
|
||||
|
||||
% signal_cell = {Opt_sig_1, Opt_sig_2};
|
||||
% Opt_sig_wdm = Optical_Multiplex("fs_in",Opt_sig.fs,"fs_out",4*Opt_sig.fs,...
|
||||
@@ -117,7 +117,7 @@ classdef Optical_Multiplex < handle
|
||||
% adapt frequency shifts to match the FFT grid! Find nearest grid point
|
||||
[glitch(o),pos] = min(abs( freqaxis-obj.df_T(o) ));
|
||||
obj.df_T(o) = freqaxis(pos);
|
||||
|
||||
|
||||
polrots = [polrots, data_in{o}.polrot];
|
||||
end
|
||||
|
||||
@@ -139,36 +139,25 @@ classdef Optical_Multiplex < handle
|
||||
|
||||
x_envelopes = NaN([blocklen_out N]);
|
||||
y_envelopes = x_envelopes;
|
||||
|
||||
% ---- OPTIMIZED: Pre-compute all LO phases as [blocklen_out × N] matrix ----
|
||||
time_idx = (0:blocklen_out-1).'; % [blocklen_out × 1]
|
||||
lo_phases = mod(2*pi * time_idx * obj.df_T / obj.fs_out, 2*pi); % [blocklen_out × N]
|
||||
lo_all = cos(lo_phases) + 1i*sin(lo_phases); % [blocklen_out × N]
|
||||
|
||||
% Collect all resampled signals first (still requires loop due to cell array)
|
||||
x_signals = zeros(blocklen_out, N);
|
||||
y_signals = zeros(blocklen_out, N);
|
||||
|
||||
|
||||
for o = 1:N
|
||||
data_in_resampled = data_in{o}.resample("fs_out", obj.fs_out);
|
||||
x_signals(:, o) = data_in_resampled.signal(:, 1);
|
||||
y_signals(:, o) = data_in_resampled.signal(:, 2);
|
||||
|
||||
pha = mod(2*pi*(0:blocklen_out-1)*obj.df_T(o)/obj.fs_out,2*pi).';
|
||||
lo = cos(pha)+1i*sin(pha);
|
||||
data_in_resampled = data_in{o}.resample("fs_out",obj.fs_out);
|
||||
|
||||
res_env = ifft(fft(data_in_resampled.signal(:,1)).*H);
|
||||
x_envelopes(:,o) = att.*res_env.*lo;
|
||||
|
||||
res_env = ifft(fft(data_in_resampled.signal(:,2)).*H);
|
||||
y_envelopes(:,o) = att.*res_env.*lo;
|
||||
|
||||
end
|
||||
|
||||
% ---- VECTORIZED: Batched FFT/IFFT for all channels ----
|
||||
% Apply filter to all channels at once
|
||||
x_filtered = ifft(fft(x_signals) .* H); % [blocklen_out × N]
|
||||
y_filtered = ifft(fft(y_signals) .* H); % [blocklen_out × N]
|
||||
|
||||
% Apply attenuation and LO shift to all channels
|
||||
x_envelopes = att .* x_filtered .* lo_all; % [blocklen_out × N]
|
||||
y_envelopes = att .* y_filtered .* lo_all; % [blocklen_out × N]
|
||||
|
||||
data_out = data_in_resampled;
|
||||
data_out.signal = [sum(x_envelopes,2), sum(y_envelopes,2)];
|
||||
data_out.lambda = obj.lambda_T;
|
||||
data_out.polrot = polrots;
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
@@ -1,72 +1,41 @@
|
||||
|
||||
function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state,useGPU,useSingle)
|
||||
function [opt_out_x,opt_out_y,state] = CNLSE_plain(opt_in_x,opt_in_y,state)
|
||||
|
||||
|
||||
% GPU auto-detection if not specified
|
||||
if nargin < 4 || isempty(useGPU)
|
||||
useGPU = canUseGPU();
|
||||
end
|
||||
if nargin < 5 || isempty(useSingle)
|
||||
useSingle = false;
|
||||
end
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% pre calculations
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% pre calculations
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
state.common_beta=struct('X',0,'Y',0);
|
||||
state.common_beta=struct('X',0,'Y',0);
|
||||
|
||||
for n=1:2
|
||||
% get current polarization name and contrary one
|
||||
curPol = state.polNames{n};
|
||||
for n=1:2
|
||||
% get current polarization name and contrary one
|
||||
curPol = state.polNames{n};
|
||||
|
||||
% extend linear transfer function depending on beta values for the current polarization
|
||||
% Was ist der Sinn dieser komischen beta notation? zB. state.beta.X = [0.3142 0 -9.1105e-28 5.1068e-41]
|
||||
for n_beta = 1:length(state.beta.(curPol))
|
||||
state.common_beta.(curPol) = state.common_beta.(curPol) + state.beta.(curPol)(n_beta) * (state.omega).^(n_beta-1) / factorial(n_beta-1);
|
||||
% extend linear transfer function depending on beta values for the current polarization
|
||||
% Was ist der Sinn dieser komischen beta notation? zB. state.beta.X = [0.3142 0 -9.1105e-28 5.1068e-41]
|
||||
for n_beta = 1:length(state.beta.(curPol))
|
||||
state.common_beta.(curPol) = state.common_beta.(curPol) + state.beta.(curPol)(n_beta) * (state.omega).^(n_beta-1) / factorial(n_beta-1);
|
||||
end
|
||||
|
||||
%opt_out_struct.(curPol)=opt_in_struct.(curPol).envelope;
|
||||
end
|
||||
|
||||
%opt_out_struct.(curPol)=opt_in_struct.(curPol).envelope;
|
||||
end
|
||||
beta_const = state.beta.('X')(1);
|
||||
beta_1 = state.beta.('X')(2);
|
||||
beta_2 = state.beta.('X')(3);
|
||||
beta_3 = state.beta.('X')(4);
|
||||
deltaomega = state.omega;
|
||||
beta_x = beta_const + beta_1 * deltaomega + 1/2 * beta_2 * deltaomega.^2 + 1/6 *beta_3 * deltaomega.^3;
|
||||
|
||||
% opt_in_x = gpuArray(opt_in_x);
|
||||
% opt_in_y = gpuArray(opt_in_y);
|
||||
|
||||
beta_const = state.beta.('X')(1);
|
||||
beta_1 = state.beta.('X')(2);
|
||||
beta_2 = state.beta.('X')(3);
|
||||
beta_3 = state.beta.('X')(4);
|
||||
deltaomega = state.omega;
|
||||
beta_x = beta_const + beta_1 * deltaomega + 1/2 * beta_2 * deltaomega.^2 + 1/6 *beta_3 * deltaomega.^3;
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% GPU Transfer (if enabled)
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
if useGPU
|
||||
% Convert to single precision if requested (faster on most GPUs)
|
||||
if useSingle
|
||||
opt_in_x = gpuArray(single(opt_in_x));
|
||||
opt_in_y = gpuArray(single(opt_in_y));
|
||||
state.common_beta.X = gpuArray(single(state.common_beta.X));
|
||||
state.common_beta.Y = gpuArray(single(state.common_beta.Y));
|
||||
state.brf.db1 = gpuArray(single(state.brf.db1));
|
||||
state.brf.db0 = gpuArray(single(state.brf.db0));
|
||||
for k = 1:numel(state.brf.matR)
|
||||
state.brf.matR{k} = gpuArray(single(state.brf.matR{k}));
|
||||
end
|
||||
else
|
||||
opt_in_x = gpuArray(opt_in_x);
|
||||
opt_in_y = gpuArray(opt_in_y);
|
||||
state.common_beta.X = gpuArray(state.common_beta.X);
|
||||
state.common_beta.Y = gpuArray(state.common_beta.Y);
|
||||
state.brf.db1 = gpuArray(state.brf.db1);
|
||||
state.brf.db0 = gpuArray(state.brf.db0);
|
||||
for k = 1:numel(state.brf.matR)
|
||||
state.brf.matR{k} = gpuArray(state.brf.matR{k});
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Split Step Method
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% Split Step Method
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
% [opt_out_x,opt_out_y] = split_step_loop(state.L,opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin,...
|
||||
% state.lin_z_test,state.corr_length,state.n_plates_done,state.missing_dz,state.brf,state.common_beta,....
|
||||
@@ -75,33 +44,35 @@ end
|
||||
% [opt_out_x,opt_out_y] = split_step_loop_mex(state.L,opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin,...
|
||||
% state.lin_z_test,state.corr_length,state.n_plates_done,state.missing_dz,state.brf,state.common_beta,....
|
||||
% state.chi,state.manakov,state.beat_len);
|
||||
|
||||
% get nonlinear step size
|
||||
[state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin);
|
||||
%[state.dz] = getNLstepsize_original(state,opt_out_struct);
|
||||
|
||||
state.n_step = 0;
|
||||
state.z_prop = 0;
|
||||
state.test_dz = [];
|
||||
state.powers = [];
|
||||
|
||||
% get nonlinear step size
|
||||
[state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin);
|
||||
%[state.dz] = getNLstepsize_original(state,opt_out_struct);
|
||||
tic
|
||||
|
||||
state.n_step = 0;
|
||||
state.z_prop = 0;
|
||||
state.test_dz = [];
|
||||
state.powers = [];
|
||||
while state.z_prop < state.L
|
||||
|
||||
while state.z_prop < state.L
|
||||
% reduce step length (dz) if we are to overshoot the fiber length
|
||||
% (L) in the next step
|
||||
if state.z_prop + state.dz > state.L
|
||||
state.dz = state.L - state.z_prop;
|
||||
end
|
||||
|
||||
% update step number (n)
|
||||
state.n_step=state.n_step+1;
|
||||
|
||||
% reduce step length (dz) if we are to overshoot the fiber length
|
||||
% (L) in the next step
|
||||
if state.z_prop + state.dz > state.L
|
||||
state.dz = state.L - state.z_prop;
|
||||
end
|
||||
% append current step length to logbook (dzs)
|
||||
state.dzs(state.n_step)=state.dz;
|
||||
|
||||
% update step number (n)
|
||||
state.n_step=state.n_step+1;
|
||||
|
||||
% append current step length to logbook (dzs)
|
||||
state.dzs(state.n_step)=state.dz;
|
||||
|
||||
|
||||
% half linear step
|
||||
[opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,...
|
||||
|
||||
% half linear step
|
||||
[opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,...
|
||||
state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,...
|
||||
state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,...
|
||||
state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X]...
|
||||
@@ -111,12 +82,12 @@ while state.z_prop < state.L
|
||||
state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,...
|
||||
state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X);
|
||||
|
||||
% complete nonlinear step
|
||||
|
||||
[opt_in_x,opt_in_y] = nl_step(opt_in_x,opt_in_y, state.dz, state.gamma, state.chi, state.manakov, state.beat_len ,state.alpha_lin.X, state.alpha_lin.Y);
|
||||
|
||||
% half linear step
|
||||
[opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,...
|
||||
% complete nonlinear step
|
||||
|
||||
[opt_in_x,opt_in_y] = nl_step(opt_in_x,opt_in_y, state.dz, state.gamma, state.chi, state.manakov, state.beat_len ,state.alpha_lin.X, state.alpha_lin.Y);
|
||||
|
||||
% half linear step
|
||||
[opt_in_x,opt_in_y,state.z_prop,state.lin_z_test,...
|
||||
state.corr_length,state.n_plates_done,state.missing_dz,state.n_step,...
|
||||
state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,...
|
||||
state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X]...
|
||||
@@ -126,34 +97,24 @@ while state.z_prop < state.L
|
||||
state.test_plates,state.test_plate_numbers,state.brf,state.common_beta.X,...
|
||||
state.common_beta.Y,state.alpha_lin.X,state.alpha_lin.X);
|
||||
|
||||
% get nonlinear step size
|
||||
[state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin);
|
||||
%[state.dz] = getNLstepsize_original(state,opt_out_struct);
|
||||
% get nonlinear step size
|
||||
[state.dz] = getNLstepsize(opt_in_x,opt_in_y,state.gamma,state.SS_dzmin,state.SS_dzmax,state.SS_dphimax,state.alpha_lin);
|
||||
%[state.dz] = getNLstepsize_original(state,opt_out_struct);
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
% GPU Gather (if enabled)
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
if useGPU
|
||||
opt_out_x = gather(opt_in_x);
|
||||
opt_out_y = gather(opt_in_y);
|
||||
|
||||
% Gather state arrays back to CPU
|
||||
state.common_beta.X = gather(state.common_beta.X);
|
||||
state.common_beta.Y = gather(state.common_beta.Y);
|
||||
state.brf.db1 = gather(state.brf.db1);
|
||||
state.brf.db0 = gather(state.brf.db0);
|
||||
for k = 1:numel(state.brf.matR)
|
||||
state.brf.matR{k} = gather(state.brf.matR{k});
|
||||
|
||||
end
|
||||
else
|
||||
opt_out_x = opt_in_x;
|
||||
opt_out_y = opt_in_y;
|
||||
end
|
||||
|
||||
toc
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
opt_out_x = (opt_in_x);
|
||||
opt_out_y = (opt_in_y);
|
||||
|
||||
% opt_out_x = gather(opt_in_x);
|
||||
% opt_out_y = gather(opt_in_y);
|
||||
|
||||
end
|
||||
@@ -1,20 +0,0 @@
|
||||
function canUse = canUseGPU()
|
||||
%CANUSEGPU Check if a compatible GPU is available for parallel computing
|
||||
% Returns true if MATLAB Parallel Computing Toolbox is available and
|
||||
% a CUDA-capable GPU is detected.
|
||||
|
||||
canUse = false;
|
||||
|
||||
% Check if Parallel Computing Toolbox is installed
|
||||
if ~license('test', 'Distrib_Computing_Toolbox')
|
||||
return;
|
||||
end
|
||||
|
||||
% Check for GPU device
|
||||
try
|
||||
gpu = gpuDevice();
|
||||
canUse = gpu.DeviceSupported;
|
||||
catch
|
||||
canUse = false;
|
||||
end
|
||||
end
|
||||
@@ -1,25 +1,25 @@
|
||||
|
||||
function [rDZ] = getNLstepsize(ux,uy,gamma,dzmin,dzmax,dphimax,alpha_lin)
|
||||
|
||||
% gather() handles gpuArray inputs - ensures scalar is on CPU for comparisons
|
||||
maxPow = gather(max(gamma.*max(real(ux).^2+imag(ux).^2+real(uy).^2+imag(uy).^2)));
|
||||
|
||||
Leff = dphimax/maxPow;
|
||||
alpha_lin = max([alpha_lin.X alpha_lin.Y]);
|
||||
nl_att_len_ratio = alpha_lin*Leff;
|
||||
|
||||
if nl_att_len_ratio >= 1
|
||||
rDZ = dzmax;
|
||||
else
|
||||
if alpha_lin == 0
|
||||
step = Leff;
|
||||
|
||||
maxPow = max(gamma.*max(real(ux).^2+imag(ux).^2+real(uy).^2+imag(uy).^2));
|
||||
|
||||
Leff = dphimax/maxPow;
|
||||
alpha_lin = max([alpha_lin.X alpha_lin.Y]);
|
||||
nl_att_len_ratio = alpha_lin*Leff;
|
||||
|
||||
if nl_att_len_ratio >= 1
|
||||
rDZ = dzmax;
|
||||
else
|
||||
%effective length?
|
||||
step = -1/alpha_lin*log(1-nl_att_len_ratio);
|
||||
end
|
||||
|
||||
rDZ = min([step dzmax]);
|
||||
rDZ = max([rDZ dzmin]);
|
||||
end
|
||||
|
||||
if alpha_lin == 0
|
||||
step = Leff;
|
||||
else
|
||||
%effective length?
|
||||
step = -1/alpha_lin*log(1-nl_att_len_ratio);
|
||||
end
|
||||
|
||||
rDZ = min([step dzmax]);
|
||||
rDZ = max([rDZ dzmin]);
|
||||
end
|
||||
|
||||
end
|
||||
@@ -22,37 +22,37 @@ n_plates_left = n_plates - n_plates_done;
|
||||
|
||||
% compute last plate size ( if it fits, it should be 0)
|
||||
if missing_dz > aStepSize
|
||||
|
||||
|
||||
last_plate = aStepSize;
|
||||
missing_dz = missing_dz-aStepSize;
|
||||
plate_sizes = last_plate;
|
||||
plate_numbers = n_plates;
|
||||
|
||||
|
||||
else
|
||||
|
||||
|
||||
last_plate = aStepSize - missing_dz - (n_plates_left-1)*corr_length;
|
||||
|
||||
|
||||
if missing_dz == 0
|
||||
missing_dz = [];
|
||||
end
|
||||
|
||||
|
||||
%build vector of plate lengths with missing plate part from prev.
|
||||
%iterartion , then some normal plates and finally a fraction of a plate
|
||||
%to fit into the step length
|
||||
plate_sizes = [missing_dz corr_length*ones(1,n_plates_left-1) last_plate];
|
||||
|
||||
|
||||
if n_plates_done == 0
|
||||
plate_numbers =[(n_plates_done+1):(n_plates-1) n_plates];
|
||||
else
|
||||
plate_numbers = [n_plates_done (n_plates_done+1):(n_plates-1) n_plates]; % not wrking yet
|
||||
end
|
||||
|
||||
%remember for next step
|
||||
|
||||
%remember for next step
|
||||
missing_dz = corr_length - last_plate;
|
||||
|
||||
|
||||
end
|
||||
|
||||
plate_steps = repmat(n_step,1,length(plate_sizes)); %#ok<NASGU>
|
||||
plate_steps = repmat(n_step,1,length(plate_sizes));
|
||||
|
||||
%
|
||||
%figure;stem(plate_sizes);
|
||||
@@ -66,39 +66,50 @@ test_plate_numbers = [test_plate_numbers, plate_numbers];
|
||||
opt_x=fft(opt_x);
|
||||
opt_y=fft(opt_y);
|
||||
|
||||
% Note: db1, db0, common_beta_x, common_beta_y are already on GPU when
|
||||
% useGPU=true (transferred in CNLSE_plain)
|
||||
db1 = brf.db1;
|
||||
db0 = brf.db0;
|
||||
% db1 = gpuArray(brf.db1);
|
||||
% db0 = gpuArray(brf.db0);
|
||||
% common_beta_x = gpuArray(common_beta_x);
|
||||
% common_beta_y = gpuArray(common_beta_y);
|
||||
|
||||
db1 = (brf.db1);
|
||||
db0 = (brf.db0);
|
||||
common_beta_x = (common_beta_x);
|
||||
common_beta_y = (common_beta_y);
|
||||
|
||||
% process every waveplate with given sizes in plate_sizes
|
||||
for n=1:length(plate_sizes)
|
||||
dz = plate_sizes(n);
|
||||
|
||||
|
||||
% figure(87);subplot(2,1,1);plot(real(x(900:1150)));subplot(2,1,2);plot(real(y(900:1150)));
|
||||
% MOV1=[MOV1 getframe(87)];
|
||||
|
||||
% extract rotation matrix from pre calculated matrices
|
||||
matR = brf.matR{plate_numbers(n)};
|
||||
|
||||
|
||||
% transform to eigenvalue of of fiber segment
|
||||
tOpt.X = conj(matR(1,1))*opt_x + conj(matR(2,1))*opt_y;
|
||||
tOpt.Y = conj(matR(1,2))*opt_x + conj(matR(2,2))*opt_y;
|
||||
|
||||
|
||||
% calculate statistical delta beta for pmd
|
||||
delta_beta = 0.5*(db1+db0(n))/corr_length;
|
||||
|
||||
% build transfer function with delta beta
|
||||
%common.beta = beta1+beta2*omega^2
|
||||
|
||||
%accumulate delta beta for log...
|
||||
brf.simdgd = brf.simdgd + (db1(length(db1)/2+1)+db0(n))/corr_length;
|
||||
|
||||
|
||||
h.X = exp(-1j*(common_beta_x-delta_beta)*dz);
|
||||
h.Y = exp(-1j*(common_beta_y+delta_beta)*dz);
|
||||
|
||||
% delta_beta has to be added to the transfer function
|
||||
|
||||
% process with transfer function
|
||||
tOpt.X = h.X.*tOpt.X ;
|
||||
tOpt.Y = h.Y.*tOpt.Y ;
|
||||
|
||||
|
||||
% rotate back
|
||||
opt_x = matR(1,1)*tOpt.X + matR(1,2)*tOpt.Y;
|
||||
opt_y = matR(2,1)*tOpt.X + matR(2,2)*tOpt.Y;
|
||||
|
||||
|
||||
end
|
||||
|
||||
lin_z_test = lin_z_test + sum(plate_sizes,2);
|
||||
@@ -106,8 +117,9 @@ lin_z_test = lin_z_test + sum(plate_sizes,2);
|
||||
%update the number of processed plates so far
|
||||
n_plates_done = n_plates_done + n_plates_left;
|
||||
|
||||
% attenuate the signal each linear step with alpha
|
||||
rOpt_x=ifft(exp(-alpha_lin_x*aStepSize/2).*opt_x);
|
||||
rOpt_y=ifft(exp(-alpha_lin_y*aStepSize/2).*opt_y);
|
||||
% attanuate the signal each linear state with alpha
|
||||
% ( 0.2dB = 4.6052e-05 )
|
||||
rOpt_x=ifft(exp(-alpha_lin_x*aStepSize/2).*opt_x); % /2 not sure why (have to find it in formulas)
|
||||
rOpt_y=ifft(exp(-alpha_lin_y*aStepSize/2).*opt_y); % but not relevant for now
|
||||
|
||||
end
|
||||
end
|
||||
@@ -10,6 +10,10 @@ classdef EQ < handle
|
||||
training_length %Number of training symbols
|
||||
training_loops %Number of loops through sequence for training mode
|
||||
ideal_dfe %Error free DFE decisions
|
||||
weighted_DFE %Weighted DFE on/off
|
||||
weighted_DFE_mode %Weighted DFE mode
|
||||
weighted_DFE_d_min %d_min threshold parameter for weighted DFE mode 1
|
||||
weighted_DFE_I_mode %[a_s, b_s, I_max]-parameters for the weighted DFE
|
||||
|
||||
DB_aim %Aim at duobinary output sequence
|
||||
|
||||
@@ -68,6 +72,10 @@ classdef EQ < handle
|
||||
options.training_length = 1024 %Number of training symbols
|
||||
options.training_loops = 1 %Number of loops through sequence for training mode
|
||||
options.ideal_dfe = 0 %Error free DFE decisions
|
||||
options.weighted_DFE = 0;
|
||||
options.weighted_DFE_mode = 'R1';
|
||||
options.weighted_DFE_d_min = 0.5;
|
||||
options.weighted_DFE_I_mode = [5,0.5,0.6];
|
||||
|
||||
options.DB_aim %Aim at duobinary output sequence
|
||||
|
||||
@@ -438,6 +446,43 @@ classdef EQ < handle
|
||||
dd_out(k) = constellation_in_(dd_idx);
|
||||
end
|
||||
|
||||
% Implementation of a weighted DFE in
|
||||
% order to prevent error propagation.
|
||||
% For further details, study [1], chapter 3.2.2 -
|
||||
% Modifications of DFE
|
||||
|
||||
if obj.weighted_DFE
|
||||
% define new constellations
|
||||
const = unique(ref_in);
|
||||
|
||||
% determine reliability factor gamma_k
|
||||
if output_vec(m) > min(const) && output_vec(m) < max(const)
|
||||
gamma_k = 1 - abs(output_vec(m) - dd_out(k));
|
||||
else
|
||||
gamma_k = 1;
|
||||
end
|
||||
|
||||
% select mode
|
||||
if strcmp(obj.weighted_DFE_mode,'R1')
|
||||
if gamma_k >= obj.weighted_DFE_d_min
|
||||
f_gamma_k = 1;
|
||||
else
|
||||
f_gamma_k = 0;
|
||||
end
|
||||
elseif strcmp(obj.weighted_DFE_mode,'R2')
|
||||
f_gamma_k = gamma_k;
|
||||
elseif strcmp(obj.weighted_DFE_mode,'I1')
|
||||
nom = 1-exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1));
|
||||
denom = 1+exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1));
|
||||
f_gamma_k = (1/2)*((nom/denom) - 1);
|
||||
elseif strcmp(obj.weighted_DFE_mode,'I2')
|
||||
nom = 1-exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1));
|
||||
denom = 1+exp(-obj.weighted_DFE_I_mode(1)*((gamma_k/obj.weighted_DFE_I_mode(2))-1));
|
||||
f_gamma_k = (obj.weighted_DFE_I_mode(3)/2)*((nom/denom) - 1);
|
||||
end
|
||||
% calculate weighted output
|
||||
output_vec(m) = f_gamma_k.*dd_out(k)+(1-f_gamma_k).*output_vec(m);
|
||||
end
|
||||
|
||||
if obj.Nb(1) > 0
|
||||
dd_DFE(2:end) = dd_DFE(1:end-1);
|
||||
@@ -739,3 +784,6 @@ classdef EQ < handle
|
||||
end
|
||||
end
|
||||
|
||||
% References
|
||||
% [1] T. J. Wettlin, “Experimental Evaluation of Advanced Digital Signal Processing for Intra-Datacenter Systems using Direct-Detection,” 2023. [Online]. Available: https://nbn-resolving.org/urn:nbn:de:gbv:8:3-2023-00703-8
|
||||
|
||||
|
||||
@@ -597,7 +597,7 @@ classdef ML_MLSE < handle
|
||||
obj.S = obj.nSym;
|
||||
obj.Nf = obj.order * obj.sps;
|
||||
obj.nStates = obj.S^obj.L;
|
||||
obj.nFeasible = obj.nStates * obj.S; %feasible state transitions
|
||||
obj.nFeasible = obj.nStates * obj.S;
|
||||
|
||||
% --- Trellis mapping
|
||||
obj.trellis_states = reshape(obj.constellation,1,[]);
|
||||
|
||||
@@ -1,47 +0,0 @@
|
||||
classdef Timing_Recovery < handle
|
||||
|
||||
properties(Access=public)
|
||||
timing_error_detector
|
||||
sps
|
||||
damping_factor
|
||||
normalized_loop_bandwidth
|
||||
detector_gain
|
||||
end
|
||||
|
||||
methods(Access=public)
|
||||
function obj = Timing_Recovery(options)
|
||||
arguments(Input)
|
||||
|
||||
options.timing_error_detector = 'Gardner';
|
||||
options.sps = 2;
|
||||
options.damping_factor = 1.0;
|
||||
options.normalized_loop_bandwidth = 0.005;
|
||||
options.detector_gain = 1;
|
||||
|
||||
end
|
||||
|
||||
fn = fieldnames(options);
|
||||
for n = 1:numel(fn)
|
||||
obj.(fn{n}) = options.(fn{n});
|
||||
end
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
||||
function [data_out,timing_error] = process(obj, data_in)
|
||||
|
||||
timing_synchronization = comm.SymbolSynchronizer( ...
|
||||
"TimingErrorDetector", obj.timing_error_detector, ...
|
||||
"SamplesPerSymbol", obj.sps, ...
|
||||
"DampingFactor", obj.damping_factor, ...
|
||||
"NormalizedLoopBandwidth", obj.normalized_loop_bandwidth, ...
|
||||
"DetectorGain", obj.detector_gain);
|
||||
|
||||
data_out = data_in;
|
||||
[data_out.signal,timing_error] = timing_synchronization(data_in.signal);
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
@@ -23,12 +23,6 @@ classdef DBHandler < handle
|
||||
arguments
|
||||
options.dataBase = ""; % Default value for pathToDB if not provided
|
||||
options.type = "mysql";
|
||||
options.server = "";
|
||||
options.port = 3306;
|
||||
options.user = "";
|
||||
options.password = "";
|
||||
|
||||
|
||||
end
|
||||
|
||||
% Assign values to class properties based on input arguments
|
||||
@@ -52,13 +46,12 @@ classdef DBHandler < handle
|
||||
|
||||
obj.conn = database( ...
|
||||
string(obj.dataBase), ... % Database name
|
||||
options.user, ... % Username
|
||||
options.password, ... % Password (or getSecret)
|
||||
"silas", ... % Username
|
||||
"silas", ... % Password (or getSecret)
|
||||
"Vendor", "MySQL", ...
|
||||
"Server", options.server, ...
|
||||
"Server", "134.245.243.254", ...
|
||||
"PortNumber", 3306, ...
|
||||
"JDBCDriverLocation", "C:\Users\Silas\Documents\mysql-connector-j-9.3.0\mysql-connector-j-9.3.0.jar");
|
||||
|
||||
end
|
||||
|
||||
catch e
|
||||
|
||||
@@ -54,6 +54,12 @@ for k = 1:numel(fn)
|
||||
end
|
||||
json_str = jsonencode(eq_small);
|
||||
|
||||
fn = fieldnames(eq_);
|
||||
for k = 1:numel(fn)
|
||||
if issparse(eq_.(fn{k}))
|
||||
eq_.(fn{k}) = full(eq_.(fn{k}));
|
||||
end
|
||||
end
|
||||
ml_mlse_results.config.eq = jsonencode(eq_);
|
||||
ml_mlse_results.config.equalizer_structure = int32(equalizer_structure.ml_mlse);
|
||||
ml_mlse_results.config.comment = 'function: ML-based MLSE';
|
||||
@@ -87,6 +93,7 @@ end
|
||||
function [bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, precode_mode, M, eth_style)
|
||||
% Calculate BER based on precoding mode
|
||||
mapper = PAMmapper(M, 0, "eth_style", eth_style);
|
||||
skip_front = 150;
|
||||
|
||||
switch precode_mode
|
||||
case db_mode.no_db
|
||||
@@ -101,11 +108,11 @@ switch precode_mode
|
||||
tx_bits_precoded = mapper.demap(tx_symbols_precoded);
|
||||
|
||||
rx_bits = mapper.demap(eq_signal_hd_precoded);
|
||||
[~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
|
||||
[~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1);
|
||||
|
||||
% B) Just determine BER
|
||||
rx_bits = mapper.demap(eq_signal_hd);
|
||||
[bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
|
||||
[bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1);
|
||||
|
||||
case db_mode.db_precoded
|
||||
% Data is precoded on TX side
|
||||
@@ -113,12 +120,12 @@ switch precode_mode
|
||||
eq_signal_hd_decoded = Duobinary().encode(eq_signal_hd, "M", M);
|
||||
eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M);
|
||||
rx_bits_decoded = mapper.demap(eq_signal_hd_decoded);
|
||||
[~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
|
||||
[~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1);
|
||||
|
||||
% B) Omit the Coding by comparing with demapped TX symbol sequence
|
||||
tx_bits_demapped = mapper.demap(tx_symbols);
|
||||
rx_bits = mapper.demap(eq_signal_hd);
|
||||
[bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
|
||||
[bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", skip_front, "skip_end", 150, "returnErrorLocation", 1);
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
@@ -1,219 +0,0 @@
|
||||
|
||||
|
||||
% Parameters
|
||||
c0 = physconst('lightspeed'); % [m/s]
|
||||
lambda0 = 1310e-9; % [m]
|
||||
omega0 = 2*pi*c0/lambda0;
|
||||
|
||||
L = 5e-3; % [m] effective phase section length (set as needed)
|
||||
n_eff = 2.2; % [-] effective index (set as needed)
|
||||
|
||||
E0 = 1; % field amplitude (arbitrary)
|
||||
Vpi = 3.2; % [V] half-wave voltage (your V_pi)
|
||||
|
||||
% Drive
|
||||
f0 = 1e9; % [Hz]
|
||||
fs = 200e9; % [Hz]
|
||||
Nper = 2; % number of periods
|
||||
Vpp = 0.6*Vpi; % [V] peak-to-peak of v_drive(t)
|
||||
|
||||
biasV = 1.1; % [V] differential bias added to v_drive
|
||||
|
||||
% Time axis + differential drive voltage v_drive(t)
|
||||
T = Nper/f0;
|
||||
t = (0:1/fs:T-1/fs).';
|
||||
|
||||
|
||||
if 1
|
||||
% SINE
|
||||
v_drive = biasV + (Vpp/2)*sin(2*pi*f0*t); % v_drive(t) (peak = Vpp/2)
|
||||
|
||||
else
|
||||
|
||||
% --- Generate PAM-4 Sequence ---
|
||||
symbols = linspace(-0.5, 0.5, 4);
|
||||
num_symbols = 12; % Increased slightly for better visual
|
||||
rng(44);
|
||||
random_data = symbols(randi(4, 1, num_symbols));
|
||||
|
||||
% Create time axis (Note: T is your period from the sine code)
|
||||
sps = round(T * fs);
|
||||
t = (0:1/fs:(num_symbols*T)-1/fs).';
|
||||
|
||||
% Upsample to rectangular waveform
|
||||
v_pam = repelem(random_data, sps).';
|
||||
|
||||
% Apply swing and bias: Resulting range is [biasV-Vpp/2, biasV+Vpp/2]
|
||||
v_drive_rect = biasV + (v_pam * Vpp);
|
||||
|
||||
% --- Round the edges ---
|
||||
filter_span = round(sps/1.5); % Increased span for smoother "rounding"
|
||||
window = gausswin(filter_span);
|
||||
window = window / sum(window);
|
||||
|
||||
% Apply filter (using 'same' to keep vector length, but be aware of edge transients)
|
||||
v_drive = conv(v_drive_rect, window, 'same');
|
||||
|
||||
end
|
||||
|
||||
|
||||
% Analytic
|
||||
v_ = linspace(-1,2, 2001);
|
||||
% Field transfer function (amplitude)
|
||||
Field_mzm_analytic = cos((pi/2)*v_);
|
||||
|
||||
% Power transfer function (intensity)
|
||||
P_mzm_analytic = Field_mzm_analytic.^2;
|
||||
|
||||
% Imbalance factor in YOUR notation:
|
||||
rho = 1;
|
||||
|
||||
% Push-pull branch voltages (consistent with v_drive = v1 - v2)
|
||||
v1 = +0.5*v_drive; % arm 1
|
||||
v2 = -0.5*v_drive; % arm 2
|
||||
|
||||
% Phases phi1, phi2
|
||||
phi1 = pi * v1 / Vpi;
|
||||
phi2 = pi * v2 / Vpi;
|
||||
|
||||
% Fields: E_in and E_out (exactly your Eq. (mzm_e_field))
|
||||
E_in = E0 .* exp(1i*omega0*t);
|
||||
|
||||
common_phase = exp(-1i * (omega0*L*n_eff/c0)); % exp(-j*omega0*L*n_eff/c0)
|
||||
|
||||
E_out = E0 .* exp(1i*omega0*t) .* common_phase .* 0.5 .* ...
|
||||
( exp(-1i*phi1) + rho .* exp(-1i*phi2) );
|
||||
|
||||
% Transfer function (numerical): E_out/E_in
|
||||
H_num = E_out ./ E_in;
|
||||
|
||||
% Power (normalized)
|
||||
Pnorm_num = abs(H_num).^2; % since |E_out/E_in|^2
|
||||
|
||||
% Ideal TF (analytic) for comparison (rho=1, push-pull)
|
||||
H_ideal = common_phase .* cos( (pi/2) * (v_drive./Vpi) );
|
||||
|
||||
Pnorm_ideal = abs(H_ideal).^2;
|
||||
Pnorm_math = cos( (pi/2) * (v_drive./Vpi) ).^2;
|
||||
|
||||
|
||||
set(groot, 'defaultLegendInterpreter', 'tex');
|
||||
set(groot, 'defaultAxesTickLabelInterpreter', 'tex');
|
||||
set(groot, 'defaultTextInterpreter', 'tex');
|
||||
|
||||
% Normalized voltage axis (multiples of Vpi)
|
||||
v_norm = v_drive./Vpi;
|
||||
|
||||
colfield = [0,0,0]; %is black
|
||||
colpow = linspecer(2);
|
||||
colpow = colpow(1,:);
|
||||
colvdrive = linspecer(2);
|
||||
colvdrive = colvdrive(2,:);
|
||||
|
||||
%% SIGNAL IN
|
||||
figure(1); clf
|
||||
plot(v_norm,t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on;
|
||||
ylabel('t [ns]'); xlabel('v_{drive}(t)/V_\pi');
|
||||
title('Drive voltage (normalized)');
|
||||
xlim([min(v_) max(v_)]);
|
||||
% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\linear_casee\mzm_input_signal.tex');
|
||||
|
||||
|
||||
%% IN/OUT (static transfer) — normalized x-axis + analytic curve
|
||||
if 0
|
||||
figure(2); clf
|
||||
plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF
|
||||
plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color',colpow); hold on;% analytic power TF
|
||||
% show input time signal
|
||||
plot(v_norm,-1+t*1e9, 'LineWidth', 1.0,'Color',colvdrive); grid on;
|
||||
% show output time signal
|
||||
plot(2+t*1e9, Pnorm_num, 'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colvdrive); hold on;
|
||||
plot(2+t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on;
|
||||
scatter(v_norm, Pnorm_num, 12, '.', 'LineWidth', 1,'MarkerEdgeColor',colvdrive);
|
||||
scatter(biasV./Vpi,(cos((pi/2)*biasV./Vpi)^2),10,'Marker','o');
|
||||
line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--');
|
||||
line([max(v_drive) max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--');
|
||||
xline([min(v_norm) max(v_norm)])
|
||||
|
||||
grid on;
|
||||
xlabel('v_{drive}(t)/V_\pi'); ylabel('|E_{out}/E_{in}|^2');
|
||||
% legend
|
||||
xlim([min(v_) max(v_)+1]);
|
||||
ylim([-1 1]);
|
||||
|
||||
% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm.tex');
|
||||
end
|
||||
%%
|
||||
|
||||
figure(3); clf
|
||||
plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF
|
||||
plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color',colpow); hold on;% analytic power TF
|
||||
|
||||
scatter(v_norm, Pnorm_num, 12, '.', 'LineWidth', 1,'MarkerEdgeColor',colvdrive);
|
||||
scatter(biasV./Vpi,(cos((pi/2)*biasV./Vpi)^2),10,'Marker','o');
|
||||
line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--');
|
||||
line([max(v_drive) max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), -2],'linewidth',0.5,'color','black','linestyle','--');
|
||||
xline([min(v_norm) max(v_norm)])
|
||||
|
||||
grid on;
|
||||
xlabel('v_{drive}(t)/V_\pi'); ylabel('|E_{out}/E_{in}|^2');
|
||||
% legend
|
||||
xlim([min(v_) max(v_)]);
|
||||
ylim([-1 1]);
|
||||
|
||||
% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\mzm_tramsfer_function_matlab.tex');
|
||||
|
||||
|
||||
%%
|
||||
|
||||
figure(4); clf
|
||||
% plot(v_, Field_mzm_analytic, 'LineWidth', 1.2,'LineStyle','--','Color',colfield); hold on;% analytic power TF
|
||||
plot(v_, P_mzm_analytic, 'LineWidth', 1.2, 'Color','black'); hold on;% analytic power TF
|
||||
input_dots = linspace(min(v_drive),max(v_drive),4)./Vpi;
|
||||
% input_dots = unique(v_drive_rect)./Vpi;
|
||||
output_dots = (cos((pi/2)*input_dots).^2);
|
||||
scatter(input_dots,output_dots,'Marker','x','LineWidth',1,'MarkerEdgeColor','black');
|
||||
scatter(input_dots,zeros(size(input_dots)),'Marker','^','LineWidth',2,'MarkerEdgeColor','black');
|
||||
% scatter(ones(size(input_dots)),output_dots,'Marker','<','LineWidth',2,'MarkerEdgeColor','black');
|
||||
|
||||
for i = 1:numel(input_dots)
|
||||
% Draw the dashed projection lines
|
||||
line([input_dots(i), input_dots(i)], [output_dots(i), 0], 'linewidth', 0.5, 'color', 'black', 'linestyle', '--', 'handlevisibility', 'off');
|
||||
line([input_dots(i), 1], [output_dots(i), output_dots(i)], 'linewidth', 0.5, 'color', 'black', 'linestyle', '--', 'handlevisibility', 'off');
|
||||
|
||||
% Add the level annotation boxes near the output (y-axis)
|
||||
% Adjust the '1.05' to move the box further right or 'output_dots(i)' for height
|
||||
j = 3-(i-1)*2;
|
||||
text(1, output_dots(i), sprintf('Level %d', j), ...
|
||||
'FontSize', 8, ...
|
||||
'EdgeColor', 'black', ...
|
||||
'BackgroundColor', 'white', ...
|
||||
'Margin', 2);
|
||||
end
|
||||
|
||||
xlim([0,1.5]);
|
||||
ylim([0,1])
|
||||
|
||||
% line([min(v_drive), min(v_drive)]./Vpi,[(cos((pi/2)*min(v_drive)./Vpi)^2), 0],'linewidth',0.5,'color','black','linestyle','--');
|
||||
% line([max(v_drive), max(v_drive)]./Vpi,[(cos((pi/2)*max(v_drive)./Vpi)^2), 0],'linewidth',0.5,'color','black','linestyle','--');
|
||||
|
||||
% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\linear_casee\mzm_tf.tex');
|
||||
% xticks(sort(input_dots));
|
||||
% yticks(sort(output_dots));
|
||||
grid off
|
||||
|
||||
|
||||
%%
|
||||
% % FIELD TF (only field here; do not mix power into this figure)
|
||||
figure(5); clf
|
||||
% plot(t*1e9, real(H_num), 'LineWidth', 1.0); hold on;
|
||||
% plot(t*1e9, real(H_ideal), '--', 'LineWidth', 1.0,'DisplayName','Field','Color',colfield); hold on;
|
||||
plot(t*1e9, Pnorm_num, 'LineWidth', 1.0,'DisplayName','Intensity', 'Color',colpow); hold on;
|
||||
grid on;
|
||||
xlabel('t [ns]'); ylabel('Re\{E_{out}/E_{in}\}');
|
||||
legend
|
||||
yticks(sort(output_dots));
|
||||
% mat2tikz_improved('C:\Users\Silas\Documents\6971e0b65b380ca6d71c837f\02_IMDD_System\tikz\linear_casee\mzm_output_signal.tex');
|
||||
|
||||
|
||||
|
||||
@@ -1,23 +0,0 @@
|
||||
function mat2tikz_improved(filename)
|
||||
arguments
|
||||
% Default to the path in your example if no argument is provided
|
||||
filename (1,1) string = 'C:\Users\Silas\Documents\Dissertation\00_Examples\tikz\textfig.tikz';
|
||||
end
|
||||
cleanfigure;
|
||||
matlab2tikz(char(filename), ...
|
||||
'width','\fwidth', ...
|
||||
'height','\fheight', ...
|
||||
'showInfo',false, ...
|
||||
'extraAxisOptions',{ ...
|
||||
'legend style={font=\footnotesize}', ...
|
||||
'xlabel style={font=\color{white!15!black},font=\small},',...
|
||||
'ylabel style={font=\color{white!15!black},font=\small},',...
|
||||
'legend columns=1', ...
|
||||
'every axis/.append style={font=\scriptsize}',...
|
||||
'legend columns=1',...
|
||||
'legend style={at={(0.02,0.98)},font=\footnotesize,draw=black!60,rounded corners=2pt,inner sep=1pt,fill=white,column sep=6pt,anchor= north west}',...
|
||||
'legend style={at={(0.02,0.98)},draw=white!0!white,font=\scriptsize,inner sep=0.1pt,fill=white,column sep=1pt,anchor= north west}',...
|
||||
'every axis/.append style={font=\scriptsize}',...
|
||||
});
|
||||
|
||||
end
|
||||
@@ -1,13 +1,10 @@
|
||||
classdef WesPalette
|
||||
% WESPALETTE Wes Anderson color palettes with auto-completion
|
||||
% Usage:
|
||||
% cmap = WesPalette.Zissou1.rgb() % full palette
|
||||
% cmap = WesPalette.Zissou1.rgb(3) % 3 colors (discrete default)
|
||||
% cmap = WesPalette.Zissou1.rgb(12,"discrete") % any n, no interpolation
|
||||
% cmap = WesPalette.Zissou1.rgb(256,"continuous") % smooth colormap (Lab interpolation)
|
||||
%
|
||||
% Requires:
|
||||
% - colorspace.m (Pascal Getreuer) on MATLAB path for "continuous" mode
|
||||
% cmap = WesPalette.Zissou1.rgb()
|
||||
% cmap = WesPalette.Zissou1.rgb(3)
|
||||
|
||||
% https://github.com/karthik/wesanderson?tab=readme-ov-file
|
||||
|
||||
enumeration
|
||||
BottleRocket1
|
||||
@@ -37,33 +34,21 @@ classdef WesPalette
|
||||
end
|
||||
|
||||
methods
|
||||
function cmap = rgb(obj, n, mode)
|
||||
function cmap = rgb(obj, n)
|
||||
% Return palette as Nx3 RGB colormap [0–1]
|
||||
%
|
||||
% n : number of requested colors (optional)
|
||||
% mode : "discrete" (default) or "continuous"
|
||||
|
||||
base_hex = obj.hex();
|
||||
base_rgb = WesPalette.hex2rgb(base_hex);
|
||||
hex = obj.hex();
|
||||
|
||||
if nargin < 2 || isempty(n)
|
||||
cmap = base_rgb;
|
||||
return;
|
||||
end
|
||||
if nargin < 3 || isempty(mode)
|
||||
mode = "discrete";
|
||||
end
|
||||
mode = lower(string(mode));
|
||||
rgb = hex2rgb(hex);
|
||||
|
||||
validateattributes(n, {'numeric'}, {'scalar','integer','positive'}, mfilename, 'n');
|
||||
if mode ~= "discrete" && mode ~= "continuous"
|
||||
error('mode must be "discrete" or "continuous".');
|
||||
end
|
||||
|
||||
if mode == "discrete"
|
||||
cmap = WesPalette.sample_discrete(base_rgb, n);
|
||||
if nargin == 2
|
||||
if n > size(rgb,1)
|
||||
error('Requested %d colors, but only %d available.', ...
|
||||
n, size(rgb,1))
|
||||
end
|
||||
cmap = rgb(1:n,:);
|
||||
else
|
||||
cmap = WesPalette.interpolate_continuous_lab(base_rgb, n);
|
||||
cmap = rgb;
|
||||
end
|
||||
end
|
||||
end
|
||||
@@ -71,116 +56,78 @@ classdef WesPalette
|
||||
methods (Access = private)
|
||||
function hex = hex(obj)
|
||||
% Internal HEX storage
|
||||
|
||||
switch obj
|
||||
case WesPalette.BottleRocket1
|
||||
hex = {'#A42820','#5F5647','#9B110E','#3F5151','#4E2A1E','#550307','#0C1707'};
|
||||
|
||||
case WesPalette.BottleRocket2
|
||||
hex = {'#FAD510','#CB2314','#273046','#354823','#1E1E1E'};
|
||||
|
||||
case {WesPalette.Rushmore1, WesPalette.Rushmore}
|
||||
hex = {'#E1BD6D','#EABE94','#0B775E','#35274A','#F2300F'};
|
||||
|
||||
case WesPalette.Royal1
|
||||
hex = {'#899DA4','#C93312','#FAEFD1','#DC863B'};
|
||||
|
||||
case WesPalette.Royal2
|
||||
hex = {'#9A8822','#F5CDB4','#F8AFA8','#FDDDA0','#74A089'};
|
||||
|
||||
case WesPalette.Zissou1
|
||||
hex = {'#3B9AB2','#78B7C5','#EBCC2A','#E1AF00','#F21A00'};
|
||||
|
||||
case WesPalette.Zissou1Continuous
|
||||
hex = {'#3A9AB2','#6FB2C1','#91BAB6','#A5C2A3','#BDC881', ...
|
||||
'#DCCB4E','#E3B710','#E79805','#EC7A05','#EF5703','#F11B00'};
|
||||
|
||||
case WesPalette.Darjeeling1
|
||||
hex = {'#FF0000','#00A08A','#F2AD00','#F98400','#5BBCD6'};
|
||||
|
||||
case WesPalette.Darjeeling2
|
||||
hex = {'#ECCBAE','#046C9A','#D69C4E','#ABDDDE','#000000'};
|
||||
|
||||
case WesPalette.Chevalier1
|
||||
hex = {'#446455','#FDD262','#D3DDDC','#C7B19C'};
|
||||
|
||||
case WesPalette.FantasticFox1
|
||||
hex = {'#DD8D29','#E2D200','#46ACC8','#E58601','#B40F20'};
|
||||
|
||||
case WesPalette.Moonrise1
|
||||
hex = {'#F3DF6C','#CEAB07','#D5D5D3','#24281A'};
|
||||
|
||||
case WesPalette.Moonrise2
|
||||
hex = {'#798E87','#C27D38','#CCC591','#29211F'};
|
||||
|
||||
case WesPalette.Moonrise3
|
||||
hex = {'#85D4E3','#F4B5BD','#9C964A','#CDC08C','#FAD77B'};
|
||||
|
||||
case WesPalette.Cavalcanti1
|
||||
hex = {'#D8B70A','#02401B','#A2A475','#81A88D','#972D15'};
|
||||
|
||||
case WesPalette.GrandBudapest1
|
||||
hex = {'#F1BB7B','#FD6467','#5B1A18','#D67236'};
|
||||
|
||||
case WesPalette.GrandBudapest2
|
||||
hex = {'#E6A0C4','#C6CDF7','#D8A499','#7294D4'};
|
||||
|
||||
case WesPalette.IsleofDogs1
|
||||
hex = {'#9986A5','#79402E','#CCBA72','#0F0D0E','#D9D0D3','#8D8680'};
|
||||
|
||||
case WesPalette.IsleofDogs2
|
||||
hex = {'#EAD3BF','#AA9486','#B6854D','#39312F','#1C1718'};
|
||||
|
||||
case WesPalette.FrenchDispatch
|
||||
hex = {'#90D4CC','#BD3027','#B0AFA2','#7FC0C6','#9D9C85'};
|
||||
|
||||
case WesPalette.AsteroidCity1
|
||||
hex = {'#0A9F9D','#CEB175','#E54E21','#6C8645','#C18748'};
|
||||
|
||||
case WesPalette.AsteroidCity2
|
||||
hex = {'#C52E19','#AC9765','#54D8B1','#B67C3B','#175149','#AF4E24'};
|
||||
|
||||
case WesPalette.AsteroidCity3
|
||||
hex = {'#FBA72A','#D3D4D8','#CB7A5C','#5785C1'};
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
methods (Static, Access = private)
|
||||
function rgb = hex2rgb(hex)
|
||||
% hex: cellstr like {'#RRGGBB', ...}
|
||||
if isstring(hex), hex = cellstr(hex); end
|
||||
n = numel(hex);
|
||||
rgb = zeros(n,3);
|
||||
for i = 1:n
|
||||
h = char(hex{i});
|
||||
if startsWith(h,'#'), h = h(2:end); end
|
||||
if numel(h) ~= 6
|
||||
error('Invalid HEX color: %s', hex{i});
|
||||
end
|
||||
rgb(i,1) = hex2dec(h(1:2))/255;
|
||||
rgb(i,2) = hex2dec(h(3:4))/255;
|
||||
rgb(i,3) = hex2dec(h(5:6))/255;
|
||||
end
|
||||
end
|
||||
|
||||
function cmap = sample_discrete(base_rgb, n)
|
||||
% No interpolation; allow any n by sampling/repeating.
|
||||
k = size(base_rgb,1);
|
||||
|
||||
if n <= k
|
||||
idx = round(linspace(1, k, n)); % spread across palette
|
||||
idx = max(1, min(k, idx));
|
||||
cmap = base_rgb(idx,:);
|
||||
else
|
||||
reps = floor(n / k);
|
||||
rmd = mod(n, k);
|
||||
cmap = [repmat(base_rgb, reps, 1); base_rgb(1:rmd,:)];
|
||||
end
|
||||
end
|
||||
|
||||
function cmap = interpolate_continuous_lab(base_rgb, n)
|
||||
% Smooth interpolation in Lab using colorspace().
|
||||
% Requires colorspace.m by Pascal Getreuer on MATLAB path.
|
||||
|
||||
k = size(base_rgb,1);
|
||||
if k == 1
|
||||
cmap = repmat(base_rgb, n, 1);
|
||||
return;
|
||||
end
|
||||
|
||||
% Convert to Lab, interpolate each channel, convert back
|
||||
lab = colorspace('Lab<-RGB', base_rgb);
|
||||
|
||||
t_base = linspace(0, 1, k);
|
||||
t_new = linspace(0, 1, n);
|
||||
|
||||
lab_new = zeros(n,3);
|
||||
for c = 1:3
|
||||
lab_new(:,c) = interp1(t_base, lab(:,c), t_new, 'linear');
|
||||
end
|
||||
|
||||
rgb_new = colorspace('RGB<-Lab', lab_new);
|
||||
|
||||
% Clamp to displayable gamut
|
||||
cmap = min(max(rgb_new, 0), 1);
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
@@ -5,8 +5,8 @@ y2 = 1e0 ./ (1 + exp(-0.4*(x-12))); % NLPN
|
||||
y3 = 1e-6 * 10.^(0.45*x); % RP on gamma
|
||||
y4 = 1e-2 * 10.^(0.18*(x-8)); % RP on beta2
|
||||
|
||||
cmap = WesPalette.AsteroidCity1;
|
||||
% cmap = linspecer(4);
|
||||
cmap = WesPalette.AsteroidCity1.rgb(4);
|
||||
cmap = linspecer(4);
|
||||
figure1=figure(202998);clf;hold on
|
||||
lw = 0.8; ms = 4;
|
||||
plot(x,y1,'LineWidth',lw,'Color',cmap(1,:),'Marker','o','MarkerEdgeColor',cmap(1,:),'MarkerFaceColor',[1,1,1],'MarkerSize',ms);
|
||||
|
||||
@@ -1,252 +0,0 @@
|
||||
|
||||
base = "C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC";
|
||||
mode = 0; %0 oder 1
|
||||
M = 2;
|
||||
|
||||
all_files = dir(fullfile(base, "**/*.mat"));
|
||||
|
||||
if M == 2
|
||||
tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\14G_PAM2\tx_info\tx_info_PAM2_14Gbd0.75RRC.mat");
|
||||
filename = fullfile(base, "14G_PAM2\M=2_Rs=1.4e10_Fs=8e10_I=265mA_RoP=46.3mW_L=31m_PS=RRC_rolloff=0.75_Mode=Rise.mat");
|
||||
elseif M == 4
|
||||
tx_data = load("C:\Users\Silas\Nextcloud4\Dokumente\02_Ablage_Office\FSO_FP_QCL_60umUTC\6G_PAM4\tx_info\tx_info_PAM4_6Gbd0.6RRC.mat");
|
||||
filename = fullfile(base, "6G_PAM4\M=4_Rs=6e9_Fs=8e10_I=255mA_RoP=42.3mW_L=31m_PS=RRC_rolloff=0.6_Mode=Rise.mat");
|
||||
end
|
||||
|
||||
if mode == 1
|
||||
[f, p] = uigetfile(fullfile(base, "**/*.mat"));
|
||||
if f~=0
|
||||
filename = fullfile(p,f);
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
datas = load(filename);
|
||||
|
||||
%%
|
||||
str = filename;
|
||||
M_ = str2double(regexp(str, 'M=([^_]+)', 'tokens', 'once'));
|
||||
assert(M==M_);
|
||||
fsym = str2double(regexp(str, 'Rs=([^_]+)', 'tokens', 'once'));
|
||||
fs = str2double(regexp(str, 'Fs=([^_]+)', 'tokens', 'once'));
|
||||
I = sscanf(char(regexp(str, 'I=([^_]+)', 'tokens', 'once')), '%f');
|
||||
rop = sscanf(char(regexp(str, 'RoP=([^_]+)', 'tokens', 'once')), '%f');
|
||||
L = sscanf(char(regexp(str, 'L=([^_]+)', 'tokens', 'once')), '%f');
|
||||
pulseshape = string( regexp(str, 'PS=([^_]+)', 'tokens', 'once'));
|
||||
rolloff = str2double(regexp(str, 'rolloff=([^_]+)', 'tokens', 'once'));
|
||||
mode = string( regexp(str, 'Mode=([^\.]+)', 'tokens', 'once'));
|
||||
|
||||
%%
|
||||
% Tx data
|
||||
|
||||
Bits = Informationsignal(tx_data.tx_data,"fs",fsym);
|
||||
Symbols = Informationsignal(real(tx_data.tx_PAM_sym),"fs",fsym);
|
||||
|
||||
mapping_style = M==4; % Pam2 is like move-it; PAM-4 is different, same mapping like ETH peopled used in Zurich... hence the "eth_style" argument here and there
|
||||
PM = PAMmapper(M,0,"eth_style",mapping_style); % one should rename "eth style" as this is simply a different mapping scheme
|
||||
|
||||
Symbols_ = PM.map(Bits) .* PM.scaling;
|
||||
assert(isequal(Symbols.signal,Symbols_.signal));
|
||||
|
||||
Bits_ = PM.demap(Symbols);
|
||||
[bits,errors,ber,errorIndice] = calc_ber(Bits_.signal,Bits.signal);
|
||||
assert(ber == 0);
|
||||
|
||||
%% For comparison, apply pulsef on Tx Symbols
|
||||
Pform = Pulseformer("fsym",fsym,"fdac",fs,"pulse","rrc","pulselength",16,"alpha",rolloff);
|
||||
Digi_sig_tx_compare = Pform.process(Symbols);
|
||||
|
||||
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff);
|
||||
Digi_sync = Pform.process(Symbols);
|
||||
|
||||
MF = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff);
|
||||
Rx_sig_compare = MF.process(Digi_sig_tx_compare);
|
||||
|
||||
%%
|
||||
|
||||
% Rx Data
|
||||
traceData = datas.tr.lastData(2).trace.ch3;
|
||||
|
||||
%FYI: Voltage=(RawData−YReference)×YIncrement+YOrigin
|
||||
scoperead_volts = (traceData.RawData - traceData.YReference) * traceData.YIncrement + traceData.YOrigin;
|
||||
demystified = isequal(traceData.YData,scoperead_volts);
|
||||
assert(demystified);
|
||||
|
||||
timesig_compare = [0:1:datas.tr.lastData(1).trace.ch3.Points-1] ./ fs;
|
||||
timesig = datas.tr.lastData(1).trace.ch3.XData;
|
||||
assert(isequal(traceData.YData,scoperead_volts));
|
||||
|
||||
Scope_sig = Electricalsignal(traceData.YData,"fs",fs);
|
||||
|
||||
%%
|
||||
% 1) matched filter
|
||||
% pulse is symmetric, hence we can use pulsef directly as matched filter.
|
||||
% It feels off (bit I think correct) that the fsym is now the output freq.!!
|
||||
% -> output 2 sps to omit timing recovery!?
|
||||
apply_matched_filter = 1;
|
||||
k = 1;
|
||||
if apply_matched_filter
|
||||
|
||||
Pform = Pulseformer("fsym",fsym,"fdac",k*fsym,"pulse","rrc","pulselength",16,"alpha",rolloff,"matched",1);
|
||||
Rx_matched = Pform.process(Scope_sig);
|
||||
else
|
||||
|
||||
Rx_matched = Filter('filtdegree',4,"f_cutoff",fsym*0.5,"fs",Scope_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scope_sig);
|
||||
Rx_matched = Rx_matched.resample("fs_out",k*fsym);
|
||||
end
|
||||
Rx_matched.spectrum();
|
||||
|
||||
|
||||
%%
|
||||
|
||||
coefficients = arburg(Rx_matched.signal,25);
|
||||
|
||||
figure()
|
||||
[h,w] = freqz(1,coefficients,Rx_matched.length,"whole",Rx_matched.fs);
|
||||
h = h/max(abs(h));
|
||||
hold on
|
||||
w_ = (w - Rx_matched.fs/2);
|
||||
plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(coefficients,2)), ' '],'LineWidth',2);
|
||||
|
||||
%% Timing Rec
|
||||
apply_timing_rec = 1;
|
||||
if apply_timing_rec
|
||||
[Rx_symbolsync, timing_error] = Timing_Recovery("timing_error_detector",'Gardner (non-data-aided)','sps',k,'damping_factor',0.1,'normalized_loop_bandwidth',0.1,'detector_gain',2.7).process(Rx_matched);
|
||||
figure();plot(timing_error);
|
||||
Rx_symbolsync.fs = fsym;
|
||||
% Tsynch
|
||||
[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_symbolsync.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0);
|
||||
Rx_synced = Rx_synced_cell{1};
|
||||
sps = 1;
|
||||
else
|
||||
% Tsynch
|
||||
[~,Rx_synced_cell,inverted,sequenceFound,sequenceStarts] = Rx_matched.tsynch("reference", Symbols, "fs_ref", fsym, "debug_plots", 0);
|
||||
Rx_synced = Rx_synced_cell{1};
|
||||
Rx_synced = Rx_synced.resample("fs_out",2*fsym);
|
||||
sps = 2;
|
||||
end
|
||||
|
||||
if M == 2
|
||||
ber_in_paper = 10^(-2.6); %fig 3a) 4 Gb/s MWIR FSO Transmission using Directly Modulated QCL and an Uncooled UTC-PD at Room-Temperature
|
||||
elseif M == 4
|
||||
ber_in_paper = 10^(-2.5);
|
||||
end
|
||||
|
||||
|
||||
%%
|
||||
|
||||
Rx_synced = Rx_synced_cell{1};
|
||||
|
||||
|
||||
% -------------------- FFE --------------------
|
||||
% requires some more digging what is going on :-)
|
||||
eq_ffe = EQ("Ne",[50, 1, 1],"Nb",[2,0,0], ...
|
||||
"training_length",512,"training_loops",5,"dd_loops",5, ...
|
||||
"K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
|
||||
"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
|
||||
|
||||
vars = logspace(-4,-3,36);
|
||||
|
||||
parfor i = 1:numel(vars)
|
||||
|
||||
len_tr = 4096;
|
||||
mu_ffe1 = 0.01;% mus(i);%0.0001;
|
||||
mu_ffe2 = 0.0008;
|
||||
mu_ffe3 = 0.001;
|
||||
mu_dc = 0.005;
|
||||
mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3];
|
||||
mu_dfe = vars(i);
|
||||
duob_mode = db_mode.no_db;
|
||||
|
||||
% requires some more digging what is going on :-)
|
||||
eq_ffe_1 = EQ("Ne",[150, 1, 0],"Nb",[50,0,0], ...
|
||||
"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
|
||||
"K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
|
||||
"FFEmu",0,"plotfinal",0,"ideal_dfe",0);
|
||||
|
||||
eq_ffe_2 = FFE("epochs_tr",1,"epochs_dd",vars(i),"len_tr",4096,"mu_dd",vars(i),"mu_tr",vars(i),"order",999,"sps",1,"decide",0, "adaption",adaption_method.nlms,"dd_mode",0);
|
||||
% eq_ffe_2 = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"ffe_mu_dd",1e-5,"dfe_mu_dd",mus(i),"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",50,"dfe_order",10,"sps",1,"decide",1);
|
||||
|
||||
|
||||
ffe_results = ffe(eq_ffe_1,M,Rx_synced,Symbols,Bits, ...
|
||||
"precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ...
|
||||
"eth_style_symbol_mapping",mapping_style);
|
||||
|
||||
ffe_results.metrics.BER
|
||||
bers(i) = ffe_results.metrics.BER;
|
||||
end
|
||||
|
||||
figure();
|
||||
plot(vars,bers);
|
||||
yline(ber_in_paper)
|
||||
beautifyBERplot();
|
||||
%
|
||||
fprintf('Paper: %.1e \n \n',ber_in_paper);
|
||||
ffe_results.metrics.print("description",'FFE');
|
||||
fprintf('FFE: %.1e \n',ffe_results.metrics.BER);
|
||||
|
||||
|
||||
%% -------------------- VNLE + MLSE --------------------
|
||||
len_tr = 4096;
|
||||
mu_ffe1 = 0.0001;% mus(i);%0.0001;
|
||||
mu_ffe2 = 0.0008;
|
||||
mu_ffe3 = 0.001;
|
||||
mu_dc = 0.005;
|
||||
mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3];
|
||||
mu_dfe = 0.0004;
|
||||
duob_mode = db_mode.no_db;
|
||||
|
||||
pf_ncoeffs = 4;
|
||||
|
||||
vars = 1:7;
|
||||
bers = zeros(size(vars));
|
||||
parfor i = 1:numel(vars)
|
||||
|
||||
eqv = EQ("Ne",[200, 1, 0],"Nb",[2, 0, 0], ...
|
||||
"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
|
||||
"K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
|
||||
"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
|
||||
pf_ncoeffs = vars(i);
|
||||
pf = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
|
||||
mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0,"eth_style",mapping_style).levels);
|
||||
[vnle_results, mlse_results] = vnle_postfilter_mlse(eqv, pf, mlse_, M, Rx_synced, Symbols, Bits, ...
|
||||
"precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", mapping_style);
|
||||
fprintf('Paper: %.1e \n \n',ber_in_paper);
|
||||
vnle_results.metrics.print("description",'VNLE');
|
||||
mlse_results.metrics.print("description",'MLSE');
|
||||
bers(i) = mlse_results.metrics.BER;
|
||||
end
|
||||
%%
|
||||
figure();hold on
|
||||
plot(vars,bers_ffe,'DisplayName','FFE [200,0,0] + PF + MLSE');
|
||||
plot(vars,bers_vnle,'DisplayName','VNLE [200,1,0] + PF + MLSE');
|
||||
plot(vars,bers_vnledfe,'DisplayName','VNLE [200,1,0] + DFE [2] + PF + MLSE');
|
||||
plot(vars,bers_vnledfe_ideal,'DisplayName','VNLE [200,1,0] + ideal DFE [2] + PF + MLSE');
|
||||
yline(ber_in_paper);
|
||||
yline([2e-2, 4.85e-3, 3.8e-3, 2,2e-4],'LineWidth',2,'Color',[0.8,0.8,0.8],'LineStyle',':','HandleVisibility','off');
|
||||
ylim([1e-5,0.1]);
|
||||
beautifyBERplot();
|
||||
|
||||
%% -------------------- DB target --------------------
|
||||
mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels);
|
||||
|
||||
eq_ = EQ("Ne",[50, 5, 5],"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
|
||||
"K",sps,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
|
||||
|
||||
dbt_results = duobinary_target(eq_,mlse_db_, M, Rx_synced, Symbols, Bits, ...
|
||||
"precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [],"eth_style_symbol_mapping",mapping_style);
|
||||
|
||||
fprintf('Paper: %.1e \n \n',ber_in_paper);
|
||||
dbt_results.metrics.print("description",'Duobinary');
|
||||
|
||||
|
||||
%%
|
||||
|
||||
%ML-based MLSE (L=2)
|
||||
mu_ml = 0.01; training_epochs = 100;
|
||||
ml_mlse_equalizer = ML_MLSE("epochs_tr",training_epochs,"epochs_dd",1, ...
|
||||
"len_tr",len_tr,"mu_dd",mu_ml,"mu_tr",mu_ml,"order",11,"sps",1, ...
|
||||
"traceback_depth",128,"L",3,"delta",4,"adaptive_mu",0);
|
||||
|
||||
[ml_mlse_results] = ml_mlse(ml_mlse_equalizer, M, Rx_synced, Symbols, Bits,"precode_mode",duob_mode);
|
||||
ml_mlse_results.metrics.print("description",'ML ');
|
||||
@@ -6,7 +6,7 @@ db = DBHandler("dataBase", "labor_highspeed", "type", database_type);
|
||||
|
||||
pam_levels = [4, 6, 8]; % three tiles
|
||||
bitrate_set = 360e9;
|
||||
fiberL = 2;
|
||||
fiberL = 10;
|
||||
|
||||
fields = [
|
||||
db.getTableFieldNames('power_state_info');
|
||||
@@ -96,7 +96,7 @@ end
|
||||
%% ============================================================
|
||||
% PLOT — 1×3 (PAM-4, PAM-6, PAM-8)
|
||||
% ============================================================
|
||||
fig = figure(9112); clf;
|
||||
fig = figure(9110); clf;
|
||||
tiledlayout(1,3,'TileSpacing','compact','Padding','compact');
|
||||
|
||||
lw = 1.8;
|
||||
@@ -223,19 +223,19 @@ ylabel('');
|
||||
|
||||
end
|
||||
|
||||
% pos = 1e3.*[2.7770 1.2017 1.4000 0.3200];
|
||||
% set(fig, 'Position', pos);
|
||||
pos = 1e3.*[2.7770 1.2017 1.4000 0.3200];
|
||||
set(fig, 'Position', pos);
|
||||
|
||||
%% === EXPORT ===
|
||||
% outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz';
|
||||
% matlab2tikz(outfile, ...
|
||||
% 'width','\fwidth', ...
|
||||
% 'height','\fheight', ...
|
||||
% 'showInfo',false, ...
|
||||
% 'extraAxisOptions',{ ...
|
||||
% 'legend style={font=\footnotesize}', ...
|
||||
% 'legend columns=1' ...
|
||||
% });
|
||||
outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\wavelength_analysis.tikz';
|
||||
matlab2tikz(outfile, ...
|
||||
'width','\fwidth', ...
|
||||
'height','\fheight', ...
|
||||
'showInfo',false, ...
|
||||
'extraAxisOptions',{ ...
|
||||
'legend style={font=\footnotesize}', ...
|
||||
'legend columns=1' ...
|
||||
});
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -1,139 +1,83 @@
|
||||
tablename = 'C:\Users\Silas\Documents\latex\JLT_400G_submission\HighSpeedExperiments_oneandonly_csv.csv';
|
||||
% Returns a Table
|
||||
data = readtable(tablename,"Delimiter",';','DecimalSeparator',',');
|
||||
|
||||
|
||||
%% ============================================================
|
||||
% PLOT
|
||||
% ============================================================
|
||||
%% 1. DATA EXTRACTION & SETUP
|
||||
%% 1. DATA EXTRACTION & SETUP
|
||||
raw_M = data.M;
|
||||
raw_baud = data.BaudRate;
|
||||
raw_net = data.NetRate;
|
||||
raw_codes = string(data.ZoteroCode);
|
||||
raw_names = string(data.Name);
|
||||
raw_band = string(data.Band);
|
||||
|
||||
% Filter Valid Data
|
||||
target_M = [2, 4, 6, 8];
|
||||
validIdx = ismember(raw_M, target_M) & ~isnan(raw_baud) & ~isnan(raw_net);
|
||||
|
||||
Mvals = raw_M(validIdx);
|
||||
baud = raw_baud(validIdx);
|
||||
netrate = raw_net(validIdx);
|
||||
codes = raw_codes(validIdx);
|
||||
names = raw_names(validIdx);
|
||||
bands = raw_band(validIdx);
|
||||
|
||||
pam_list = target_M;
|
||||
colors = flip(cbrewer2('SET1',4));
|
||||
|
||||
%% 2. PLOT (For Visual Check only)
|
||||
figure; hold on;
|
||||
ms = 20;
|
||||
lw = 0.5;
|
||||
ms = 32; % scatter size
|
||||
lw = 0.8; % line width
|
||||
|
||||
for k = 1:4 % PAM-2/4/6/8
|
||||
|
||||
for k = 1:length(pam_list)
|
||||
M = pam_list(k);
|
||||
idxPam = (Mvals == M);
|
||||
|
||||
|
||||
% Extract for this PAM
|
||||
x = baud(idxPam);
|
||||
y = netrate(idxPam);
|
||||
b = bands(idxPam);
|
||||
n = names(idxPam);
|
||||
|
||||
% Get color for this PAM format
|
||||
col = colors(k,:);
|
||||
|
||||
|
||||
% ----- LEGEND FLAG (only add one entry per PAM) -----
|
||||
firstLegend = true;
|
||||
|
||||
|
||||
% ---- PLOT ALL POINTS (marker based on publication) ----
|
||||
for i = 1:sum(idxPam)
|
||||
|
||||
% Marker Logic
|
||||
ms = 20;
|
||||
if strcmpi(b(i), 'O')
|
||||
marker = 'o';
|
||||
elseif strcmpi(b(i), 'C')
|
||||
marker = 'd';
|
||||
else
|
||||
marker = 's';
|
||||
end
|
||||
if strcmpi(n(i), 'THIS WORK')
|
||||
marker = 'pentagram';
|
||||
ms = 100;
|
||||
end
|
||||
|
||||
% Plot Scatter
|
||||
% marker selection by publication
|
||||
pubIdx = find(pub_list == n(i), 1);
|
||||
marker = markerlist{mod(pubIdx-1, nMarkers) + 1};
|
||||
|
||||
if firstLegend
|
||||
scatter(x(i), y(i), ms, 'Marker', marker, ...
|
||||
'MarkerEdgeColor', col, 'MarkerFaceColor', col, ...
|
||||
h = scatter(x(i), y(i), ms, ...
|
||||
'Marker', marker, ...
|
||||
'MarkerEdgeColor', col, ...
|
||||
'MarkerFaceColor', col, ...
|
||||
'DisplayName', sprintf('PAM-%d', M));
|
||||
firstLegend = false;
|
||||
else
|
||||
scatter(x(i), y(i), ms, 'Marker', marker, ...
|
||||
'MarkerEdgeColor', col, 'MarkerFaceColor', col, ...
|
||||
h = scatter(x(i), y(i), ms, ...
|
||||
'Marker', marker, ...
|
||||
'MarkerEdgeColor', col, ...
|
||||
'MarkerFaceColor', col, ...
|
||||
'HandleVisibility','off');
|
||||
end
|
||||
|
||||
% ====== CUSTOM DATATIP CONTENT ======
|
||||
dt = h.DataTipTemplate;
|
||||
dt.DataTipRows(1).Label = 'Baud rate';
|
||||
dt.DataTipRows(2).Label = 'Net rate';
|
||||
|
||||
% Add publication name
|
||||
dt.DataTipRows(end+1) = dataTipTextRow('Publication', n(i));
|
||||
|
||||
|
||||
end
|
||||
|
||||
% Fit lines
|
||||
if length(x) >= 3
|
||||
[p, S, mu] = polyfit(x, y, 2);
|
||||
xfit = linspace(min(x), max(x), 200);
|
||||
yfit = polyval(p, xfit, S, mu);
|
||||
plot(xfit, yfit, '-', 'LineWidth', lw, 'Color', col, 'HandleVisibility', 'off');
|
||||
|
||||
% ---- Fit (PAM-specific) ----
|
||||
valid = ~isnan(x) & ~isnan(y);
|
||||
if sum(valid) >= 3
|
||||
p = polyfit(x(valid), y(valid), 2);
|
||||
xfit = linspace(min(x(valid)), max(x(valid)), 200);
|
||||
yfit = polyval(p, xfit);
|
||||
|
||||
plot(xfit, yfit, ':', ...
|
||||
'LineWidth', lw, ...
|
||||
'Color', col, ...
|
||||
'HandleVisibility', 'off'); % do NOT add to legend
|
||||
end
|
||||
end
|
||||
|
||||
grid on; box on;
|
||||
grid on;
|
||||
xlabel('Baud rate [GBd]');
|
||||
ylabel('Net rate [Gb/s]');
|
||||
% title('Check Command Window for TikZ Code');
|
||||
% legend('Location','northwest');
|
||||
|
||||
%% 3. GENERATE TIKZ ANNOTATION CODE
|
||||
% This prints the manual \draw commands to the console
|
||||
legend('Location','northwest');
|
||||
set(gca,'FontSize',11);
|
||||
|
||||
%% GENERATE TIKZ ANNOTATION CODE
|
||||
% This prints the manual \draw commands to the console
|
||||
|
||||
%% GENERATE TIKZ ANNOTATION CODE (Colored Borders + Tiny Font)
|
||||
%% GENERATE TIKZ ANNOTATION CODE (No Arrow, Close Text)
|
||||
fprintf('\n\n%% ===========================================================\n');
|
||||
fprintf('%% COPY THE FOLLOWING LINES INTO YOUR .TEX FILE \n');
|
||||
fprintf('%% (Paste them just before \\end{axis})\n');
|
||||
fprintf('%% ===========================================================\n\n');
|
||||
|
||||
for i = 1:length(baud)
|
||||
bx = baud(i);
|
||||
by = netrate(i);
|
||||
key = codes(i);
|
||||
M_val = Mvals(i);
|
||||
|
||||
% --- PLACEMENT LOGIC ---
|
||||
if M_val == 8
|
||||
% PAM-8: Place Top-Left
|
||||
% 'south east' anchor means the text's bottom-right corner touches the coordinate
|
||||
% shift moves it slightly up and left to clear the marker
|
||||
anchorStr = 'south east';
|
||||
shiftStr = 'shift={(-3pt, 3pt)}';
|
||||
else
|
||||
% Others: Place Bottom-Right
|
||||
% 'north west' anchor means the text's top-left corner touches the coordinate
|
||||
% shift moves it slightly down and right
|
||||
anchorStr = 'north west';
|
||||
shiftStr = 'shift={(3pt, -3pt)}';
|
||||
end
|
||||
|
||||
% --- PRINT COMMAND ---
|
||||
% Uses \node directly at the coordinate (axis cs:...)
|
||||
fprintf('\\node[anchor=%s, %s, font=\\tiny, fill=white, inner sep=1pt] at (axis cs:%.2f, %.2f) {\\cite{%s}};\n', ...
|
||||
anchorStr, shiftStr, bx, by, key);
|
||||
end
|
||||
fprintf('\n')
|
||||
|
||||
%% === EXPORT ===
|
||||
outfile = 'C:\Users\Silas\Documents\latex\JLT_400G_submission\media\matlab2tikz\highspeedresults_test.tikz';
|
||||
|
||||
outfile = 'C:\Users\Silas\Documents\latex\JLT_400G copy\media\matlab2tikz\highspeedresults.tikz';
|
||||
matlab2tikz(outfile, ...
|
||||
'width','\fwidth', ...
|
||||
'height','\fheight', ...
|
||||
|
||||
@@ -10,11 +10,10 @@ if dsp_options.mode == "load_run_id"
|
||||
|
||||
if experiment == "highspeed_2024"
|
||||
|
||||
dsp_options.database_type = "mysql";
|
||||
dsp_options.database_type = 'mysql';
|
||||
dsp_options.dataBase = 'labor_highspeed';%'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\silas_labor_newdsp_newstructure.db';
|
||||
dsp_options.storage_path = 'Z:\2024\sioe_labor\';
|
||||
db = DBHandler("dataBase", [dsp_options.dataBase],...
|
||||
"type", dsp_options.database_type,"server","192.168.178.192","user","silas","password","silas");
|
||||
db = DBHandler("dataBase", [dsp_options.dataBase], "type", dsp_options.database_type);
|
||||
|
||||
elseif experiment == "mpi_ecoc_2025"
|
||||
|
||||
|
||||
@@ -5,9 +5,9 @@ M = 4;
|
||||
randkey = 1;
|
||||
|
||||
% --- Parameter sweep
|
||||
order_range = 5:5:50; % FFE order
|
||||
delta_range = 0:5:20; % delta
|
||||
SNR_dB = 30;
|
||||
order_range = 2:3:11; % FFE order
|
||||
delta_range = 0:2:4; % delta
|
||||
SNR_dB = 20;
|
||||
|
||||
% --- Prepare bit sequence
|
||||
order_bits = 19;
|
||||
@@ -21,14 +21,10 @@ Symbols = PAMmapper(M,0).map(Bits);
|
||||
Symbols.fs = 200e9;
|
||||
|
||||
% --- Channel (minimal ISI + AWGN)
|
||||
h = abs([0.3 0.9 0.3]); h = h/norm(h);
|
||||
|
||||
% h = [1 -1.67085330039878 1.17918163282514 -0.805210559745616 0.571564213123367 -0.296337147529674 0.00649773445209780 0.0854177610195952 -0.0576009020965258 0.0520994427061551 -0.0624586034913656 0.0553280962699552 -0.00705582559925755 -0.0336399056707792 0.0706903719452810 -0.0334124287931977 0.0131699455037966 0.0587431373842994 -0.0515902976066452 0.00647904355473619 0.0137506750904990 -0.0547974515885928 0.00994735499340592 -0.0135513582534086 -0.00463322575007739 0.0277311946101940];
|
||||
% h = h/norm(h);
|
||||
h = [0.3 0.9 0.3]; h = h/norm(h);
|
||||
symbols_filt = Symbols.filter(h,1);
|
||||
symbols_noi = symbols_filt;
|
||||
symbols_noi.signal = awgn(symbols_filt.signal,SNR_dB,'measured');
|
||||
symbols_noi.spectrum();
|
||||
|
||||
% --- Generate all parameter pairs
|
||||
[O,D] = ndgrid(order_range, delta_range);
|
||||
@@ -42,7 +38,7 @@ ce_vec = nan(size(pairs,1),1);
|
||||
ce_training = nan(size(pairs,1),training_len);
|
||||
|
||||
% --- Parallel loop over parameter pairs
|
||||
for k = 1:size(pairs,1)
|
||||
parfor k = 1:size(pairs,1)
|
||||
order_k = pairs(k,1);
|
||||
delta_k = pairs(k,2);
|
||||
|
||||
@@ -93,7 +89,7 @@ end
|
||||
beautifyBERplot
|
||||
ylabel('BER'); xlabel('Filter Order [N]');
|
||||
title('BER vs. Filter order');
|
||||
% ylim([1e-4, 0.1]);
|
||||
ylim([1e-4, 0.1]);
|
||||
yline(3.8e-3,'HandleVisibility','off');
|
||||
yline(2.2e-4,'HandleVisibility','off');
|
||||
|
||||
|
||||
@@ -48,36 +48,36 @@ for i = 1:numel(SNR_dB)
|
||||
symbols_noi = symbols_filt;
|
||||
symbols_noi.signal = awgn(symbols_filt.signal, SNR_dB(i), 'measured'); % AWGN with given SNR
|
||||
|
||||
% Sequence Est L=5
|
||||
mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',h);
|
||||
mlse_.DIR = h;
|
||||
[y_mlse] = mlse_.process(symbols_noi,Symbols);
|
||||
mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse);
|
||||
[~, ~, ber_mlse_l5(i), ~] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
|
||||
fprintf('MLSE L5: %.2e \n',ber_mlse_l5(i));
|
||||
|
||||
% 2nd Approach
|
||||
mu_lms = 0.0005;
|
||||
pf_ncoeffs = 1;
|
||||
eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",16,"sps",1,"dd_mode",1,"adaption_technique","lms");
|
||||
pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
|
||||
|
||||
% FFE
|
||||
[y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols);
|
||||
|
||||
Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe);
|
||||
[~, ~, ber_ffe(i), ~] = calc_ber(Eq_bits.signal, Bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1);
|
||||
fprintf('FFE: %.2e \n',ber_ffe(i));
|
||||
|
||||
% Postfilter
|
||||
[y_white,~] = pf_.process(y_ffe, ffe_noise);
|
||||
|
||||
% Sequence Est
|
||||
mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',pf_.coefficients);
|
||||
[y_mlse] = mlse_.process(y_white,Symbols);
|
||||
mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse);
|
||||
[~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
|
||||
fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i));
|
||||
% % Sequence Est L=5
|
||||
% mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',h);
|
||||
% mlse_.DIR = h;
|
||||
% [y_mlse] = mlse_.process(symbols_noi,Symbols);
|
||||
% mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse);
|
||||
% [~, ~, ber_mlse_l5(i), ~] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
|
||||
% fprintf('MLSE L5: %.2e \n',ber_mlse_l5(i));
|
||||
%
|
||||
% % 2nd Approach
|
||||
% mu_lms = 0.0005;
|
||||
% pf_ncoeffs = 1;
|
||||
% eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",2^13,"mu_dd",mu_lms,"mu_tr",mu_lms,"order",16,"sps",1,"dd_mode",1,"adaption_technique","lms");
|
||||
% pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);
|
||||
%
|
||||
% % FFE
|
||||
% [y_ffe, ffe_noise] = eq_.process(symbols_noi, Symbols);
|
||||
%
|
||||
% Eq_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_ffe);
|
||||
% [~, ~, ber_ffe(i), ~] = calc_ber(Eq_bits.signal, Bits.signal, "skip_front", 0, "skip_end", 0, "returnErrorLocation", 1);
|
||||
% fprintf('FFE: %.2e \n',ber_ffe(i));
|
||||
%
|
||||
% % Postfilter
|
||||
% [y_white,~] = pf_.process(y_ffe, ffe_noise);
|
||||
%
|
||||
% % Sequence Est
|
||||
% mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels,'scale_mode',0,'trellis_exclusion',0,'trellis_state_mode',2,'debug',0,'DIR',pf_.coefficients);
|
||||
% [y_mlse] = mlse_.process(y_white,Symbols);
|
||||
% mlse_bits = PAMmapper(M, 0, "eth_style", 0).demap(y_mlse);
|
||||
% [~, errors, ber_nwf_mlse_l2(i), errpos] = calc_ber(mlse_bits.signal, Bits.signal, "skip_front", 10, "skip_end", 10, "returnErrorLocation", 1);
|
||||
% fprintf('MLSE: %.2e \n',ber_nwf_mlse_l2(i));
|
||||
|
||||
% ML-base MLSE L=2
|
||||
adaptive_mu = 0;
|
||||
|
||||
@@ -1,60 +0,0 @@
|
||||
L = 4; % Oversampling factor
|
||||
fsym = 1e9/4;
|
||||
rollOff = 0.5; % Pulse shaping roll-off factor
|
||||
htx = rcosdesign(rollOff,16, L, 'sqrt');
|
||||
hrx = conj(fliplr(htx));
|
||||
|
||||
signal = zeros(1000,1);
|
||||
signal(500) = 1;
|
||||
Digi_sig = Informationsignal(signal,"fs",fsym);
|
||||
|
||||
% Digi_sig = Digi_sig.resample("fs_out",L*fsym);
|
||||
%
|
||||
% Digi_sig_tx = Digi_sig.filter(htx,1);
|
||||
% Digi_sig_rx = Digi_sig_tx.filter(hrx,1);
|
||||
%
|
||||
% figure;plot(Digi_sig_tx.signal);hold on;plot(Digi_sig_rx.signal),plot(Digi_sig.signal);
|
||||
%
|
||||
%
|
||||
%
|
||||
%
|
||||
pulsef = Pulseformer("alpha",1,"matched",0,"fdac",Digi_sig.fs*L,"fsym",fsym,"pulse","rrc","pulselength",16);
|
||||
Digi_sig = pulsef.process(Digi_sig);
|
||||
|
||||
pulsef = Pulseformer("alpha",1,"matched",1,"fdac",Digi_sig.fs,"fsym",fsym,"pulse","rrc","pulselength",16);
|
||||
Digi_sig_matched = pulsef.process(Digi_sig);
|
||||
|
||||
% Filter:
|
||||
htx = rcosdesign(rollOff,16, L, 'sqrt');
|
||||
% Note half of the target delay is used, because when combined
|
||||
% to the matched filter, the total delay will be achieved.
|
||||
hrx = conj(fliplr(htx));
|
||||
|
||||
figure
|
||||
plot(htx)
|
||||
title('Transmit Filter')
|
||||
xlabel('Index')
|
||||
ylabel('Amplitude')
|
||||
|
||||
figure
|
||||
plot(hrx)
|
||||
title('Rx Filter (Matched Filter)')
|
||||
xlabel('Index')
|
||||
ylabel('Amplitude')
|
||||
|
||||
p = conv(htx,hrx);
|
||||
|
||||
figure
|
||||
plot(p)
|
||||
title('Combined Tx-Rx = Raised Cosine')
|
||||
xlabel('Index')
|
||||
ylabel('Amplitude')
|
||||
|
||||
% And let's highlight the zero-crossings
|
||||
zeroCrossings = NaN*ones(size(p));
|
||||
zeroCrossings(1:L:end) = 0;
|
||||
zeroCrossings((rcDelay)*L + 1) = NaN; % Except for the central index
|
||||
hold on
|
||||
plot(zeroCrossings, 'o')
|
||||
legend('RC Pulse', 'Zero Crossings')
|
||||
hold off
|
||||
@@ -312,7 +312,7 @@ for realiz = 1:s.num_realiz
|
||||
Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ...
|
||||
"beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
|
||||
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
|
||||
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",1,"useSingle",1).process(Opt_sig_wdm_fib);
|
||||
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib);
|
||||
|
||||
% -------- Evaluate at intermediate distance (enqueue jobs) --------
|
||||
if eval_ptr <= nEval && seg == eval_seg(eval_ptr)
|
||||
|
||||
@@ -1,219 +0,0 @@
|
||||
%% Bayesian Optimization for FFE Parameter Tuning
|
||||
% This script uses bayesopt to find optimal mu_dd and mu_tr values
|
||||
% that minimize BER for the FFE equalizer.
|
||||
|
||||
clear; clc;
|
||||
|
||||
%% Setup - Same as gpu_processing_dpfiber.m
|
||||
s.wavelengthplan = calcWavelengthPlan(4, 400e9, 1310);
|
||||
link_length = 10;
|
||||
s.pmd = 0.1;
|
||||
s.gamma = 0.0023;
|
||||
|
||||
s.M = 4;
|
||||
fsym = 112e9;
|
||||
fdac = 2*fsym;
|
||||
fadc = 120000000000;
|
||||
s.random_key = 1;
|
||||
|
||||
% Laser / Modulator
|
||||
vbias_rel = 0.5;
|
||||
u_pi = 4.6;
|
||||
vbias = -vbias_rel*u_pi;
|
||||
laser_linewidth = 0e6;
|
||||
|
||||
duob_mode = db_mode.no_db;
|
||||
rcalpha = 0.05;
|
||||
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
|
||||
|
||||
s.chirpalpha = 0;
|
||||
s.p_launch = 3;
|
||||
s.p = "co";
|
||||
|
||||
N = numel(s.wavelengthplan);
|
||||
|
||||
switch s.p
|
||||
case "co"
|
||||
pol_rot = 100.*ones(1,N);
|
||||
d_local = 0;
|
||||
end
|
||||
|
||||
f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
|
||||
margin = 25e12;
|
||||
f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
|
||||
f_nyq = f_span/2;
|
||||
|
||||
kover = 4;
|
||||
upsample_required = f_nyq./(fdac*kover/2);
|
||||
upsample_pow = 2^nextpow2(upsample_required);
|
||||
|
||||
s.f_opt = fdac*kover*upsample_pow;
|
||||
s.f_opt_nyq = s.f_opt/2;
|
||||
|
||||
s.rop = -8; % Fixed ROP for optimization
|
||||
|
||||
%% Generate TX signals (run once)
|
||||
fprintf('Generating TX signals...\n');
|
||||
for l = 1:N
|
||||
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
|
||||
"fsym",fsym,"M",s.M,"order",15,"useprbs",0, ...
|
||||
"fs_out",fdac, ...
|
||||
"applyclipping",0,"clipfactor",1.5, ...
|
||||
"applypulseform",1,"pulseformer",Pform, ...
|
||||
"randkey",s.random_key+l, ...
|
||||
"mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ...
|
||||
).process();
|
||||
|
||||
Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ...
|
||||
"filterType",filtertypes.gaussian,"active",true);
|
||||
|
||||
El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ...
|
||||
"bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ...
|
||||
"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
|
||||
|
||||
El_sig = El_sig.normalize("mode","oneone");
|
||||
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
|
||||
El_sig = El_sig .* scaling;
|
||||
|
||||
Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ...
|
||||
"lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ...
|
||||
"linewidth",laser_linewidth,"randomkey",s.random_key+l,"alpha",s.chirpalpha).process(El_sig);
|
||||
|
||||
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
|
||||
end
|
||||
|
||||
%% WDM mux + launch
|
||||
Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ...
|
||||
"lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);
|
||||
|
||||
Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
|
||||
"amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm);
|
||||
|
||||
%% Fiber propagation
|
||||
segment_length = 1;
|
||||
nSegments = link_length/segment_length;
|
||||
nSegments = round(nSegments);
|
||||
|
||||
zdw = 1310;
|
||||
randomize_D = true;
|
||||
Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key);
|
||||
|
||||
Opt_sig_wdm_fib = Opt_sig_wdm;
|
||||
fprintf('Running fiber propagation...\n');
|
||||
for seg = 1:nSegments
|
||||
fprintf('Segment %d/%d\n', seg, nSegments);
|
||||
Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ...
|
||||
"beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
|
||||
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
|
||||
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",true).process(Opt_sig_wdm_fib);
|
||||
end
|
||||
|
||||
%% Pre-process to get Rx_sig (do demux once)
|
||||
fprintf('Pre-processing receiver chain...\n');
|
||||
l = 1; % Use channel 1 for optimization
|
||||
|
||||
Opt_sig_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ...
|
||||
"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);
|
||||
|
||||
Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
|
||||
"amplification_db",s.rop).process(Opt_sig_demux{l});
|
||||
|
||||
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ...
|
||||
"nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx);
|
||||
|
||||
rx_bwl = 100e9;
|
||||
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ...
|
||||
"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
|
||||
|
||||
Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
|
||||
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc, ...
|
||||
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth, ...
|
||||
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0, ...
|
||||
"adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig);
|
||||
|
||||
Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
|
||||
|
||||
[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0);
|
||||
Rx_sig = Scpe_cell{1};
|
||||
Rx_sig = Rx_sig.normalize("mode","rms");
|
||||
|
||||
fprintf('Receiver pre-processing complete. Ready for optimization.\n\n');
|
||||
|
||||
%% Define the objective function for bayesopt
|
||||
function ber = ffe_objective(params, Rx_sig, Symbols_l, Tx_bits_l, M, duob_mode)
|
||||
mu_dd = params.mu_dd;
|
||||
mu_tr = params.mu_tr;
|
||||
|
||||
try
|
||||
eq_ffe = FFE("epochs_tr", 5, "epochs_dd", 2, "len_tr", 2^13, ...
|
||||
"mu_dd", mu_dd, "mu_tr", mu_tr, ...
|
||||
"order", 50, "sps", 2, "decide", 0, ...
|
||||
"adaption", adaption_method.nlms, "dd_mode", 1);
|
||||
|
||||
ffe_results = ffe(eq_ffe, M, Rx_sig, Symbols_l, Tx_bits_l, ...
|
||||
"precode_mode", duob_mode, ...
|
||||
'showAnalysis', 0, ...
|
||||
"postFFE", [], ...
|
||||
"eth_style_symbol_mapping", 0);
|
||||
|
||||
ber = ffe_results.metrics.BER;
|
||||
|
||||
if ber == 0
|
||||
ber = 1e-10;
|
||||
end
|
||||
|
||||
if ~isfinite(ber)
|
||||
ber = 0.5;
|
||||
end
|
||||
|
||||
fprintf(' mu_dd=%.4e, mu_tr=%.4e -> BER=%.4e\n', mu_dd, mu_tr, ber);
|
||||
|
||||
catch ME
|
||||
fprintf(' mu_dd=%.4e, mu_tr=%.4e -> FAILED (%s)\n', mu_dd, mu_tr, ME.message);
|
||||
ber = 0.5;
|
||||
end
|
||||
end
|
||||
|
||||
%% Define optimizable variables
|
||||
mu_dd_var = optimizableVariable('mu_dd', [1e-5, 0.1], 'Transform', 'log');
|
||||
mu_tr_var = optimizableVariable('mu_tr', [1e-5, 0.1], 'Transform', 'log');
|
||||
|
||||
%% Run Bayesian Optimization
|
||||
fprintf('========== Starting Bayesian Optimization ==========\n');
|
||||
fprintf('Optimizing mu_dd and mu_tr to minimize BER\n');
|
||||
fprintf('Search range: mu_dd=[1e-5, 0.1], mu_tr=[1e-5, 0.1]\n\n');
|
||||
|
||||
objective_fn = @(params) ffe_objective(params, Rx_sig, Symbols{l}, Tx_bits{l}, s.M, duob_mode);
|
||||
|
||||
results = bayesopt(objective_fn, [mu_dd_var, mu_tr_var], ...
|
||||
'MaxObjectiveEvaluations', 30, ...
|
||||
'AcquisitionFunctionName', 'expected-improvement-plus', ...
|
||||
'IsObjectiveDeterministic', false, ...
|
||||
'ExplorationRatio', 0.5, ...
|
||||
'Verbose', 1, ...
|
||||
'PlotFcn', []);
|
||||
|
||||
%% Display Results
|
||||
fprintf('\n========== FFE Optimization Complete ==========\n');
|
||||
fprintf('Best FFE parameters found:\n');
|
||||
fprintf(' mu_dd = %.6e\n', results.XAtMinObjective.mu_dd);
|
||||
fprintf(' mu_tr = %.6e\n', results.XAtMinObjective.mu_tr);
|
||||
fprintf(' BER = %.6e\n', results.MinObjective);
|
||||
|
||||
%% Verify with optimal parameters
|
||||
fprintf('\nVerifying optimal FFE parameters...\n');
|
||||
best_mu_dd = results.XAtMinObjective.mu_dd;
|
||||
best_mu_tr = results.XAtMinObjective.mu_tr;
|
||||
|
||||
eq_ffe_best = FFE("epochs_tr", 5, "epochs_dd", 2, "len_tr", 2^13, ...
|
||||
"mu_dd", best_mu_dd, "mu_tr", best_mu_tr, ...
|
||||
"order", 50, "sps", 2, "decide", 0, ...
|
||||
"adaption", adaption_method.nlms, "dd_mode", 1);
|
||||
|
||||
ffe_results_best = ffe(eq_ffe_best, s.M, Rx_sig, Symbols{l}, Tx_bits{l}, ...
|
||||
"precode_mode", duob_mode, ...
|
||||
'showAnalysis', 1, ...
|
||||
"postFFE", [], ...
|
||||
"eth_style_symbol_mapping", 0);
|
||||
|
||||
fprintf('\nFinal FFE BER with optimal parameters: %.6e\n', ffe_results_best.metrics.BER);
|
||||
@@ -1,43 +0,0 @@
|
||||
|
||||
|
||||
%%%%% SETTINGS %%%%%%
|
||||
useprbs = 1;
|
||||
M = 8;
|
||||
randkey = 1;
|
||||
fsym = 112e9;
|
||||
viewresults = 0;
|
||||
|
||||
%%%%% Mapping %%%%%
|
||||
M = 6;
|
||||
data = 0:M-1;
|
||||
bitpersymbol = log2(M);
|
||||
|
||||
s = RandStream('twister','Seed',1);
|
||||
bitpattern = randi(s,[0 1], 2^18, 1);
|
||||
bitpattern = bitpattern(1:end-mod(length(bitpattern),5));
|
||||
|
||||
bits_tx = Informationsignal(bitpattern);
|
||||
|
||||
symbols = PAMmapper(M,0).map(bits);
|
||||
|
||||
pam6transitions = combvec(PAMmapper(M,0).levels,PAMmapper(M,0).levels)';
|
||||
|
||||
pam6bits = PAMmapper(6,0,"eth_style",0).demap(reshape(pam6transitions',[],1)./sqrt(10));
|
||||
symbols_rx = PAMmapper(M,0).map(pam6bits).*sqrt(10);
|
||||
|
||||
pam6bits = reshape(pam6bits',5,[])';
|
||||
|
||||
figure; hold on
|
||||
scatter(pam6transitions(:,1), pam6transitions(:,2), 'x', 'LineWidth', 1);
|
||||
n = size(pam6transitions,1);
|
||||
labels = cellstr(char(pam6bits + '0')); % -> N x 1 cell array of char rows
|
||||
text(pam6transitions(:,1), pam6transitions(:,2), labels, ...
|
||||
'HorizontalAlignment','left', 'VerticalAlignment','bottom');
|
||||
|
||||
|
||||
|
||||
bits_rx = PAMmapper(M,0).demap(symbols);
|
||||
|
||||
[~,error_num,ber,error_pos] = calc_ber(bits_tx.signal,bits_rx.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1);
|
||||
|
||||
PAMmapper(8,0).showBitMapping
|
||||
@@ -29,6 +29,9 @@ para.skip =0;
|
||||
para.bruijn = 0;
|
||||
para.reset_prms = 0;
|
||||
para.method = 1;
|
||||
para.pcs = 0;
|
||||
para.shape_para = 0;
|
||||
para.rng_num = 0;
|
||||
|
||||
data_in = [];
|
||||
global loop;
|
||||
|
||||
@@ -1,349 +0,0 @@
|
||||
%% GPU vs CPU Comparison Test for DP_Fiber
|
||||
% This script runs the fiber simulation with and without GPU acceleration
|
||||
% and compares the numerical results.
|
||||
|
||||
% clear; clc;
|
||||
|
||||
%% Setup (same as gpu_processing_dpfiber.m but simplified)
|
||||
s.wavelengthplan = calcWavelengthPlan(4, 400e9, 1310);
|
||||
link_length = 10;
|
||||
s.pmd = 0.1;
|
||||
s.gamma = 0.0023;
|
||||
|
||||
s.M = 4;
|
||||
fsym = 112e9;
|
||||
fdac = 2*fsym;
|
||||
fadc = 120000000000;
|
||||
s.random_key = 1;
|
||||
|
||||
% Laser / Modulator
|
||||
vbias_rel = 0.5;
|
||||
u_pi = 4.6;
|
||||
vbias = -vbias_rel*u_pi;
|
||||
laser_linewidth = 0e6;
|
||||
|
||||
% DB Stuff
|
||||
duob_mode = db_mode.no_db;
|
||||
|
||||
rcalpha = 0.05;
|
||||
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
|
||||
|
||||
s.chirpalpha = 0;
|
||||
s.p_launch = 3;
|
||||
s.p = "co";
|
||||
|
||||
N = numel(s.wavelengthplan);
|
||||
|
||||
switch s.p
|
||||
case "co"
|
||||
pol_rot = 100.*ones(1,N);
|
||||
d_local = 0;
|
||||
case "pair"
|
||||
pol_rot = repmat([100,100,0,0],1,N/4);
|
||||
d_local = 0;
|
||||
case "alt"
|
||||
pol_rot = repmat([100,0,100,0],1,N/4);
|
||||
d_local = 0;
|
||||
case "seg"
|
||||
pol_rot = 100.*ones(1,N);
|
||||
d_local = 3;
|
||||
otherwise
|
||||
error('Unknown fwm_mitigation_technique: %s', string(s.p));
|
||||
end
|
||||
|
||||
f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
|
||||
margin = 5e12;
|
||||
f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
|
||||
f_nyq = f_span/2;
|
||||
|
||||
kover = 4;
|
||||
upsample_required = f_nyq./(fdac*kover/2);
|
||||
upsample_pow = 2^nextpow2(upsample_required);
|
||||
|
||||
s.f_opt = fdac*kover*upsample_pow;
|
||||
s.f_opt_nyq = s.f_opt/2;
|
||||
|
||||
%% TX per channel
|
||||
for l = 1:N
|
||||
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
|
||||
"fsym",fsym,"M",s.M,"order",15,"useprbs",0, ...
|
||||
"fs_out",fdac, ...
|
||||
"applyclipping",0,"clipfactor",1.5, ...
|
||||
"applypulseform",1,"pulseformer",Pform, ...
|
||||
"randkey",s.random_key+l, ...
|
||||
"mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ...
|
||||
).process();
|
||||
|
||||
Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ...
|
||||
"filterType",filtertypes.gaussian,"active",true);
|
||||
|
||||
El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ...
|
||||
"bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ...
|
||||
"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
|
||||
|
||||
clear Digi_sig
|
||||
|
||||
El_sig = El_sig.normalize("mode","oneone");
|
||||
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
|
||||
El_sig = El_sig .* scaling;
|
||||
|
||||
Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ...
|
||||
"lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ...
|
||||
"linewidth",laser_linewidth,"randomkey",s.random_key+l,"alpha",s.chirpalpha).process(El_sig);
|
||||
|
||||
clear El_sig
|
||||
|
||||
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
|
||||
|
||||
clear Eml_out Lp_awg
|
||||
end
|
||||
|
||||
disp('Signal generated for all channels.');
|
||||
|
||||
%% WDM mux + launch
|
||||
Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ...
|
||||
"lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);
|
||||
|
||||
Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
|
||||
"amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm);
|
||||
|
||||
% Save input for both runs
|
||||
Opt_sig_input = Opt_sig_wdm;
|
||||
|
||||
segment_length = 1;
|
||||
nSegments = link_length/segment_length;
|
||||
if abs(nSegments - round(nSegments)) > 1e-12
|
||||
error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length);
|
||||
end
|
||||
nSegments = round(nSegments);
|
||||
|
||||
zdw = 1310;
|
||||
randomize_D = true;
|
||||
|
||||
if nSegments > 0
|
||||
Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key);
|
||||
else
|
||||
Dvec = [];
|
||||
end
|
||||
|
||||
%% Run WITHOUT GPU
|
||||
fprintf('\n========== Running WITHOUT GPU (CPU) ==========\n');
|
||||
Opt_sig_cpu = Opt_sig_input;
|
||||
|
||||
tic;
|
||||
for seg = 1:nSegments
|
||||
fprintf('CPU Segment %d/%d \n',seg, nSegments);
|
||||
|
||||
Opt_sig_cpu = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ...
|
||||
"beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
|
||||
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
|
||||
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",false).process(Opt_sig_cpu);
|
||||
end
|
||||
time_cpu = toc;
|
||||
fprintf('CPU Time: %.3f seconds\n', time_cpu);
|
||||
|
||||
%% Run WITH GPU (Double Precision)
|
||||
fprintf('\n========== Running WITH GPU (Double Precision) ==========\n');
|
||||
Opt_sig_gpu_double = Opt_sig_input;
|
||||
|
||||
tic;
|
||||
for seg = 1:nSegments
|
||||
fprintf('GPU-Double Segment %d/%d \n',seg, nSegments);
|
||||
|
||||
Opt_sig_gpu_double = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ...
|
||||
"beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
|
||||
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
|
||||
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",false).process(Opt_sig_gpu_double);
|
||||
end
|
||||
time_gpu_double = toc;
|
||||
fprintf('GPU Double Time: %.3f seconds\n', time_gpu_double);
|
||||
|
||||
%% Run WITH GPU (Single Precision)
|
||||
fprintf('\n========== Running WITH GPU (Single Precision) ==========\n');
|
||||
Opt_sig_gpu_single = Opt_sig_input;
|
||||
|
||||
tic;
|
||||
for seg = 1:nSegments
|
||||
fprintf('GPU-Single Segment %d/%d \n',seg, nSegments);
|
||||
|
||||
Opt_sig_gpu_single = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ...
|
||||
"beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
|
||||
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
|
||||
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",true).process(Opt_sig_gpu_single);
|
||||
end
|
||||
time_gpu_single = toc;
|
||||
fprintf('GPU Single Time: %.3f seconds\n', time_gpu_single);
|
||||
|
||||
%% Compare Results
|
||||
fprintf('\n========== Numerical Comparison ==========\n');
|
||||
|
||||
sig_cpu = Opt_sig_cpu.signal;
|
||||
sig_gpu_double = Opt_sig_gpu_double.signal;
|
||||
sig_gpu_single = Opt_sig_gpu_single.signal;
|
||||
|
||||
% Check dimensions
|
||||
fprintf('CPU signal size: [%d x %d]\n', size(sig_cpu,1), size(sig_cpu,2));
|
||||
fprintf('GPU Double signal size: [%d x %d]\n', size(sig_gpu_double,1), size(sig_gpu_double,2));
|
||||
fprintf('GPU Single signal size: [%d x %d]\n', size(sig_gpu_single,1), size(sig_gpu_single,2));
|
||||
|
||||
% CPU vs GPU Double
|
||||
fprintf('\n--- CPU vs GPU Double ---\n');
|
||||
diff_cpu_double = abs(sig_cpu - sig_gpu_double);
|
||||
max_diff_cpu_double = max(diff_cpu_double(:));
|
||||
mean_diff_cpu_double = mean(diff_cpu_double(:));
|
||||
rel_diff_cpu_double = max_diff_cpu_double / max(abs(sig_cpu(:)));
|
||||
fprintf('Max absolute difference: %.6e\n', max_diff_cpu_double);
|
||||
fprintf('Mean absolute difference: %.6e\n', mean_diff_cpu_double);
|
||||
fprintf('Max relative difference: %.6e\n', rel_diff_cpu_double);
|
||||
|
||||
% CPU vs GPU Single
|
||||
fprintf('\n--- CPU vs GPU Single ---\n');
|
||||
diff_cpu_single = abs(sig_cpu - sig_gpu_single);
|
||||
max_diff_cpu_single = max(diff_cpu_single(:));
|
||||
mean_diff_cpu_single = mean(diff_cpu_single(:));
|
||||
rel_diff_cpu_single = max_diff_cpu_single / max(abs(sig_cpu(:)));
|
||||
fprintf('Max absolute difference: %.6e\n', max_diff_cpu_single);
|
||||
fprintf('Mean absolute difference: %.6e\n', mean_diff_cpu_single);
|
||||
fprintf('Max relative difference: %.6e\n', rel_diff_cpu_single);
|
||||
|
||||
% GPU Double vs GPU Single
|
||||
fprintf('\n--- GPU Double vs GPU Single ---\n');
|
||||
diff_double_single = abs(sig_gpu_double - sig_gpu_single);
|
||||
max_diff_double_single = max(diff_double_single(:));
|
||||
mean_diff_double_single = mean(diff_double_single(:));
|
||||
rel_diff_double_single = max_diff_double_single / max(abs(sig_gpu_double(:)));
|
||||
fprintf('Max absolute difference: %.6e\n', max_diff_double_single);
|
||||
fprintf('Mean absolute difference: %.6e\n', mean_diff_double_single);
|
||||
fprintf('Max relative difference: %.6e\n', rel_diff_double_single);
|
||||
|
||||
% Check tolerances
|
||||
fprintf('\n--- Tolerance Check ---\n');
|
||||
tol_double = 1e-10;
|
||||
tol_single = 1e-5; % Single precision has ~7 significant digits
|
||||
|
||||
if max_diff_cpu_double < tol_double
|
||||
fprintf('✓ CPU vs GPU Double: EQUIVALENT (diff < %.0e)\n', tol_double);
|
||||
else
|
||||
fprintf('✗ CPU vs GPU Double: DIFFER beyond tolerance (%.0e)\n', tol_double);
|
||||
end
|
||||
|
||||
if max_diff_cpu_single < tol_single
|
||||
fprintf('✓ CPU vs GPU Single: ACCEPTABLE (diff < %.0e)\n', tol_single);
|
||||
else
|
||||
fprintf('⚠ CPU vs GPU Single: Precision loss detected (diff = %.2e, tol = %.0e)\n', max_diff_cpu_single, tol_single);
|
||||
end
|
||||
|
||||
% Performance comparison
|
||||
fprintf('\n========== Performance Summary ==========\n');
|
||||
fprintf('CPU Time: %.3f s\n', time_cpu);
|
||||
fprintf('GPU Double Time: %.3f s\n', time_gpu_double);
|
||||
fprintf('GPU Single Time: %.3f s\n', time_gpu_single);
|
||||
fprintf('\n');
|
||||
fprintf('Speedup (GPU Double vs CPU): %.2fx\n', time_cpu/time_gpu_double);
|
||||
fprintf('Speedup (GPU Single vs CPU): %.2fx\n', time_cpu/time_gpu_single);
|
||||
fprintf('Speedup (GPU Single vs GPU Double): %.2fx\n', time_gpu_double/time_gpu_single);
|
||||
|
||||
%% ========== BER Comparison ==========
|
||||
% Process each fiber output through simplified receiver to check if
|
||||
% single-precision affects actual BER performance
|
||||
|
||||
fprintf('\n========== BER Comparison ==========\n');
|
||||
fprintf('Processing signals through receiver chain...\n');
|
||||
|
||||
% Receiver parameters
|
||||
rop = -7; % Received optical power [dBm]
|
||||
len_tr = 4096; % Training length
|
||||
mu_dc = 0.005;
|
||||
mu_ffe = [0.0001 0.0008 0.001];
|
||||
mu_dfe = 0.0004;
|
||||
|
||||
% Helper function to process through receiver and get BER
|
||||
function ber = process_receiver(Opt_sig_fib, l, Symbols, Tx_bits, ...
|
||||
fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode)
|
||||
|
||||
% Demux single channel
|
||||
Opt_sig_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ...
|
||||
"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_fib);
|
||||
|
||||
% ROP amplifier
|
||||
Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
|
||||
"amplification_db",rop).process(Opt_sig_demux{l});
|
||||
|
||||
% Photodiode
|
||||
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ...
|
||||
"nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx);
|
||||
|
||||
% Low-pass filter
|
||||
rx_bwl = 100e9;
|
||||
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ...
|
||||
"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
|
||||
|
||||
% Scope
|
||||
Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
|
||||
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc, ...
|
||||
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth, ...
|
||||
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0, ...
|
||||
"adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig);
|
||||
|
||||
% Resample to 2 sps
|
||||
Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
|
||||
|
||||
% Time sync
|
||||
[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0);
|
||||
Rx_sig = Scpe_cell{1};
|
||||
Rx_sig = Rx_sig.normalize("mode","rms");
|
||||
|
||||
% FFE Equalizer
|
||||
ffe_order = [50, 0, 0];
|
||||
eq_ffe = EQ("Ne",ffe_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
|
||||
"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0);
|
||||
|
||||
ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols{l},Tx_bits{l}, ...
|
||||
"precode_mode",duob_mode, ...
|
||||
'showAnalysis',0, ...
|
||||
"postFFE",[], ...
|
||||
"eth_style_symbol_mapping",0);
|
||||
|
||||
ber = ffe_results.metrics.BER;
|
||||
end
|
||||
|
||||
% Process each mode for channel 1
|
||||
l = 4; % Use first channel for comparison
|
||||
|
||||
fprintf('Processing CPU result...\n');
|
||||
ber_cpu = process_receiver(Opt_sig_cpu, l, Symbols, Tx_bits, ...
|
||||
fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode);
|
||||
|
||||
fprintf('Processing GPU Double result...\n');
|
||||
ber_gpu_double = process_receiver(Opt_sig_gpu_double, l, Symbols, Tx_bits, ...
|
||||
fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode);
|
||||
|
||||
fprintf('Processing GPU Single result...\n');
|
||||
ber_gpu_single = process_receiver(Opt_sig_gpu_single, l, Symbols, Tx_bits, ...
|
||||
fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode);
|
||||
|
||||
% Display BER results
|
||||
fprintf('\n========== BER Results (Channel %d, ROP = %d dBm) ==========\n', l, rop);
|
||||
fprintf('CPU: BER = %.4e\n', ber_cpu);
|
||||
fprintf('GPU Double: BER = %.4e\n', ber_gpu_double);
|
||||
fprintf('GPU Single: BER = %.4e\n', ber_gpu_single);
|
||||
|
||||
fprintf('\n--- BER Comparison ---\n');
|
||||
if ber_cpu == 0 && ber_gpu_double == 0 && ber_gpu_single == 0
|
||||
fprintf('✓ All BERs are zero (no errors detected)\n');
|
||||
else
|
||||
ber_diff_double = abs(ber_cpu - ber_gpu_double);
|
||||
ber_diff_single = abs(ber_cpu - ber_gpu_single);
|
||||
fprintf('|BER_cpu - BER_gpu_double| = %.4e\n', ber_diff_double);
|
||||
fprintf('|BER_cpu - BER_gpu_single| = %.4e\n', ber_diff_single);
|
||||
|
||||
if ber_diff_double < 1e-6 && ber_diff_single < 1e-6
|
||||
fprintf('✓ BER differences are negligible\n');
|
||||
elseif ber_diff_single > ber_diff_double * 10
|
||||
fprintf('⚠ Single precision shows measurable BER impact\n');
|
||||
else
|
||||
fprintf('✓ BER differences within acceptable range\n');
|
||||
end
|
||||
end
|
||||
|
||||
fprintf('\n========== Test Complete ==========\n');
|
||||
@@ -1,265 +0,0 @@
|
||||
|
||||
|
||||
s.wavelengthplan = calcWavelengthPlan(16, 400e9, 1310);
|
||||
N = numel(s.wavelengthplan);
|
||||
link_length = 10;
|
||||
s.pmd = 0.1;%0.1;
|
||||
s.gamma = 0.0023;
|
||||
|
||||
s.M = 4;
|
||||
fsym = 112e9;
|
||||
fdac = 2*fsym;
|
||||
fadc = 120000000000;
|
||||
s.random_key = 1;
|
||||
|
||||
% Laser / s.Modulator
|
||||
vbias_rel = 0.5;
|
||||
u_pi = 4.6;
|
||||
vbias = -vbias_rel*u_pi;
|
||||
laser_linewidth = 0e6;
|
||||
|
||||
% DB Stuff
|
||||
duob_mode = db_mode.no_db;
|
||||
|
||||
rcalpha = 0.05;
|
||||
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);
|
||||
|
||||
s.chirpalpha = 0;
|
||||
|
||||
s.p_launch = 3;
|
||||
s.p = "co";
|
||||
|
||||
switch s.p
|
||||
case "co"
|
||||
pol_rot = 100.*ones(1,N);
|
||||
d_local = 0;
|
||||
case "pair"
|
||||
pol_rot = repmat([100,100,0,0],1,N/4);
|
||||
d_local = 0;
|
||||
case "alt"
|
||||
pol_rot = repmat([100,0,100,0],1,N/4);
|
||||
d_local = 0;
|
||||
case "seg"
|
||||
pol_rot = 100.*ones(1,N);
|
||||
d_local = 3;
|
||||
otherwise
|
||||
error('Unknown fwm_mitigation_technique: %s', string(s.p));
|
||||
end
|
||||
|
||||
f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
|
||||
margin = 25e12; % some THz left and right
|
||||
f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
|
||||
f_nyq = f_span/2;
|
||||
|
||||
kover = 4;
|
||||
upsample_required = f_nyq./(fdac*kover/2);
|
||||
upsample_pow = 2^nextpow2(upsample_required);
|
||||
|
||||
s.f_opt = fdac*kover*upsample_pow;
|
||||
s.f_opt_nyq = s.f_opt/2;
|
||||
|
||||
s.rop = -10:1:0;
|
||||
|
||||
profile on
|
||||
|
||||
%% ---------- TX per channel ----------
|
||||
for l = 1:N
|
||||
|
||||
[Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
|
||||
"fsym",fsym,"M",s.M,"order",17,"useprbs",0, ...
|
||||
"fs_out",fdac, ...
|
||||
"applyclipping",0,"clipfactor",1.5, ...
|
||||
"applypulseform",1,"pulseformer",Pform, ...
|
||||
"randkey",s.random_key+l, ...
|
||||
"mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ...
|
||||
).process();
|
||||
|
||||
Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ...
|
||||
"filterType",filtertypes.gaussian,"active",true);
|
||||
|
||||
El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ...
|
||||
"bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ...
|
||||
"H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);
|
||||
|
||||
% Digi_sig not needed after AWG
|
||||
clear Digi_sig
|
||||
|
||||
% Electrical Driver Amplifier
|
||||
El_sig = El_sig.normalize("mode","oneone");
|
||||
scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
|
||||
El_sig = El_sig .* scaling;
|
||||
|
||||
% E/O Conversion
|
||||
Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ...
|
||||
"lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ...
|
||||
"linewidth",laser_linewidth,"randomkey",s.random_key+l,"alpha",s.chirpalpha).process(El_sig);
|
||||
|
||||
% El_sig not needed after EML
|
||||
clear El_sig
|
||||
|
||||
signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);
|
||||
|
||||
% Eml_out not needed after pol controller
|
||||
clear Eml_out Lp_awg
|
||||
end
|
||||
|
||||
disp('Signal generated for all channels.');
|
||||
|
||||
%% ---------- WDM mux + launch ----------
|
||||
Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ...
|
||||
"lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);
|
||||
|
||||
Opt_sig_wdm.spectrum();
|
||||
|
||||
Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
|
||||
"amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm);
|
||||
|
||||
Opt_sig_wdm_fib = Opt_sig_wdm;
|
||||
|
||||
segment_length = 1;
|
||||
nSegments = link_length/segment_length;
|
||||
if abs(nSegments - round(nSegments)) > 1e-12
|
||||
error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length);
|
||||
end
|
||||
nSegments = round(nSegments);
|
||||
|
||||
zdw = 1310;
|
||||
randomize_D = true;
|
||||
|
||||
% Guard for 0 km: avoid calling getDispersionVector(0,...) if it doesn't support it
|
||||
if nSegments > 0
|
||||
Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key);
|
||||
else
|
||||
Dvec = [];
|
||||
end
|
||||
|
||||
for seg = 1:nSegments
|
||||
|
||||
fprintf('Segment %d/%d \n',seg, nSegments);
|
||||
|
||||
Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ...
|
||||
"beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
|
||||
"gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
|
||||
"SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1,"useGPU",true,"useSingle",1).process(Opt_sig_wdm_fib);
|
||||
|
||||
end
|
||||
|
||||
profile off
|
||||
profile viewer
|
||||
|
||||
Opt_sig_wdm_fib.spectrum();
|
||||
|
||||
%% ========== BER Evaluation ==========
|
||||
fprintf('\n========== BER Evaluation ==========\n');
|
||||
fprintf('Processing signals through receiver chain...\n');
|
||||
|
||||
% Receiver parameters
|
||||
len_tr = 4096; % Training length
|
||||
mu_dc = 0.005;
|
||||
mu_ffe = [0.0001 0.0008 0.001];
|
||||
mu_dfe = 0.0004;
|
||||
|
||||
% Helper function to process through receiver and get BER
|
||||
function ber = process_receiver(Opt_sig_fib, l, Symbols, Tx_bits, ...
|
||||
fdac, kover, upsample_pow, fsym, fadc, rop, s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode)
|
||||
|
||||
% Demux single channel
|
||||
Opt_sig_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ...
|
||||
"fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_fib);
|
||||
|
||||
% ROP amplifier
|
||||
Opt_sig_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
|
||||
"amplification_db",rop).process(Opt_sig_demux{l});
|
||||
|
||||
% Photodiode
|
||||
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ...
|
||||
"nep",1.8e-11,"randomkey",s.random_key+l).process(Opt_sig_rx);
|
||||
|
||||
% Low-pass filter
|
||||
rx_bwl = 100e9;
|
||||
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ...
|
||||
"filterType",filtertypes.butterworth,"active",true).process(PD_sig);
|
||||
|
||||
% Scope
|
||||
Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true);
|
||||
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc, ...
|
||||
"delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth, ...
|
||||
"samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0, ...
|
||||
"adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(PD_sig);
|
||||
|
||||
% Resample to 2 sps
|
||||
Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);
|
||||
|
||||
% Time sync
|
||||
[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols{l}, "fs_ref", fsym, "debug_plots", 0);
|
||||
Rx_sig = Scpe_cell{1};
|
||||
Rx_sig = Rx_sig.normalize("mode","rms");
|
||||
|
||||
% FFE Equalizer
|
||||
ffe_order = [50, 0, 0];
|
||||
eq_ffe = EQ("Ne",ffe_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5, ...
|
||||
"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0);
|
||||
eq_ffe = FFE("epochs_tr",5,"epochs_dd",2,"len_tr",2^13,"mu_dd",6.624e-05,"mu_tr",0.058136,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1);
|
||||
|
||||
ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols{l},Tx_bits{l}, ...
|
||||
"precode_mode",duob_mode, ...
|
||||
'showAnalysis',0, ...
|
||||
"postFFE",[], ...
|
||||
"eth_style_symbol_mapping",0);
|
||||
|
||||
ber = ffe_results.metrics.BER;
|
||||
end
|
||||
|
||||
% Process each ROP value for selected channels using parfor
|
||||
ber_results = zeros(length(s.rop), N);
|
||||
|
||||
% Flatten loop for parfor: iterate over all (ROP, channel) combinations
|
||||
num_rop = length(s.rop);
|
||||
rop_vals = s.rop;
|
||||
ber_flat = zeros(num_rop * N, 1);
|
||||
|
||||
parfor idx = 1:(num_rop * N)
|
||||
% Convert linear index to (ri, l) subscripts
|
||||
ri = ceil(idx / N);
|
||||
l = mod(idx - 1, N) + 1;
|
||||
|
||||
fprintf('ROP %d dBm, Channel %d/%d\n', rop_vals(ri), l, N);
|
||||
ber_flat(idx) = process_receiver(Opt_sig_wdm_fib, l, Symbols, Tx_bits, ...
|
||||
fdac, kover, upsample_pow, fsym, fadc, rop_vals(ri), s, len_tr, mu_dc, mu_ffe, mu_dfe, duob_mode);
|
||||
end
|
||||
|
||||
% Reshape back to [num_rop × N] matrix
|
||||
ber_results = reshape(ber_flat, [N, num_rop]).';
|
||||
|
||||
|
||||
%% Display BER Results
|
||||
fprintf('\n========== BER Results ==========\n');
|
||||
fprintf('ROP [dBm] | ');
|
||||
for l = 1:N
|
||||
fprintf('Ch%d | ', l);
|
||||
end
|
||||
fprintf('\n');
|
||||
|
||||
for ri = 1:length(s.rop)
|
||||
fprintf('%8d | ', s.rop(ri));
|
||||
for l = 1:N
|
||||
fprintf('%.2e | ', ber_results(ri, l));
|
||||
end
|
||||
fprintf('\n');
|
||||
end
|
||||
|
||||
% Plot BER vs ROP
|
||||
figure;
|
||||
semilogy(s.rop, mean(ber_results, 2), '-o', 'LineWidth', 2);
|
||||
hold on;
|
||||
for l = 1:N
|
||||
semilogy(s.rop, ber_results(:, l), '--', 'LineWidth', 1);
|
||||
end
|
||||
hold off;
|
||||
xlabel('ROP [dBm]');
|
||||
ylabel('BER');
|
||||
title('BER vs Received Optical Power (GPU Single Precision)');
|
||||
legend(['Mean', arrayfun(@(x) sprintf('Ch%d', x), 1:N, 'UniformOutput', false)]);
|
||||
grid on;
|
||||
|
||||
fprintf('\n========== Test Complete ==========\n');
|
||||
@@ -1,70 +0,0 @@
|
||||
% 1. Setup Data OUTSIDE the timer
|
||||
d = gpuDevice;
|
||||
N = 10000;
|
||||
|
||||
fprintf('Preparing data...\n');
|
||||
A_cpu_double = rand(N, N); % Create double on CPU
|
||||
A_cpu_single = single(A_cpu_double); % Create single on CPU
|
||||
|
||||
% Warmup run (wakes up the GPU from idle state)
|
||||
A_warm = gpuArray.rand(1000, 1000, 'single');
|
||||
B_warm = A_warm * A_warm;
|
||||
wait(d);
|
||||
|
||||
fprintf('------------------------------------------------\n');
|
||||
|
||||
% TEST 1: Double Precision (The "Slow" way)
|
||||
% We move data to GPU first so we only measure calculation time
|
||||
A_gpu_double = gpuArray(A_cpu_double);
|
||||
wait(d); % Ensure transfer is done before starting timer
|
||||
|
||||
fprintf('Running DOUBLE precision test... ');
|
||||
tic;
|
||||
B_gpu = A_gpu_double * A_gpu_double;
|
||||
wait(d); % FORCE MATLAB TO WAIT FOR GPU
|
||||
time_double = toc;
|
||||
fprintf('Done.\n');
|
||||
fprintf('Double Precision Time: %.4f seconds\n', time_double);
|
||||
|
||||
% TEST 2: Single Precision (The "Fast" way)
|
||||
A_gpu_single = gpuArray(A_cpu_single);
|
||||
wait(d); % Ensure transfer is done
|
||||
|
||||
fprintf('Running SINGLE precision test... ');
|
||||
tic;
|
||||
B_gpu = A_gpu_single * A_gpu_single;
|
||||
wait(d); % FORCE MATLAB TO WAIT FOR GPU
|
||||
time_single = toc;
|
||||
fprintf('Done.\n');
|
||||
fprintf('Single Precision Time: %.4f seconds\n', time_single);
|
||||
|
||||
% Calculate Speedup
|
||||
fprintf('------------------------------------------------\n');
|
||||
fprintf('Speedup Factor using single precision: %.2fx\n', time_double / time_single);
|
||||
|
||||
%
|
||||
% Create 10,000 small matrices (10x10) stacked in a 3D array
|
||||
A_stack = gpuArray.rand(10, 10, 10000, 'single');
|
||||
B_stack = gpuArray.rand(10, 10, 10000, 'single');
|
||||
|
||||
% BAD: Looping (GPU overhead kills you)
|
||||
tic;
|
||||
for i=1:10000
|
||||
C(:,:,i) = A_stack(:,:,i) * B_stack(:,:,i);
|
||||
end
|
||||
wait(d);
|
||||
loop_time = toc;
|
||||
|
||||
% GOOD: Pagefun (Executes all 10,000 mults simultaneously)
|
||||
tic;
|
||||
C_stack = pagefun(@mtimes, A_stack, B_stack);
|
||||
wait(d);
|
||||
pagefun_time = toc;
|
||||
|
||||
fprintf('------------------------------------------------\n');
|
||||
fprintf('Speedup Factor using pagefun: %.2fx\n', loop_time / pagefun_time);
|
||||
|
||||
|
||||
fprintf('Total VRAM: %.2f GB\n', d.TotalMemory / 1e9);
|
||||
fprintf('Available VRAM: %.2f GB\n', d.AvailableMemory / 1e9);
|
||||
fprintf('Usage: %.1f%%\n', 100 * (1 - d.AvailableMemory / d.TotalMemory));
|
||||
Reference in New Issue
Block a user