From 6b0f9de118e830d0da4e81897341fc0cd7cffff2 Mon Sep 17 00:00:00 2001 From: Silas Labor Zizou Date: Thu, 10 Oct 2024 11:07:59 +0200 Subject: [PATCH 1/2] Lab changes --- Classes/00_signals/Signal.m | 129 +++++++-- Classes/01_transmit/ChannelFreqResp.m | 9 +- Classes/01_transmit/PAMmapper.m | 7 +- Classes/01_transmit/PAMsource.m | 1 + Classes/02_etc/Amplifier.m | 4 +- Classes/04_DSP/CIC_filter.m | 59 ++++ .../04_DSP/Equalizer/FFE_adaptive_decision.m | 24 +- Classes/05_Lab/Awg2Scope.m | 109 ++++++++ Classes/05_Lab/AwgKeysight.m | 29 +- Classes/05_Lab/DC_supply.m | 5 +- Classes/05_Lab/OptAtten.m | 5 +- Classes/05_Lab/ScopeKeysight.m | 12 +- Classes/Warehouse_class/classes/DataStorage.m | 8 + .../Warehouse_class/classes/DataStorage2.m | 114 ++++++++ Classes/Warehouse_class/classes/Parameter2.m | 51 ++++ .../classes/minimalExample_gen2.m | 45 ++++ Datatypes/scope_fadc.m | 2 +- .../400G_FTN_setups/imdd_dsp_approaches.m | 7 +- projects/Lab_2024/lab_db_precode.m | 253 ++++++++++++++++++ projects/Lab_2024/lab_minimal.m | 96 ++++--- projects/Lab_2024/runEQ.m | 32 ++- 21 files changed, 901 insertions(+), 100 deletions(-) create mode 100644 Classes/04_DSP/CIC_filter.m create mode 100644 Classes/05_Lab/Awg2Scope.m create mode 100644 Classes/Warehouse_class/classes/DataStorage2.m create mode 100644 Classes/Warehouse_class/classes/Parameter2.m create mode 100644 Classes/Warehouse_class/classes/minimalExample_gen2.m create mode 100644 projects/Lab_2024/lab_db_precode.m diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index e955679..8224f4b 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -6,6 +6,10 @@ classdef Signal signal logbook fs + + gitSHA + gitStatus + gitPatch end methods @@ -17,18 +21,26 @@ classdef Signal options.fs = []; end - obj.signal = signal; - obj.signal = obj.signal; - obj.fs = options.fs; - + obj.signal = signal; + obj.signal = obj.signal; + obj.fs = options.fs; + + [~,obj.gitSHA] = system('git rev-parse HEAD'); + [~,obj.gitStatus] = system('git status --porcelain'); + [~,obj.gitPatch] = system('git diff'); + + %%% Stuff for Logbook %%% SignalType = []; TimeStamp = []; Length = []; SignalPower = []; Nase = []; + SignalCopy = []; + ModifierName = []; + ModifierCopy= {}; Description = []; - obj.logbook = table(SignalType,TimeStamp,Length,SignalPower,Nase,Description); + obj.logbook = table(SignalType,TimeStamp,Length,SignalPower,Nase,SignalCopy,ModifierName, ModifierCopy, Description); end @@ -131,9 +143,13 @@ classdef Signal options.fignum options.displayname = []; options.timeframe = 0; + options.clear = 0; end figure(options.fignum); % If figure does not exist, create new figure + if options.clear + clf + end % 2) Plot into the figure handle found or created in one t = (0:length(obj.signal)-1) / obj.fs; % time vector @@ -227,24 +243,58 @@ classdef Signal %% Write Logbook Entry function obj = logbookentry(obj,varargin) - if nargin > 1 + if nargin == 2 Description = varargin{1}; + CallingModifier = evalin('caller','obj'); + elseif nargin == 3 + Description = varargin{1}; + CallingModifier = varargin{2}; else Description = ""; + CallingModifier = evalin('caller','obj'); end + CallingModifierStruct = obj.objToStructFilteredRecursive(CallingModifier); + SignalType = [string(class(obj))]; TimeStamp = [(datetime('now','TimeZone','local','Format','HH:mm:ss'))]; Length = num2str(obj.length, ['%' sprintf('.%df', 0)]);%[obj.length]; SignalPower = [obj.power]; Nase = [0]; + SignalCopy = obj.signal; + ModifierName = class(CallingModifier); + ModifierCopy = {CallingModifierStruct}; - cell = {SignalType , TimeStamp , Length , SignalPower(1) , Nase, Description}; + cell = {SignalType , TimeStamp , Length , SignalPower(1) , Nase, SignalCopy, ModifierName, ModifierCopy, Description}; - obj.logbook = [obj.logbook;cell]; + obj.logbook = [obj.logbook; cell]; end + function s = objToStructFilteredRecursive(~,obj) + % Convert the object to a structure using 'struct' and catch warnings + warnState = warning('off', 'MATLAB:structOnObject'); + s = struct(obj); % Convert to struct + warning(warnState); % Restore previous warning state + + % Get all field names of the struct + fields = fieldnames(s); + + % Loop over each field and handle filtering + for i = 1:numel(fields) + fieldData = s.(fields{i}); + + if isstruct(fieldData) % If the field is a struct, call recursively + s.(fields{i}) = obj.objToStructFilteredRecursive(fieldData); + elseif numel(fieldData) > 1000 % Remove field if it has more than 1000 elements + %s = rmfield(s, fields{i}); + s.(fields{i}) = []; + elseif isa(fieldData,'table') + s = rmfield(s, fields{i}); + end + end + end + %% Resample Signal function obj = resample(obj,options) @@ -260,13 +310,23 @@ classdef Signal warning('The signals fs is different from the given fs_in while it should be the same.'); end - obj.signal = resample(obj.signal,options.fs_out,options.fs_in,options.n,options.beta); + if options.fs_in == options.fs_out - desc = ['resample signal from ', num2str(options.fs_in*1e-9), ' GHz to ', num2str(options.fs_out*1e-9), ' GHz' ]; + desc = ['No need to resample signal from ', num2str(options.fs_in*1e-9), ' GHz to ', num2str(options.fs_out*1e-9), ' GHz' ]; + + obj = obj.logbookentry(desc,obj); - obj = obj.logbookentry(desc); + else - obj.fs = options.fs_out; + obj.signal = resample(obj.signal,options.fs_out,options.fs_in,options.n,options.beta); + + desc = ['resample signal from ', num2str(options.fs_in*1e-9), ' GHz to ', num2str(options.fs_out*1e-9), ' GHz' ]; + + obj = obj.logbookentry(desc,obj); + + obj.fs = options.fs_out; + + end end @@ -284,7 +344,7 @@ classdef Signal N = 2^(nextpow2(length(obj.signal))-8); [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean"); - normalize = 1; + normalize = 0; if normalize p_lin = p_lin./ max(p_lin); p_dbm = 10*log10(p_lin); %dB to dBm in case of "power" @@ -295,6 +355,7 @@ classdef Signal end figure(options.fignum); % If figure does not exist, create new figure + ax = gca; hold on plot(w.*1e-9,p_dbm,'DisplayName',options.displayname,'LineWidth',1); xlabel("Frequency in GHz"); @@ -304,11 +365,11 @@ classdef Signal edgetick = 2^(nextpow2(obj.fs*1e-9)); % xticks([-edgetick:16:edgetick]); xlim([100*round( min(w.*1e-9)/100,1)-10,100*round( max(w.*1e-9)/100,1)+10]) - ylim([100*round( min(p_dbm)/100,1)-3,100*round( max(p_dbm)/100,1)+3]); + ylim([min(floor( min(p_dbm))-3 , ax.YLim(1)), max(ceil( max(p_dbm) )+(3), ax.YLim(2))]); yticks([-200:10:10]); grid on grid minor - legend + legend('Interpreter','none'); end @@ -451,13 +512,21 @@ classdef Signal [pks,pkpos] = findpeaks(co./max(co),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum); shifts = lags(pkpos); - %Cut occurences of ref signal from signal + %Cut occurences of ref signal from signal (only positive shifts) S = {}; - for c = shifts + for c = shifts(shifts>0) sig = obj.delay(-c,'mode','samples'); sig.signal = sig.signal(1:length(b)); S{end+1,1} = sig; end + + %return/keep the sinal with the highest correlation (only within positive shifts) + [~,idx]=max(pks(shifts>0)); + obj.signal = S{idx}.signal; + + for c = 1:numel(shifts(shifts>0)) + S{c}.logbook = []; + end %plot all synced signals and the ref signal debug = 0; @@ -469,14 +538,16 @@ classdef Signal end end - %return the sinal with the highest correlation... - [~,idx]=max(pks); - obj.signal = S{idx}.signal; + end function obj = filter(obj,a,b) + + lbdesc = ['Filtering signal with H = a: ',num2str(a),' / b: ',num2str(b)]; + obj = obj.logbookentry(lbdesc,obj); + obj.signal = filter(a,b,obj.signal); end @@ -560,7 +631,15 @@ classdef Signal end - function eye(obj,fsym,M) + function eye(obj,fsym,M,options) + + arguments + obj + fsym + M + options.fignum = 100; + options.displayname = ""; + end mode = 1; @@ -587,7 +666,7 @@ classdef Signal eye_mat = reshape(x(1:end-mod(length(x),histpoints_horizontal)),histpoints_horizontal,floor(length(x)/histpoints_horizontal)); %% reshape signal into 256 rows each row has the histogram(eye data of all symbols) - figure(922) + figure(options.fignum) clf if mode == 2 % generate "intuitive eye diagram" by drawing lines on top over @@ -633,19 +712,19 @@ classdef Signal colormap(cbrewer2("Blues",4096)); if isa(obj,'Opticalsignal') - title("Optical Eye") + title(['Optical Eye ',options.displayname]) ylabel("Power in mW"); y_tickstring = string(linspace(maxA.*1e3,minA.*1e3,16)); min_ = min(abs(obj.signal(100:end-100)).^2); max_ = abs(max(obj.signal(100:end-100)).^2); elseif isa(obj,'Electricalsignal') - title("Electrical Eye") + title(['Electrical Eye ',options.displayname]) ylabel("Voltage in V"); y_tickstring = string(linspace(maxA,minA,16)); min_ = min(obj.signal(100:end-100)); max_ = abs(max(obj.signal(100:end-100))); else - title("Digital Eye") + title(['Digital Eye ',options.displayname]) ylabel("Digital Signal Amplitude"); y_tickstring = string(linspace(maxA,minA,16)); min_ = min(obj.signal(100:end-100)); diff --git a/Classes/01_transmit/ChannelFreqResp.m b/Classes/01_transmit/ChannelFreqResp.m index b8da77c..75960f3 100644 --- a/Classes/01_transmit/ChannelFreqResp.m +++ b/Classes/01_transmit/ChannelFreqResp.m @@ -184,9 +184,12 @@ classdef ChannelFreqResp < handle end % iH(1) is DC ---> iH(end) is High Freq. - H_inv = [iH(1) iH fliplr(conj(iH)) conj(iH(1))]; - H_inv = [iH(1) iH 0 fliplr(conj(iH))]; - + + if mod(length(Target.signal),2) %ungerade + H_inv = [iH(1) iH iH(end) fliplr(conj(iH)) conj(iH(1))]; + else + H_inv = [iH(1) iH 0 fliplr(conj(iH))]; + end % H_inv = [iH(1) iH iH(end) fliplr(conj(iH)) conj(iH(1))]; obj.H_apply = H_inv; diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index 9414d2e..64ef02f 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -32,8 +32,8 @@ classdef PAMmapper signal_in.signal = obj.map_(signal_in.signal); % signal_in = signal_in.normalize("mode","rms"); - - signal_in = signal_in.logbookentry(); + lbdesc = ['Map bat stream to PAM ',num2str(obj.M),' symbols']; + signal_in = signal_in.logbookentry(lbdesc,obj); out = signal_in; else out = signal_in; @@ -43,7 +43,8 @@ classdef PAMmapper function signalclass_out = demap(obj,signalclass_in) signalclass_in.signal = obj.demap_(signalclass_in.signal); - signalclass_in = signalclass_in.logbookentry(); + lbdesc = ['Demap PAM ',num2str(obj.M),' symbols to bit stream']; + signalclass_in = signalclass_in.logbookentry(lbdesc,obj); signalclass_out = signalclass_in; end diff --git a/Classes/01_transmit/PAMsource.m b/Classes/01_transmit/PAMsource.m index a4e1723..099b6be 100644 --- a/Classes/01_transmit/PAMsource.m +++ b/Classes/01_transmit/PAMsource.m @@ -93,6 +93,7 @@ classdef PAMsource end bits = Informationsignal(bitpattern); + bits = bits.logbookentry(['Generate bit stream with size: ', num2str(size(bitpattern))]); symbols = PAMmapper(obj.M,0).map(bits); symbols.fs = obj.fsym; diff --git a/Classes/02_etc/Amplifier.m b/Classes/02_etc/Amplifier.m index 4289061..31324e2 100644 --- a/Classes/02_etc/Amplifier.m +++ b/Classes/02_etc/Amplifier.m @@ -40,8 +40,8 @@ classdef Amplifier signalclass_in = obj.process_(signalclass_in); % append to logbook - lbdesc = ['Amp ']; - signalclass_in = signalclass_in.logbookentry(lbdesc); + lbdesc = ['Optical Amplifier ']; + signalclass_in = signalclass_in.logbookentry(lbdesc,obj); % write to output signalclass_out = signalclass_in; diff --git a/Classes/04_DSP/CIC_filter.m b/Classes/04_DSP/CIC_filter.m new file mode 100644 index 0000000..97c4db0 --- /dev/null +++ b/Classes/04_DSP/CIC_filter.m @@ -0,0 +1,59 @@ +% Moving Average filter +N = 7; +xn = sin(2*pi*[0:.1:10]); +hn = ones(1,N); +y1n = conv(xn,hn) .* 1/N; + +% transfer function of Moving Average filter +figure() +hF = fft(hn,1024); +plot([-512:511]/1024, abs(fftshift(hF))); +xlabel('Normalized frequency') +ylabel('Amplitude') +title('frequency response of Moving average filter') + +% Implementing Cascaded Integrator Comb filter with the +% comb section following the integrator stage +N = 10; +delayBuffer = zeros(1,N); +intOut = 0; +xn = sin(2*pi*[0:.1:10]); +for ii = 1:length(xn) + % comb section + combOut = xn(ii) - delayBuffer(end); + delayBuffer(2:end) = delayBuffer(1:end-1); + delayBuffer(1) = xn(ii); + + % integrator + intOut = intOut + combOut; + y2n(ii) = intOut; +end + +err12 = y1n(1:length(xn)) - y2n; +err12dB = 10*log10(err12*err12'/length(err12)); % identical outputs + + +% Implementing Cascaded Integrator Comb filter with the +% integrator section following the comb stage + +N = 10; +delayBuffer = zeros(1,N); +intOut = 0; +xn = sin(2*pi*[0:.1:10]); +for ii = 1:length(xn) + % integrator + intOut = intOut + xn(ii); + % comb section + combOut = intOut - delayBuffer(end); + delayBuffer(2:end) = delayBuffer(1:end-1); + delayBuffer(1) = intOut; + y3n(ii) = combOut; + +end +err13 = y1n(1:length(xn)) - y3n; +err13dB = 10*log10(err13*err13'/length(err13)); % identical outputs + +figure() +hold on +plot(xn) +plot(y1n) \ No newline at end of file diff --git a/Classes/04_DSP/Equalizer/FFE_adaptive_decision.m b/Classes/04_DSP/Equalizer/FFE_adaptive_decision.m index f1f9c9d..f1d9332 100644 --- a/Classes/04_DSP/Equalizer/FFE_adaptive_decision.m +++ b/Classes/04_DSP/Equalizer/FFE_adaptive_decision.m @@ -82,7 +82,7 @@ classdef FFE_adaptive_decision < handle end X.fs = D.fs; %change sampling frequency of outgoing signal from fdac e.g. 2 sps to symbol spaced = fsym lbdesc = [num2str(obj.order),' tap FFE']; - X = X.logbookentry(lbdesc); % append to logbook + X = X.logbookentry(lbdesc,obj); % append to logbook end @@ -106,6 +106,8 @@ classdef FFE_adaptive_decision < handle symbol = 0; % y_buffer = zeros(numel(obj.constellation),500); y_buffer = repmat(obj.constellation,1,obj.buffer_length); + adap_constellation = zeros(length(obj.constellation),length(x)); + cnt=0; for sample = 1 : obj.sps : N symbol = symbol+1; @@ -119,16 +121,24 @@ classdef FFE_adaptive_decision < handle d_hat(symbol,1) = d(symbol); else y_buffer(y_buffer==0) = NaN; - adap_constellation = mean(y_buffer,2,"omitnan"); - [~,symbol_idx] = min(abs(y(symbol) - adap_constellation)); % decision for closest constellation point + + adap_constellation(:,symbol) = mean(y_buffer,2,"omitnan"); + + [~,symbol_idx] = min(abs(y(symbol) - adap_constellation(:,symbol))); % decision for closest constellation point + d_hat(symbol,1) = obj.constellation(symbol_idx); end - y_buffer(symbol_idx,1) = y(symbol); y_buffer(symbol_idx,:) = circshift(y_buffer(symbol_idx,:),1); + + %y_(symbol) = y(symbol) - (adap_constellation(symbol_idx,symbol)-d_hat(symbol,1)); - err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + if training + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + else + err(symbol) = y(symbol) - d_hat(symbol); % Instantaneous error + end if mio ~= 0 obj.e = obj.e - (mio * err(symbol) * U) ; % Weight update rule of LMS @@ -136,13 +146,13 @@ classdef FFE_adaptive_decision < handle normalizationfactor = (U.' * U); obj.e = obj.e - err(symbol) * U / normalizationfactor; % Weight update rule of NLMS end - obj.error(epoch,symbol) = err(symbol) * err(symbol)'; % Instantaneous square error end end - + + %figure(1234);scatter(1:length(d_hat),y,1,'.');hold on;scatter(1:length(d_hat),y_,1,'.');plot(adap_constellation(1,1:length(d_hat)));plot(adap_constellation(2,1:length(d_hat)));plot(adap_constellation(3,1:length(d_hat)));plot(adap_constellation(4,1:length(d_hat))) end end diff --git a/Classes/05_Lab/Awg2Scope.m b/Classes/05_Lab/Awg2Scope.m new file mode 100644 index 0000000..765546c --- /dev/null +++ b/Classes/05_Lab/Awg2Scope.m @@ -0,0 +1,109 @@ +classdef Awg2Scope + %NAME Summary of this class goes here + % Detailed explanation goes here + + properties(Access=public) + Awg + Scope + + mapping; + + end + + methods (Access=public) + function obj = Awg2Scope(Awg,Scope,mapping) + %Simple class to call the Awg and Scope and map the signals + %accordingly in the correct formats with correct l + % ogbook + %entries... + + arguments + Awg + Scope + mapping + end + + obj.Awg = Awg; + obj.Scope = Scope; + + obj.mapping = mapping; % AWG CH [1,2,3,4] -> Scope CH [0,0,0,1] + + end + + function [S1,S2,S3,S4] = process(obj,channels) + + arguments + obj + % leave this as it is! Important for further handling/ parsing + channels.signal1 Informationsignal = Informationsignal([]) + channels.signal2 Informationsignal = Informationsignal([]) + channels.signal3 Informationsignal = Informationsignal([]) + channels.signal4 Informationsignal = Informationsignal([]) + + % add new optional arguments here + end + + [S1,S2,S3,S4]=obj.Awg.upload("signal1",channels.signal1,... + "signal2",channels.signal2,... + "signal3",channels.signal3,... + "signal4",channels.signal4... + ); + + scpe_sig_cell = obj.Scope.read(); + + % Map Scope measurement to output signal + % mapping index is the AWG chanel and mapping number is the + % respective Scope channel + % mapping=[1 2 3 4] means that AWG chann 1 is mapped to Scope ch 1 and so on + % mapping=[0 0 2 1] means that + % AWG chann 3 is mapped to Scope ch 2 + % AWG chann 4 is mapped to Scope ch 1 + + lbdesc = ['Scope Record']; + try + S1 = Electricalsignal(S1,"fs",S1.fs,"logbook",S1.logbook); + S1.signal = scpe_sig_cell{obj.mapping(1)}.signal; + S1.fs = scpe_sig_cell{obj.mapping(1)}.fs; + S1 = S1.logbookentry(lbdesc,obj); + catch + % S1.signal = scpe_sig_cell{1}; + end + + try + S2 = Electricalsignal(S2,"fs",S2.fs,"logbook",S2.logbook); + S2.signal = scpe_sig_cell{obj.mapping(2)}.signal; + S2.fs = scpe_sig_cell{obj.mapping(2)}.fs; + S2 = S2.logbookentry(lbdesc,obj); + catch + % S2.signal = scpe_sig_cell{2}; + S2 = S2.logbookentry(lbdesc,obj); + end + + try + S3 = Electricalsignal(S3,"fs",S3.fs,"logbook",S3.logbook); + S3.signal = scpe_sig_cell{obj.mapping(3)}.signal; + S3.fs = scpe_sig_cell{obj.mapping(3)}.fs; + S3 = S3.logbookentry(lbdesc,obj); + catch + % S3.signal = scpe_sig_cell{3}; + S3 = S3.logbookentry(lbdesc,obj); + end + + try + S4 = Electricalsignal(S4,"fs",S4.fs,"logbook",S4.logbook); + S4.signal = scpe_sig_cell{obj.mapping(4)}.signal; + S4.fs = scpe_sig_cell{obj.mapping(4)}.fs; + S4 = S4.logbookentry(lbdesc,obj); + catch + % S4.signal = scpe_sig_cell{4}; + S4 = S4.logbookentry(lbdesc,obj); + end + + + + end + + + end + +end diff --git a/Classes/05_Lab/AwgKeysight.m b/Classes/05_Lab/AwgKeysight.m index f02a7c5..0270ed4 100644 --- a/Classes/05_Lab/AwgKeysight.m +++ b/Classes/05_Lab/AwgKeysight.m @@ -24,9 +24,9 @@ classdef AwgKeysight properties(Access=protected) model awg_model - fdac double + %fdac double skews double - voltages double + %voltages double scaletodac logical numChannels double @@ -36,6 +36,11 @@ classdef AwgKeysight numProvidedSignals double end + properties(Access=public) + fdac double + voltages double + end + methods (Access=public) function obj = AwgKeysight(options) @@ -78,7 +83,7 @@ classdef AwgKeysight end - function upload(obj,channels,options) + function [signal1,signal2,signal3,signal4] = upload(obj,channels,options) arguments obj @@ -116,7 +121,18 @@ classdef AwgKeysight success = obj.upload_(unpackedSignals); + for s = 1:obj.numChannels + lbdesc = ['Upload to Awg']; + obj = channels.(fn{s}).logbookentry(lbdesc,obj); + end + + signal1 = channels.signal1; + signal2 = channels.signal2; + signal3 = channels.signal3; + signal4 = channels.signal4; + end + end @@ -203,8 +219,11 @@ classdef AwgKeysight v = visadev('TCPIP0::localhost::hislip0::INSTR'); end - disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); - + debug = 0; + if debug + disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); + end + %check if channel config matches writeline(v,'*opt?'); aw=readline(v); diff --git a/Classes/05_Lab/DC_supply.m b/Classes/05_Lab/DC_supply.m index 5d79681..a04cd3e 100644 --- a/Classes/05_Lab/DC_supply.m +++ b/Classes/05_Lab/DC_supply.m @@ -47,7 +47,10 @@ classdef DC_supply %connect to device v = visadev("GPIB1::19::INSTR"); - disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); + debug = 0; + if debug + disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); + end cmd = 'INST:SEL?'; writeline(v, cmd); diff --git a/Classes/05_Lab/OptAtten.m b/Classes/05_Lab/OptAtten.m index 5e0e320..84e98b5 100644 --- a/Classes/05_Lab/OptAtten.m +++ b/Classes/05_Lab/OptAtten.m @@ -55,7 +55,10 @@ classdef OptAtten < handle %connect to device v = visadev('TCPIP::134.245.243.248::INSTR'); - disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); + debug = 0; + if debug + disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); + end %Keysight Technologies, N7764A, MY49A00696, 1.13.1 diff --git a/Classes/05_Lab/ScopeKeysight.m b/Classes/05_Lab/ScopeKeysight.m index 9b33923..5578a32 100644 --- a/Classes/05_Lab/ScopeKeysight.m +++ b/Classes/05_Lab/ScopeKeysight.m @@ -79,8 +79,11 @@ classdef ScopeKeysight end end - disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); - + debug = 0; + if debug + disp(['Connected to Instrument: ',char(v.Vendor),' ',char(v.Model),' SerNo:',char(v.SerialNumber)]); + end + % Check if Scope is ready to go opdone = 0; acqdone = 0; @@ -202,6 +205,11 @@ classdef ScopeKeysight procdone = sscanf(obj.writeReceiveCheck(v,'pder?'), '%f'); pause(0.005); end + + % arrange to Electrcalsinal class as output %% + for ch = 1:numel(obj.channel) + recordedSignals{ch} = Electricalsignal(recordedSignals{ch},"fs",obj.fadc.getValue);%fs is a enum and requires get function + end end diff --git a/Classes/Warehouse_class/classes/DataStorage.m b/Classes/Warehouse_class/classes/DataStorage.m index 5111ade..6fd1e97 100644 --- a/Classes/Warehouse_class/classes/DataStorage.m +++ b/Classes/Warehouse_class/classes/DataStorage.m @@ -48,6 +48,14 @@ classdef DataStorage < handle end + function save(obj,path) + try + save(path,"obj"); + catch + + end + end + function showInfo(obj) disp("Data Structure with fields:"); fprintf('%-12s', 'Name'); fprintf('%1s', '| '); fprintf('%0s ', 'Dimension'); fprintf('%4s', '| '); fprintf('%0s ', 'Physical Values'); fprintf('\n'); diff --git a/Classes/Warehouse_class/classes/DataStorage2.m b/Classes/Warehouse_class/classes/DataStorage2.m new file mode 100644 index 0000000..018f51c --- /dev/null +++ b/Classes/Warehouse_class/classes/DataStorage2.m @@ -0,0 +1,114 @@ +classdef DataStorage2 < handle + % DATASTORAGE: Stores data with physical parameter mappings + + properties + inputParams = struct; + parameter = struct; + fn = []; + dim = []; + sto = struct; + end + + methods + function obj = DataStorage2(inputParams) + % Constructor to initialize the DataStorage object + if nargin > 0 + obj.inputParams = inputParams; + obj.fn = string(fieldnames(inputParams)); + obj = obj.buildParameter(); + obj.dim = obj.getDimension(); + obj.sto = struct; + else + error('Input parameters are required.'); + end + + end + + function showInfo(obj) + % Displays information about the storage and its dimensions + disp("Data Structure with fields:"); + fprintf('%-12s | %-8s | %-12s\n', 'Name', 'Dimension', 'Physical Values'); + disp('-------------------------------------------------------'); + for i = 1:numel(obj.fn) + fprintf('%-12s | %-8d | %-12s\n', ... + char(obj.fn(i)), obj.dim(i), ... + strjoin(string(obj.parameter.(obj.fn(i)).values), ', ')); + end + disp('-------------------------------------------------------'); + end + + function dim = getDimension(obj) + % Get the dimensions based on the length of parameters + dim = zeros(1, numel(obj.fn)); + for p = 1:numel(obj.fn) + dim(p) = obj.parameter.(obj.fn(p)).length; + end + end + + function obj = buildParameter(obj) + % Build the Parameter objects for each input parameter + for p = 1:numel(obj.fn) + name = obj.fn(p); + values = obj.inputParams.(name); + obj.parameter.(name) = Parameter2(name, values); + end + end + + function addStorage(obj, varName) + % Create an empty storage for a specific variable name + obj.sto.(string(varName)) = cell(obj.dim); + end + + function addValueToStorage(obj, valueToStore, storageVarName, varargin) + % Add a value to the storage at the specified indices + if nargin - 3 == numel(obj.fn) + lin_idx = obj.getIndicesByPhys(varargin); + obj.sto.(storageVarName){lin_idx} = valueToStore; + else + error('Please provide all indices for the storage.'); + end + end + + function value = getStoValue(obj, storageVarName, varargin) + % Retrieve a value from storage based on physical parameters + if nargin - 2 == numel(obj.fn) + lin_idx = obj.getIndicesByPhys(varargin); + value = cell(1, numel(lin_idx)); + for i = 1:numel(lin_idx) + value{i} = obj.sto.(storageVarName){lin_idx(i)}; + end + value = value(~cellfun('isempty', value)); % Remove empty entries + else + error('Please provide all physical parameters.'); + end + end + + function lin_idx = getIndicesByPhys(obj, varargin) + % Unpack nested cell array if needed + if numel(varargin) == 1 && iscell(varargin{1}) + varargin = varargin{1}; % Unpack if single cell array is passed + end + + indices = cell(1, numel(obj.fn)); + + % Loop through each parameter (e.g., L, D) + for p = 1:numel(obj.fn) + % Unwrap if it's a cell + if iscell(varargin{p}) + physVal = varargin{p}{1}; % Extract scalar from cell + else + physVal = varargin{p}; % It's already a scalar + end + + paramName = obj.fn(p); % Get the parameter name (e.g., 'L' or 'D') + + % Call getIndexByPhys on the corresponding Parameter2 object + indices{p} = obj.parameter.(paramName).getIndexByPhys(physVal); + end + + % Convert subscript indices to a linear index + lin_idx = sub2ind(obj.dim, indices{:}); + end + + end +end diff --git a/Classes/Warehouse_class/classes/Parameter2.m b/Classes/Warehouse_class/classes/Parameter2.m new file mode 100644 index 0000000..84852ff --- /dev/null +++ b/Classes/Warehouse_class/classes/Parameter2.m @@ -0,0 +1,51 @@ +classdef Parameter2 < handle + % PARAMETER2: Represents a physical parameter with mappings between values and indices + + properties + name + values + length + physToIndexMap % Rename this from 'getPhysForIndex' + indexToPhysMap % Rename this from 'getIndexForPhys' + end + + methods + function obj = Parameter2(name, values) + % Constructor to initialize the Parameter2 object + obj.name = name; + obj.values = values; + obj.length = numel(values); + + % Initialize the mappings + obj.physToIndexMap = containers.Map('KeyType', 'double', 'ValueType', 'any'); + obj.indexToPhysMap = containers.Map('KeyType', 'double', 'ValueType', 'any'); + obj = obj.buildMappings(); + end + + function obj = buildMappings(obj) + % Build mappings between physical values and indices + for idx = 1:obj.length + obj.indexToPhysMap(idx) = obj.values(idx); + obj.physToIndexMap(obj.values(idx)) = idx; + end + end + + function physVal = getPhysForIndex(obj, idx) + % Return the physical value corresponding to the index + if isKey(obj.indexToPhysMap, idx) + physVal = obj.indexToPhysMap(idx); + else + error('Index out of range for parameter %s', obj.name); + end + end + + function idx = getIndexByPhys(obj, physVal) + % Return the index corresponding to the physical value + if isKey(obj.physToIndexMap, physVal) + idx = obj.physToIndexMap(physVal); + else + error('Physical value %g not found in parameter %s', physVal, obj.name); + end + end + end +end diff --git a/Classes/Warehouse_class/classes/minimalExample_gen2.m b/Classes/Warehouse_class/classes/minimalExample_gen2.m new file mode 100644 index 0000000..c348197 --- /dev/null +++ b/Classes/Warehouse_class/classes/minimalExample_gen2.m @@ -0,0 +1,45 @@ +% Define input parameters for the DataStorage2 +inputParams.L = [1, 2, 10, 80]; % Length in kilometers +inputParams.D = [16, 17, 18]; % Diameter in millimeters + +% Create a DataStorage2 instance with the input parameters +dataStorage = DataStorage2(inputParams); % Using DataStorage2 class + +% Display the current information about the data storage structure +dataStorage.showInfo(); + +% Add a storage variable named 'testStorage' +dataStorage.addStorage('testStorage'); + +% Add a value (e.g., 100) to the storage at specific physical parameter values +% For example, we store the value 100 at L = 10 km and D = 17 mm +dataStorage.addValueToStorage(100, 'testStorage', 10, 16); +dataStorage.addValueToStorage(100, 'testStorage', 10, 17); +dataStorage.addValueToStorage(100, 'testStorage', 10, 18); + +% Retrieve the value from the storage at the same physical parameter values +storedValue = dataStorage.getStoValue('testStorage', 10, 16:18); +disp('Retrieved value from storage:'); +disp(storedValue); + +% Retrieve another value at a non-existent location (L = 2 km, D = 8 mm) +% This will show how the function handles empty storage entries +nonExistentValue = dataStorage.getStoValue('testStorage', 2, 8); +disp('Retrieved value from empty location:'); +disp(nonExistentValue); + +% Use the internal mappings to check how physical values map to indices +% Get the linear index for physical values L = 10 km and D = 17 mm +lin_idx = dataStorage.getIndicesByPhys(10, 17); +disp('Linear index for L=10 km and D=17 mm:'); +disp(lin_idx); + +% Check the reverse mapping: physical value for index 2 of parameter L +physValForIndex = dataStorage.parameter.L.getPhysForIndex(2); +disp('Physical value for index 2 of parameter L:'); +disp(physValForIndex); + +% Check the mapping: index for physical value D = 21 mm +indexForPhys = dataStorage.parameter.D.getIndexByPhys(21); +disp('Index for physical value D=21 mm:'); +disp(indexForPhys); diff --git a/Datatypes/scope_fadc.m b/Datatypes/scope_fadc.m index cc9c3c6..e4ab2c6 100644 --- a/Datatypes/scope_fadc.m +++ b/Datatypes/scope_fadc.m @@ -15,7 +15,7 @@ classdef scope_fadc < int32 methods function val = getValue(obj) - val = double(obj); % Convert int32 to double to get the numeric value + val = double(obj).*1e9; % Convert int32 to double and from GHz to Hz to get the numeric value end end diff --git a/projects/400G_FTN_setups/imdd_dsp_approaches.m b/projects/400G_FTN_setups/imdd_dsp_approaches.m index 8192c39..8ebf223 100644 --- a/projects/400G_FTN_setups/imdd_dsp_approaches.m +++ b/projects/400G_FTN_setups/imdd_dsp_approaches.m @@ -27,6 +27,7 @@ name = ['wh_',strrep(num2str(now),'.','')]; wh = DataStorage(params); +wh.addStorage("Rx_Bits"); wh.addStorage("ber_ffe"); %% Init Params @@ -200,7 +201,7 @@ for M = wh.parameter.M.values end end - Rx_bits = PAMmapper(M,0).demap(EQ_sig); + Rx_bits(i) = PAMmapper(M,0).demap(EQ_sig); [~,errors_bm,ber_ffe(i),errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); disp(['BER: ',sprintf('%.1E',ber_ffe(i)),' - - ROP: ',num2str(patten(i)),'dBm - - PAM-',num2str(M),' - - ',num2str(fsym*1e-9),' GBd']); @@ -210,12 +211,14 @@ for M = wh.parameter.M.values rop=wh.parameter.rop.values(i); wh.addValueToStorage(ber_ffe(i),'ber_ffe',M,datarate,rop); + wh.addValueToStorage(Rx_bits(i),'Rx_Bits',M,datarate,rop); end toc - % wh.save('C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\MPI_August\auswertung\') + filename = 'bla3'; + wh.save(['C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\MPI_August\auswertung\',filename]); end end diff --git a/projects/Lab_2024/lab_db_precode.m b/projects/Lab_2024/lab_db_precode.m new file mode 100644 index 0000000..1cdcd1a --- /dev/null +++ b/projects/Lab_2024/lab_db_precode.m @@ -0,0 +1,253 @@ + +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\mpi_ofc_2024\'; +experiment_name = '10km_db_transmit_no_mpi_'; + +only_dsp = 0; + +ffe_only = 0; +postfilter_approach = 1; +db_channel_approach = 0; +db_coding_approach = 0; + +db_precode = db_coding_approach || db_channel_approach; + +%%% SIR Sweep for MPI Experiment %%% +params = struct; +params.i_atten = [40]; + +wh = DataStorage(params); + +wh.addStorage("ber"); +wh.addStorage("sir"); +wh.addStorage("pd_in"); +wh.addStorage("signals"); + +precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active +precomp_amp_max = 3; + +M = 4; +pn_key = 2; +usemrds = 0; +fdac = 92e9; +fsym = 92e9; +fadc = 160e9; +rrcalpha = 0.05; +v_bias = 2.25; +i_atten = params.i_atten(1); +pd_in_desired = 7; + +if ~only_dsp + + disp(['Start Measurement of ',num2str(prod(wh.dim)),' loops...']) + + %%%%% SET Volatges %%%%%% + dcs = DC_supply("active",[1,1],"voltage",[v_bias, 9]); + dcs.set("voltage",[v_bias, 9]); + + voa = OptAtten("active",[1,2,1,1],"value",[0,pd_in_desired,0,i_atten],"wavelength",[1310,1310,1310,1310]); + % voa.set('active',[1,2,1,1],'value',[0,pd_in_desired,0,i_atten]); + % voa.readvals(); + + 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,0.62]); + A2S = Awg2Scope(AWG,SCP,[0,0,0,1]); + + + if 1 + precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\standard_system\"; + precomp_fn = "lab_mpi_setup_2"; + else + precomp_path = "C:\Users\sioe\Documents\MATLAB\model-collection\sioe_models\Labor_2024\Lab_PAM4\"; + precomp_fn = "precomp_bla__loop1_1"; + end + + + %%%%% Symbol Generation %%%%%% + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... + "fs_out",fdac,"applyclipping",0,"clipfactor",1.7,... + "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... + "db_precode",db_precode,... + "mrds_code",usemrds,"mrds_blocklength",512,"db_encode",db_coding_approach).process(); + + if precomp_mode == 1 % measure channel + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.buildOFDM(); + elseif precomp_mode == 2 % apply precomp + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',precomp_amp_max,'loadPath',precomp_path,'fileName',precomp_fn); + end + + Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); + + Digi_sig.spectrum("displayname","Normal Tx","fignum",10); + + % Digi_sig = Filter('filtdegree',1,"f_cutoff",45e9,"fs",Digi_sig.fs,"filterType",filtertypes.butterworth,"active",true).process(Digi_sig); + + % Digi_sig.spectrum("displayname","Lowpass Tx","fignum",999); + % Digi_sig.eye(fsym,M); + + + save([folderpath,[experiment_name,'bits']],"Bits"); + save([folderpath,[experiment_name,'symbols']],"Symbols"); + +end + +for i_atten = wh.parameter.i_atten.values + + if ~only_dsp + + %%%%% SET ATTENUATOR %%%%%% + voa.set('active',[1,2,1,1],'value',[0,7,0,i_atten]); + + %%%%% AWG --> Scope %%%%%% + [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); + + Scpe_sig.spectrum("displayname","Scope PSD","fignum",20); + + + + % Scpe_sig.spectrum("displayname",'Rx Signal','fignum',10); + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*fsym); + + if precomp_mode == 1 + freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + freqresp.plot(); + end + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + save([folderpath,experiment_name,'rx_signal_iatten_',num2str(i_atten),''],"S"); + + average_signals = 0; + if average_signals + scope_mean = zeros(size(S{1}.signal)); + for n=1:numel(S) + scope_mean = scope_mean + S{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig.signal = scope_mean; + end + + Scpe_sig.plot("displayname","Scope PSD","fignum",30); + Scpe_sig.eye(fsym,M,"fignum",40,"displayname",' after Scope'); + + + sir = voa.power_state(3)-voa.power_state(4); + pd_in = voa.power_state(2); + + end + + %%%%% 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); + + if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + 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); + + disp(['FFE: ',sprintf('%.1E',ber),'| SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + Noi = EQ_sig-Symbols; + + 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); + + nc = 2; + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig = EQ_sig.filter(burg_coeff,1); + + if 0 + 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 + + 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),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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); + + EQ_sig.plot("fignum",50,"displayname",'After EQ'); + + disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + end + + + + wh.addValueToStorage(ber,'ber',i_atten); + wh.addValueToStorage(sir,'sir',i_atten); + wh.addValueToStorage(pd_in,'pd_in',i_atten); + wh.addValueToStorage(Rx_bits,'signals',i_atten); + +end + +wh.save([folderpath,experiment_name,'_wh']); + +cols = linspecer(8); + +sir_vals = wh.parameter.i_atten.values; +bers = wh.getStoValue('ber',sir_vals); + +figure(90); +a = gca; +hold on; % Retain the plot so new points can be added without complete redraw + +plot(sir_vals,bers,"LineWidth",0.5,"LineStyle","-","Marker",".","MarkerSize",15,"DisplayName",[experiment_name]); +yline(3.8e-3,'DisplayName','HD-FEC','LineStyle','--','HandleVisibility','off'); +xlabel('Received Optical Power (dBm)'); +ylabel('Bit Error Rate (BER)'); +title('Bit Error Rate vs. ROP'); +set(gca,'yscale','log'); +set(gca,'Box','on'); +grid on; +grid minor +legend('Interpreter','none') + +autoArrangeFigures(2,3,2) + diff --git a/projects/Lab_2024/lab_minimal.m b/projects/Lab_2024/lab_minimal.m index dbb690f..fe61cfa 100644 --- a/projects/Lab_2024/lab_minimal.m +++ b/projects/Lab_2024/lab_minimal.m @@ -1,14 +1,18 @@ + + +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\mpi_ofc_2024\'; + %%% SIR Sweep for MPI Experiment %%% params = struct; -params.i_atten = 10; +params.i_atten = [40]; wh = DataStorage(params); wh.addStorage("ber"); +wh.addStorage("ber_pf"); wh.addStorage("sir"); wh.addStorage("pd_in"); - -playandrecord = 0; +wh.addStorage("signals"); precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active M = 4; @@ -18,8 +22,8 @@ fdac = 92e9; fsym = 92e9; fadc = 160e9; rrcalpha = 0.05; -v_bias = 2.27; -i_atten = 40; +v_bias = 2.25; +i_atten = params.i_atten(1); pd_in_desired = 7; disp(['Start Measurement of ',num2str(prod(wh.dim)),' loops...']) @@ -28,10 +32,13 @@ disp(['Start Measurement of ',num2str(prod(wh.dim)),' loops...']) dcs = DC_supply("active",[1,1],"voltage",[v_bias, 9]); dcs.set("voltage",[v_bias, 9]); -voa = OptAtten("active",[0,2,1,1],"value",[0,pd_in_desired,0,i_atten],"wavelength",[1310,1310,1310,1310]); -voa.set('active',[0,2,1,1],'value',[0,pd_in_desired,0,i_atten]); -voa.readvals(); -%%%%% SET Volatges %%%%%% +voa = OptAtten("active",[1,2,1,1],"value",[0,pd_in_desired,0,i_atten],"wavelength",[1310,1310,1310,1310]); +% voa.set('active',[1,2,1,1],'value',[0,pd_in_desired,0,i_atten]); +% voa.readvals(); + +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,0.62]); +A2S = Awg2Scope(AWG,SCP,[0,0,0,1]); if 1 @@ -45,6 +52,7 @@ end %%%%% Symbol Generation %%%%%% Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... "fs_out",fdac,"applyclipping",0,"clipfactor",1.7,... "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... @@ -59,24 +67,23 @@ elseif precomp_mode == 2 % apply precomp Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',3,'loadPath',precomp_path,'fileName',precomp_fn); end -%%%%% AWG %%%%%% -AWG = AwgKeysight("model","M8196A","fdac",fdac,"scaletodac",[1,1,1,1],"skews",[0,0,0,0],"voltages",[0,0,0,0.62]); -AWG.upload("signal4",Digi_sig); +Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); + +save([folderpath,'bits'],"Bits"); +save([folderpath,'symbols'],"Symbols"); for i_atten = wh.parameter.i_atten.values %%%%% SET ATTENUATOR %%%%%% - voa.set('active',[0,2,1,1],'value',[0,7,0,i_atten]); + voa.set('active',[1,2,1,1],'value',[0,7,0,i_atten]); - %%%%% Scope %%%%%% - SCP = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc",'GSa_160',"channel",[1],"recordLen",2000000,"removeDC",1); - Scpe_sig = SCP.read(); - Scpe_sig = Electricalsignal(Scpe_sig{1},"fs",fadc); + %%%%% AWG --> Scope %%%%%% + [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); % Scpe_sig.spectrum("displayname",'Rx Signal','fignum',10); %%%%%% Sample to 2x fsym %%%%%% - Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*92e9); + Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*fsym); if precomp_mode == 1 freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); @@ -85,44 +92,57 @@ for i_atten = wh.parameter.i_atten.values %%%%%% Sync Rx signal with reference %%%%%% [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + save([folderpath,'rx_iatten_',num2str(i_atten),''],"S"); Scpe_sig.eye(fsym,M); + Scpe_sig.spectrum("displayname","Scope PSD","fignum",1996); %%%%% EQUALIZE %%%%%% + + % Scpe_sig.signal = Scpe_sig.signal - movmean(Scpe_sig.signal,[100,0]); - Eq_ffe_only = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0); - ber_ffe_only = runEQ(Eq_ffe_only,Scpe_sig,Symbols,Bits,M); + Eq_ffe_only = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",50,"sps",2,"decide",0); + [ber_ffe_only,EQ_sig,Noi] = runEQ(Eq_ffe_only,Scpe_sig,Symbols,Bits,M); - Eq_move_it = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.05,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); - ber_move_it = runEQ(Eq_move_it,Scpe_sig,Symbols,Bits,M); + nc = 2; + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig = EQ_sig.filter(burg_coeff,1); + + EQ_sig.plot("displayname",'After EQ','fignum',1112); + + if 0 + 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 + + % EQ_sig.spectrum("displayname","Signal Spectrum after Postfilter","fignum",1234); - Eq_nonlin = 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); - %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",0); - ber_nonlin = runEQ(Eq_nonlin,Scpe_sig,Symbols,Bits,M); - - Eq_adaptive_decision = FFE_adaptive_decision("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"buffer_length",80); - ber_adaptive_decision = runEQ(Eq_adaptive_decision,Scpe_sig,Symbols,Bits,M); - - FFE_FFDCAVG("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"mu_buff",128); - ber_ff_dcavg = runEQ(FFE_FFDCAVG,Scpe_sig,Symbols,Bits,M); - - Eq_feedback_removal = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",25,"sps",2,"decide",0,"mu_dc",0.05,"dc_buffer_len",1); - ber_feedback_removal = runEQ(Eq_feedback_removal,Scpe_sig,Symbols,Bits,M); - + 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_pf,errors] = calc_ber(Rx_bits.signal,Bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); sir = voa.power_state(3)-voa.power_state(4); pd_in = voa.power_state(2); - disp(['BER: ',sprintf('%.1E',ber_eq),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + disp(['FFE: ',sprintf('%.1E',ber_ffe_only),' -> PF -> MLSE: ',sprintf('%.1E',ber_pf),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); - wh.addValueToStorage(ber_eq,'ber',i_atten); + wh.addValueToStorage(ber_ffe_only,'ber',i_atten); + wh.addValueToStorage(ber_pf,'ber_pf',i_atten); wh.addValueToStorage(sir,'sir',i_atten); wh.addValueToStorage(pd_in,'pd_in',i_atten); + wh.addValueToStorage(Rx_bits,'signals',i_atten); end - +name = ['mpidata_92GBd']; +wh.save([folderpath,name]); cols = linspecer(8); diff --git a/projects/Lab_2024/runEQ.m b/projects/Lab_2024/runEQ.m index 1fe66a8..da1661f 100644 --- a/projects/Lab_2024/runEQ.m +++ b/projects/Lab_2024/runEQ.m @@ -1,13 +1,25 @@ -function BER = runEQ(Eq_class,Scpe_sig,Symbols,Bits,M) - -[EQ_sig] = Eq_class.process(Scpe_sig,Symbols); - -error = EQ_sig-Symbols; -error.spectrum("displayname",['Error PSD after: ',inputname(1),' '],'fignum',564); - -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); +function [BER,EQ_sig,Noi] = runEQ(Eq_class,Scpe_sig,Symbols,Bits,M) + [EQ_sig] = Eq_class.process(Scpe_sig,Symbols); + + Noi = EQ_sig-Symbols; + + + 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); + + if 0 + figure(); + 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)); + hold on + histogram(received(lvl,:),1000,"EdgeAlpha",0); + end + + Noi.spectrum("displayname",['Error PSD after: ',inputname(1),' '],'fignum',564); + end end \ No newline at end of file From bf94e3dc2f6bf9f6276a2a35b5d96091652eaffc Mon Sep 17 00:00:00 2001 From: Silas Labor Zizou Date: Tue, 15 Oct 2024 08:43:27 +0200 Subject: [PATCH 2/2] Work on lab PC - many new automations - scripts to record and save MPI and bias optimizations... --- Classes/04_DSP/Coding/Duobinary.m | 2 +- Classes/05_Lab/Exfo_laser.m | 74 ++++ Classes/Warehouse_class/classes/DataStorage.m | 16 +- Functions/showCurrentMeasurement.m | 136 +++++++ Functions/updateWaitbar.m | 45 +++ projects/Lab_2024/bias_sweep_evaluation.m | 145 +++++++ projects/Lab_2024/bias_sweep_evaluation_2.m | 124 ++++++ projects/Lab_2024/lab_baudrate_sweep.m | 274 +++++++++++++ projects/Lab_2024/lab_bias_sweep.m | 269 +++++++++++++ projects/Lab_2024/lab_bias_sweep_gigantisch.m | 320 +++++++++++++++ projects/Lab_2024/lab_db_precode.m | 197 +++++----- projects/Lab_2024/lab_modulator_tf_sweep.m | 68 ++++ projects/Lab_2024/lab_precompensation_sweep.m | 265 +++++++++++++ projects/Lab_2024/lab_rop_sweep.m | 261 +++++++++++++ projects/Lab_2024/lab_sir_sweep.m | 368 ++++++++++++++++++ projects/Lab_2024/sweep_laser_vs_power.m | 88 +++++ 16 files changed, 2548 insertions(+), 104 deletions(-) create mode 100644 Classes/05_Lab/Exfo_laser.m create mode 100644 Functions/showCurrentMeasurement.m create mode 100644 Functions/updateWaitbar.m create mode 100644 projects/Lab_2024/bias_sweep_evaluation.m create mode 100644 projects/Lab_2024/bias_sweep_evaluation_2.m create mode 100644 projects/Lab_2024/lab_baudrate_sweep.m create mode 100644 projects/Lab_2024/lab_bias_sweep.m create mode 100644 projects/Lab_2024/lab_bias_sweep_gigantisch.m create mode 100644 projects/Lab_2024/lab_modulator_tf_sweep.m create mode 100644 projects/Lab_2024/lab_precompensation_sweep.m create mode 100644 projects/Lab_2024/lab_rop_sweep.m create mode 100644 projects/Lab_2024/lab_sir_sweep.m create mode 100644 projects/Lab_2024/sweep_laser_vs_power.m diff --git a/Classes/04_DSP/Coding/Duobinary.m b/Classes/04_DSP/Coding/Duobinary.m index e9888b9..259e04d 100644 --- a/Classes/04_DSP/Coding/Duobinary.m +++ b/Classes/04_DSP/Coding/Duobinary.m @@ -108,7 +108,7 @@ classdef Duobinary data = data - b; data = data ./ 2; - assert(isequal((0:M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar + % assert(isequal((0:M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar % duobinary coding (1+D) % coeff = [1,1]; diff --git a/Classes/05_Lab/Exfo_laser.m b/Classes/05_Lab/Exfo_laser.m new file mode 100644 index 0000000..1fa6acf --- /dev/null +++ b/Classes/05_Lab/Exfo_laser.m @@ -0,0 +1,74 @@ +classdef Exfo_laser + + properties(Access=public) + wavelength + power + end + + methods (Access=public) + + function obj = Exfo_laser(options) + + + arguments + options.wavelength = 1310; %dbm + options.power = -10; %dbm + end + + % + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end + + end + + function success = set(obj,options) + + + arguments + obj + options.wavelength = obj.wavelength; %dbm + options.power = obj.power; %dbm + end + + % Connect to the laser + o = serialport("COM8", 9600); + configureTerminator(o, "CR"); % Set the terminator to carriage return (CR) + writeline(o, "*IDN?"); + pause(1); + if o.NumBytesAvailable ~= 0 + disp(['Laser Mainframe: ', readline(o)]); + else + error('No connection to the mainframe'); + clear o; + end + + + + end + + % Function to set the wavelength of the laser + function setLaserWavelength(~,serialObj, channel, wavelength) + command = ['CH', num2str(channel), ':L=', num2str(wavelength)]; + writeline(serialObj, command); + pause(0.5); % Allow time for the wavelength to change + writeline(serialObj, ['CH', num2str(channel), ':L?']); % Query current wavelength + current_wavelen = readline(serialObj); + disp(['Current Wavelength: ', current_wavelen]); + end + + % Function to set the laser power + function setLaserPower(~,serialObj, channel, power_dBm) + command = ['CH', num2str(channel), ':P=', num2str(power_dBm)]; + writeline(serialObj, command); + pause(0.2); % Allow time for power to adjust + writeline(serialObj, ['CH', num2str(channel), ':P?']); % Query current power + current_power = readline(serialObj); + disp(['Current Power: ', current_power, ' dBm']); + end + + end +end diff --git a/Classes/Warehouse_class/classes/DataStorage.m b/Classes/Warehouse_class/classes/DataStorage.m index 6fd1e97..dd9eb4b 100644 --- a/Classes/Warehouse_class/classes/DataStorage.m +++ b/Classes/Warehouse_class/classes/DataStorage.m @@ -128,7 +128,21 @@ classdef DataStorage < handle try tmp = obj.sto.(storageVarName){lin_idx(i)}; if ~isempty(tmp) - value(i,:) = tmp ; + if isa(tmp,'Signal') + if i == 1 + value = {}; + end + value{i} = tmp ; + elseif isa(tmp,'cell') + if isa(tmp{1},'Signal') + if i == 1 + value = {}; + end + value{i} = tmp{1} ; + end + else + value(i,:) = tmp ; + end else errcnt = errcnt+1; diff --git a/Functions/showCurrentMeasurement.m b/Functions/showCurrentMeasurement.m new file mode 100644 index 0000000..b9710b8 --- /dev/null +++ b/Functions/showCurrentMeasurement.m @@ -0,0 +1,136 @@ +function showCurrentMeasurement(varargin) + % showCurrentMeasurement displays measurement data in a figure with variable names + % as column headers and values listed below. Calling the function multiple times + % with the same variable names but different values adds more data points. + % + % Usage: + % showCurrentMeasurement('VariableName1', VariableValue1, 'VariableName2', VariableValue2, ...) + % + % Example: + % % First measurement + % voltage = 5.12; + % ber = 3.86e-5; + % power = voltage * 0.85; + % showCurrentMeasurement('Voltage', voltage, 'ber', ber, 'Power', power); + % + % % Second measurement + % voltage = 5.15; + % ber = 2.54e-5; + % power = voltage * 0.90; + % showCurrentMeasurement('Voltage', voltage, 'ber', ber, 'Power', power); + + % Validate that inputs are in name-value pairs + if mod(nargin, 2) ~= 0 + error('Inputs must be provided as name-value pairs.'); + end + + % Extract variable names and values + numPairs = nargin / 2; + names = varargin(1:2:end); + values = varargin(2:2:end); + + % Ensure variable names are strings + for i = 1:length(names) + if ~ischar(names{i}) && ~isstring(names{i}) + error('Variable names must be strings.'); + end + names{i} = char(names{i}); + end + + % Convert values to strings for display + for i = 1:length(values) + varNameLower = lower(names{i}); % Convert variable name to lowercase for case-insensitive comparison + if isnumeric(values{i}) + if strcmp(varNameLower, 'ber') + % Format 'ber' values in exponential notation with two decimal places + values{i} = sprintf('%.2e', values{i}); + else + values{i} = num2str(values{i}); + end + else + values{i} = char(values{i}); + end + end + + % Define a unique tag for the figure to locate it later + figTag = 'CurrentMeasurementsFigure'; + + % Try to find an existing figure with the specified tag + hFig = findobj('Type', 'figure', 'Tag', figTag); + + if isempty(hFig) + % Create a new figure and table + hFig = figure('Name', 'Current Measurements', 'NumberTitle', 'off', ... + 'MenuBar', 'none', 'ToolBar', 'none', 'Resize', 'on', ... + 'Tag', figTag); + + % Initialize data and column names + data = values; + columnNames = names; + + % Create the uitable + hTable = uitable('Parent', hFig, 'Data', data, ... + 'ColumnName', columnNames, ... + 'FontSize', 14, ... + 'RowName', [], ... + 'Units', 'normalized', ... + 'Position', [0, 0, 1, 1]); + % Adjust column widths + setColumnWidths(hTable); + + % Store the table handle for future use + setappdata(hFig, 'DataTable', hTable); + else + % Retrieve the existing table handle + hTable = getappdata(hFig, 'DataTable'); + + % Get current data and column names + currentData = get(hTable, 'Data'); + columnNames = get(hTable, 'ColumnName'); + + % Ensure that variable names are consistent + if ~isequal(columnNames, names') + error('Variable names must be consistent with previous calls.'); + end + + % Append new data to the existing data + updatedData = [currentData; values]; + + % Update the table data + set(hTable, 'Data', updatedData); + + % Adjust column widths + setColumnWidths(hTable); + end + + % Adjust the figure size to fit the table content without changing its position + drawnow; + tableExtent = get(hTable, 'Extent'); + % Get the current figure position + figPosition = get(hFig, 'Position'); + % Update the figure size while preserving the position + figPosition(3) = max(figPosition(3), tableExtent(3) + 20); % Width + figPosition(4) = max(figPosition(4), tableExtent(4) + 20); % Height + set(hFig, 'Position', figPosition); +end + +function setColumnWidths(hTable) + % Helper function to adjust column widths based on content + data = get(hTable, 'Data'); + columnNames = get(hTable, 'ColumnName'); + numColumns = length(columnNames); + columnWidths = cell(1, numColumns); + + % Calculate the maximum width needed for each column + for col = 1:numColumns + maxContentLength = max(cellfun(@length, data(:, col))); + headerLength = length(columnNames{col}); + maxLength = max(maxContentLength, headerLength); + + % Estimate pixel width (approximate, adjust as needed) + pixelWidth = maxLength * 14; % 8 pixels per character as an estimate + columnWidths{col} = pixelWidth; + end + + set(hTable, 'ColumnWidth', columnWidths); +end diff --git a/Functions/updateWaitbar.m b/Functions/updateWaitbar.m new file mode 100644 index 0000000..94accd9 --- /dev/null +++ b/Functions/updateWaitbar.m @@ -0,0 +1,45 @@ +function updateWaitbar(currentIteration, totalIterations) + +if currentIteration == 1 + + % Check if the waitbar already exists using its unique Tag + hWaitbar = findobj('Tag', 'MyUniqueWaitbar'); + + if isempty(hWaitbar) || ~ishandle(hWaitbar) + % Create a waitbar with a unique Tag if it doesn't exist + hWaitbar = waitbar(0, 'Starting process...', 'Name', 'Processing Progress', 'Tag', 'MyUniqueWaitbar'); + else + % Waitbar exists, reset the progress bar + waitbar(0, hWaitbar, 'Resuming process...'); + end + +elseif currentIteration == totalIterations+1 + % Check if the waitbar already exists using its unique Tag + hWaitbar = findobj('Tag', 'MyUniqueWaitbar'); + + % Close the waitbar after the loop is completed + if ishandle(hWaitbar) + close(hWaitbar); + end + +else + + % Check if the waitbar already exists using its unique Tag + hWaitbar = findobj('Tag', 'MyUniqueWaitbar'); + + % Calculate the progress fraction + progressFraction = currentIteration / totalIterations; + + % Update the waitbar's progress and message + if ishandle(hWaitbar) + waitbar(progressFraction, hWaitbar, ... + sprintf('Progress: %d/%d', currentIteration, totalIterations)); + else + % % If the waitbar was closed, recreate it + % hWaitbar = waitbar(progressFraction, 'Resuming process...', 'Name', 'Processing Progress', 'Tag', 'MyUniqueWaitbar'); + end + +end + + + diff --git a/projects/Lab_2024/bias_sweep_evaluation.m b/projects/Lab_2024/bias_sweep_evaluation.m new file mode 100644 index 0000000..08ffcac --- /dev/null +++ b/projects/Lab_2024/bias_sweep_evaluation.m @@ -0,0 +1,145 @@ + + +filename = "C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\bias_sweep\PAM6_10km_ffe__wh.mat"; +a = load(filename); +wh2 = a.obj; + +m = wh2.getStoValue('m',wh2.parameter.vbias.values(1),wh2.parameter.awg_vpp.values(1)); + +v_bias_vals = wh2.parameter.vbias.values; +awg_vpp_vals = wh2.parameter.awg_vpp.values; + +bers = []; +rop_measured = []; +cnt = 0; +for awg_vpp_cur = awg_vpp_vals + cnt = cnt+1; + bers(cnt,:) = wh2.getStoValue('ber',v_bias_vals,awg_vpp_cur); + rop_measured(cnt,:) = wh2.getStoValue('rop',v_bias_vals,awg_vpp_cur); +end + +[bestber,bestindex] = min(bers,[],'all'); +[awg_pos,v_bias_pos]=ind2sub(size(bers),bestindex); +bestawgvpp=awg_vpp_vals(awg_pos); +bestvbias=v_bias_vals(v_bias_pos); + +disp(['Best Vpp: ',num2str(bestvbias),' V; Best Vpp AWG: ',num2str(bestawgvpp),' V' ]) + +figure(); +sgtitle(['PAM ', num2str(m)]) +subplot1 = subplot(1,2,1); + +% Compute the logarithm of BER data +% Adding a small epsilon to avoid log(0) +epsilon = 1e-12; +log_bers = log10(bers + epsilon); + +% Set limits for z-data scaling in log scale +zmin = log10(1e-4 + epsilon); +zmax = log10(0.5 + epsilon); + +% Plot the filled contour plot with log-scaled z-data +contourf_handle = contourf(v_bias_vals, awg_vpp_vals, log_bers, 'Parent', subplot1, "ShowText",true,"LabelFormat", @mylabelfun); + +% Set x and y labels with subscripts for clarity +xlabel('V_{bias}'); +ylabel('V_{pp} AWG'); +title('BER (Mind used Equalizer!)'); + +% Adjust the grid to display white lines +grid on; +set(subplot1, 'GridColor', [1 1 1]); % Set grid color to white + +% Set limits for z-data scaling +clim([zmin zmax]); + +% Adjust the colormap +colormap(flipud(cbrewer2('RdBu',64))); + +% Add a colorbar and adjust its ticks to represent actual BER values +c = colorbar; +% Set colorbar ticks at log-spaced intervals +tick_values = [1e-4 1e-3 1e-2 1e-1 0.5]; +tick_positions = log10(tick_values + epsilon); +set(c, 'Ticks', tick_positions, 'TickLabels', arrayfun(@num2str, tick_values, 'UniformOutput', false)); + +% Store variables in the figure's application data for use in the data tip function +setappdata(gcf, 'v_bias_vals', v_bias_vals); +setappdata(gcf, 'awg_vpp_vals', awg_vpp_vals); +setappdata(gcf, 'bers', bers); +setappdata(gcf, 'power', rop_measured); % Store the Power data + +% Set up the data cursor mode to display custom data tips +dcm_obj = datacursormode(gcf); +set(dcm_obj, 'UpdateFcn', @customDataTip); +hold on +scatter(bestvbias,bestawgvpp,100,"red",'Marker','x','LineWidth',2); + + + +subplot2 = subplot(1,2,2); + +% Plot the filled contour plot +contourf_handle = contourf(v_bias_vals, awg_vpp_vals, rop_measured, 'Parent', subplot2); + +% Set x and y labels +xlabel('V_{bias}'); +ylabel('V_{pp} AWG'); +title('Power at Rx'); + +% Adjust the grid to display white lines +grid on; +set(subplot2, 'GridColor', [1 1 1]); % Set grid color to white + +% Set limits for z-data scaling +clim([-5 -3]); + +% Adjust the colormap +colormap(flipud(cbrewer2('RdBu',64))); + +% Add a colorbar +colorbar; + +% Store variables in the figure's application data for use in the data tip function +setappdata(gcf, 'v_bias_vals', v_bias_vals); +setappdata(gcf, 'awg_vpp_vals', awg_vpp_vals); +setappdata(gcf, 'bers', bers); +setappdata(gcf, 'power', rop_measured); % Store the Power data + +% Set up the data cursor mode to display custom data tips +dcm_obj = datacursormode(gcf); +set(dcm_obj, 'UpdateFcn', @customDataTip); + +function labels = mylabelfun(vals) + lab = 10.^vals; + labels = arrayfun(@(x) num2str(x, '%.1e'), lab, 'UniformOutput', false); +end + + +% Define the custom data tip function +function txt = customDataTip(~, event_obj) + % Retrieve stored variables + v_bias_vals = getappdata(gcf, 'v_bias_vals'); + awg_vpp_vals = getappdata(gcf, 'awg_vpp_vals'); + bers = getappdata(gcf, 'bers'); + power = getappdata(gcf, 'power'); % Retrieve the Power data + + % Get the position of the data cursor + pos = event_obj.Position; + xdata = pos(1); + ydata = pos(2); + + % Find the nearest indices in the data arrays + [~, xInd] = min(abs(v_bias_vals - xdata)); + [~, yInd] = min(abs(awg_vpp_vals - ydata)); + + % Get the corresponding BER and Power values + berValue = bers(yInd, xInd); + powerValue = power(yInd, xInd); % Get the Power value + + % Format the text for the data tip + txt = {['V_{bias} = ', num2str(xdata)], ... + ['V_{pp} AWG = ', num2str(ydata)], ... + ['BER = ', num2str(berValue, '%.1e')], ... + ['Power = ', num2str(powerValue)]}; +end \ No newline at end of file diff --git a/projects/Lab_2024/bias_sweep_evaluation_2.m b/projects/Lab_2024/bias_sweep_evaluation_2.m new file mode 100644 index 0000000..c310afd --- /dev/null +++ b/projects/Lab_2024/bias_sweep_evaluation_2.m @@ -0,0 +1,124 @@ + + +filename = "C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\bias_sweep_gigantisch\wh_pam4.mat"; +a = load(filename); +wh2 = a.wh; + + +v_bias_vals = wh2.parameter.vbias.values; +awg_vpp_vals = wh2.parameter.awg_vpp.values; +eq_mode_vals = wh2.parameter.eq_mode.values; +eq_mode_show = eq_mode_vals(2); +eq_modes = ["FFE","FFE+MLSE","DB precoded","DB encoded"]; + +precomp_amp_max_vals = wh2.parameter.precomp_amp_max.values; +precomp_amp_max_show = precomp_amp_max_vals(2); + +m = wh2.getStoValue('m',wh2.parameter.vbias.values(1),wh2.parameter.awg_vpp.values(1),wh2.parameter.eq_mode.values(1),wh2.parameter.precomp_amp_max.values(1)); + +figure(); +sgtitle(['PAM ', num2str(m),' | EQ: ', char(eq_modes(eq_mode_show))]) + + +for p = 1:numel(precomp_amp_max_vals) + precomp_amp_max_show = precomp_amp_max_vals(p); + subplot1 = subplot(2,3,p); + + bers = []; + rop_measured = []; + cnt = 0; + for awg_vpp_cur = awg_vpp_vals + cnt = cnt+1; + bers(cnt,:) = wh2.getStoValue('ber',v_bias_vals,awg_vpp_cur,eq_mode_show,precomp_amp_max_show); + rop_measured(cnt,:) = wh2.getStoValue('rop',v_bias_vals,awg_vpp_cur,eq_mode_show,precomp_amp_max_show); + end + + [bestber,bestindex] = min(bers,[],'all'); + [awg_pos,v_bias_pos]=ind2sub(size(bers),bestindex); + bestawgvpp=awg_vpp_vals(awg_pos); + bestvbias=v_bias_vals(v_bias_pos); + + disp(['Best Vpp: ',num2str(bestvbias),' V; Best Vpp AWG: ',num2str(bestawgvpp),' V' ]) + + % Compute the logarithm of BER data + % Adding a small epsilon to avoid log(0) + epsilon = 1e-12; + log_bers = log10(bers + epsilon); + + % Set limits for z-data scaling in log scale + zmin = log10(1e-4 + epsilon); + zmax = log10(0.5 + epsilon); + + % Plot the filled contour plot with log-scaled z-data + contourf_handle = contourf(v_bias_vals, awg_vpp_vals, log_bers, 'Parent', subplot1, "ShowText",true,"LabelFormat", @mylabelfun); + + % Set x and y labels with subscripts for clarity + xlabel('V_{bias}'); + ylabel('V_{pp} AWG'); + title(['Prec. Ampl.: ',num2str(precomp_amp_max_show), 'dB']); + + % Adjust the grid to display white lines + grid on; + set(subplot1, 'GridColor', [1 1 1]); % Set grid color to white + + % Set limits for z-data scaling + clim([zmin zmax]); + + % Adjust the colormap + colormap(flipud(cbrewer2('RdBu',64))); + + % Add a colorbar and adjust its ticks to represent actual BER values + c = colorbar; + % Set colorbar ticks at log-spaced intervals + tick_values = [1e-4 1e-3 1e-2 1e-1 0.5]; + tick_positions = log10(tick_values + epsilon); + set(c, 'Ticks', tick_positions, 'TickLabels', arrayfun(@num2str, tick_values, 'UniformOutput', false)); + + % Store variables in the figure's application data for use in the data tip function + setappdata(gcf, 'v_bias_vals', v_bias_vals); + setappdata(gcf, 'awg_vpp_vals', awg_vpp_vals); + setappdata(gcf, 'bers', bers); + setappdata(gcf, 'power', rop_measured); % Store the Power data + + % Set up the data cursor mode to display custom data tips + dcm_obj = datacursormode(gcf); + set(dcm_obj, 'UpdateFcn', @customDataTip); + hold on + scatter(bestvbias,bestawgvpp,100,"red",'Marker','x','LineWidth',2); + +end + + +function labels = mylabelfun(vals) + lab = 10.^vals; + labels = arrayfun(@(x) num2str(x, '%.1e'), lab, 'UniformOutput', false); +end + + +% Define the custom data tip function +function txt = customDataTip(~, event_obj) + % Retrieve stored variables + v_bias_vals = getappdata(gcf, 'v_bias_vals'); + awg_vpp_vals = getappdata(gcf, 'awg_vpp_vals'); + bers = getappdata(gcf, 'bers'); + power = getappdata(gcf, 'power'); % Retrieve the Power data + + % Get the position of the data cursor + pos = event_obj.Position; + xdata = pos(1); + ydata = pos(2); + + % Find the nearest indices in the data arrays + [~, xInd] = min(abs(v_bias_vals - xdata)); + [~, yInd] = min(abs(awg_vpp_vals - ydata)); + + % Get the corresponding BER and Power values + berValue = bers(yInd, xInd); + powerValue = power(yInd, xInd); % Get the Power value + + % Format the text for the data tip + txt = {['V_{bias} = ', num2str(xdata)], ... + ['V_{pp} AWG = ', num2str(ydata)], ... + ['BER = ', num2str(berValue, '%.1e')], ... + ['Power = ', num2str(powerValue)]}; +end \ No newline at end of file diff --git a/projects/Lab_2024/lab_baudrate_sweep.m b/projects/Lab_2024/lab_baudrate_sweep.m new file mode 100644 index 0000000..43b88a2 --- /dev/null +++ b/projects/Lab_2024/lab_baudrate_sweep.m @@ -0,0 +1,274 @@ + +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\baudrate_sweep\'; +experiment_name = 'PAM4_10km_ffe_'; + +ffe_only = 0; +postfilter_approach = 0; +db_channel_approach = 1; +db_coding_approach = 0; + +db_precode = db_coding_approach || db_channel_approach; + +%%% SIR Sweep for MPI Experiment %%% +params = struct; +params.fsym = [56,68,80,92].*1e9; +params.fsym = [92].*1e9; + +wh = DataStorage(params); + +wh.addStorage("ber"); +wh.addStorage("pd_in"); +wh.addStorage("rop"); +wh.addStorage("m"); +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 = 2; + +M = 4; +pn_key = 2; +usemrds = 0; +fdac = 92e9; +awg_vpp = 0.15; +fadc = 160e9; +rrcalpha = 0.05; +v_bias = 2.25; +pd_in_set = 4; +rop_atten = 0; + +disp(['Start Measurement of ',num2str(prod(wh.dim)),' loops...']) + +for fsym = wh.parameter.fsym.values + + loop_name = ['_fsym_',num2str(fsym)]; + + %%%%% SET Voltages %%%%%% + dcs = DC_supply("active",[1,1],"voltage",[v_bias, 9]); + dcs.set("voltage",[v_bias, 9]); + + %%%%% SET Attenuator %%%%%% + voa = OptAtten("active",[1,2,1,1],"value",[rop_atten,pd_in_set,0,0],"wavelength",[1310,1310,1310,1310]); + voa.set('active',[1,2,1,1],'value',[rop_atten,pd_in_set,0,0]); + % voa.readvals(); + + %%%%% 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]); + + + %%%%% Symbol Generation %%%%%% + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... + "fs_out",fdac,"applyclipping",0,"clipfactor",1.7,... + "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... + "db_precode",db_precode,... + "mrds_code",usemrds,"mrds_blocklength",512,"db_encode",db_coding_approach).process(); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 % measure channel + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.buildOFDM(); + elseif precomp_mode == 2 % apply precomp + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',precomp_amp_max,'loadPath',precomp_path,'fileName',precomp_fn); + end + + %%%%% Resample to DAC rate %%%%%% + Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); + + + %%%%% Plot and Save Routine 1 %%%%%%%%%%%%%%%%%%%%%%%%% + Digi_sig.spectrum("displayname","Normal Tx","fignum",10); + + save([folderpath,[experiment_name,'bits'],loop_name],"Bits"); + save([folderpath,[experiment_name,'symbols'],loop_name],"Symbols"); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%%% AWG --> Scope %%%%%% + [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Scpe_sig.spectrum("displayname","Scope PSD","fignum",20); + + Scpe_sig.plot("displayname","Scope raw signal","fignum",25); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*fsym); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 + freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + freqresp.plot(); + end + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + + %%%%%% SNR CHEAT - Avg. the measured signal occurences %%%%%% + average_signals = 0; + if average_signals + scope_mean = zeros(size(S{1}.signal)); + for n=1:numel(S) + scope_mean = scope_mean + S{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig.signal = scope_mean; + end + + %%%%% Plot and Save Routine 2 %%%%%%%%%%%%%%%%%%%%%%%%% + save([folderpath,experiment_name,'rx_signal',loop_name],"S"); + Scpe_sig.eye(fsym,M,"fignum",40,"displayname",' after Scope'); + + voa.readvals(); + rop = voa.power_state(1); + pd_in = voa.power_state(2); + + + %%%%% 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 = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.05,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + 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); + + disp(['FFE: ',sprintf('%.1E',ber),'| ROP: ',num2str(rop),' dB | PD_in: ',num2str(pd_in),' dBm']); + + if 1 + 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)); + hold on + histogram(received(lvl,:),1000,"EdgeAlpha",0); + end + end + + elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + Noi = EQ_sig-Symbols; + + 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); + + nc = 2; + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig = EQ_sig.filter(burg_coeff,1); + + if 0 + 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 + + 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),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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); + + EQ_sig.plot("fignum",50,"displayname",'After EQ'); + + disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + end + + wh.addValueToStorage(ber,'ber',fsym); + wh.addValueToStorage(rop,'rop',fsym); + wh.addValueToStorage(pd_in,'pd_in',fsym); + wh.addValueToStorage(Rx_bits,'signals',fsym); + wh.addValueToStorage(M,'m',fsym); + + showCurrentMeasurement('BER', ber, 'ROP', rop, 'PD in', pd_in, 'PAM',M, 'Vbias', v_bias, 'AWG Vpp', Awg_vpp); + autoArrangeFigures(3,3,2); +end + +wh.save([folderpath,experiment_name,'_wh']); + +cols = linspecer(8); + +fsym_vals = wh.parameter.fsym.values; + +bers = wh.getStoValue('ber',fsym_vals); +rop_measured = wh.getStoValue('rop',fsym_vals); +pd_in_measured = wh.getStoValue('pd_in',fsym_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(fsym_vals.*1e-9, bers, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", experiment_name); + +% Store pd_in_measured in the ZData property +hLine.ZData = pd_in_measured; + +% 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" + + +% Add a new data tip row for PDin +pdinRow = dataTipTextRow('PDin', 'ZData'); +hLine.DataTipTemplate.DataTipRows(3) = pdinRow; + +% Continue with the rest of your plot settings +yline(3.8e-3, 'DisplayName', 'HD-FEC', 'LineStyle', '--', 'HandleVisibility', 'off'); +xlabel('Symbol Rate'); +ylabel('Bit Error Rate (BER)'); +title('Bit Error Rate vs. ROP'); +set(gca, 'yscale', 'log'); +set(gca, 'Box', 'on'); +grid on; +grid minor; +legend('Interpreter', 'none'); + +autoArrangeFigures(3,3,2) + diff --git a/projects/Lab_2024/lab_bias_sweep.m b/projects/Lab_2024/lab_bias_sweep.m new file mode 100644 index 0000000..b22f415 --- /dev/null +++ b/projects/Lab_2024/lab_bias_sweep.m @@ -0,0 +1,269 @@ + +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\bias_sweep\'; +experiment_name = 'PAM4_DB_precoded_10km_ffe_'; + +% a = load([folderpath,experiment_name,'_wh']); +% wh2 = a.obj; + + +ffe_only = 0; +postfilter_approach = 0; +db_channel_approach = 1; +db_coding_approach = 0; + +db_precode = db_coding_approach || db_channel_approach; + +%%% SIR Sweep for MPI Experiment %%% +params = struct; + +params.vbias = [2.1:0.05:2.4]; +params.awg_vpp = [0.15:0.05:0.6]; + +wh = DataStorage(params); + +wh.addStorage("ber"); +wh.addStorage("pd_in"); +wh.addStorage("rop"); +wh.addStorage("m"); +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 = -2; + +M = 4; +pn_key = 2; +usemrds = 0; +fsym = 92e9; +fdac = 92e9; +Awg_vpp = 0.35; +fadc = 160e9; +rrcalpha = 0.05; +v_bias = 2.25; +pd_in_set = 6; +rop_atten = 0; + +looptatal = prod(wh.dim); +disp(['Start Measurement of ',num2str(looptatal),' loops...']) +hWaitbar = waitbar(0, 'Starting measurement...', 'Name', 'Processing Progress'); +loopcnt = 0; + +for v_bias = wh.parameter.vbias.values + for awg_vpp = wh.parameter.awg_vpp.values + + loopcnt = loopcnt+1; + progressFraction = loopcnt / looptatal; + waitbar(progressFraction, hWaitbar, ... + sprintf('Progress: %d/%d', loopcnt, looptatal)); + + try + a + % ber = wh2.getStoValue('ber', v_bias,awg_vpp); + % rop = wh2.getStoValue('rop', v_bias,awg_vpp); + % pd_in = wh2.getStoValue('pd_in', v_bias,awg_vpp); + % Rx_bits = wh2.getStoValue('signals', v_bias,awg_vpp); + % Rx_bits = Rx_bits{1}; + % M = wh2.getStoValue('m', v_bias,awg_vpp); + + + catch + + loop_name = ['_fsym_',num2str(fsym)]; + + %%%%% SET Voltages %%%%%% + dcs = DC_supply("active",[1,1],"voltage",[v_bias, 9]); + dcs.set("voltage",[v_bias, 9]); + + %%%%% SET Attenuator %%%%%% + voa = OptAtten("active",[1,2,1,1],"value",[rop_atten,pd_in_set,0,0],"wavelength",[1310,1310,1310,1310]); + voa.set('active',[1,2,1,1],'value',[rop_atten,pd_in_set,0,0]); + % voa.readvals(); + + %%%%% Construct AWG and Scope Modules %%%%%% + SCP = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc",'GSa_160',"channel",[1,0],"recordLen",1000000,"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]); + + + %%%%% Symbol Generation %%%%%% + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... + "fs_out",fdac,"applyclipping",0,"clipfactor",1.7,... + "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... + "db_precode",db_precode,... + "mrds_code",usemrds,"mrds_blocklength",512,"db_encode",db_coding_approach).process(); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 % measure channel + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.buildOFDM(); + elseif precomp_mode == 2 % apply precomp + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',precomp_amp_max,'loadPath',precomp_path,'fileName',precomp_fn); + end + + %%%%% Resample to DAC rate %%%%%% + Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); + + + %%%%% Plot and Save Routine 1 %%%%%%%%%%%%%%%%%%%%%%%%% + % Digi_sig.spectrum("displayname","Normal Tx","fignum",10); + + % save([folderpath,[experiment_name,'bits'],loop_name],"Bits"); + % save([folderpath,[experiment_name,'symbols'],loop_name],"Symbols"); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%%% AWG --> Scope %%%%%% + [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Scpe_sig.spectrum("displayname","Scope PSD","fignum",20); + + % Scpe_sig.plot("displayname","Scope raw signal","fignum",25); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*fsym); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 + freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + freqresp.plot(); + end + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + + %%%%%% SNR CHEAT - Avg. the measured signal occurences %%%%%% + average_signals = 0; + if average_signals + scope_mean = zeros(size(S{1}.signal)); + for n=1:numel(S) + scope_mean = scope_mean + S{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig.signal = scope_mean; + end + + %%%%% Plot and Save Routine 2 %%%%%%%%%%%%%%%%%%%%%%%%% + save([folderpath,experiment_name,'rx_signal',loop_name],"S"); + Scpe_sig.eye(fsym,M,"fignum",40,"displayname",' after Scope'); + + voa.readvals(); + rop = voa.power_state(1); + pd_in = voa.power_state(2); + + + %%%%% 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 = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.05,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + 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); + + disp(['FFE: ',sprintf('%.1E',ber),'| ROP: ',num2str(rop),' dB | PD_in: ',num2str(pd_in),' dBm']); + + if 1 + 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)); + hold on + histogram(received(lvl,:),1000,"EdgeAlpha",0); + end + end + + elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + Noi = EQ_sig-Symbols; + + 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); + + nc = 2; + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig = EQ_sig.filter(burg_coeff,1); + + if 0 + 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 + + 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']); + + elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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); + + EQ_sig.plot("fignum",50,"displayname",'After EQ'); + + disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | PD_in: ',num2str(pd_in),' dBm']); + + end + + end + + wh.addValueToStorage(ber,'ber',v_bias,awg_vpp); + wh.addValueToStorage(rop,'rop',v_bias,awg_vpp); + wh.addValueToStorage(pd_in,'pd_in',v_bias,awg_vpp); + wh.addValueToStorage(Rx_bits,'signals',v_bias,awg_vpp); + wh.addValueToStorage(M,'m',v_bias,awg_vpp); + + showCurrentMeasurement('BER', ber, 'ROP', rop, 'PD in', pd_in, 'PAM',M, 'Vbias', v_bias, 'AWG Vpp', awg_vpp, 'Precomp MaxAmp',precomp_amp_max); + autoArrangeFigures(3,3,2); + + + end +end + +close(hWaitbar); + +wh.save([folderpath,experiment_name,'_wh']); + +autoArrangeFigures(3,3,2) + diff --git a/projects/Lab_2024/lab_bias_sweep_gigantisch.m b/projects/Lab_2024/lab_bias_sweep_gigantisch.m new file mode 100644 index 0000000..d00c0ca --- /dev/null +++ b/projects/Lab_2024/lab_bias_sweep_gigantisch.m @@ -0,0 +1,320 @@ + +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\bias_sweep_4db_pdin\'; +experiment_name = 'PAM6_alles_10km_ffe_'; + +wh2 = obj; + +ffe_only = 0; +postfilter_approach = 0; +db_channel_approach = 1; +db_coding_approach = 0; + +db_precode = db_coding_approach || db_channel_approach; + +%%% SIR Sweep for MPI Experiment %%% +params = struct; + +params.vbias = [2.2:0.05:2.5]; +params.awg_vpp = [0.15:0.05:0.6]; +params.eq_mode = [1]; +params.precomp_amp_max = [0:2:5]; + + +wh = DataStorage(params); + +wh.addStorage("ber"); +wh.addStorage("pd_in"); +wh.addStorage("rop"); +wh.addStorage("m"); +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 = -2; + +M = 6; +pn_key = 2; +usemrds = 0; +fsym = 70e9; +fdac = 92e9; +awg_vpp = 0.35; +fadc = 160e9; +rrcalpha = 0.05; +v_bias = 2.25; +pd_in_set = 4; +rop_atten = 0; + +looptatal = prod(wh.dim); +iterationTimes = zeros(looptatal, 1); % Preallocate for speed + +disp(['Start Measurement of ',num2str(looptatal),' loops...']) +hWaitbar = waitbar(0, 'Starting measurement...', 'Name', 'Processing Progress'); +loopcnt = 0; + +estimatedTimeRemaining = 0; +estimatedTotalTime = 0; +for eq_mode = wh.parameter.eq_mode.values + for precomp_amp_max = wh.parameter.precomp_amp_max.values + for v_bias = wh.parameter.vbias.values + for awg_vpp = wh.parameter.awg_vpp.values + + iterationStartTime = tic; + experiment_name = ['PAM6_alles_10km_eq',num2str(eq_mode),'maxamp',num2str(precomp_amp_max),'vbias',num2str(v_bias),'awgvpp',num2str(awg_vpp)]; + experiment_name = strrep(experiment_name,'.','_'); + + loopcnt = loopcnt+1; + progressFraction = loopcnt / looptatal; + waitbar(progressFraction, hWaitbar, ... + sprintf('Progress: %d/%d\nEstimated time remaining: %.2f hours\nEstimated time remaining: %.2f hours', ... + loopcnt, looptatal, estimatedTimeRemaining/60/60, estimatedTotalTime/60/60)); + + + loop_name = ['_fsym_',num2str(fsym)]; + + try + + a + 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 + + case 2 + ffe_only = 0; + postfilter_approach = 1; + db_channel_approach = 0; + db_coding_approach = 0; + + db_precode = db_coding_approach || db_channel_approach; + case 3 + ffe_only = 0; + postfilter_approach = 0; + db_channel_approach = 1; + db_coding_approach = 0; + + db_precode = db_coding_approach || db_channel_approach; + case 4 + ffe_only = 0; + postfilter_approach = 0; + db_channel_approach = 0; + db_coding_approach = 1; + + db_precode = db_coding_approach || db_channel_approach; + end + + + + %%%%% SET Voltages %%%%%% + dcs = DC_supply("active",[1,1],"voltage",[v_bias, 9]); + dcs.set("voltage",[v_bias, 9]); + + %%%%% SET Attenuator %%%%%% + voa = OptAtten("active",[1,2,1,1],"value",[rop_atten,pd_in_set,0,0],"wavelength",[1310,1310,1310,1310]); + voa.set('active',[1,2,1,1],'value',[rop_atten,pd_in_set,0,0]); + % voa.readvals(); + + %%%%% Construct AWG and Scope Modules %%%%%% + SCP = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc",'GSa_160',"channel",[1,0],"recordLen",1000000,"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]); + + %%%%% Symbol Generation %%%%%% + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... + "fs_out",fdac,"applyclipping",0,"clipfactor",1.7,... + "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... + "db_precode",db_precode,... + "mrds_code",usemrds,"mrds_blocklength",512,"db_encode",db_coding_approach).process(); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 % measure channel + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.buildOFDM(); + elseif precomp_mode == 2 % apply precomp + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',precomp_amp_max,'loadPath',precomp_path,'fileName',precomp_fn); + end + + %%%%% Resample to DAC rate %%%%%% + Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); + + + %%%%% Plot and Save Routine 1 %%%%%%%%%%%%%%%%%%%%%%%%% + % Digi_sig.spectrum("displayname","Normal Tx","fignum",10); + + % save([folderpath,[experiment_name,'bits'],loop_name],"Bits"); + % save([folderpath,[experiment_name,'symbols'],loop_name],"Symbols"); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%%% AWG --> Scope %%%%%% + [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Scpe_sig.spectrum("displayname","Scope PSD","fignum",20); + + % Scpe_sig.plot("displayname","Scope raw signal","fignum",25); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*fsym); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 + freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + freqresp.plot(); + end + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + + %%%%%% SNR CHEAT - Avg. the measured signal occurences %%%%%% + average_signals = 0; + if average_signals + scope_mean = zeros(size(S{1}.signal)); + for n=1:numel(S) + scope_mean = scope_mean + S{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig.signal = scope_mean; + end + + %%%%% Plot and Save Routine 2 %%%%%%%%%%%%%%%%%%%%%%%%% + save([folderpath,experiment_name,'rx_signal',loop_name],"S"); + % Scpe_sig.eye(fsym,M,"fignum",40,"displayname",' after Scope'); + + voa.readvals(); + rop = voa.power_state(1); + pd_in = voa.power_state(2); + + + %%%%% 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 = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.05,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + % 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); + + disp(['FFE: ',sprintf('%.1E',ber),'| ROP: ',num2str(rop),' dB | PD_in: ',num2str(pd_in),' dBm']); + + 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)); + hold on + histogram(received(lvl,:),1000,"EdgeAlpha",0); + end + end + + elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + % EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + Noi = EQ_sig-Symbols; + + 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); + + nc = 2; + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig = EQ_sig.filter(burg_coeff,1); + + if 0 + 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 + + 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']); + + elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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); + + % EQ_sig.plot("fignum",50,"displayname",'After EQ'); + + disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | PD_in: ',num2str(pd_in),' dBm']); + + end + + end + + wh.addValueToStorage(ber,'ber',v_bias,awg_vpp,eq_mode,precomp_amp_max); + wh.addValueToStorage(rop,'rop',v_bias,awg_vpp,eq_mode,precomp_amp_max); + wh.addValueToStorage(pd_in,'pd_in',v_bias,awg_vpp,eq_mode,precomp_amp_max); + %wh.addValueToStorage(Rx_bits,'signals',v_bias,awg_vpp); + wh.addValueToStorage(M,'m',v_bias,awg_vpp,eq_mode,precomp_amp_max); + + showCurrentMeasurement('BER', ber, 'ROP', rop, '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 * looptatal; + estimatedTimeRemaining = estimatedTotalTime - sum(iterationTimes(1:loopcnt)); + %autoArrangeFigures(3,3,2); + + end + end + end +end + +close(hWaitbar); + +wh.save([folderpath,experiment_name,'wh']); + +autoArrangeFigures(3,3,2) + +disp("measurement done") \ No newline at end of file diff --git a/projects/Lab_2024/lab_db_precode.m b/projects/Lab_2024/lab_db_precode.m index 1cdcd1a..91900e7 100644 --- a/projects/Lab_2024/lab_db_precode.m +++ b/projects/Lab_2024/lab_db_precode.m @@ -1,11 +1,11 @@ -folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\mpi_ofc_2024\'; -experiment_name = '10km_db_transmit_no_mpi_'; +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\no_mpi_2024\'; +experiment_name = '10km_ffe_no_mpi_'; only_dsp = 0; -ffe_only = 0; -postfilter_approach = 1; +ffe_only = 1; +postfilter_approach = 0; db_channel_approach = 0; db_coding_approach = 0; @@ -25,34 +25,35 @@ wh.addStorage("signals"); precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active precomp_amp_max = 3; -M = 4; +M = 6; pn_key = 2; usemrds = 0; fdac = 92e9; -fsym = 92e9; +fsym = 56e9; fadc = 160e9; rrcalpha = 0.05; v_bias = 2.25; i_atten = params.i_atten(1); -pd_in_desired = 7; +pd_in_desired = 6; -if ~only_dsp - disp(['Start Measurement of ',num2str(prod(wh.dim)),' loops...']) - +disp(['Start Measurement of ',num2str(prod(wh.dim)),' loops...']) + +for i = 1 + %%%%% SET Volatges %%%%%% dcs = DC_supply("active",[1,1],"voltage",[v_bias, 9]); dcs.set("voltage",[v_bias, 9]); - + voa = OptAtten("active",[1,2,1,1],"value",[0,pd_in_desired,0,i_atten],"wavelength",[1310,1310,1310,1310]); - % voa.set('active',[1,2,1,1],'value',[0,pd_in_desired,0,i_atten]); + %voa.set('active',[1,2,1,1],'value',[0,pd_in_desired,0,i_atten]); % voa.readvals(); - + 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,0.62]); A2S = Awg2Scope(AWG,SCP,[0,0,0,1]); - - + + if 1 precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\standard_system\"; precomp_fn = "lab_mpi_setup_2"; @@ -60,17 +61,17 @@ if ~only_dsp precomp_path = "C:\Users\sioe\Documents\MATLAB\model-collection\sioe_models\Labor_2024\Lab_PAM4\"; precomp_fn = "precomp_bla__loop1_1"; end - - + + %%%%% Symbol Generation %%%%%% Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); - + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... "fs_out",fdac,"applyclipping",0,"clipfactor",1.7,... "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... "db_precode",db_precode,... "mrds_code",usemrds,"mrds_blocklength",512,"db_encode",db_coding_approach).process(); - + if precomp_mode == 1 % measure channel freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); Digi_sig = freqresp.buildOFDM(); @@ -78,50 +79,42 @@ if ~only_dsp freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',precomp_amp_max,'loadPath',precomp_path,'fileName',precomp_fn); end - + Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); Digi_sig.spectrum("displayname","Normal Tx","fignum",10); - % Digi_sig = Filter('filtdegree',1,"f_cutoff",45e9,"fs",Digi_sig.fs,"filterType",filtertypes.butterworth,"active",true).process(Digi_sig); - % Digi_sig.spectrum("displayname","Lowpass Tx","fignum",999); - % Digi_sig.eye(fsym,M); - - save([folderpath,[experiment_name,'bits']],"Bits"); save([folderpath,[experiment_name,'symbols']],"Symbols"); -end -for i_atten = wh.parameter.i_atten.values - - if ~only_dsp + for i_atten = wh.parameter.i_atten.values %%%%% SET ATTENUATOR %%%%%% voa.set('active',[1,2,1,1],'value',[0,7,0,i_atten]); - + %%%%% AWG --> Scope %%%%%% [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); Scpe_sig.spectrum("displayname","Scope PSD","fignum",20); - - + + % Scpe_sig.spectrum("displayname",'Rx Signal','fignum',10); - + %%%%%% Sample to 2x fsym %%%%%% Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*fsym); - + if precomp_mode == 1 freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); freqresp.plot(); end - + %%%%%% Sync Rx signal with reference %%%%%% [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); save([folderpath,experiment_name,'rx_signal_iatten_',num2str(i_atten),''],"S"); - + average_signals = 0; if average_signals scope_mean = zeros(size(S{1}.signal)); @@ -135,96 +128,96 @@ for i_atten = wh.parameter.i_atten.values Scpe_sig.plot("displayname","Scope PSD","fignum",30); Scpe_sig.eye(fsym,M,"fignum",40,"displayname",' after Scope'); - sir = voa.power_state(3)-voa.power_state(4); pd_in = voa.power_state(2); - end - %%%%% 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); + %%%%% 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 = 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); - if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% + if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% - [EQ_sig] = Eq.process(Scpe_sig,Symbols); + [EQ_sig] = Eq.process(Scpe_sig,Symbols); - EQ_sig.plot("fignum",50,"displayname",'After EQ'); + 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); + 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),'| SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + disp(['FFE: ',sprintf('%.1E',ber),'| SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); - elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% - [EQ_sig] = Eq.process(Scpe_sig,Symbols); + [EQ_sig] = Eq.process(Scpe_sig,Symbols); - EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); - Noi = EQ_sig-Symbols; + Noi = EQ_sig-Symbols; - 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); + 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); - nc = 2; - burg_coeff = arburg(Noi.signal,nc); + nc = 2; + burg_coeff = arburg(Noi.signal,nc); - EQ_sig = EQ_sig.filter(burg_coeff,1); + EQ_sig = EQ_sig.filter(burg_coeff,1); + + if 0 + 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 + + 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),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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); + + EQ_sig.plot("fignum",50,"displayname",'After EQ'); + + disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); - if 0 - 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 - 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),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); - - elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% - - [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); - - EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); - - EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); - EQ_sig = Duobinary().decode(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([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); - - elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% - - [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); - EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); - EQ_sig = Duobinary().decode(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); - - EQ_sig.plot("fignum",50,"displayname",'After EQ'); - - disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + wh.addValueToStorage(ber,'ber',i_atten); + wh.addValueToStorage(sir,'sir',i_atten); + wh.addValueToStorage(pd_in,'pd_in',i_atten); + wh.addValueToStorage(Rx_bits,'signals',i_atten); end - - - - wh.addValueToStorage(ber,'ber',i_atten); - wh.addValueToStorage(sir,'sir',i_atten); - wh.addValueToStorage(pd_in,'pd_in',i_atten); - wh.addValueToStorage(Rx_bits,'signals',i_atten); - end wh.save([folderpath,experiment_name,'_wh']); diff --git a/projects/Lab_2024/lab_modulator_tf_sweep.m b/projects/Lab_2024/lab_modulator_tf_sweep.m new file mode 100644 index 0000000..6c41e4a --- /dev/null +++ b/projects/Lab_2024/lab_modulator_tf_sweep.m @@ -0,0 +1,68 @@ +% 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; + +% Create a DataStorage object with the defined parameters +wh = DataStorage(params); + +% Add a storage field for output power measurements +wh.addStorage("p_out"); + +% 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 + "value", [0, 0, 0, 0], ... % Set attenuation values to 0 dB + "wavelength", [1310, 1310, 1310, 1310]); % Set the wavelength for each channel + +% 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 + +% 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 + + % Retrieve all stored output power measurements up to the current point + p_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 + grid on; % Enable grid for better readability +end + +% After completing the measurements, retrieve all output power data +p_out = wh.getStoValue('p_out', wh.parameter.v_bias.values); + +% Plot the output power converted from dBm to linear scale (Watts) +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 +grid on; % Enable grid for better readability diff --git a/projects/Lab_2024/lab_precompensation_sweep.m b/projects/Lab_2024/lab_precompensation_sweep.m new file mode 100644 index 0000000..1786ebc --- /dev/null +++ b/projects/Lab_2024/lab_precompensation_sweep.m @@ -0,0 +1,265 @@ + +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\precompensation_sweep\'; +experiment_name = 'PAM4_DBencode_10km_'; + +% a = load([folderpath,experiment_name,'_wh']); +% wh2 = a.obj; + + +ffe_only = 0; +postfilter_approach = 0; +db_channel_approach = 0; +db_coding_approach = 1; + +db_precode = db_coding_approach || db_channel_approach; + +%%% SIR Sweep for MPI Experiment %%% +params = struct; + +params.precomp_amp_max = [-6:6]; + +wh = DataStorage(params); + +wh.addStorage("ber"); +wh.addStorage("pd_in"); +wh.addStorage("rop"); +wh.addStorage("m"); +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 = -2; + +M = 4; +pn_key = 2; +usemrds = 0; +fsym = 92e9; +fdac = 92e9; +awg_vpp = 0.2; +fadc = 160e9; +rrcalpha = 0.05; +v_bias = 2.35; +pd_in_set = 6; +rop_atten = 0; + +looptatal = prod(wh.dim); +disp(['Start Measurement of ',num2str(looptatal),' loops...']) +hWaitbar = waitbar(0, 'Starting measurement...', 'Name', 'Processing Progress'); +loopcnt = 0; + +for precomp_amp_max = wh.parameter.precomp_amp_max.values + + loopcnt = loopcnt+1; + progressFraction = loopcnt / looptatal; + waitbar(progressFraction, hWaitbar, ... + sprintf('Progress: %d/%d', loopcnt, looptatal)); + + + loop_name = ['_fsym_',num2str(fsym)]; + + %%%%% SET Voltages %%%%%% + dcs = DC_supply("active",[1,1],"voltage",[v_bias, 9]); + dcs.set("voltage",[v_bias, 9]); + + %%%%% SET Attenuator %%%%%% + voa = OptAtten("active",[1,2,1,1],"value",[rop_atten,pd_in_set,0,0],"wavelength",[1310,1310,1310,1310]); + voa.set('active',[1,2,1,1],'value',[rop_atten,pd_in_set,0,0]); + % voa.readvals(); + + %%%%% Construct AWG and Scope Modules %%%%%% + SCP = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc",'GSa_160',"channel",[1,0],"recordLen",1000000,"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]); + + + %%%%% Symbol Generation %%%%%% + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... + "fs_out",fdac,"applyclipping",0,"clipfactor",1.7,... + "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... + "db_precode",db_precode,... + "mrds_code",usemrds,"mrds_blocklength",512,"db_encode",db_coding_approach).process(); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 % measure channel + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.buildOFDM(); + elseif precomp_mode == 2 % apply precomp + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',precomp_amp_max,'loadPath',precomp_path,'fileName',precomp_fn); + end + + %%%%% Resample to DAC rate %%%%%% + Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); + + + %%%%% Plot and Save Routine 1 %%%%%%%%%%%%%%%%%%%%%%%%% + % Digi_sig.spectrum("displayname","Normal Tx","fignum",10); + + % save([folderpath,[experiment_name,'bits'],loop_name],"Bits"); + % save([folderpath,[experiment_name,'symbols'],loop_name],"Symbols"); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%%% AWG --> Scope %%%%%% + [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Scpe_sig.spectrum("displayname","Scope PSD","fignum",20); + + % Scpe_sig.plot("displayname","Scope raw signal","fignum",25); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*fsym); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 + freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + freqresp.plot(); + end + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + + %%%%%% SNR CHEAT - Avg. the measured signal occurences %%%%%% + average_signals = 0; + if average_signals + scope_mean = zeros(size(S{1}.signal)); + for n=1:numel(S) + scope_mean = scope_mean + S{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig.signal = scope_mean; + end + + %%%%% Plot and Save Routine 2 %%%%%%%%%%%%%%%%%%%%%%%%% + save([folderpath,experiment_name,'rx_signal',loop_name],"S"); + Scpe_sig.eye(fsym,M,"fignum",40,"displayname",' after Scope'); + + voa.readvals(); + rop = voa.power_state(1); + pd_in = voa.power_state(2); + + + %%%%% 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 = EQ("Ne",[50,7,7],"Nb",[0,0,0],"training_length",4096*2,"training_loops",5,"dd_loops",5,"K",2,"DCmu",0.05,"DDmu",[0.0004 0.0004 0.0004 0.0004 ],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + + if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + 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); + + disp(['FFE: ',sprintf('%.1E',ber),'| ROP: ',num2str(rop),' dB | PD_in: ',num2str(pd_in),' dBm']); + + if 1 + 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)); + hold on + histogram(received(lvl,:),1000,"EdgeAlpha",0); + end + end + + elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + Noi = EQ_sig-Symbols; + + 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); + + nc = 2; + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig = EQ_sig.filter(burg_coeff,1); + + if 0 + 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 + + 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']); + + elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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); + + EQ_sig.plot("fignum",50,"displayname",'After EQ'); + + disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | PD_in: ',num2str(pd_in),' dBm']); + + end + + + wh.addValueToStorage(ber,'ber',precomp_amp_max); + wh.addValueToStorage(rop,'rop',precomp_amp_max); + wh.addValueToStorage(pd_in,'pd_in',precomp_amp_max); + wh.addValueToStorage(Rx_bits,'signals',precomp_amp_max); + wh.addValueToStorage(M,'m',precomp_amp_max); + + showCurrentMeasurement('BER', ber, 'ROP', rop, 'PD in', pd_in, 'PAM',M, 'Vbias', v_bias, 'AWG Vpp', awg_vpp, 'Precomp MaxAmp',precomp_amp_max); + autoArrangeFigures(3,3,2); + + +end + +close(hWaitbar); + +wh.save([folderpath,experiment_name,'_wh']); + +autoArrangeFigures(3,3,2) + +figure(90) +amp_vals = wh.parameter.precomp_amp_max.values; +ber = wh.getStoValue('ber',wh.parameter.precomp_amp_max.values); +plot(amp_vals,ber,'DisplayName',['PAM ',num2str(M)]); +xlabel('Precompensation Max Amp'); +ylabel('BER'); +grid on +grid minor +title('Bit Error Rate vs. ROP'); +set(gca, 'yscale', 'log'); +legend diff --git a/projects/Lab_2024/lab_rop_sweep.m b/projects/Lab_2024/lab_rop_sweep.m new file mode 100644 index 0000000..0e52e1a --- /dev/null +++ b/projects/Lab_2024/lab_rop_sweep.m @@ -0,0 +1,261 @@ + +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\no_mpi_2024\'; +experiment_name = '10km_ffe_no_mpi_'; + +ffe_only = 1; +postfilter_approach = 0; +db_channel_approach = 0; +db_coding_approach = 0; + +db_precode = db_coding_approach || db_channel_approach; + +%%% SIR Sweep for MPI Experiment %%% +params = struct; +params.rop_atten = [7:-1:0]; % high atten to low atten to make sure there is no sudden opening of VOA +params.pd_in_set = [6]; % desired P_out (outp. power mode=2) + +wh = DataStorage(params); + +wh.addStorage("ber"); +wh.addStorage("pd_in"); +wh.addStorage("rop"); +wh.addStorage("signals"); + +precomp_mode = 2; %0=do nothing ; 1= measure; 2=precomp active +precomp_amp_max = 3; + +M = 4; +pn_key = 2; +usemrds = 0; +fdac = 92e9; +fsym = 92e9; +fadc = 160e9; +rrcalpha = 0.05; +v_bias = 2.25; + +disp(['Start Measurement of ',num2str(prod(wh.dim)),' loops...']) + +for rop_atten = wh.parameter.rop_atten.values + for pd_in_set = wh.parameter.pd_in_set.values + + loop_name = ['_ropatten_',num2str(rop_atten),'_pdin_',num2str(pd_in_set)]; + + %%%%% SET Voltages %%%%%% + dcs = DC_supply("active",[1,1],"voltage",[v_bias, 9]); + dcs.set("voltage",[v_bias, 9]); + + voa = OptAtten("active",[1,2,1,1],"value",[rop_atten,pd_in_set,0,0],"wavelength",[1310,1310,1310,1310]); + voa.set('active',[1,2,1,1],'value',[rop_atten,pd_in_set,0,0]); + % voa.readvals(); + + SCP = ScopeKeysight("model","DSAZ634A",'autoscale',1,"fadc",'GSa_160',"channel",[1,0],"recordLen",1000000,"removeDC",1); + AWG = AwgKeysight("model","M8196A","fdac",fdac,"scaletodac",[1,1,1,1],"skews",[0,0,0,0],"voltages",[0,0,0,0.62]); + A2S = Awg2Scope(AWG,SCP,[0,0,0,1]); + + if 1 + precomp_path = "C:\Users\sioe\Documents\MATLAB\imdd_simulation\projects\standard_system\"; + precomp_fn = "lab_mpi_setup_2"; + else + precomp_path = "C:\Users\sioe\Documents\MATLAB\model-collection\sioe_models\Labor_2024\Lab_PAM4\"; + precomp_fn = "precomp_bla__loop1_1"; + end + + + %%%%% Symbol Generation %%%%%% + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",18,"useprbs",1,... + "fs_out",fdac,"applyclipping",0,"clipfactor",1.7,... + "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... + "db_precode",db_precode,... + "mrds_code",usemrds,"mrds_blocklength",512,"db_encode",db_coding_approach).process(); + + if precomp_mode == 1 % measure channel + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.buildOFDM(); + elseif precomp_mode == 2 % apply precomp + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',precomp_amp_max,'loadPath',precomp_path,'fileName',precomp_fn); + end + + Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); + + Digi_sig.spectrum("displayname","Normal Tx","fignum",10); + + + save([folderpath,[experiment_name,'bits'],loop_name],"Bits"); + save([folderpath,[experiment_name,'symbols'],loop_name],"Symbols"); + + %%%%% AWG --> Scope %%%%%% + [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); + + Scpe_sig.spectrum("displayname","Scope PSD","fignum",20); + + + + % Scpe_sig.spectrum("displayname",'Rx Signal','fignum',10); + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*fsym); + + if precomp_mode == 1 + freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + freqresp.plot(); + end + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + save([folderpath,experiment_name,'rx_signal',loop_name],"S"); + + average_signals = 0; + if average_signals + scope_mean = zeros(size(S{1}.signal)); + for n=1:numel(S) + scope_mean = scope_mean + S{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig.signal = scope_mean; + end + + Scpe_sig.plot("displayname","Scope PSD","fignum",30); + Scpe_sig.eye(fsym,M,"fignum",40,"displayname",' after Scope'); + + voa.readvals(); + rop = voa.power_state(1); + pd_in = voa.power_state(2); + + + %%%%% 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 = 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); + + if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + 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); + + disp(['FFE: ',sprintf('%.1E',ber),'| ROP: ',num2str(rop),' dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + Noi = EQ_sig-Symbols; + + 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); + + nc = 2; + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig = EQ_sig.filter(burg_coeff,1); + + if 0 + 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 + + 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),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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); + + EQ_sig.plot("fignum",50,"displayname",'After EQ'); + + disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + end + + wh.addValueToStorage(ber,'ber',rop_atten,pd_in_set); + wh.addValueToStorage(rop,'rop',rop_atten,pd_in_set); + wh.addValueToStorage(pd_in,'pd_in',rop_atten,pd_in_set); + wh.addValueToStorage(Rx_bits,'signals',rop_atten,pd_in_set); + + showCurrentMeasurement('BER', ber, 'ROP', rop, 'PD in', pd_in); + autoArrangeFigures(3,3,2); + end +end + +wh.save([folderpath,experiment_name,'_wh']); + +cols = linspecer(8); + +rop_vals = wh.parameter.rop_atten.values; +pd_in_set = wh.parameter.pd_in_set.values(1); + +bers = wh.getStoValue('ber',rop_vals,pd_in_set); +rop_measured = wh.getStoValue('rop',rop_vals,pd_in_set); +pd_in_measured = wh.getStoValue('pd_in',rop_vals,pd_in_set); + +s = wh.getStoValue('signals',rop_vals,pd_in_set); + +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(rop_measured, bers, "LineWidth", 0.5, "LineStyle", "-", "Marker", ".", "MarkerSize", 15, "DisplayName", experiment_name); + +% Store pd_in_measured in the ZData property +hLine.ZData = pd_in_measured; + +% Customize the data tips +% Set labels for existing data tip rows +hLine.DataTipTemplate.DataTipRows(1).Label = 'ROP'; +hLine.DataTipTemplate.DataTipRows(2).Label = 'BER'; +hLine.DataTipTemplate.DataTipRows(2).Format = '%.2e'; % Format BER as "3e-4" + + +% Add a new data tip row for PDin +pdinRow = dataTipTextRow('PDin', 'ZData'); +hLine.DataTipTemplate.DataTipRows(3) = pdinRow; + +% Continue with the rest of your plot settings +yline(3.8e-3, 'DisplayName', 'HD-FEC', 'LineStyle', '--', 'HandleVisibility', 'off'); +xlabel('Received Optical Power (dBm)'); +ylabel('Bit Error Rate (BER)'); +title('Bit Error Rate vs. ROP'); +set(gca, 'yscale', 'log'); +set(gca, 'Box', 'on'); +grid on; +grid minor; +legend('Interpreter', 'none'); + +autoArrangeFigures(3,3,2) + diff --git a/projects/Lab_2024/lab_sir_sweep.m b/projects/Lab_2024/lab_sir_sweep.m new file mode 100644 index 0000000..2bb3efc --- /dev/null +++ b/projects/Lab_2024/lab_sir_sweep.m @@ -0,0 +1,368 @@ + +folderpath = 'C:\Users\sioe\Nextcloud\Dokumente\02_Ablage_Office\Lab_Data_24\sir_sweep_pam4\'; +experiment_name = 'PAM4_DB_encoded_10km_'; + +%%% SIR Sweep for MPI Experiment %%% +params = struct; + +params.vbias = [2.45]; +params.awg_vpp = [0.25]; +params.eq_mode = [4]; +params.i_atten = [0:4:40]; + +wh = DataStorage(params); + +wh.addStorage("ber"); +wh.addStorage("pd_in"); +wh.addStorage("m"); +wh.addStorage("sir"); +wh.addStorage("s_pow"); +wh.addStorage("i_pow"); +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; + +M = 4; +pn_key = 2; +usemrds = 0; +fsym = 92e9; +fdac = 92e9; +awg_vpp = 0.35; +fadc = 160e9; +rrcalpha = 0.05; +v_bias = 2.25; +pd_in_set = 6; + +looptotal = prod(wh.dim); +iterationTimes = zeros(looptotal, 1); % Preallocate for speed + +disp(['Start Measurement of ',num2str(looptotal),' loops...']) +hWaitbar = waitbar(0, 'Starting measurement...', 'Name', 'Processing Progress'); +loopcnt = 0; + +estimatedTimeRemaining = 0; +estimatedTotalTime = 0; +for eq_mode = wh.parameter.eq_mode.values + for i_atten = wh.parameter.i_atten.values + for v_bias = wh.parameter.vbias.values + for awg_vpp = wh.parameter.awg_vpp.values + + iterationStartTime = tic; + + loop_name = ['_iatten_',num2str(i_atten)]; + + loopcnt = loopcnt+1; + progressFraction = loopcnt / looptotal; + waitbar(progressFraction, hWaitbar, ... + 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 + case 1 + ffe_only = 1; + postfilter_approach = 0; + db_channel_approach = 0; + db_coding_approach = 0; + + db_precode = db_coding_approach || db_channel_approach; + + case 2 + ffe_only = 0; + postfilter_approach = 1; + db_channel_approach = 0; + db_coding_approach = 0; + + db_precode = db_coding_approach || db_channel_approach; + case 3 + ffe_only = 0; + postfilter_approach = 0; + db_channel_approach = 1; + db_coding_approach = 0; + + db_precode = db_coding_approach || db_channel_approach; + case 4 + ffe_only = 0; + postfilter_approach = 0; + db_channel_approach = 0; + db_coding_approach = 1; + + db_precode = db_coding_approach || db_channel_approach; + end + + + + %%%%% SET Voltages %%%%%% + dcs = DC_supply("active",[1,1],"voltage",[v_bias, 9]); + dcs.set("voltage",[v_bias, 9]); + + %%%%% SET Attenuator %%%%%% + voa = OptAtten("active",[1,2,1,1],"value",[0,pd_in_set,0,i_atten],"wavelength",[1310,1310,1310,1310]); + voa.set('active',[1,2,1,1],'value',[0,pd_in_set,0,i_atten]); + % voa.readvals(); + + %%%%% 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]); + + %%%%% Symbol Generation %%%%%% + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rrcalpha); + + [Digi_sig,Symbols,Bits] = PAMsource("fsym",fsym,"M",M,"order",19,"useprbs",1,... + "fs_out",fdac,"applyclipping",0,"clipfactor",1.7,... + "applypulseform",0,"pulseformer",Pform,"randkey",pn_key,... + "db_precode",db_precode,... + "mrds_code",usemrds,"mrds_blocklength",512,"db_encode",db_coding_approach).process(); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 % measure channel + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.buildOFDM(); + elseif precomp_mode == 2 % apply precomp + freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',Digi_sig.fs); + Digi_sig = freqresp.precomp(Digi_sig,'maxampdb',precomp_amp_max,'loadPath',precomp_path,'fileName',precomp_fn); + end + + %%%%% Resample to DAC rate %%%%%% + Digi_sig = Digi_sig.resample("fs_out",AWG.fdac); + + + %%%%% 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"); + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%%% AWG --> Scope %%%%%% + [~,~,~,Scpe_sig] = A2S.process("signal4",Digi_sig); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Scpe_sig.spectrum("displayname","Scope PSD","fignum",20); + + % Scpe_sig.plot("displayname","Scope raw signal","fignum",25); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_in",160e9,"fs_out",2*fsym); + + %%%%% Precompensation Routine %%%%%% + if precomp_mode == 1 + freqresp.estimate(Scpe_sig,"save",true,"savePath",precomp_path,"fileName",precomp_fn); + freqresp.plot(); + end + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + + %%%%%% SNR CHEAT - Avg. the measured signal occurences %%%%%% + average_signals = 0; + if average_signals + scope_mean = zeros(size(S{1}.signal)); + for n=1:numel(S) + scope_mean = scope_mean + S{n}.signal; + end + scope_mean = scope_mean ./ n; + Scpe_sig.signal = scope_mean; + end + + %%%%% Plot and Save Routine 2 %%%%%%%%%%%%%%%%%%%%%%%%% + save([folderpath,experiment_name,'rx_signal',loop_name],"S"); + % 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 = 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); + + if ffe_only %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + 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); + + disp(['FFE: ',sprintf('%.1E',ber),'| SIR: ',num2str(sir),' dB | PD_in: ',num2str(pd_in),' dBm']); + + 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)); + hold on + histogram(received(lvl,:),1000,"EdgeAlpha",0); + end + end + + elseif postfilter_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig] = Eq.process(Scpe_sig,Symbols); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + Noi = EQ_sig-Symbols; + + 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); + + nc = 2; + burg_coeff = arburg(Noi.signal,nc); + + EQ_sig = EQ_sig.filter(burg_coeff,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 + + 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)); + hold on + histogram(received(lvl,:),1000,"EdgeAlpha",0); + end + 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']); + + elseif db_channel_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Duobinary().encode(Symbols)); + + EQ_sig.plot("fignum",50,"displayname",'After EQ','clear',1); + + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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([' DB Precode -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | PD_in: ',num2str(pd_in),' dBm']); + + elseif db_coding_approach %%%%%%%%%%%%%%%%%%%%%%%%%%% + + [EQ_sig, Noi] = Eq.process(Scpe_sig,Symbols); + EQ_sig = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels).process(EQ_sig); + EQ_sig = Duobinary().decode(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); + + % EQ_sig.plot("fignum",50,"displayname",'After EQ'); + + disp([' DB Precode -> DB Code -> Channel -> FFE -> Decode/ Mod ',sprintf('%.1E',ber),' | PD_in: ',num2str(pd_in),' dBm']); + + end + + end + + 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); + wh.addValueToStorage(sir,'sir',v_bias,awg_vpp,eq_mode,i_atten); + wh.addValueToStorage(s_pow,'s_pow',v_bias,awg_vpp,eq_mode,i_atten); + wh.addValueToStorage(i_pow,'i_pow',v_bias,awg_vpp,eq_mode,i_atten); + 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 +end + +close(hWaitbar); + +wh.save([folderpath,experiment_name,'wh']); + + +cols = linspecer(8); + +i_atten_vals = wh.parameter.i_atten.values; +v_bias = wh.parameter.vbias.values(1); +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); + +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" + + +% Continue with the rest of your plot settings +yline(3.8e-3, 'DisplayName', 'HD-FEC', 'LineStyle', '--', 'HandleVisibility', 'off'); +xlabel('Signal to Interference Ratio in dB'); +ylabel('Bit Error Rate (BER)'); +title('Bit Error Rate vs. SIR'); +set(gca, 'yscale', 'log'); +set(gca, 'Box', 'on'); +grid on; +grid minor; +legend('Interpreter', 'none'); + + +autoArrangeFigures(3,3,2) + +disp("measurement done") \ No newline at end of file diff --git a/projects/Lab_2024/sweep_laser_vs_power.m b/projects/Lab_2024/sweep_laser_vs_power.m new file mode 100644 index 0000000..60df6e1 --- /dev/null +++ b/projects/Lab_2024/sweep_laser_vs_power.m @@ -0,0 +1,88 @@ + +% 1) Establish connection to laser +o = serialport("COM9",9600); %per USB angeschlossen +configureTerminator(o,"CR") +writeline(o,"*IDN?"); +wait(1) +if o.NumBytesAvailable ~= 0 + disp(['Laser Mainframe: ',readline(o)]); +else + error('Keine Verbindung zum Mainframe mglich?') + clear o + +end + +% 2) connect to device +v = visa('keysight', 'TCPIP0::134.245.243.248::inst0::INSTR'); +fopen(v); +fprintf(v, '*IDN?;'); +disp(['Powermeter: ' fscanf(v)]); + +% Define channels +laser_channel = 7; +powermeter_slot = 1; +l = 1300:0.1:1320; + +% turn on the laser +writeline(o,['CH',num2str(laser_channel),':ENABLE']); +wait(0.2) +current_wavelen = readline(o); + +clear power +clear lambda + +for n = 1:length(l) + + % 2 change wavelength in laser slot + command = string(['CH',num2str(laser_channel),':L=',num2str(l(n))]); + writeline(o,command); + wait(0.5); + readline(o); + + % query and check wavelength + writeline(o,['CH',num2str(laser_channel),':L?']); + wait(0.2) + current_wavelen = readline(o); + current_wavelen = str2double(strrep(regexp(current_wavelen,'([CH7:L=])+([\d]*)+([.])+([\d]*)','match'),'CH7:L=','')); + + if l(n) ~= current_wavelen + clear o + fclose(v); + delete(v); + clear v + error('Wellenlnge wurde nicht bernommen'); + end + + wait(0.75); + + % get current power in slot + slot = 1; + fprintf(v, [':READ' num2str(slot) ':POW?']); + power(n) = sscanf(fscanf(v),'%f'); + lambda(n) = current_wavelen; + + if mod(n,10)==1 + disp(['Measured ',num2str(power(n)),' dBm at ',num2str(lambda(n)),' nm']) + end +end + +figure(2) +hold on +plot(lambda,power,'Marker','*'); +xlabel('Wavelngth in nm') +ylabel('Power in dBm') +grid minor + +figure(211) +hold on +plot(lambda,10.^(power/10),'Marker','*'); +xlabel('Wavelngth in nm') +ylabel('Power in mW') +grid minor + +%close the serial connection + +clear o +fclose(v); +delete(v); +clear v