diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index edccd25..872b9d3 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -703,7 +703,8 @@ classdef Signal end if options.debug_plots - figure() + figure(121212);clf + subplot(1,2,1); findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend') end @@ -711,7 +712,7 @@ classdef Signal shifts = lags(pkpos); - shifts = shifts(shifts>=0); + % shifts = shifts(shifts>=0); if numel(shifts) > 0 @@ -720,25 +721,21 @@ classdef Signal isFlipped = 1; end - for c = shifts(shifts>=0) + for c = shifts sig = obj.delay(-c,'mode','samples'); sig.signal = sig.signal(1:length(b)).*-isFlipped; S{end+1,1} = sig; end - % - - - %return/keep the sinal with the highest correlation (only within positive shifts) - [~,idx]=max(pks(shifts>=0)); + [~,idx]=max(pks); obj.signal = S{idx}.signal; %put signal with highest corr. to first index in S array swap = S{1}; S{1} = S{idx}; S{idx} = swap; - for c = 1:numel(shifts(shifts>0)) + for c = 1:numel(shifts) S{c}.logbook = []; end @@ -750,11 +747,14 @@ classdef Signal %plot all synced signals and the ref signal if options.debug_plots - figure;hold on; + figure(121212);hold on; + subplot(1,2,2); for i = 1:size(S,1) + hold on plot(S{i}.normalize('mode','oneone').signal(1000:1100),'LineWidth',0.1,'Color',[0.2157 0.4941 0.7216]); plot(b(1000:1100),'LineWidth',1); end + drawnow; end @@ -913,14 +913,12 @@ classdef Signal elseif mode == 1 % generate eye diagram using histogram - maxA = max(sig(100:end-100)); - minA = min(sig(100:end-100)); + maxA = max(sig(100:end-100))*2; + minA = min(sig(100:end-100))*2; % maxA = 0.0015; % minA = 0; - - difference= maxA-minA; data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1; @@ -940,11 +938,9 @@ classdef Signal % beautify cm=flip(cbrewer2("RdYlBu",4096)); % cm=flip(cbrewer2("RdBu",4096)); - cm(1,:) = [1,1,1]; - cm(1,:) = [0,0,0]; + cm(1,:) = [1,1,1]; % set zeros to white => clean background colormap(cm); % colormap('turbo'); - ax.Colormap(1,:) = [1,1,1]; % ax.CLim = [0 50]; if isa(obj,'Opticalsignal') title(['Optical Eye ',options.displayname]) diff --git a/Classes/01_transmit/ChannelFreqResp.m b/Classes/01_transmit/ChannelFreqResp.m index f896f16..0c9fda8 100644 --- a/Classes/01_transmit/ChannelFreqResp.m +++ b/Classes/01_transmit/ChannelFreqResp.m @@ -38,6 +38,9 @@ classdef ChannelFreqResp < handle Nfft df faxis + + fileName + path end methods (Access=public) @@ -194,7 +197,9 @@ classdef ChannelFreqResp < handle obj.H_apply = H_inv; - Target.signal = real((ifft( ( fft(real( Target.signal' )) .* H_inv ) ))); + H_inv = reshape(H_inv,size(Target.signal)); + + Target.signal = real((ifft( ( fft(real( Target.signal )) .* H_inv ) ))); Target.signal = Target.signal'; end @@ -250,8 +255,18 @@ classdef ChannelFreqResp < handle % Havg_smooth = smooth(Havg,10); symaxis = (obj.faxis-(obj.f_ref/2))/1e9; Havg = fftshift(Havg); - %Havg = smooth(Havg); - plot(symaxis, 20*log10(abs(Havg)),'LineWidth',0.5); + Havg = smooth(Havg); + Havg = Havg./max(Havg); + + if ~isempty(obj.fileName) + dpname = obj.fileName; + else + dpname = 'Measured System response'; + end + + plot(symaxis, 20*log10(abs(Havg)),'LineWidth',0.5,'DisplayName',dpname); + + legend('Interpreter','latex') grid on; @@ -368,6 +383,8 @@ classdef ChannelFreqResp < handle end fprintf('Frequency response information successfully loaded from %s\n', fullFileName); + obj.fileName = options.fileName; + obj.path = options.loadPath; end end diff --git a/Classes/01_transmit/PAMsource.m b/Classes/01_transmit/PAMsource.m index ce15418..f1f8e5a 100644 --- a/Classes/01_transmit/PAMsource.m +++ b/Classes/01_transmit/PAMsource.m @@ -11,6 +11,7 @@ classdef PAMsource db_precode db_encode + duobinary_mode mrds_code mrds_blocklength @@ -38,6 +39,7 @@ classdef PAMsource options.db_precode = 0; options.db_encode = 0; + options.duobinary_mode = db_mode.no_db; options.mrds_code = 0; options.mrds_blocklength = 512; @@ -102,8 +104,7 @@ classdef PAMsource para.reset_prms = 0; para.method = 1; - - + data_in = []; global loop; loop = 0; @@ -134,43 +135,52 @@ classdef PAMsource symbols.fs = obj.fsym; %%%%%% Duobinary %%%%%%%%%%% - if obj.db_precode - symbols = Duobinary().precode(symbols); + + switch obj.duobinary_mode + + case db_mode.no_db + + case db_mode.db_precoded + symbols = Duobinary().precode(symbols); + case db_mode.db_encoded + symbols = Duobinary().precode(symbols); + symbols = Duobinary().encode(symbols); end - if obj.db_encode - symbols = Duobinary().encode(symbols); - end - % figure(12);hold on;histogram(symbols.signal,'Normalization','probability'); + %%%%%% MRDS %%%%%%%%%%% if obj.mrds_code symbols = MRDS_coding("blocklength",obj.mrds_blocklength).encode(symbols); end + %%%%%% SAVE MIN-MAX FOR CLIP %%%%%%%%%%% + if obj.applyclipping sym_min = min(symbols.signal); sym_max = max(symbols.signal); end - - %%%%% Pulse-forming %%%%%% if obj.applypulseform %%% MY CODE - digi_sig = obj.pulseformer.process(symbols); + if 0 + digi_sig = obj.pulseformer.process(symbols); - %%% MOVEIT WRAPPER BETA - % symbols.spectrum("fignum",111,"displayname","1) RAW SIGNAL"); - % - % pulsf = Moveit_wrapper('pulsef'); - % pulsf.para.alpharacos = obj.pulseformer.alpha; - % pulsf.para.f_sym = obj.fsym; - % pulsf.para.fs = obj.fsym; - % pulsf.para.pulse = 2; - % - % digi_sig = pulsf.process(symbols); + else + %%% MOVEIT WRAPPER BETA + % symbols.spectrum("fignum",111,"displayname","1) RAW SIGNAL"); + pulsf = Moveit_wrapper('pulsef'); + pulsf.para.alpharacos = obj.pulseformer.alpha; + pulsf.para.f_sym = obj.pulseformer.fsym; + pulsf.para.fs = obj.pulseformer.fdac; + pulsf.para.pulse = 1; + + digi_sig = pulsf.process(symbols); + digi_sig.fs = obj.pulseformer.fdac; + + end % Design raised cosine filter with given order in symbols % digi_sig = symbols; % sps = 4; diff --git a/Classes/01_transmit/Pulseformer.m b/Classes/01_transmit/Pulseformer.m index cbf62f3..4efe336 100644 --- a/Classes/01_transmit/Pulseformer.m +++ b/Classes/01_transmit/Pulseformer.m @@ -6,10 +6,8 @@ classdef Pulseformer end properties(Access=public) - fs + fdac fsym - matched_sps - output_sps pulse pulselength alpha @@ -22,10 +20,8 @@ classdef Pulseformer % Detailed explanation goes here arguments - options.fs double + options.fdac double options.fsym double - options.matched_sps double - options.output_sps double options.pulse pulseform = pulseform.rc options.pulselength double {mustBeInteger} = 32 options.alpha double = 0.05 @@ -54,11 +50,7 @@ classdef Pulseformer signalclass_in = signalclass_in.logbookentry(lbdesc); % write fs to signal - if obj.matched - signalclass_in.fs = obj.fsym .* obj.output_sps;%.* (obj.fdac./obj.fsym); - else - signalclass_in.fs = obj.fs;%.* (obj.fdac./obj.fsym); - end + signalclass_in.fs = obj.fdac;%.* (obj.fdac./obj.fsym); % write to output signalclass_out = signalclass_in; @@ -85,19 +77,15 @@ classdef Pulseformer data_out end - if ~isempty(obj.output_sps) - obj.fsym = obj.fsym.*obj.output_sps; - end - - if ~rem(obj.fs,obj.fsym) + if ~rem(obj.fdac,obj.fsym) %ist ein Vielfaches - sps = obj.fs / obj.fsym; + sps = obj.fdac / obj.fsym; p = sps; q = 1; else %ist kein Vielfaches - p = obj.fsym / gcd(obj.fs, obj.fsym); %upsampling p->->-> - q = obj.fs/ gcd(obj.fs, obj.fsym); %downsampling <-q + p = obj.fsym / gcd(obj.fdac, obj.fsym); %upsampling p->->-> + q = obj.fdac/ gcd(obj.fdac, obj.fsym); %downsampling <-q sps= q; %sps während dem pulse shaping end @@ -117,7 +105,6 @@ classdef Pulseformer end manual_cyclic_convolution = 0; - upsample_filter = 0; upfirdn_convolution = 1; if manual_cyclic_convolution @@ -135,28 +122,16 @@ classdef Pulseformer data_out=ifft( fft(symbolov.') .* repmat( H,size(data_in,1),1 ).' ).'; data_out = circshift(data_out,[0 -(obj.pulselength*sps)]); - if rem(obj.fs,obj.fsym) + if rem(obj.fdac,obj.fsym) data_out = data_out(1:q:end); end end - if upsample_filter - data_out_ = upsample(data_in,p); - - mfOutput = filter(h, 1, data_out_); % Matched filter output - - figure() - hold on - stem(mfOutput(1:1000),'Marker','o','MarkerSize',1,'LineStyle','-','LineWidth',1); - stem(data_out_(1:1000),'Marker','o','MarkerSize',1,'LineStyle','-','LineWidth',1); - - - end - if upfirdn_convolution %Apply Filter using Matlab build in fctn. + data_out_ = upfirdn(data_in,h,p,q); %cut signal, which is longer due to fir filter @@ -183,7 +158,7 @@ classdef Pulseformer % scale = max(max([abs(real(data_out)) abs(imag(data_out))])); %find max value from real and imag part % data_out = data_out./scale; - data_out = data_out'; + % data_out = data_out'; %Check output integrity if abs(round(p/q * length(data_in)) - length(data_out)) > 4 diff --git a/Classes/05_Lab/AwgKeysight.m b/Classes/05_Lab/AwgKeysight.m index e0ef03a..ec82430 100644 --- a/Classes/05_Lab/AwgKeysight.m +++ b/Classes/05_Lab/AwgKeysight.m @@ -69,16 +69,16 @@ classdef AwgKeysight obj.sampleMemorySize = 512000; case awg_model.M8199A obj.numChannels = 4; - obj.waveformGranularity = 512; - obj.sampleMemorySize = 1024000; %to check + obj.waveformGranularity = 256; + obj.sampleMemorySize = 512000; %see documentation "strg+f: granularity" case awg_model.M8199A_ILV obj.numChannels = 2; - obj.waveformGranularity = 512; - obj.sampleMemorySize = 1024000; %to check + obj.waveformGranularity = 512; %see documentation "strg+f: granularity" + obj.sampleMemorySize = 1024000; case awg_model.M8199B obj.numChannels = 2; - obj.waveformGranularity = 512; - obj.sampleMemorySize = 1024000; %to check + obj.waveformGranularity = 512; %see documentation "strg+f: granularity" + obj.sampleMemorySize = 1024000; end end @@ -175,7 +175,10 @@ classdef AwgKeysight % The waveform granularity is e.g. 128. This means that the waveform length must be a multiple of this granularity. The minimum waveform length is 128 samples. if rem(length(signal{s}),obj.waveformGranularity) ~= 0 - signal{s} = [signal{s} signal{s}(:,1:(obj.waveformGranularity-rem(length(signal{s}),obj.waveformGranularity)))]; + + fill_length = (obj.waveformGranularity-rem(length(signal{s}),obj.waveformGranularity)); + segment__fill_samples = signal{s}(:,1:fill_length); + signal{s} = [signal{s} segment__fill_samples]; debug = 0; if debug @@ -196,7 +199,7 @@ classdef AwgKeysight siglen = min(length(signal{s}),obj.sampleMemorySize); signal{s} = signal{s}(1:siglen); if siglen==obj.sampleMemorySize - warning(['Signal ',num2str(s),' was truncated to ',num2str(obj.sampleMemorySize),' to fit into AWG sample memory']) + warning(['Signal ',num2str(s),' too long was truncated to ',num2str(obj.sampleMemorySize),' to fit into AWG sample memory! Reduce Num Symbols or Sampling Rate!']) end @@ -362,6 +365,7 @@ classdef AwgKeysight writeline(visaobj,char(command)); writeline(visaobj,':syst:err?'); + err_aw=readline(visaobj); assert(strfind(err_aw,"No error"),['AWG System Error? after command: ', char(command)]) diff --git a/Classes/05_Lab/OptAtten.m b/Classes/05_Lab/OptAtten.m index 32cf1b9..c9051df 100644 --- a/Classes/05_Lab/OptAtten.m +++ b/Classes/05_Lab/OptAtten.m @@ -53,7 +53,29 @@ classdef OptAtten < handle try %connect to device - v = visadev('TCPIP::134.245.243.248::INSTR'); + ipAddress = 'TCPIP::134.245.243.248::INSTR'; + maxRetries = 5; + pauseTime = 2; % seconds between retries + + v = []; + success = false; + + for attempt = 1:maxRetries + try + v = visadev(ipAddress); + success = true; + + break; % Exit loop if successful + catch ME + fprintf('Connection attempt %d failed: %s\n', attempt, ME.message); + pause(pauseTime); + end + end + + if ~success + error(' Could not connect to device after %d attempts.', maxRetries); + end + debug = 0; if debug @@ -315,7 +337,8 @@ classdef OptAtten < handle p3 = obj.power_state(3); p4 = obj.power_state(4); - if options.silent == 0 + delete(v); + if options.silent ~= 0 fprintf("\n CH1: %.2f dB | CH2: %.2f dB | CH3: %.2f dB | CH4: %.2f dB\n \n",p1, p2, p3, p4); end end diff --git a/Classes/05_Lab/Thor_PDFA.m b/Classes/05_Lab/Thor_PDFA.m index b0912c3..fefad8f 100644 --- a/Classes/05_Lab/Thor_PDFA.m +++ b/Classes/05_Lab/Thor_PDFA.m @@ -32,14 +32,30 @@ classdef Thor_PDFA < handle end end + obj.getStatus(); + end function o = connectSerial(obj) - % Connect to the PDFA - o = serialport(obj.serialport_number, 9600); - configureTerminator(o, "CR", "CR"); % Set the terminator to carriage return (CR) - flush(o); + try + % Connect to the PDFA + o = serialport(obj.serialport_number, 9600); + configureTerminator(o, "CR", "CR"); % Set the terminator to carriage return (CR) + flush(o); + catch ME + disp(['Error connecting to PDFA: ', ME.message]); + % List available serial ports + availablePorts = serialportlist("available"); + + if isempty(availablePorts) + disp('No available serial ports found.'); + else + disp('Available serial ports:'); + disp(availablePorts); + end + + end end diff --git a/Classes/DataBaseHandler/DBHandler.m b/Classes/DataBaseHandler/DBHandler.m index d3bc066..6f90ca7 100644 --- a/Classes/DataBaseHandler/DBHandler.m +++ b/Classes/DataBaseHandler/DBHandler.m @@ -217,12 +217,31 @@ classdef DBHandler < handle end - % Append the new row to the database table - try - sqlwrite(obj.conn, tableName, newRow); - % disp(['Successfully appended new row to the table ', tableName]); - catch e - error('Failed to append to the table %s: %s', tableName, e.message); + % Parameters for retry logic + maxRetries = 5; + basePause = 0.05; % seconds + + attempt = 0; + success = false; + + while ~success && attempt <= maxRetries + try + sqlwrite(obj.conn, tableName, newRow); + success = true; % Write successful + catch e + if contains(e.message, 'database is locked') || contains(e.message, 'cannot rollback transaction') + attempt = attempt + 1; + pauseTime = basePause * (1 + rand()); % Add random jitter to reduce collision chance + fprintf('Database locked. Retry %d/%d after %.3f seconds...\n', attempt, maxRetries, pauseTime); + pause(pauseTime); + else + error('Failed to append to the table %s: %s', tableName, e.message); + end + end + end + + if ~success + error('Failed to append to table %s after %d retries due to database lock.', tableName, maxRetries); end % Retrieve the measurement_id of the newly inserted row for linking other tables @@ -539,27 +558,38 @@ classdef DBHandler < handle end end + % --- Adaptive FROM Clause --- - % Use "Runs" as the main table and add LEFT JOINs for every other table in obj.tables - % (except "sqlite_sequence") that has a run_id field. mainTable = 'Runs'; fromClause = ['FROM ', mainTable, ' ']; + + % Prepare join buffers + normalJoins = ''; + equalizerJoin = ''; + tableNamesAll = fieldnames(obj.tables); for t = 1:numel(tableNamesAll) tableName = tableNamesAll{t}; + if strcmpi(tableName, mainTable) || strcmpi(tableName, 'sqlite_sequence') continue; end if isfield(obj.tables.(tableName), 'run_id') - % most tables are directly linked to runs table - fromClause = [fromClause, 'LEFT JOIN ', tableName, ' ON ', mainTable, '.run_id = ', tableName, '.run_id ']; + % Direct join to Runs + normalJoins = [normalJoins, 'LEFT JOIN ', tableName, ' ON ', mainTable, '.run_id = ', tableName, '.run_id ']; elseif isfield(obj.tables.(tableName), 'eq_id') - % equalizer is only linked to results table - fromClause = [fromClause, 'LEFT JOIN ', tableName, ' ON ', 'Results', '.eqParam_id = ', tableName, '.eq_id ']; + % Equalizer depends on Results, collect this join separately + equalizerJoin = ['LEFT JOIN ', tableName, ' ON Results.eqParam_id = ', tableName, '.eq_id ']; end end + % Combine joins: normal joins first, Equalizer last + fromClause = [fromClause, normalJoins, equalizerJoin]; + + + + % --- WHERE Clause Construction --- baseQuery = [selectClause, ' ', fromClause, 'WHERE ']; filterClauses = []; diff --git a/Classes/Moveit_wrapper.m b/Classes/Moveit_wrapper.m index 96dce6a..621ab67 100644 --- a/Classes/Moveit_wrapper.m +++ b/Classes/Moveit_wrapper.m @@ -49,29 +49,30 @@ classdef Moveit_wrapper < handle isSignalClass = 0; if isa(signal_in,"Signal") - signal_out = signal_in; + signalclass_out = signal_in; signal_in = signal_in.signal; isSignalClass = 1; end - % signal_in = signal_in.'; + signal_in = signal_in.'; % INIT MOVEIT global loop; loop = 0; - signal_out = obj.moveit_init(signal_in); + signal_init = obj.moveit_init(signal_in); % RUN MOVEIT global loop; loop = 1; - signal_out = obj.moveit_init(signal_out); + signal_out = obj.moveit_simulation_computation(signal_in); if isSignalClass % append to logbook + signalclass_out.signal = signal_out; + signal_out=signalclass_out; lbdesc = ['Logbookentry']; signal_out = signal_out.logbookentry(lbdesc); - signal_out.signal = signal_out; end end diff --git a/Functions/EQ_structures/vnle_postfilter_mlse.m b/Functions/EQ_structures/vnle_postfilter_mlse.m index 35ea356..70c95af 100644 --- a/Functions/EQ_structures/vnle_postfilter_mlse.m +++ b/Functions/EQ_structures/vnle_postfilter_mlse.m @@ -19,7 +19,9 @@ end [eq_signal_sd,eq_noise] = eq_.process(rx_signal,tx_symbols); if ~isempty(options.postFFE) + tic [eq_signal_sd,eq_noise] = options.postFFE.process(eq_signal_sd,tx_symbols); + toc end eq_signal_hd = PAMmapper(M,0).quantize(eq_signal_sd); @@ -237,29 +239,32 @@ if options.showAnalysis % fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl); - fprintf('VNLE BER: %.2e \n',ber_vnle); - - fprintf('MLSE BER: %.2e \n',ber_mlse); + % fprintf('VNLE BER: %.2e \n',ber_vnle); + % + % fprintf('MLSE BER: %.2e \n',ber_mlse); + figure(336);clf; showEQNoisePSD(eq_noise,"fignum",336,"displayname",'Residual Noise after VNLE','postfilter_taps',pf_.coefficients); + % figure(337);clf; + % rx_signal.spectrum("normalizeTo0dB",1,"fignum",337,"displayname",'Rx Signal'); - rx_signal.spectrum("normalizeTo0dB",1,"fignum",337,"displayname",'Rx Signal'); + % figure(338);clf; + % showEQcoefficients('n1',eq_.e,'n2',eq_.e2,'n3',eq_.e3,"displayname",'Coefficients','fignum',338); - tx_symbols.spectrum("normalizeTo0dB",1,"fignum",337,'displayname','Tx Signal'); + if ~isempty(options.postFFE) + showEQcoefficients('n1',options.postFFE.e,"displayname",'Coefficients','fignum',338); + end - showLevelHistogram(eq_signal_sd,tx_symbols) - % showLevelHistogram(mlse_sig_sd,tx_symbols) + showEQfilter(eq_.e,eq_signal_sd.fs.*2); - showEQcoefficients('n1',eq_.e,'n2',eq_.e2,'n3',eq_.e3,"displayname",'Coefficients'); + figure(340);clf; + eq_signal_sd.eye(eq_signal_sd.fs,M,"fignum",340); - showEQNoiseSNR(tx_symbols,eq_noise,"displayname",'vnle snr','fignum',101); + figure(341);clf; + showLevelHistogram(eq_signal_sd,tx_symbols,"fignum",341); - %%% EQ SNR Spectrum %230 - %snr - - - % showErrorBurstCount(eq_signal_sd,tx_symbols) + % autoArrangeFigures(3,3,2) end diff --git a/Functions/EQ_visuals/showEQcoefficients.m b/Functions/EQ_visuals/showEQcoefficients.m index cb747b8..a844155 100644 --- a/Functions/EQ_visuals/showEQcoefficients.m +++ b/Functions/EQ_visuals/showEQcoefficients.m @@ -53,8 +53,8 @@ function showEQcoefficients(options) for i = 1:numSubplots subplot(1, numSubplots, i); - stem(coeffs{i}, 'Color', options.color, 'LineWidth', 1, ... - 'Marker', '.', 'MarkerSize', 10); + stem(coeffs{i}, 'Color', options.color, 'LineWidth', 0.1, ... + 'Marker', '.', 'MarkerSize', 1); title(titles{i}); ylim([-1, 1]); % Set y-axis limits to [-1, 1] grid on; diff --git a/Functions/EQ_visuals/showEQfilter.m b/Functions/EQ_visuals/showEQfilter.m index 3ef6732..1d13022 100644 --- a/Functions/EQ_visuals/showEQfilter.m +++ b/Functions/EQ_visuals/showEQfilter.m @@ -14,7 +14,7 @@ function showEQfilter(coefficients,fs) H = H(1:half_nfft); % Plot the magnitude and phase responses. - figure; + figure(339); % Magnitude response (in dB) subplot(2,1,1); diff --git a/Functions/Lab_helper/holdAndShowValue.m b/Functions/Lab_helper/holdAndShowValue.m index e45c588..1740aaa 100644 --- a/Functions/Lab_helper/holdAndShowValue.m +++ b/Functions/Lab_helper/holdAndShowValue.m @@ -15,7 +15,7 @@ function holdAndShowValue() 'ButtonPushedFcn', @(src, event) closeWindow()); % Initialize the timer - updateTimer = timer('ExecutionMode', 'fixedRate', 'Period', 1, ... + updateTimer = timer('ExecutionMode', 'fixedRate', 'Period', 0.1, ... 'TimerFcn', @(src, event) updatePowerValues()); % Start the timer diff --git a/Functions/Lab_helper/loadFreqResp.m b/Functions/Lab_helper/loadFreqResp.m index 108dc5d..019ea07 100644 --- a/Functions/Lab_helper/loadFreqResp.m +++ b/Functions/Lab_helper/loadFreqResp.m @@ -1,9 +1,34 @@ +% Define the precomp path +precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\"; -precomp_path = "D:\kiel_dsp\imdd_simulation\projects\HighSpeed_ETH"; -precomp_filename = "lab_high_speed_nachtschichtmzm"; -% precomp_filename = "lab_high_speed"; -% precomp_filename = "lab_high_speed"; +% Step 1: Find all valid files (assume .mat files for ChannelFreqResp) +fileList = dir(fullfile(precomp_path, '*.mat')); +fileNames = {fileList.name}; -freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',92e9); -freqresp.load('loadPath',precomp_path,'fileName',precomp_filename); -freqresp.plot(); +% Check if any files found +if isempty(fileNames) + error('No valid .mat files found in the specified folder.'); +end + +% Step 2: Let user select which files to plot +[selection, ok] = listdlg( ... + 'PromptString', 'Select Transfer Functions to Plot:', ... + 'SelectionMode', 'multiple', ... + 'ListString', fileNames); + +% Check if user pressed OK +if ~ok || isempty(selection) + disp('No files selected. Exiting.'); + return; +end + +% Step 3: Loop over selected files and plot +for idx = selection + precomp_filename = erase(fileNames{idx}, '.mat'); % Remove .mat extension + + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',92e9); + freqresp.load('loadPath', precomp_path, 'fileName', precomp_filename); + + fprintf('Plotting: %s\n', precomp_filename); + freqresp.plot(); +end diff --git a/projects/ECOC_2025/dsp_standalone.m b/projects/ECOC_2025/dsp_standalone.m index 371ac40..6424498 100644 --- a/projects/ECOC_2025/dsp_standalone.m +++ b/projects/ECOC_2025/dsp_standalone.m @@ -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; diff --git a/projects/ECOC_2025/estimate_linewidth.m b/projects/ECOC_2025/estimate_linewidth.m new file mode 100644 index 0000000..d97f3c7 --- /dev/null +++ b/projects/ECOC_2025/estimate_linewidth.m @@ -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; diff --git a/projects/ECOC_2025/measure_run_script.m b/projects/ECOC_2025/measure_run_script.m index be6ef25..15405f8 100644 --- a/projects/ECOC_2025/measure_run_script.m +++ b/projects/ECOC_2025/measure_run_script.m @@ -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 diff --git a/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m index 9747fa2..e737228 100644 --- a/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m +++ b/projects/ECOC_2025/plots_from_database/plot_mpi_trial.m @@ -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 + 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 diff --git a/projects/run_loop.m b/projects/run_loop.m new file mode 100644 index 0000000..35a7178 --- /dev/null +++ b/projects/run_loop.m @@ -0,0 +1,101 @@ +% loop_records.m + +% === Prepare base configuration === +% databasePath = 'C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\ECOC_2025\'; +databasePath = 'Z:\2025\ECOC Silas\'; +database_name = 'ecoc2025_loops.db'; +db = DBHandler("pathToDB", [databasePath, database_name]); +conf = db.tables.Configurations; +confPrev = db.tables.Configurations; + +awg_upload_required = 0; + +% Define vector parameters directly in conf: +conf.unique_elab_id = "20250408-842b0a3078172b48dd032795226fbe683190afc4"; +conf.symbolrate = 112e9; +conf.bitrate = conf.symbolrate .* floor(log2(6)*10)/10; +conf.pam_level = 4; +conf.db_mode = db_mode.no_db; +conf.pulsef_alpha = 0.2; +conf.v_bias = 2.5; +conf.v_awg = 0.95; +conf.precomp_amp = -64; %-64 +conf.wavelength = 1310; +conf.laser_power = 11; +conf.fiber_length = 0; +conf.pd_in_desired = 9; +conf.is_mpi = 1; +conf.signal_attenuation = 0; +conf.interference_path_length = 1000; %m +conf.interference_attenuation = [40,20:-1:0]; +conf.pam_source = []; + +if conf.is_mpi + current_folder = ['mpi_',num2str(conf.interference_path_length),'m_pam',num2str(conf.pam_level)]; +else + current_folder = 'testing'; +end + +% Init empty loop id +loop_id = []; + +% === Detect vector fields automatically === +paramNames = fieldnames(conf); +vectorParams = struct(); +scalarConf = conf; +parallel_dsp = false; +max_occurences = 1; +for i = 1:numel(paramNames) + thisName = paramNames{i}; + thisValue = conf.(thisName); + if isnumeric(thisValue) && numel(thisValue) > 1 + vectorParams.(thisName) = thisValue; + parallel_dsp = true; + max_occurences = 10; + else + scalarConf.(thisName) = thisValue; + end +end + +vectorFieldNames = fieldnames(vectorParams); + +if isempty(vectorFieldNames) + % No vector parameters, single run + fprintf('Running single measurement (no parameter sweep)\n'); + confRequest = scalarConf; + confRequest.bitrate = confRequest.symbolrate * floor(log2(confRequest.pam_level) * 10) / 10; + measure_run_script +else + % Vector parameters found, run parameter sweep + vectorValues = struct2cell(vectorParams); + nParams = numel(vectorValues); + grid = cell(1, nParams); + [grid{:}] = ndgrid(vectorValues{:}); + totalRuns = numel(grid{1}); + fprintf('Total runs: %d\n', totalRuns); + + for runIdx = 1:totalRuns + confRequest = scalarConf; + paramDisplayCell = cell(1, nParams); + + for vIdx = 1:nParams + fieldName = vectorFieldNames{vIdx}; + paramValue = grid{vIdx}(runIdx); + confRequest.(fieldName) = paramValue; + if isnumeric(paramValue) + paramDisplayCell{vIdx} = sprintf('%s = %.4g', fieldName, paramValue); + else + paramDisplayCell{vIdx} = sprintf('%s = %s', fieldName, char(paramValue)); + end + end + + fprintf('Running measurement: %s\n', strjoin(paramDisplayCell, ', ')); + confRequest.bitrate = confRequest.symbolrate * floor(log2(confRequest.pam_level) * 10) / 10; + % measure_run("config",confRun,"parallel",parallel_dsp,"max_occurences",max_occurences); + measure_run_script; + end +end + + +disp('- LOOP FINISHED - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') +fprintf('\n')