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

Merge commit '1b554f2d252c0bae558330a96a11c4d8ffa652ac'

Browse Source
This commit is contained in:
sioe
2024-10-20 17:28:37 +02:00
parent 70fd0a44be 1b554f2d25
commit d6beff670c
6 changed files with 279 additions and 125 deletions
Download Patch File Download Diff File Expand all files Collapse all files

134
projects/Lab_2024/lab_modulator_tf_sweep.m
View File

@@ -1,68 +1,136 @@
measurement_name = 'MZM_SD40';
% Initialize a structure to hold parameters
params = struct;
% Define the bias voltage range from 1.2V to 2.8V with 0.01V increments
params.v_bias = 1.2:0.01:2.8;
params.v_bias = 0:0.1:6;
% Create a DataStorage object with the defined parameters
wh = DataStorage(params);
% Add a storage field for output power measurements
wh.addStorage("p_out");
wh.addStorage("i_draw");
wh.addStorage("v_set");
% Display the total number of measurement loops to be executed
disp(['Start Measurement of ', num2str(prod(wh.dim)), ' loops...']);
% Initialize the Optical Attenuator (VOA) with specified settings
voa = OptAtten(...
"active", [1, 0, 0, 0], ... % Activate only the first channel
"active", [0, 1, 0, 0], ... % Activate only the first channel
"value", [0, 0, 0, 0], ... % Set attenuation values to 0 dB
"wavelength", [1310, 1310, 1310, 1310]); % Set the wavelength for each channel
% Set the VOA active channels and attenuation values
voa.set('active', [0, 1, 0, 0], 'value', [0, 0, 0, 0]);
% Initialize the DC Power Supply with specified settings
dcs = DC_supply(...
"active", [1, 1], ... % Activate the first two channels
"voltage", [v_bias, 9]); % Set initial voltages for channels
v_bias_init = wh.parameter.v_bias.values(1);
dcs = DC_supply(...
"active", [1, 0], ... % Activate the first two channels
"voltage", [v_bias_init, 0]); % Set initial voltages for channels
% Set the VOA active channels and attenuation values
voa.set('active', [1, 0, 0, 0], 'value', [0, 0, 0, 0]);
% Loop over each bias voltage value to perform measurements
for v_bias = wh.parameter.v_bias.values
try
% Try to retrieve existing output power measurement to avoid repetition
p_out = wh2.getStoValue('p_out', v_bias);
wh.addValueToStorage(p_out, 'p_out', v_bias);
catch
% If no existing measurement, proceed with the measurement
%%%%% SET Voltages %%%%%%
% Update the DC supply voltage for the current bias voltage
dcs.set("voltage", [v_bias, 9]);
%%%%% Measure VOA %%%%%%
% Read the current values from the VOA
voa.readvals();
% Store the measured output power in the DataStorage object
wh.addValueToStorage(voa.power_state(1), 'p_out', v_bias);
end
%%%%% SET Voltages %%%%%%
% Update the DC supply voltage for the current bias voltage
dcs.set("voltage", [v_bias, 0]);
%%%%% Measure Voltages %%%%%%
[v,c] = dcs.readVals();
% assert(v(1)==v_bias,['Desired voltage is ', num2str(v_bias),'V but currently is ', num2str(v(1)), ' V ']);
%%%%% Measure VOA %%%%%%
% Read the current values from the VOA
voa.readvals();
% Store the measured output power in the DataStorage object
wh.addValueToStorage(voa.power_state(2), 'p_out', v_bias);
wh.addValueToStorage(c(1), 'i_draw', v_bias);
wh.addValueToStorage(v(1), 'v_set', v_bias);
% Retrieve all stored output power measurements up to the current point
p_out = wh.getStoValue('p_out', wh.parameter.v_bias.values);
p_opt_out = wh.getStoValue('p_out', wh.parameter.v_bias.values);
% Plot the measured output power in dBm versus the negative bias voltage
figure(90);
plot(-wh.parameter.v_bias.values(1:numel(p_out)), p_out, 'DisplayName', 'Measured Output Power in dBm');
xlim([-max(params.v_bias), -min(params.v_bias)]); % Set x-axis limits
% Plot the output power converted from dBm to linear scale (Watts)
figure(91);
clf
hold on
plot(wh.parameter.v_bias.values(1:numel(p_opt_out)), db2pow(p_opt_out), 'DisplayName', 'Measured Output Power in Watts','LineWidth',2);
xlim([min(params.v_bias), max(params.v_bias)]); % Set x-axis limits
grid on; % Enable grid for better readability
xlabel('applied voltage in V')
ylabel('measured optical power in dB')
end
save(['C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\HighSpeedMeasurement_2024\Messreihe_TFLN_MZM\',measurement_name],'wh');
p_opt_in = 0; %dbm
% After completing the measurements, retrieve all output power data
p_out = wh.getStoValue('p_out', wh.parameter.v_bias.values);
p_opt_out = wh.getStoValue('p_out', wh.parameter.v_bias.values);
transmission = p_opt_out;
i_draw = wh.getStoValue('i_draw', wh.parameter.v_bias.values);
v_set = wh.getStoValue('v_set', wh.parameter.v_bias.values);
power_draw = i_draw .* v_set .* 1e3; %watt
resistance = v_set ./ i_draw;
i_draw = i_draw.* 1e3;
% Plot the output power converted from dBm to linear scale (Watts)
cols = linspecer(4);
c = 1;
linstyle = '-';
figure(91);
plot(-wh.parameter.v_bias.values(1:numel(p_out)), db2pow(p_out), 'DisplayName', 'Measured Output Power in Watts');
xlim([-max(params.v_bias), -min(params.v_bias)]); % Set x-axis limits
hold on
subplot(2,2,1)
hold on
plot(power_draw, transmission, 'DisplayName', measurement_name,'LineWidth',1,'Marker','none','Color',cols(c,:),'LineStyle',linstyle);
title('MZM Attenuation')
xlim([min(power_draw), max(power_draw)]); % Set x-axis limits
grid on; % Enable grid for better readability
xlabel('DUT Power (mW)');
ylabel('Transmission in dB');
legend
subplot(2,2,2)
hold on
plot(v_set, i_draw, 'DisplayName', measurement_name,'LineWidth',1,'Marker','none','Color',cols(c,:),'LineStyle',linstyle);
title('Voltage vs. Current')
xlim([min(v_set), max(v_set)]); % Set x-axis limits
grid on; % Enable grid for better readability
xlabel('Voltage (mA)');
ylabel('Current (mA)');
legend
subplot(2,2,3)
hold on
plot(power_draw, db2pow(transmission), 'DisplayName', measurement_name,'LineWidth',1,'Marker','none','Color',cols(c,:),'LineStyle',linstyle);
title('Measured Output Power in Watt')
xlim([min(power_draw), max(power_draw)]); % Set x-axis limits
grid on; % Enable grid for better readability
xlabel('DUT Power (mW)');
ylabel('Trabnsmission in %');
legend
subplot(2,2,4)
hold on
plot(power_draw, resistance, 'DisplayName', measurement_name,'LineWidth',1,'Marker','none','Color',cols(c,:),'LineStyle',linstyle);
title('Heater Resistance = U/I')
xlim([min(power_draw), max(power_draw)]); % Set x-axis limits
grid on; % Enable grid for better readability
xlabel('DUT Power (mW)');
ylabel('Resistance in Ohm');
ylim([455, 470])
legend
disp('Measurement Complete!')

164
projects/Lab_2024/lab_sir_sweep.m
View File

@@ -1,13 +1,13 @@
folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\sir_sweep_pam4\';
experiment_name = 'PAM4_DB_encoded_10km_';
folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\sir_sweep_sd40\';
experiment_name = 'PAM4_56_v2';
%%% SIR Sweep for MPI Experiment %%%
params = struct;
params.vbias = [2.45];
params.awg_vpp = [0.25];
params.eq_mode = [4];
params.vbias = [2.5];
params.awg_vpp = [0.35];
params.eq_mode = [2];
params.i_atten = [0:4:40];
wh = DataStorage(params);
@@ -23,13 +23,13 @@ wh.addStorage("signals");
precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\standard_system\";
precomp_fn = "lab_mpi_setup_2";
precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active
precomp_amp_max = 5;
precomp_amp_max = 2;
M = 4;
pn_key = 2;
usemrds = 0;
fsym = 92e9;
fdac = 92e9;
fsym = 68e9;
fdac = 256e9;
awg_vpp = 0.35;
fadc = 160e9;
rrcalpha = 0.05;
@@ -60,21 +60,6 @@ for eq_mode = wh.parameter.eq_mode.values
sprintf('Progress: %d/%d\nEstimated time remaining: %.2f hours\nEstimated time remaining: %.2f hours', ...
loopcnt, looptotal, estimatedTimeRemaining/60/60, estimatedTotalTime/60/60));
try
z
ber = wh2.getStoValue('ber',v_bias,awg_vpp,eq_mode,precomp_amp_max);
assert(~isempty(ber),'err')
rop = wh2.getStoValue('rop',v_bias,awg_vpp,eq_mode,precomp_amp_max);
assert(~isempty(rop))
pd_in = wh2.getStoValue('pd_in',v_bias,awg_vpp,eq_mode,precomp_amp_max);
assert(~isempty(pd_in))
M = wh2.getStoValue('m',v_bias,awg_vpp,eq_mode,precomp_amp_max);
assert(~isempty(M))
catch
switch eq_mode
@@ -93,6 +78,7 @@ for eq_mode = wh.parameter.eq_mode.values
db_coding_approach = 0;
db_precode = db_coding_approach || db_channel_approach;
case 3
ffe_only = 0;
postfilter_approach = 0;
@@ -122,8 +108,9 @@ for eq_mode = wh.parameter.eq_mode.values
%%%%% Construct AWG and Scope Modules %%%%%%
SCP = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc",'GSa_160',"channel",[1,0],"recordLen",2000000,"removeDC",1);
AWG = AwgKeysight("model","M8196A","fdac",fdac,"scaletodac",[1,1,1,1],"skews",[0,0,0,0],"voltages",[0,0,0,awg_vpp]);
A2S = Awg2Scope(AWG,SCP,[0,0,0,1]);
%AWG = AwgKeysight("model","M8196A","fdac",fdac,"scaletodac",[1,1,1,1],"skews",[0,0,0,0],"voltages",[0,0,0,awg_vpp]);
AWG = AwgKeysight("model","M8199B","fdac",fdac,"scaletodac",[1,1],"skews",[0,0],"voltages",[0,awg_vpp]);
A2S = Awg2Scope(AWG,SCP,[0,1,0,0]);
%%%%% Symbol Generation %%%%%%
Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha);
@@ -146,10 +133,12 @@ for eq_mode = wh.parameter.eq_mode.values
%%%%% Resample to DAC rate %%%%%%
Digi_sig = Digi_sig.resample("fs_out",AWG.fdac);
Digi_sig = Filter('filtdegree',5,"f_cutoff",1.1*(fsym/log2(M)),"fs",Digi_sig.fs,"filterType",filtertypes.gaussian).process(Digi_sig);
%%%%% Plot and Save Routine 1 %%%%%%%%%%%%%%%%%%%%%%%%%
Digi_sig.spectrum("displayname","Normal Tx","fignum",10);
if loopcnt == 1
save([folderpath,[experiment_name,'bits'],loop_name],"Bits");
save([folderpath,[experiment_name,'symbols'],loop_name],"Symbols");
@@ -157,7 +146,7 @@ for eq_mode = wh.parameter.eq_mode.values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% AWG --> Scope %%%%%%
[~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig);
[~,Scpe_sig,~,~] = A2S.process("signal2",Digi_sig);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Scpe_sig.spectrum("displayname","Scope PSD","fignum",20);
@@ -194,17 +183,43 @@ for eq_mode = wh.parameter.eq_mode.values
% Scpe_sig.eye(fsym,M,"fignum",40,"displayname",' after Scope');
voa.readvals();
pd_in = voa.power_state(2);
s_pow = voa.power_state(3);
i_pow = voa.power_state(4);
sir = s_pow- i_pow;
%%%%% EQUALIZE %%%%%%
Eq = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",50,"sps",2,"decide",0);
Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[50,7,7],"sps",2,"decide",1);
% Eq = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",[50,7,7],"sps",2,"decide",1);
Eq = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
Eq = EQ("Ne",[50,0,0],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.0,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1);
S = Scpe_sig.signal;
% recursion
N1 = 201;
% Initialize the moving sum for the first window
half_window = (N1 - 1) / 2;
moving_sum = sum(S(1:N1));
% Calculate the first element of R1
cic_filtrd(1:half_window+1) = S(1:half_window+1) - (moving_sum / N1);
% Loop over the signal and apply the recursive moving average subtraction
for n = (half_window+2):(length(S)-half_window)
% Update the moving sum by subtracting the oldest value and adding the new one
moving_sum = moving_sum - S(n-half_window-1) + S(n+half_window);
% Calculate the new value of R1
cic_filtrd(n) = S(n) - (moving_sum / N1);
end
cic_filtrd(n+1:length(S)) = S(n+1:length(S)) - (moving_sum / N1);
Scpe_sig_ = Scpe_sig;
Scpe_sig_.signal = cic_filtrd';
if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -212,7 +227,6 @@ for eq_mode = wh.parameter.eq_mode.values
EQ_sig.plot("fignum",50,"displayname",'After EQ');
Rx_bits = PAMmapper(M,0).demap(EQ_sig);
[~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
@@ -231,46 +245,62 @@ for eq_mode = wh.parameter.eq_mode.values
elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%%
[EQ_sig] = Eq.process(Scpe_sig,Symbols);
[EQ_sig, Noi] = Eq.process(Scpe_sig_,Symbols);
EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1);
% EQ_sig.plot("displayname",'After VNLE','fignum',90,'clear',1);
Noi = EQ_sig-Symbols;
% Quantization is too far from orig. symbols ->
% error psd is quite different
% Sym_ = PAMmapper(M,0).quantize(EQ_sig);
% Noi_ = Sym_-EQ_sig;
% Noi_.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:));
Rx_bits = PAMmapper(M,0).demap(EQ_sig);
[~,errors_bm,ber_ffe_only,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
[~,~,ber_vnle,~] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
nc = 2;
burg_coeff = arburg(Noi.signal,nc);
EQ_sig.normalize('mode','rms').spectrum('displayname','EQ Out','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc,:));
EQ_sig = EQ_sig.filter(burg_coeff,1);
Noi.normalize('mode','rms').spectrum('displayname','Noise PSD optimal','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+1,:));
if 1
Noi.spectrum('displayname','Noise PSD','fignum',123)
[h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs);
h = h/max(abs(h));
hold on
w_ = (w - Noi.fs/2);
plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']);
end
for nc = 1:3
if 0
figure(53);
constellation = unique(Symbols.signal);
received = NaN(numel(constellation),length(Symbols));
for lvl = 1:numel(constellation)
received(lvl,Symbols.signal==constellation(lvl)) = EQ_sig.signal(Symbols.signal==constellation(lvl));
burg_coeff = arburg(Noi.signal,nc);
EQ_sig_filt = EQ_sig.filter(burg_coeff,1);
% EQ_sig.spectrum("displayname","Signal Spectrum after Postfilter","fignum",1234);
EQ_sig_mlse = MLSE("DIR",burg_coeff,"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig_filt);
% EQ_sig.spectrum("displayname","Signal Spectrum after MLSE","fignum",1234);
if 0
cols = linspecer(12);
EQ_sig_filt.normalize('mode','rms').spectrum('displayname','Noise PSD','fignum',1234,'normalizeTo0dB',1,'normalizeToNyquist',1,'color',cols(nc+2,:));
% [h,w] = freqz(1,burg_coeff,length(Noi),"whole",Noi.fs);
% h = h/max(abs(h));
% hold on
% w_ = (w - Noi.fs/2);
% figure(123)
% plot(w_.*1e-9,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']);
[h,w] = freqz(1,burg_coeff,length(Noi),"whole");
h = h/max(abs(h));
hold on
histogram(received(lvl,:),1000,"EdgeAlpha",0);
w_ = (w - pi);
plot(w_,20*log10(fftshift(h)),'DisplayName',['', num2str(nc), ' coefficients for burg alg.']);
end
Rx_bits = PAMmapper(M,0).demap(EQ_sig_mlse);
[~,errors_bm,ber(nc),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% disp(['BER: ',sprintf('%.1E',ber_mlse(i,j)),' - - ROP: ',num2str(patten(i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']);
end
EQ_sig = MLSE("DIR",burg_coeff,"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig);
Rx_bits = PAMmapper(M,0).demap(EQ_sig);
[~,errors_bm,ber,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
disp(['FFE: ',sprintf('%.1E',ber_ffe_only),' -> PF -> MLSE: ',sprintf('%.1E',ber),' dB | PD_in: ',num2str(pd_in),' dBm']);
disp(['FFE: ',sprintf('%.1E',ber_vnle),' -> PF -> MLSE: ',sprintf('%.1E',ber(nc)),' dB | PD_in: ',num2str(pd_in),' dBm']);
elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -301,8 +331,8 @@ for eq_mode = wh.parameter.eq_mode.values
end
end
wh.addValueToStorage(ber_vnle,'vnle',v_bias,awg_vpp,eq_mode,i_atten);
wh.addValueToStorage(ber,'ber',v_bias,awg_vpp,eq_mode,i_atten);
wh.addValueToStorage(pd_in,'pd_in',v_bias,awg_vpp,eq_mode,i_atten);
wh.addValueToStorage(Rx_bits,'signals',v_bias,awg_vpp,eq_mode,i_atten);
@@ -312,13 +342,13 @@ for eq_mode = wh.parameter.eq_mode.values
wh.addValueToStorage(M,'m',v_bias,awg_vpp,eq_mode,i_atten);
showCurrentMeasurement('BER', ber, 'PD in', pd_in, 'PAM',M, 'Vbias', v_bias, 'AWG Vpp', awg_vpp, 'Precomp MaxAmp',precomp_amp_max);
iterationTimes(loopcnt) = toc(iterationStartTime);
averageTimePerIteration = mean(iterationTimes(1:loopcnt));
estimatedTotalTime = averageTimePerIteration * looptotal;
estimatedTimeRemaining = estimatedTotalTime - sum(iterationTimes(1:loopcnt));
%autoArrangeFigures(3,3,2);
end
end
end
@@ -337,18 +367,12 @@ awg_vpp = wh.parameter.awg_vpp.values(1);
eq_mode = wh.parameter.eq_mode.values(1);
bers = wh.getStoValue('ber',v_bias,awg_vpp,eq_mode,i_atten_vals);
sirs = wh.getStoValue('sir',v_bias,awg_vpp,eq_mode,i_atten_vals);
figure(90);
hold on; % Retain the plot so new points can be added without complete redraw
% Plot the data and get the line handle
hLine = plot(i_atten_vals, bers, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", experiment_name);
% Customize the data tips
% Set labels for existing data tip rows
hLine.DataTipTemplate.DataTipRows(1).Label = 'Fsym';
hLine.DataTipTemplate.DataTipRows(2).Label = 'BER';
hLine.DataTipTemplate.DataTipRows(2).Format = '%.2e'; % Format BER as "3e-4"
plot(sirs, bers, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", experiment_name);
% Continue with the rest of your plot settings