ecoc measurements

This commit is contained in:
Silas Labor Zizou
2025-04-13 16:29:12 +02:00
parent 00f1c557c0
commit 0236103b13
19 changed files with 709 additions and 359 deletions

View File

@@ -1,9 +1,9 @@
basePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\';
database_name='ecoc2025.db';
savePath = 'Z:\2025\ECOC Silas\ecoc_2025\';
databasePath = 'Z:\2025\ECOC Silas\';
database_name = 'ecoc2025_loops.db';
db = DBHandler("pathToDB", [databasePath, database_name]);
num_occ = 1;
run_par = false;
run_id = 562;

View File

@@ -0,0 +1,44 @@
% Parameters
reclen = 133169152;
Scp_long = ScopeKeysight("model","DSAZ634A",'autoscale',0,"fadc","GSa_160","channel",[0,0,1,0],"recordLen",reclen,"removeDC",1,"extRef",1);
Scpe_sig_raw = Scp_long.read("channel",[0,0,1,0]);
Scpe_sig_raw = Scpe_sig_raw{3};
signal = Scpe_sig_raw.signal;
Fs = Scpe_sig_raw.fs; % Sampling rate [Hz]
N = length(signal); % Number of samples
f_center = 20e9; % Frequency of interest (~20 GHz)
% Create time vector
t = (0:N-1)' / Fs;
% Mix down to baseband
signal_mix = signal .* exp(-1j * 2 * pi * f_center * t);
% Low-pass filter (optional but recommended)
% Design a filter with bandwidth ~ linewidth * safety margin
BW = 1000e6; % 100 MHz bandwidth for example
% Decimate (reduce sampling rate)
decimation_factor = 10;%floor(Fs / (2 * BW));
signal_decim = downsample(signal_mix, decimation_factor);
Fs_decim = Fs / decimation_factor;
% FFT of downmixed signal
Nfft = 2^nextpow2(length(signal_decim) * 4); % Zero-pad
spectrum = fft(signal_decim .* blackmanharris(length(signal_decim)), Nfft);
freq = (0:Nfft-1) * (Fs_decim / Nfft);
% Only positive frequencies
halfIdx = 1:floor(Nfft/2);
spectrum_half = abs(spectrum(halfIdx));
freq_half = freq(halfIdx);
% Plot
figure;
plot(freq_half / 1e6, 20*log10(spectrum_half));
xlabel('Frequency (MHz)');
ylabel('Magnitude (dB)');
title('Zoomed Spectrum around 20 GHz');
grid on;

View File

@@ -2,16 +2,11 @@
savePath = 'Z:\2025\ECOC Silas\ecoc_2025\';
current_folder = 'precomp_optimization_friday';
fullFolderPath = fullfile(savePath, current_folder);
if ~exist(fullFolderPath, 'dir')
mkdir(fullFolderPath);
end
basePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\';
database_name='ecoc2025.db';
db = DBHandler("pathToDB",[basePath,database_name]);
precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\";
precomp_fn = "precomp_1km_1mm_cable_70ghz_pd_shf_t850_2p8_bias_shot2";
precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active
@@ -24,9 +19,6 @@ referenceFields = fieldnames(db.tables.Configurations);
assert(isequal(fieldnames(conf), referenceFields), 'Fieldnames do not match the reference structure!');
currentTime = datetime('now', 'Format', 'yyyyMMdd_HHmmss');
timeStr = char(currentTime);
filename = sprintf('%s_PAM_%d_R_%d',timeStr,conf.pam_level,conf.symbolrate.*1e-9);
%% Lab Automation
@@ -38,7 +30,6 @@ if run_lab_automation
dcs.set("voltage",[conf.v_bias, 5]);
[v_bias_meas,i_bias_meas]=dcs.readVals();
end
% Laser
if confPrev.laser_power ~= conf.laser_power
laser = Exfo_laser("mainframe_channel",1,"safety_mode",0,"connection_id",'10','lab_interface','gpib');
@@ -55,7 +46,7 @@ if run_lab_automation
pdfa = Thor_PDFA("safety_mode",0,"serialport_number",'COM15');
% Construct AWG and Scope Modules %%%%%%
fdac = 256e9;
fdac = 224e9;
fadc = 160e9;
Scp = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc","GSa_160","channel",[0,0,1,0],"recordLen",5000000,"removeDC",1,"extRef",1);
Awg = AwgKeysight("model","M8199A_ILV","fdac",fdac,"scaletodac",[1,1],"skews",[0,0],"voltages",[0,conf.v_awg]);
@@ -65,108 +56,151 @@ end
%% Signal generation and transmission
if 0% awg_upload_required || isempty(loop_id)
%%%%% Symbol Generation %%%%%%
Pform = Pulseformer("fsym",conf.symbolrate,"fdac",8*conf.symbolrate,"pulse","rc","pulselength",16,"alpha",conf.pulsef_alpha);
Pamsource = PAMsource(...
"fsym",conf.symbolrate,"M",conf.pam_level,"order",19,"useprbs",1,...
"fs_out",fdac,...
"applyclipping",0,...
"clipfactor",4,...
"applypulseform",1,...
"pulseformer",Pform,...
"randkey",1,...
"duobinary_mode", conf.db_mode);
conf.pam_source = Pamsource;
confPrev = conf;
[Digi_sig,Symbols,Bits] = Pamsource.process();
%%%%% Precompensation Routine %%%%%%
if precomp_mode == 1 % measure channel
precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',fdac);
Digi_sig = precomp_est.buildOFDM();
elseif precomp_mode == 2 % apply precomp
precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs);
Digi_sig = precomp_est.precomp(Digi_sig,'maxampdb',conf.precomp_amp,'loadPath',precomp_path,'fileName',precomp_fn);
end
% Digi_sig.spectrum("displayname",'Tx Spectrum','fignum',10,'normalizeTo0dB',0);
%%%%% Resample to DAC rate %%%%%%
Digi_sig = Digi_sig.resample("fs_out",Awg.fdac);
% [~,~,Scpe_sig_raw,~] = A2S.process("signal2",Digi_sig);
%%%%% Symbol Generation %%%%%%
Pform = Pulseformer("fsym",conf.symbolrate,"fdac",8*conf.symbolrate,"pulse","rc","pulselength",16,"alpha",conf.pulsef_alpha);
Pamsource = PAMsource(...
"fsym",conf.symbolrate,"M",conf.pam_level,"order",19,"useprbs",1,...
"fs_out",fdac,...
"applyclipping",0,...
"clipfactor",4,...
"applypulseform",1,...
"pulseformer",Pform,...
"randkey",1,...
"duobinary_mode", conf.db_mode);
conf.pam_source = Pamsource;
upload_required = ~isequal(Pamsource, confPrev.pam_source);
confPrev = conf;
[Digi_sig,Symbols,Bits] = Pamsource.process();
%%%%% Precompensation Routine %%%%%%
if precomp_mode == 1 % measure channel
precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',fdac);
Digi_sig = precomp_est.buildOFDM();
upload_required = 1;
elseif precomp_mode == 2 % apply precomp
precomp_est = ChannelFreqResp("Nacq",2048,"Navg",100,"Ncp",63,'f_ref',Digi_sig.fs);
Digi_sig = precomp_est.precomp(Digi_sig,'maxampdb',conf.precomp_amp,'loadPath',precomp_path,'fileName',precomp_fn);
end
Digi_sig.spectrum("displayname",'Tx Spectrum','fignum',10,'normalizeTo0dB',0);
%%%%% Resample to DAC rate %%%%%%
Digi_sig = Digi_sig.resample("fs_out",Awg.fdac);
% [~,~,Scpe_sig_raw,~] = A2S.process("signal2",Digi_sig);
if upload_required || confPrev.v_awg ~= conf.v_awg
Awg.upload("signal2",Digi_sig);
end
Scpe_sig_raw = Scp.read("channel",[0,0,1,0]);
Scpe_sig_raw = Scpe_sig_raw{[0,0,1,0]==1};
%%%%% Precompensation Routine %%%%%%
% Precompensation Routine (optional)
if precomp_mode == 1
% Scpe_sig_resampled = Scpe_sig_raw.resample("fs_in",fadc,"fs_out",2*fsym);
precomp_est.estimate(Scpe_sig_raw,"save",true,"savePath",precomp_path,"fileName",precomp_fn);
Scpe_sig_raw = Scp.read("channel",[0,0,1,0]);
Scpe_sig_raw = Scpe_sig_raw{[0,0,1,0]==1};
precomp_est.estimate(Scpe_sig_raw, "save", true, "savePath", precomp_path, "fileName", precomp_fn);
precomp_est.plot();
return;
return; % End routine if precomp mode is active
end
Scpe_sig_raw.spectrum("displayname",'Rx Spectrum','fignum',10,'normalizeTo0dB',0);
Scpe_sig_raw.plot("clear",1,"displayname",'Rx Signal','fignum',11);
%% Save static TX data (only once)
currentTime = datetime('now', 'Format', 'yyyyMMdd_HHmmss');
timeStr = char(currentTime);
base_filename = sprintf('%s_PAM_%d_R_%d', timeStr, conf.pam_level, round(conf.symbolrate.*1e-9));
%% %%%%% Save Routine %%%%%%%%%%%%%%%%%%%%%%%%%
tx_bits_path = [filesep, current_folder, filesep, base_filename, '_bits'];
save([savePath, tx_bits_path], "Bits");
tx_bits_path=[filesep,current_folder,filesep,filename,'_bits'];
save([savePath,tx_bits_path],"Bits");
tx_symbols_path=[filesep,current_folder,filesep,filename,'_symbols'];
save([savePath,tx_symbols_path],"Symbols");
tx_signal_path=[filesep,current_folder,filesep,filename,'_signal'];
save([savePath,tx_signal_path],"Digi_sig");
rx_raw_path = [filesep,current_folder,filesep,filename,'_rx_signal_raw'];
save([savePath,rx_raw_path],"Scpe_sig_raw");
tx_symbols_path = [filesep, current_folder, filesep, base_filename, '_symbols'];
save([savePath, tx_symbols_path], "Symbols");
tx_signal_path = [filesep, current_folder, filesep, base_filename, '_signal'];
save([savePath, tx_signal_path], "Digi_sig");
n_recording = 2; % Or any number you want
for recIdx = 1:n_recording
fprintf('Recording %d / %d\n', recIdx, n_recording);
%% Acquire Scope Signal
Scpe_sig_raw = Scp.read("channel",[0,0,1,0]);
Scpe_sig_raw = Scpe_sig_raw{[0,0,1,0]==1};
% Precompensation Routine (optional)
if precomp_mode == 1
precomp_est.estimate(Scpe_sig_raw, "save", true, "savePath", precomp_path, "fileName", precomp_fn);
precomp_est.plot();
return; % End routine if precomp mode is active
end
% Plot spectrum and time domain (optional)
% Scpe_sig_raw.spectrum("displayname", 'Rx Spectrum', 'fignum', 10, 'normalizeTo0dB', 0);
Scpe_sig_raw.plot("clear", 1, "displayname", 'Rx Signal', 'fignum', 11);
%% Save Routine
currentTime = datetime('now', 'Format', 'yyyyMMdd_HHmmss');
timeStr = char(currentTime);
filename = sprintf('%s_PAM_%d_R_%d_rec%02d', timeStr, conf.pam_level, round(conf.symbolrate.*1e-9), recIdx); % Added recIdx
rx_raw_path = [filesep, current_folder, filesep, filename, '_rx_signal_raw'];
save([savePath, rx_raw_path], "Scpe_sig_raw");
% Table 1: Append to Runs
newRun = db.tables.Runs;
newRun.run_id = NaN;
newRun.date_of_run = datetime(currentTime, 'InputFormat', 'yyyyMMdd_HHmmss');
newRun.tx_bits_path = tx_bits_path;
newRun.tx_symbols_path = tx_symbols_path;
newRun.rx_sync_path = [];
newRun.rx_raw_path = rx_raw_path;
newRun.filename = filename;
newRun.tx_signal_path = tx_signal_path;
run_id = db.appendToTable('Runs', newRun);
% Check if loop_id is already created
if isempty(loop_id)
newLoop = db.tables.LoopControl;
newLoop.loop_id = NaN;
newLoop.date_of_loop = datetime('now', 'Format', 'yyyy-MM-dd HH:mm:ss');
newLoop.description = 'Automated parameter sweep';
loop_id = db.appendToTable('LoopControl', newLoop);
% Table 2: Append to Configurations
conf.configuration_id = NaN;
conf.run_id = run_id;
db.appendToTable('Configurations', conf);
fprintf('LoopControl entry created: loop_id = %d\n', loop_id);
end
% Table 3: Append to Measurements
meas = db.tables.Measurements;
meas.measurement_id = NaN; % Auto-increment, leave empty
meas.run_id = run_id;
meas.power_laser = laser.cur_power;
meas.power_rop = [];
meas.power_pd_in = voa.power_state(2);
meas.power_mpi_interference = voa.power_state(4);
meas.power_mpi_signal = voa.power_state(3);
meas.voa_class = voa;
meas.pdfa_class = pdfa;
meas.laser_class = laser;
assert(isequal(fieldnames(meas), fieldnames(db.tables.Measurements)), 'Fieldnames of "meas" do not match the reference structure!');
db.appendToTable('Measurements', meas);
%% Append to Runs Table
newRun = db.tables.Runs;
newRun.run_id = NaN;
newRun.loop_id = loop_id;
newRun.date_of_run = datetime(currentTime, 'InputFormat', 'yyyyMMdd_HHmmss');
newRun.tx_bits_path = tx_bits_path;
newRun.tx_symbols_path = tx_symbols_path;
newRun.rx_sync_path = []; % Leave empty for now
newRun.rx_raw_path = rx_raw_path;
newRun.filename = filename;
newRun.tx_signal_path = tx_signal_path;
%% %%%%% DSP Routine %%%%%%%%%%%%%%%%%%%%%%%%%
run_id = db.appendToTable('Runs', newRun);
[out, ~] = submit_dsp(run_id, basePath, database_name, savePath,"parallel",parallel_dsp,'max_occurences',max_occurences);
%% Append to Configurations Table
conf.configuration_id = NaN;
conf.run_id = run_id;
db.appendToTable('Configurations', conf);
%% Append to Measurements Table
meas = db.tables.Measurements;
meas.measurement_id = NaN;
meas.run_id = run_id;
meas.power_laser = laser.cur_power;
meas.power_rop = [];
meas.power_pd_in = voa.power_state(2);
meas.power_mpi_interference = voa.power_state(4);
meas.power_mpi_signal = voa.power_state(3);
meas.voa_class = voa;
meas.pdfa_class = pdfa;
meas.laser_class = laser;
% Ensure consistency
assert(isequal(fieldnames(meas), fieldnames(db.tables.Measurements)), 'Fieldnames of "meas" do not match the reference structure!');
db.appendToTable('Measurements', meas);
%% Submit DSP Routine
% [out, ~] = submit_dsp(run_id, basePath, database_name, savePath, ...
% "parallel", parallel_dsp, "max_occurences", max_occurences);
end

View File

@@ -1,7 +1,7 @@
basePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\';
database_name = 'ecoc2025.db';
database_name = 'ecoc2025_loops.db';
database = DBHandler("pathToDB", [basePath, database_name]);
filterParams = database.tables;
@@ -10,14 +10,14 @@ filterParams.Configurations = struct( ...
'fiber_length', 0, ...
'db_mode', '"no_db"', ...
'interference_attenuation', [], ...
'interference_path_length', 0.06, ...
'interference_path_length', 50, ...
'is_mpi', 1, ...
'pam_level', 4, ...
'wavelength', 1310, ...
'precomp_amp', [], ...
'signal_attenuation', [], ...
'v_awg', 0.9, ...
'v_bias', 2.8 ...
'v_awg', 0.95, ...
'v_bias', 2.5 ...
);
% filterParams.EqualizerParameters.diff_precode = int32(db_mode.db_encoded);
@@ -29,12 +29,12 @@ selectedFields = {'Configurations.run_id' 'Runs.date_of_run' 'Runs.rx_raw_path'
'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'};
[dataTable,sql_query] = database.queryDB(filterParams, selectedFields);
dataTable.SIR = dataTable.power_mpi_signal - dataTable.power_mpi_interference;
dataTable.SIR = round(-6 - dataTable.power_mpi_interference);
dataTable = cleanUpTable(dataTable);
% Filter by time
startTime = datetime('2025-04-11 14:40:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
stopTime = datetime('2025-04-11 15:30:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
startTime = datetime('2025-04-12 18:00:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
stopTime = datetime('2025-04-13 19:30:00', 'InputFormat', 'yyyy-MM-dd HH:mm:ss');
dataTable.date_of_run = datetime(dataTable.date_of_run, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
dataTable.date_of_processing = datetime(dataTable.date_of_processing, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
dataTable = dataTable(dataTable.date_of_processing > startTime, :);
@@ -46,21 +46,23 @@ x_var = 'SIR';
loop_var = 'eq_id';
fixedVars = {'eq_id',x_var};
% dataTableGrpd_mean = groupIt(fixedVars,dataTable);
[dataTable, outliersTable] = removeGroupOutliers(dataTable, fixedVars, y_var);
dataTableGrpd_mean = groupIt(fixedVars, dataTable, @mean);
dataTableGrpd_min = groupIt(fixedVars, dataTable, @min);
dataTableGrpd_max = groupIt(fixedVars, dataTable, @max);
plotRealizations = 1;
% Create a new figure
mkr = '*';
mkr = '.';
figure();
hold on
unique_loop_var = unique(dataTable.(loop_var));
cols = linspecer(numel(unique_loop_var)); % Ensure color count matches
for i = 1:numel(unique_loop_var)
for i = 1:2%1:numel(unique_loop_var)
% Prepare filtered data for this loop variable
loopValue = unique_loop_var(i);
@@ -139,7 +141,7 @@ ylabel(y_var);
title([x_var, ' vs. ', y_var]);
if y_var == 'BER'
yline(3.8e-3, 'LineWidth', 1, 'LineStyle', '--', 'HandleVisibility', 'off');
yline(4e-4, 'LineWidth', 1, 'LineStyle', '--', 'HandleVisibility', 'off');
ylim([1e-5, 0.1]);
end
@@ -163,57 +165,58 @@ function resultTable = groupIt(fixedVars, dataTable, aggregationFunction)
% Output:
% resultTable - Grouped and aggregated table
% Group data
[G, groupKeys] = findgroups(dataTable(:, fixedVars));
% Prepare aggregation
varNames = dataTable.Properties.VariableNames;
nVars = numel(varNames);
aggData = cell(height(groupKeys), nVars);
groupCount = zeros(height(groupKeys), 1); % Store number of rows in each group
% Group data
[G, groupKeys] = findgroups(dataTable(:, fixedVars));
% Loop over groups
for i = 1:height(groupKeys)
idx = (G == i); % Logical index for group i
groupCount(i) = sum(idx); % Count rows in group
% Prepare aggregation
varNames = dataTable.Properties.VariableNames;
nVars = numel(varNames);
aggData = cell(height(groupKeys), nVars);
groupCount = zeros(height(groupKeys), 1); % Store number of rows in each group
% Loop over each variable
for j = 1:nVars
colData = dataTable.(varNames{j});
% Loop over groups
for i = 1:height(groupKeys)
idx = (G == i); % Logical index for group i
groupCount(i) = sum(idx); % Count rows in group
if isnumeric(colData)
% Numeric: apply aggregation function (skip empty groups safely)
if any(idx)
aggData{i, j} = aggregationFunction(colData(idx));
% Loop over each variable
for j = 1:nVars
colData = dataTable.(varNames{j});
if isnumeric(colData)
% Numeric: apply aggregation function (skip empty groups safely)
if any(idx)
% aggData{i, j} = rmoutliers(double(colData(idx)));
aggData{i, j} = aggregationFunction(colData(idx));
else
aggData{i, j} = NaN;
end
else
% Non-numeric: take first non-empty value
if iscell(colData)
nonEmptyIdx = find(idx & ~cellfun(@isempty, colData), 1);
if ~isempty(nonEmptyIdx)
aggData{i, j} = colData{nonEmptyIdx};
else
aggData{i, j} = NaN;
aggData{i, j} = [];
end
else
% Non-numeric: take first non-empty value
if iscell(colData)
nonEmptyIdx = find(idx & ~cellfun(@isempty, colData), 1);
if ~isempty(nonEmptyIdx)
aggData{i, j} = colData{nonEmptyIdx};
else
aggData{i, j} = [];
end
nonEmptyIdx = find(idx, 1);
if ~isempty(nonEmptyIdx)
aggData{i, j} = colData(nonEmptyIdx);
else
nonEmptyIdx = find(idx, 1);
if ~isempty(nonEmptyIdx)
aggData{i, j} = colData(nonEmptyIdx);
else
aggData{i, j} = [];
end
aggData{i, j} = [];
end
end
end
end
end
% Convert aggregated data to table
resultTable = cell2table(aggData, 'VariableNames', varNames);
% Convert aggregated data to table
resultTable = cell2table(aggData, 'VariableNames', varNames);
% Add group size as new column
resultTable.nRows = groupCount;
% Add group size as new column
resultTable.nRows = groupCount;
end
@@ -221,7 +224,7 @@ end
function addDatatips(sc, varargin)
% addDatatips Adds custom data tip rows to a scatter plot.
%
% addDatatips(sc, pair1, pair2, ...) adds one or more custom rows to the
% addDatatips(sc, pair1, pair2, ...) adds one or more custom rows to the
% data tip display of the scatter plot identified by sc.
%
% Each pair should be provided as a 1x2 cell array: {label, value}.
@@ -233,27 +236,27 @@ function addDatatips(sc, varargin)
% pair_one = {'Attenuation', attenuationVector};
% addDatatips(sc, pair_one);
numPoints = numel(sc.XData);
for k = 1:length(varargin)
pair = varargin{k};
if ~iscell(pair) || numel(pair) ~= 2
error('Each pair must be a 1x2 cell array: {label, value}.');
end
label = pair{1};
value = pair{2};
% If value is a vector, ensure its length is either 1 or equal to the number of scatter points.
if isvector(value) && numel(value) ~= 1 && numel(value) ~= numPoints
error('The vector for "%s" must be a scalar or have %d elements matching the scatter data points.', label, numPoints);
end
% Create a new data tip row using the provided label and vector.
newRow = dataTipTextRow(label, value);
sc.DataTipTemplate.DataTipRows(end+1) = newRow;
numPoints = numel(sc.XData);
for k = 1:length(varargin)
pair = varargin{k};
if ~iscell(pair) || numel(pair) ~= 2
error('Each pair must be a 1x2 cell array: {label, value}.');
end
label = pair{1};
value = pair{2};
% If value is a vector, ensure its length is either 1 or equal to the number of scatter points.
if isvector(value) && numel(value) ~= 1 && numel(value) ~= numPoints
error('The vector for "%s" must be a scalar or have %d elements matching the scatter data points.', label, numPoints);
end
% Create a new data tip row using the provided label and vector.
newRow = dataTipTextRow(label, value);
sc.DataTipTemplate.DataTipRows(end+1) = newRow;
end
end
@@ -274,69 +277,135 @@ function cleanedTable = cleanUpTable(inputTable)
% Output:
% cleanedTable - Cleaned MATLAB table with proper numeric types
cleanedTable = inputTable;
varNames = cleanedTable.Properties.VariableNames;
for i = 1:numel(varNames)
col = cleanedTable.(varNames{i});
% Case 1: If it's a cell array (likely mixed strings/struct)
if iscell(col)
% Convert struct 'NaN' entries to string 'NaN'
col = cellfun(@(x) convertStructToString(x), col, 'UniformOutput', false);
% Try to convert string numbers to actual numbers
numericCol = str2double(col);
if all(isnan(numericCol) == strcmpi(col, 'NaN') | cellfun(@isempty, col))
% If conversion is successful (NaNs correspond to 'NaN' strings), use it
cleanedTable.(varNames{i}) = numericCol;
else
% Else, try to convert to datetime
try
cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
catch
% If it fails, leave as cell array of strings
cleanedTable.(varNames{i}) = string(col);
end
end
% Case 2: If it's already a string array
elseif isstring(col)
numericCol = str2double(col);
if all(isnan(numericCol) == strcmpi(col, "NaN"))
cleanedTable.(varNames{i}) = numericCol;
else
% Try convert to datetime
try
cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
catch
% Leave as string
end
end
% Case 3: If it's already numeric, keep as is
elseif isnumeric(col)
continue;
% Case 4: If it's datetime, keep as is
elseif isdatetime(col)
continue;
cleanedTable = inputTable;
varNames = cleanedTable.Properties.VariableNames;
for i = 1:numel(varNames)
col = cleanedTable.(varNames{i});
% Case 1: If it's a cell array (likely mixed strings/struct)
if iscell(col)
% Convert struct 'NaN' entries to string 'NaN'
col = cellfun(@(x) convertStructToString(x), col, 'UniformOutput', false);
% Try to convert string numbers to actual numbers
numericCol = str2double(col);
if all(isnan(numericCol) == strcmpi(col, 'NaN') | cellfun(@isempty, col))
% If conversion is successful (NaNs correspond to 'NaN' strings), use it
cleanedTable.(varNames{i}) = numericCol;
else
% Catch-all for unexpected types, convert to string
cleanedTable.(varNames{i}) = string(col);
% Else, try to convert to datetime
try
cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
catch
% If it fails, leave as cell array of strings
cleanedTable.(varNames{i}) = string(col);
end
end
% Case 2: If it's already a string array
elseif isstring(col)
numericCol = str2double(col);
if all(isnan(numericCol) == strcmpi(col, "NaN"))
cleanedTable.(varNames{i}) = numericCol;
else
% Try convert to datetime
try
cleanedTable.(varNames{i}) = datetime(col, 'InputFormat', 'yyyy-MM-dd HH:mm:ss.SSSSSS');
catch
% Leave as string
end
end
% Case 3: If it's already numeric, keep as is
elseif isnumeric(col)
continue;
% Case 4: If it's datetime, keep as is
elseif isdatetime(col)
continue;
else
% Catch-all for unexpected types, convert to string
cleanedTable.(varNames{i}) = string(col);
end
end
end
function out = convertStructToString(x)
% Helper function to convert struct NaN to string 'NaN'
if isstruct(x)
out = "NaN";
elseif isstring(x) || ischar(x)
out = string(x);
else
out = x;
end
if isstruct(x)
out = "NaN";
elseif isstring(x) || ischar(x)
out = string(x);
else
out = x;
end
end
function [cleanedTable, outliersTable] = removeGroupOutliers(dataTable, fixedVars, y_var)
% Group the data
[G, groupKeys] = findgroups(dataTable(:, fixedVars));
% Initialize logical index to keep rows
keepIdx = true(height(dataTable), 1);
% Prepare storage for outliers
outlierRecords = [];
% Loop over each group
for groupIdx = 1:height(groupKeys)
% Find indices of current group
groupRows = (G == groupIdx);
% Extract y-values of this group
y_values = dataTable.(y_var)(groupRows);
% Skip groups with fewer than 3 points (optional)
if sum(groupRows) < 3
continue;
end
% Detect outliers
outlierMask = isoutlier(y_values, 'quartiles');
% If any outliers found, collect their data
if any(outlierMask)
groupData = dataTable(groupRows, :);
% Prepare table for current group outliers
outlierGroupTable = groupData(outlierMask, :);
% Add group key values for traceability
for k = 1:numel(fixedVars)
outlierGroupTable.(['Group_', fixedVars{k}]) = repmat(groupKeys{groupIdx, k}, height(outlierGroupTable), 1);
end
% Append to collection
outlierRecords = [outlierRecords; outlierGroupTable]; %#ok<AGROW>
end
% Mark outliers for removal
groupRowIdx = find(groupRows);
keepIdx(groupRowIdx(outlierMask)) = false;
end
% Apply mask to dataTable
cleanedTable = dataTable(keepIdx, :);
% Prepare output: if no outliers, return empty table
if isempty(outlierRecords)
outliersTable = table();
else
outliersTable = outlierRecords;
end
nRemoved = sum(~keepIdx);
nTotalOriginal = height(dataTable) + nRemoved;
percentageRemoved = (nRemoved / nTotalOriginal) * 100;
fprintf('Removed %d outliers from the data table (%.2f%% of total %d entries).\n', ...
nRemoved, percentageRemoved, nTotalOriginal);
end