diff --git a/.gitignore b/.gitignore index 43b6a18..de1fa6c 100644 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,5 @@ sccprj/ # Matlab code generation folders codegen/ + +.mat diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 97ccbac..edccd25 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -359,10 +359,10 @@ classdef Signal end if options.normalizeToNyquist == 0 - [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,obj.fs,"centered","psd","mean"); + [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,obj.fs,"centered","power","mean"); w = w.*1e-9; else - [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,"centered","psd","mean"); + [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,"centered","power","mean"); end if options.normalizeTo0dB @@ -411,8 +411,9 @@ classdef Signal try ylim([max(min(floor(min(p_dbm))-3, ax.YLim(1)),-40), min(max(ceil(max(p_dbm))+3, ax.YLim(2)),10)]); catch - ylim([min(floor(min(p_dbm))-3, ax.YLim(1)), max(ceil(max(p_dbm))+3, ax.YLim(2))]); + ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]); end + ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]); yticks(-200:10:10); grid on; grid minor; legend @@ -655,13 +656,22 @@ classdef Signal end %% - function [obj,S,isFlipped] = tsynch(obj,options) + function [obj,S,isFlipped,sequenceFound] = tsynch(obj,options) % time sync and cut arguments obj Signal options.reference Signal options.fs_ref = 0; + options.debug_plots = 0; end + + + + S = {}; + isFlipped=0; + sequenceFound = 0; + + %normalize the signal a = obj.normalize("mode","oneone").signal; @@ -679,29 +689,46 @@ classdef Signal %estimate start pos of signal maxpeaknum = floor(length(a)/length(b)); - [pks,pkpos] = findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend'); + try + [pks,pkpos,w,p] = findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend'); + catch + warning(['Error in findpeaks, ususally the seuqnece is too short. Max peak num: ', num2str(maxpeaknum)]); + return + end + + if mean(w) > 10 || mean(p) > 10 + return + else + sequenceFound = 1; + end + + if options.debug_plots + figure() + findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend') + end + + + shifts = lags(pkpos); shifts = shifts(shifts>=0); - S = {}; - isFlipped=0; - if numel(shifts) > 0 %Cut occurences of ref signal from signal (only positive shifts) + if all(sign(co(pkpos))) + isFlipped = 1; + end for c = shifts(shifts>=0) sig = obj.delay(-c,'mode','samples'); - sig.signal = sig.signal(1:length(b)); + sig.signal = sig.signal(1:length(b)).*-isFlipped; S{end+1,1} = sig; end % - if all(sign(co(pkpos))) - isFlipped = 1; - end + %return/keep the sinal with the highest correlation (only within positive shifts) [~,idx]=max(pks(shifts>=0)); @@ -722,8 +749,7 @@ classdef Signal end %plot all synced signals and the ref signal - debug = 0; - if debug + if options.debug_plots figure;hold on; for i = 1:size(S,1) plot(S{i}.normalize('mode','oneone').signal(1000:1100),'LineWidth',0.1,'Color',[0.2157 0.4941 0.7216]); @@ -735,6 +761,12 @@ classdef Signal end + + + + + + %% function obj = filter(obj,a,b) @@ -850,7 +882,10 @@ classdef Signal sig = obj.signal; end - x = (sig); %% make input signal rea)l + startpos = floor(0.1*length(sig)); + endpos = floor(0.9*length(sig)); + endpos = min(endpos,startpos+200000); + x = sig(startpos:endpos); %% make input signal rea)l x = resample(x,fsym*histpoints_horizontal/2,obj.fs); %% up sample to original fsym rate diff --git a/Classes/01_transmit/ChannelFreqResp.m b/Classes/01_transmit/ChannelFreqResp.m index 9a1be19..f896f16 100644 --- a/Classes/01_transmit/ChannelFreqResp.m +++ b/Classes/01_transmit/ChannelFreqResp.m @@ -241,10 +241,13 @@ classdef ChannelFreqResp < handle xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on; %%% plot for publication - figure(1234);hold all;box on;title('Magnitude Freq. Response'); - %xlim([0.2 .5*max(obj.faxis)*1e-9]); - %ylim([-40, 2]); - Havg_smooth = smooth(Havg,50); + figure(1234); + hold all; + box on; + title('Magnitude Freq. Response'); + xlim([0.2 .5*max(obj.faxis)*1e-9]); + ylim([-40, 2]); + % Havg_smooth = smooth(Havg,10); symaxis = (obj.faxis-(obj.f_ref/2))/1e9; Havg = fftshift(Havg); %Havg = smooth(Havg); diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index 00ed0ca..6ba8558 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -8,28 +8,38 @@ classdef PAMmapper thresholds levels scaling + eth_style end methods - function obj = PAMmapper(M, unipolar) + function obj = PAMmapper(M, unipolar, options) %PAMMAPPER Construct an instance of this class % Detailed explanation goes here + arguments + M + unipolar + options.eth_style = 0; + end + obj.M = M; + obj.unipolar = unipolar; - + obj.thresholds = obj.get_demodulation_thresholds(); obj.levels = obj.get_levels(); obj.scaling = rms(obj.get_levels()); + obj.eth_style = options.eth_style; + end function out = map(obj,signal_in) - + if isa(signal_in,'Signal') - signal_in.signal = obj.map_(signal_in.signal); - + signal_in.signal = obj.map_(signal_in.signal); + % signal_in = signal_in.normalize("mode","rms"); lbdesc = ['Map bat stream to PAM ',num2str(obj.M),' symbols']; signal_in = signal_in.logbookentry(lbdesc,obj); @@ -37,19 +47,19 @@ classdef PAMmapper else out = signal_in; end - + end function signal_out = demap(obj,signal_in) issignalclass = 0; if isa(signal_in,'Signal') - signalclass = signal_in; - signal_in = signal_in.signal; - issignalclass = 1; + signalclass = signal_in; + signal_in = signal_in.signal; + issignalclass = 1; end - signal_out = obj.demap_(signal_in); + signal_out = obj.demap_(signal_in); if issignalclass lbdesc = ['Demap PAM ',num2str(obj.M),' symbols to bit stream']; @@ -58,7 +68,7 @@ classdef PAMmapper signal_in.signal = signal_out; signal_out = signal_in; end - + end function pam_sig = map_(obj,bitpattern) @@ -66,50 +76,114 @@ classdef PAMmapper switch obj.M case 2 % 2-ASK: BPSK / OOK - pam_sig=bitpattern(:,1); - - if obj.unipolar==0 - pam_sig=2*pam_sig-1; + if ~obj.eth_style + pam_sig = bitpattern(:,1); + if obj.unipolar==0 + pam_sig=2*pam_sig-1; + end + else + pam_sig = -2*bitpattern(:,1) + 1; end case 4 % 4-ASK: - pam_sig=2*bitpattern(:,1)+(bitpattern(:,1)==bitpattern(:,2)); + if ~obj.eth_style + pam_sig=2*bitpattern(:,1)+(bitpattern(:,1)==bitpattern(:,2)); - if obj.unipolar==0 - pam_sig=2*pam_sig-3; + if obj.unipolar==0 + pam_sig=2*pam_sig-3; + end + else + pam_sig = (2*bitpattern(:,1)-1).*(-2*bitpattern(:,2)+3); end pam_sig = pam_sig/sqrt(5); case 6 - - m = 1; - if size(bitpattern,2)>size(bitpattern,1) - bitpattern = bitpattern'; %vector aufrecht stellen - end - % LUT based mapping - for k = 1:5:fix(length(bitpattern)/5)*5 - pam_sig(m:m+1,1) = obj.thresholds(bin2dec(int2str(bitpattern(k:k+4)'))+1,:); - m = m+2; + + if ~obj.eth_style + + m = 1; + if size(bitpattern,2)>size(bitpattern,1) + bitpattern = bitpattern'; %vector aufrecht stellen + end + + % LUT based mapping + for k = 1:5:fix(length(bitpattern)/5)*5 + pam_sig(m:m+1,1) = obj.thresholds(bin2dec(int2str(bitpattern(k:k+4)'))+1,:); + m = m+2; + end + + else + + bitsPerSymbol = reshape(bitpattern,5,[]).'; % reorder 5 bits per symbol + normFactor = 1; + %====================32 QAM===================% + %=============================================% + % Coding % + % 01000 01001 |11001 11000 % + % | % + % 01010 01110 01100 |11100 11110 11010 % + % | % + % 01011 01111 01101 |11101 11111 11011 % + % --------------------|------------------- % + % 00011 00111 00101 |10101 10111 10011 % + % | % + % 00010 00110 00100 |10100 10110 10010 % + % | % + % 00000 00001 |10001 10000 % + %=============================================% + % modulate three LSB first in first Quadrant + % first bit inverts real part if 0 + % second bit inverts imaginary part if 0 + + LSB_symbols = normFactor*(... + +(1+1i) *(bitsPerSymbol(:,3)==1 & bitsPerSymbol(:,4)==0 & bitsPerSymbol(:,5)==1)... + +(3+1i) *(bitsPerSymbol(:,3)==1 & bitsPerSymbol(:,4)==1 & bitsPerSymbol(:,5)==1)... + +(5+1i) *(bitsPerSymbol(:,3)==0 & bitsPerSymbol(:,4)==1 & bitsPerSymbol(:,5)==1)... + +(1+3i) *(bitsPerSymbol(:,3)==1 & bitsPerSymbol(:,4)==0 & bitsPerSymbol(:,5)==0)... + +(3+3i) *(bitsPerSymbol(:,3)==1 & bitsPerSymbol(:,4)==1 & bitsPerSymbol(:,5)==0)... + +(5+3i) *(bitsPerSymbol(:,3)==0 & bitsPerSymbol(:,4)==1 & bitsPerSymbol(:,5)==0)... + +(1+5i) *(bitsPerSymbol(:,3)==0 & bitsPerSymbol(:,4)==0 & bitsPerSymbol(:,5)==1)... + +(3+5i) *(bitsPerSymbol(:,3)==0 & bitsPerSymbol(:,4)==0 & bitsPerSymbol(:,5)==0)); + + Re = real(LSB_symbols); + Im = imag(LSB_symbols); + % if first bit== 0 => invert real part + % if second bit== 0 => invert imag part + modData2 = Re.*(bitsPerSymbol(:,1)*2-1) + 1i*Im.*(bitsPerSymbol(:,2)*2-1); + Re = real(modData2(1:end/2)); + Im = imag(modData2(1:end/2)); + pam_sig = zeros(length(Re)*2,1); + pam_sig(1:2:length(Re)*2) = Re; + pam_sig(2:2:length(Im)*2) = Im; + + end pam_sig = pam_sig/sqrt(10); case 8 + % 8-ASK: - x1 = bitpattern(:,1); - x2 = (bitpattern(:,1)==bitpattern(:,3)); - x3 = x2~=bitpattern(:,2); + if ~obj.eth_style + x1 = bitpattern(:,1); + x2 = (bitpattern(:,1)==bitpattern(:,3)); + x3 = x2~=bitpattern(:,2); - pam_sig = 4*x1 + 2*x2 + x3; + pam_sig = 4*x1 + 2*x2 + x3; + + if obj.unipolar==0 + pam_sig=2*pam_sig-7; + end + else + + pam_sig = (bitpattern(:,1)*2-1).*(4+(2*bitpattern(:,2)-1).*(-2*bitpattern(:,3)+3)); - if obj.unipolar==0 - pam_sig=2*pam_sig-7; end pam_sig = pam_sig/sqrt(21); - + case 16 % 16-ASK: x1 = bitpattern(:,1); @@ -143,7 +217,7 @@ classdef PAMmapper % case 16 % thres = [-10 -8 -6 -4 -2 0 2 4 6 8 10]; % end - + switch obj.M @@ -166,12 +240,12 @@ classdef PAMmapper thres = thres .* 1/sqrt(5); - case 6 %PAM 6 + case 6 %PAM 6 thres = [-3 5;-1 5;-3 -5;-1 -5;-5 3;-5 1;-5 -3;-5 -1;-1 3;-1 1;-1 -3;-1 -1;-3 3;-3 1;-3 -3;-3 -1;3 5;1 5;3 -5;1 -5;5 3;5 1;5 -3;5 -1;1 3;1 1;1 -3;1 -1;3 3;3 1;3 -3;3 -1]; % thres = [-4 -2 0 2 4]; % thres = thres ./ sqrt(10); - + case 8 % 8-ASK if obj.unipolar==0 @@ -179,7 +253,7 @@ classdef PAMmapper elseif obj.unipolar==1 thres=0.5:6.5; end - + thres=thres./sqrt(21); case 16 @@ -210,65 +284,115 @@ classdef PAMmapper end function [data_out] = demap_(obj,data_in) + data_in= data_in'; + if obj.M ~= 6 - + % create output if ~isempty(obj.thresholds) a = squeeze(repmat(real(data_in),[1 1 length(obj.thresholds)])); %Eingangssignal in 3 spalten b = squeeze(repmat(reshape(obj.thresholds(:).',[1 1 length(obj.thresholds)]),[1 length(data_in) 1])); %Threshold in 3 Spalten comp_real = a > b; %check for each symbol/ sampling if it exeeds the obj.thresholdseshold 1, 2 or 3 - + comp_real=repmat(real(data_in),[1 1 length(obj.thresholds)]) > repmat(reshape(obj.thresholds(:).',[1 1 length(obj.thresholds)]),[1 length(data_in) 1]); - + else comp_real=[]; end - + s1=size(comp_real,1); s2=size(comp_real,2); + end switch obj.M case 2 % 2-ASK - - data_out=comp_real(:,:,1); + if ~obj.eth_style + data_out=comp_real(:,:,1); + else + data_out=abs(comp_real(:,:,1)-1); + end case 4 % 4-ASK - - data_out=[comp_real(:,:,2); ones(s1,s2) - comp_real(:,:,1) + comp_real(:,:,3)]; + if ~obj.eth_style + data_out=[comp_real(:,:,2); ones(s1,s2) - comp_real(:,:,1) + comp_real(:,:,3)]; + else + data_out= [(data_in>=0); (abs(data_in)<=1)]; + end case 6 - data_in = data_in/(sqrt(mean(abs(data_in).^2))); - data_in = data_in*sqrt(10); + if ~obj.eth_style + + data_in = data_in/(sqrt(mean(abs(data_in).^2))); + data_in = data_in*sqrt(10); + + if size(data_in,2) > 1 + data_in = data_in.'; + end + + if length(data_in)/2 ~= round(length(data_in)/2) + data_in = [data_in;0]; + end + + m = 1; + for n = 1:2:length(data_in) + dist = sqrt((data_in(n)-obj.thresholds(:,1)).^2+(data_in(n+1)-obj.thresholds(:,2)).^2); + [~,dd_idx] = min(dist); + % dec_out(n:n+1) = LUT(dd_idx,:); + data_out(m:m+4) = bitget(dd_idx-1,5:-1:1); + m = m+5; + end + + else + + data_in= data_in'; + rxSym = data_in(1:2:end) + 1i*data_in(2:2:end); % 16×1 + + % Decode the sign bits (bits 1 & 2). + rxBit1 = double(real(rxSym) > 0); + rxBit2 = double(imag(rxSym) > 0); + + % Undo the quadrant inversion and normalization. + rxSym_corr = abs(real(rxSym)) + 1i*abs(imag(rxSym)); + normFactor = 1/sqrt(10); + rxSym_unscaled = rxSym_corr / normFactor; + + + cand = [1+1i, 3+1i, 5+1i, 1+3i, 3+3i, 5+3i, 1+5i, 3+5i]; + candBits = [1 0 1; + 1 1 1; + 0 1 1; + 1 0 0; + 1 1 0; + 0 1 0; + 0 0 1; + 0 0 0]; + d = abs(rxSym_unscaled - cand).^2; % 32×8 distances + [~, idx] = min(d, [], 2); + rxLSB = candBits(idx,:); % 32×3 + + + decodedSymbols = [rxBit1, rxBit2, rxLSB]; + data_out = reshape(decodedSymbols.', [], 1).'; - if size(data_in,2) > 1 - data_in = data_in.'; - end - - if length(data_in)/2 ~= round(length(data_in)/2) - data_in = [data_in;0]; - end - - m = 1; - for n = 1:2:length(data_in) - dist = sqrt((data_in(n)-obj.thresholds(:,1)).^2+(data_in(n+1)-obj.thresholds(:,2)).^2); - [~,dd_idx] = min(dist); - % dec_out(n:n+1) = LUT(dd_idx,:); - data_out(m:m+4) = bitget(dd_idx-1,5:-1:1); - m = m+5; end + case 8 % 8-ASK - data_out=[comp_real(:,:,4); - comp_real(:,:,1)-comp_real(:,:,3)+comp_real(:,:,5)-comp_real(:,:,7); - 1-comp_real(:,:,2)+comp_real(:,:,6)]; + if ~obj.eth_style + data_out=[comp_real(:,:,4); + comp_real(:,:,1)-comp_real(:,:,3)+comp_real(:,:,5)-comp_real(:,:,7); + 1-comp_real(:,:,2)+comp_real(:,:,6)]; + else + data_out = [(data_in>=0); (abs(data_in)>(4/sqrt(21))); (abs(data_in)>=(2/sqrt(21)))&(abs(data_in)<=(6/sqrt(21)))]; + end case 16 % 16-ASK @@ -289,21 +413,21 @@ classdef PAMmapper data_in Signal options.symbol_levels = [] end - + %A) normally return the preproduct of the decision a = squeeze(repmat(real(data_in.signal),[1 1 length(obj.thresholds)])); %Eingangssignal in 3 spalten b = squeeze(repmat(reshape(obj.thresholds(:).',[1 1 length(obj.thresholds)]),[1 length(data_in.signal) 1])); %Threshold in 3 Spalten comp_real = a > b; %check for each symbol/ sampling if it exeeds the obj.thresholdseshold 1, 2 or 3 - + data_out = data_in; data_out.signal = sum(comp_real,2); - + %Option: return the actual level values/ just map onto given %symbol levels if ~isempty(options.symbol_levels) data_out.signal = options.symbol_levels(data_out.signal+1); end - + end function [out] = separate_pamlevels(obj,data_in) @@ -312,7 +436,7 @@ classdef PAMmapper a = squeeze(repmat(real(data_in.signal),[1 1 length(obj.thresholds)])); %Eingangssignal in 3 spalten b = squeeze(repmat(reshape(obj.thresholds(:).',[1 1 length(obj.thresholds)]),[1 length(data_in.signal) 1])); %Threshold in 3 Spalten comp_real = a > b; %check for each symbol/ sampling if it exeeds the obj.thresholdseshold 1, 2 or 3 - comp_real_sum = sum(comp_real,2); + comp_real_sum = sum(comp_real,2); out = NaN(length(data_in),length(obj.thresholds)+1); @@ -323,7 +447,7 @@ classdef PAMmapper end function [Signal_out] = quantize(obj,Signal_in) - + constellation = obj.get_levels(); constellation = constellation ./ rms(constellation); @@ -331,9 +455,9 @@ classdef PAMmapper if isa(Signal_in,'Signal') issignalclass = 1; Sig_class = Signal_in; - Signal_in = Signal_in.signal; + Signal_in = Signal_in.signal; end - + [~,high_dim_sig] = max(size(Signal_in)); [~,high_dim_const] = max(size(constellation)); @@ -346,7 +470,7 @@ classdef PAMmapper Signal_out = constellation(symbol_idx); Signal_out = reshape(Signal_out,size(Signal_in)); - if issignalclass + if issignalclass Sig_class.signal = Signal_out; Signal_out = Sig_class; end @@ -355,6 +479,10 @@ classdef PAMmapper end + function bitmap = showBitMapping(obj) + bitmap = obj.demap([obj.levels ./ obj.scaling]'); + end + end end diff --git a/Classes/01_transmit/PAMsource.m b/Classes/01_transmit/PAMsource.m index eb94a04..ce15418 100644 --- a/Classes/01_transmit/PAMsource.m +++ b/Classes/01_transmit/PAMsource.m @@ -71,18 +71,17 @@ classdef PAMsource function [digi_sig,symbols,bits] = process(obj) %%%%% PRBS Generation in correct shape for Modulation Format %%%%%% - O = obj.order; %order of prbs - N = 2^(O-1); %length of prbs - [~,seed] = prbs(O,1); %initialize first seed of prbs bitpattern=[]; if obj.useprbs - % for i = 1:log2(obj.M) - % [bitpattern(:,i),seed] = prbs(O,N,seed); - % end + % O = obj.order; %order of prbs + % N = 2^(O-1); %length of prbs + % [~,seed] = prbs(O,1); %initialize first seed of prbs + % % for i = 1:log2(obj.M) + % % [bitpattern(:,i),seed] = prbs(O,N,seed); + % % end %%%%% MOVE-IT PRMS %%%% - state = struct(); para = struct(); @@ -103,6 +102,8 @@ classdef PAMsource para.reset_prms = 0; para.method = 1; + + data_in = []; global loop; loop = 0; @@ -116,6 +117,7 @@ classdef PAMsource else s = RandStream('twister','Seed',obj.randkey); for i = 1:log2(obj.M) + N = 2^(obj.order-1); %length of prbs bitpattern(:,i) = randi(s,[0 1], N, 1); end end @@ -151,23 +153,52 @@ classdef PAMsource sym_max = max(symbols.signal); end - % symbols.move_it_spectrum("fignum",222,"displayname","Symbols only"); - % symbols.spectrum("fignum",222,"displayname","Symbols only","normalizeTo0dB",1); - - %%%%% Pulse-forming %%%%%% if obj.applypulseform + %%% MY CODE digi_sig = obj.pulseformer.process(symbols); + + %%% MOVEIT WRAPPER BETA + % symbols.spectrum("fignum",111,"displayname","1) RAW SIGNAL"); + % + % pulsf = Moveit_wrapper('pulsef'); + % pulsf.para.alpharacos = obj.pulseformer.alpha; + % pulsf.para.f_sym = obj.fsym; + % pulsf.para.fs = obj.fsym; + % pulsf.para.pulse = 2; + % + % digi_sig = pulsf.process(symbols); + + % Design raised cosine filter with given order in symbols + % digi_sig = symbols; + % sps = 4; + % nsym = 128; + % rctFilt3 = comm.RaisedCosineTransmitFilter(... + % Shape='Square root', ... + % RolloffFactor=1, ... + % FilterSpanInSymbols=nsym, ... + % OutputSamplesPerSymbol=sps); + % + % digi_sig.signal = rctFilt3([symbols.signal; zeros(nsym/2,1)]); + % % Correct for propagation delay by removing filter transients + % fltDelay = nsym / (2*obj.fsym); + % digi_sig.signal = digi_sig.signal(fltDelay*sps*obj.fsym+1:end); + % + % digi_sig.fs = sps .* symbols.fs; + + else digi_sig = symbols; end - %%%%% Re-sample to f DAC %%%%%% - digi_sig = digi_sig.resample("fs_in",digi_sig.fs,"fs_out",obj.fs_out,"n",10,"beta",5); - % digi_sig.spectrum("fignum",111,"displayname","after pulseforming"); + %%%%% Re-sample to f DAC %%%%%% + n = 10; + digi_sig = digi_sig.resample("fs_in",digi_sig.fs,"fs_out",obj.fs_out,"n",n,"beta",5); + % digi_sig.spectrum("fignum",111,"displayname",['3) Shaped + Resampled; n: ',num2str(n)]); + %%%%% Hard clip digital signal to PAM range before DAC %%%%%% if obj.applyclipping diff --git a/Classes/01_transmit/Pulseformer.m b/Classes/01_transmit/Pulseformer.m index c8e8537..cbf62f3 100644 --- a/Classes/01_transmit/Pulseformer.m +++ b/Classes/01_transmit/Pulseformer.m @@ -3,13 +3,17 @@ classdef Pulseformer % Detailed explanation goes here properties(Access=public) - fdac + end - properties(Access=private) + properties(Access=public) + fs fsym + matched_sps + output_sps pulse pulselength - rrcalpha + alpha + matched end methods (Access=public) @@ -18,11 +22,14 @@ classdef Pulseformer % Detailed explanation goes here arguments - options.fdac double + options.fs double options.fsym double - options.pulse pulseform = pulseform.rrc + options.matched_sps double + options.output_sps double + options.pulse pulseform = pulseform.rc options.pulselength double {mustBeInteger} = 32 - options.rrcalpha double = 0.05 + options.alpha double = 0.05 + options.matched = 0; end % @@ -47,7 +54,11 @@ classdef Pulseformer signalclass_in = signalclass_in.logbookentry(lbdesc); % write fs to signal - signalclass_in.fs = obj.fdac; + if obj.matched + signalclass_in.fs = obj.fsym .* obj.output_sps;%.* (obj.fdac./obj.fsym); + else + signalclass_in.fs = obj.fs;%.* (obj.fdac./obj.fsym); + end % write to output signalclass_out = signalclass_in; @@ -74,60 +85,99 @@ classdef Pulseformer data_out end - if ~rem(obj.fdac,obj.fsym) + if ~isempty(obj.output_sps) + obj.fsym = obj.fsym.*obj.output_sps; + end + + if ~rem(obj.fs,obj.fsym) %ist ein Vielfaches - sps = obj.fdac / obj.fsym; - up = sps; - dn = 1; + sps = obj.fs / obj.fsym; + p = sps; + q = 1; else %ist kein Vielfaches - up = obj.fdac / gcd(obj.fdac, obj.fsym); - dn = obj.fsym / gcd(obj.fdac, obj.fsym); - sps= up; + p = obj.fsym / gcd(obj.fs, obj.fsym); %upsampling p->->-> + q = obj.fs/ gcd(obj.fs, obj.fsym); %downsampling <-q + sps= q; %sps während dem pulse shaping end - if obj.pulse == pulseform.rrc - %Bau das Filter (hier rrc) - racos_len = obj.pulselength*2; - alpha = obj.rrcalpha; - h = rcosdesign(alpha,racos_len,sps,"normal"); - h = h./ max(h); + if obj.pulse == pulseform.rc + filtertype = 'normal'; + elseif pulseform.rrc + filtertype = 'sqrt'; end - % Apply filter the long way (from move_it) - % block length in samples - - data_in = data_in'; - blen = length(data_in)*sps; - - % oversample symbol sequence - - symbolov=zeros(size(data_in,1),blen); - symbolov(:,1:sps:blen-sps+1)=data_in; + %Bau das Filter (hier rc) + racos_len = obj.pulselength*2; + h = rcosdesign(obj.alpha,racos_len,sps,filtertype); + % h = h./ max(h); - H=fft(h,blen); - - % Convolution of Bit sequence with impulse response + if obj.matched + h = conj(fliplr(h)); + end + + manual_cyclic_convolution = 0; + upsample_filter = 0; + upfirdn_convolution = 1; - % cyclic convolution - data_out=ifft( fft(symbolov.') .* repmat( H,size(data_in,1),1 ).' ).'; + if manual_cyclic_convolution - data_out = circshift(data_out,[0 -(obj.pulselength*sps)]); + % Apply filter the long way (from move_it) + data_in = data_in'; + blen = length(data_in)*sps; + + % oversample symbol sequence + symbolov=zeros(size(data_in,1),blen); + symbolov(:,1:sps:blen-sps+1)=data_in; + H=fft(h,blen); + + % Convolution of Bit sequence with impulse response + data_out=ifft( fft(symbolov.') .* repmat( H,size(data_in,1),1 ).' ).'; + data_out = circshift(data_out,[0 -(obj.pulselength*sps)]); + + if rem(obj.fs,obj.fsym) + data_out = data_out(1:q:end); + end - if rem(obj.fdac,obj.fsym) - data_out = data_out(1:dn:end); %! end + if upsample_filter + data_out_ = upsample(data_in,p); -% %Apply Filter using Matlab build in fctn. -% h = rcosdesign(alpha,racos_len,sps); -% h = h./ max(h); - %data_out_ = upfirdn(data_in,h,up,dn); -% -% %cut signal, which is longer due to fir filter - % st = round(up/dn*racos_len/2); %we need to cut y_out - % en = round(st + (length(data_in)*up/dn) -1); - % data_out = data_out(st:en); + mfOutput = filter(h, 1, data_out_); % Matched filter output + + figure() + hold on + stem(mfOutput(1:1000),'Marker','o','MarkerSize',1,'LineStyle','-','LineWidth',1); + stem(data_out_(1:1000),'Marker','o','MarkerSize',1,'LineStyle','-','LineWidth',1); + + + end + + if upfirdn_convolution + + %Apply Filter using Matlab build in fctn. + data_out_ = upfirdn(data_in,h,p,q); + + %cut signal, which is longer due to fir filter + st = round(p/q*racos_len/2); %we need to cut y_out + en = round(st + (length(data_in)*p/q)); + data_out = data_out_(st:en); + + end + + if upfirdn_convolution && manual_cyclic_convolution + figure() + subplot(2,1,1) + title("Convolution vs. Upfirdn and Cut") + hold on + % plot(data_out_(1:200),'DisplayName','Matlab upfirdn'); + plot(data_out(1:200),'DisplayName','By Hand cyclic convolution') + subplot(2,1,2) + hold on + plot(data_out(1:2000),'DisplayName','OUT'); + plot(data_in(1:2000),'DisplayName','IN'); + end %scaling?! see pulsef module line 696 % scale = max(max([abs(real(data_out)) abs(imag(data_out))])); %find max value from real and imag part @@ -136,7 +186,7 @@ classdef Pulseformer data_out = data_out'; %Check output integrity - if abs(round(up/dn * length(data_in)) - length(data_out)) > 4 + if abs(round(p/q * length(data_in)) - length(data_out)) > 4 warning('Check signal length after pulse shaping'); %disp('Check signal length after pulse shaping'); end diff --git a/Classes/04_DSP/Coding/Duobinary.m b/Classes/04_DSP/Coding/Duobinary.m index e93c594..6f3460b 100644 --- a/Classes/04_DSP/Coding/Duobinary.m +++ b/Classes/04_DSP/Coding/Duobinary.m @@ -26,7 +26,9 @@ classdef Duobinary M = numel(u); %make unipolar - if M == 4 + if M == 2 + data = data; + elseif M == 4 data = data .* sqrt(5); elseif M == 6 data = data .* sqrt(10); @@ -59,7 +61,9 @@ classdef Duobinary bk = bk .* 2; bk = bk + b; - if M == 4 + if M == 2 + bk = bk; + elseif M == 4 bk = bk ./ sqrt(5); elseif M == 6 bk = bk ./ sqrt(10); @@ -100,7 +104,9 @@ classdef Duobinary end %make unipolar - if options.M == 4 + if options.M == 2 + data = data; + elseif options.M == 4 data = data .* sqrt(5); elseif options.M == 6 data = data .* sqrt(10); @@ -140,7 +146,9 @@ classdef Duobinary mean_power = sum((unique_points .^ 2) .* probabilities); scaling_factor = sqrt(mean_power); - if options.M == 4 + if options.M == 2 + data = data ./ sqrt(0.5); + elseif options.M == 4 data = data ./ sqrt(2.5); % 7-level constellation weighted with probability after DB code i.e. mean([-3 3 -2 -2 2 2 -1 -1 -1 1 1 1 0 0 0 0].^2) = 2.5 --> sqrt(2.5) == rms(constellation) elseif options.M == 6 data = data ./ sqrt(5.8); @@ -186,7 +194,9 @@ classdef Duobinary end %make unipolar - if I == 7 || I == 6 + if I == 3 + data = data .* 0.5; + elseif I == 7 || I == 6 data = data .* sqrt(2.5); elseif I == 11 data = data .* sqrt(5.8); @@ -206,7 +216,9 @@ classdef Duobinary data = data .* 2; data = data - round(mean(data)); - if M == 4 + if M == 2 + data = data; + elseif M == 4 data = data ./ sqrt(5); elseif M == 6 data = data ./ sqrt(10); diff --git a/Classes/04_DSP/Equalizer/FFE.m b/Classes/04_DSP/Equalizer/FFE.m index 2f414ec..54d0b57 100644 --- a/Classes/04_DSP/Equalizer/FFE.m +++ b/Classes/04_DSP/Equalizer/FFE.m @@ -101,6 +101,17 @@ classdef FFE < handle end x = [zeros(floor(obj.order/2),1); x; zeros(obj.order,1)]; + + if training + mask = ones(obj.order,1); + else + mask = zeros(obj.order,1); + mask(900:end) = 1; + mask(ceil(length(obj.e)/2)) = 1; + end + + mask = ones(obj.order,1); + for epoch = 1 : epochs symbol = 0; @@ -110,7 +121,7 @@ classdef FFE < handle U = x(obj.order+sample-1:-1:sample); - y(symbol,1) = obj.e.' * U; % Calculating output of LMS __ * | + y(symbol,1) = (obj.e.*mask).' * U; % Calculating output of LMS __ * | if training d_hat(symbol,1) = d(symbol); @@ -129,7 +140,6 @@ classdef FFE < 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 diff --git a/Classes/04_DSP/Equalizer/Postfilter.m b/Classes/04_DSP/Equalizer/Postfilter.m index 1fed009..8d09ef1 100644 --- a/Classes/04_DSP/Equalizer/Postfilter.m +++ b/Classes/04_DSP/Equalizer/Postfilter.m @@ -43,12 +43,12 @@ classdef Postfilter < handle if ~isnan(options.useBurg) && options.useBurg - disp('using burg alg') + % disp('using burg alg') obj.coefficients = arburg(noiseclass_in.signal,obj.ncoeff); elseif ~isempty(options.coefficients) - disp('using given taps') + % disp('using given taps') obj.coefficients = options.coefficients; obj.useBurg = 0; @@ -67,22 +67,29 @@ classdef Postfilter < handle end - function showFilter(obj,noiseclass_in,options) + function showFilter(obj,options) arguments obj - noiseclass_in + options.noiseclass_in =[] options.fignum = 121 options.color = [] end % noiseclass_in.spectrum('displayname','Noise PSD shifted to 0dBm','fignum',options.fignum,'normalizeTo0dB',1); - + if ~isempty(options.noiseclass_in) + len = length(options.noiseclass_in); + fs = options.noiseclass_in.fs; + else + len = 1024; + fs = 1; + end + figure(options.fignum) - [h,w] = freqz(1,obj.coefficients,length(noiseclass_in),"whole",noiseclass_in.fs); + [h,w] = freqz(1,obj.coefficients,len,"whole",fs); h = h/max(abs(h)); hold on - w_ = (w - noiseclass_in.fs/2); + w_ = (w - fs/2); if isempty(options.color) plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(obj.coefficients,2)), ' '],'LineWidth',2); diff --git a/Classes/04_DSP/Sequence Detection/MLSE.m b/Classes/04_DSP/Sequence Detection/MLSE.m index 59b9c33..9576524 100644 --- a/Classes/04_DSP/Sequence Detection/MLSE.m +++ b/Classes/04_DSP/Sequence Detection/MLSE.m @@ -70,7 +70,7 @@ classdef MLSE < handle %%%% Separate the equalized signal into the respective levels based on the actually transmitted level constellation = unique(data_ref); decisionLevels = (constellation(1:end-1) + constellation(2:end)) / 2; - tx_bits = PAMmapper(numel(constellation),0).demap(data_ref); + tx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(data_ref); % impulse respnse i.e. [0.5, 1.0000] @@ -188,10 +188,50 @@ classdef MLSE < handle end + % Compute soft output PAM4 stream from the metric_sym + soft_output = zeros(length(data_in),1); % Expected symbol value per stage + symbol_prob = zeros(length(data_in), length(states)); % Store full probability distribution (optional) + llp = llp'; + for n = 1:length(data_in) + metrics = llp(n, :); % A posteriori metric for each PAM4 candidate + % For numerical stability, subtract the maximum metric before exponentiating + maxMetric = min(metrics); + expMetrics = exp(metrics - maxMetric); + probs = expMetrics / sum(expMetrics); % Normalize to get probabilities + symbol_prob(n, :) = probs; % (Optional) store distribution for analysis + % Compute the soft output as the expected value of the PAM4 symbols + soft_output(n) = sum(probs .* constellation'); + end + + % Number of symbols and bits per symbol + num_symbols = constellation; + num_bits = 2; % 2 bits per symbol + + bit_mapping = PAMmapper(4,0,"eth_style",1).showBitMapping; % Each row corresponds to the symbol above + + % Initialize LLR storage + llr = zeros(num_bits, length(data_in)); + + % Compute bit-wise LLRs + for bit_idx = 1:num_bits + % Find indices where bit is 0 and where it is 1 + idx_bit_0 = find(bit_mapping(:,bit_idx) == 0); + idx_bit_1 = find(bit_mapping(:,bit_idx) == 1); + + % Sum over log-probabilities (Max-Log approximation: using min instead of sum) + llr(:,bit_idx) = min(llp(:,idx_bit_1), [], 1) - min(llp(:,idx_bit_0), [], 1); + + end + + % Convert LLR values to a hard-decision bit stream + bit_stream = llr < 0; + [~,~,ber_llr,~] = calc_ber(bit_stream',tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('LLR BER : %.2e \n',ber_llr); + % directly decide based on lowest LLP index [~,llp_based_state_seq]=min(llp); LLP_EST(1:length(data_in)) = constellation(llp_based_state_seq); - rx_bits = PAMmapper(numel(constellation),0).demap(LLP_EST'); + rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(LLP_EST'); [~,~,ber_llp,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('LLP BER : %.2e \n',ber_llp); % [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); @@ -200,7 +240,7 @@ classdef MLSE < handle % fprintf('LLP BER -1: %.2e \n',ber_llp); %%%% DECIDE based on Viterbi traceback - rx_bits = PAMmapper(numel(constellation),0).demap(VITERBI_ESTIMATION_SYMBOLS'); + rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(VITERBI_ESTIMATION_SYMBOLS'); [~,~,ber_viterbi,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('Viterbi BER: %.2e \n',ber_viterbi); % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); @@ -209,14 +249,14 @@ classdef MLSE < handle % directly decide based on the FW path metrics [~,fw_direct_state_seq]=min(pm_survivor_fw); FW_EST(1:length(data_in)) = constellation(fw_direct_state_seq); - rx_bits = PAMmapper(numel(constellation),0).demap(FW_EST'); + rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(FW_EST'); [~,~,ber_fw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('FW BER: %.2e \n',ber_fw); % directly decide based on the BW path metrics [~,bw_direct_state_seq]=min(pm_survivor_bw); BW_EST(1:length(data_in)) = constellation(bw_direct_state_seq(2:end)); - rx_bits = PAMmapper(numel(constellation),0).demap(BW_EST'); + rx_bits = PAMmapper(numel(constellation),0,"eth_style",1).demap(BW_EST'); [~,~,ber_bw,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); fprintf('BW BER: %.2e \n',ber_bw); % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); @@ -224,6 +264,8 @@ classdef MLSE < handle % [~,~,ber_viterbi,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); % fprintf('BW BER: %.2e \n',ber_viterbi); + PAMmapper(4,0,"eth_style",1).showBitMapping + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tx_symbolpos = zeros(numel(constellation),length(data_ref)); diff --git a/Classes/04_DSP/TransmissionPerformance.m b/Classes/04_DSP/TransmissionPerformance.m index 49ca0c7..1550e40 100644 --- a/Classes/04_DSP/TransmissionPerformance.m +++ b/Classes/04_DSP/TransmissionPerformance.m @@ -59,8 +59,19 @@ classdef TransmissionPerformance 1.03e-2, 9.29e-3, 8.33e-3, 7.54e-3, 7.04e-3, 4.70e-3]; %% LUT for KP4-FEC and Inner Code https://grouper.ieee.org/groups/802/3/dj/public/23_03/patra_3dj_01b_2303.pdf + CODE_RATE_KP4_AND_INNER = [0.885799]; BERTHRESHOLDS_KP4_AND_INNER = 4.85e-3; + + % https://www.ieee802.org/3/bs/public/14_11/parthasarathy_3bs_01a_1114.pdf + CODE_RATE_KP4 = [1/(1+0.052)];%Reed-Solomon RS(544, 514) == KP4 + BERTHRESHOLDS_KP4 = 2.2e-4; + + CODE_RATE_HDFEC = [1/(1+0.067)]; %Beyond 300 Gbps Short-Reach Links Using TFLN MZMs With 500 mVpp and Linear Equalization + BERTHRESHOLDS_HDFEC = 3.8e-3; + + CODE_RATE_O_FEC = [1/(1+0.153)]; %Stefano im Meeting + BERTHRESHOLDS_O_FEC = 2e-2; @@ -151,6 +162,11 @@ classdef TransmissionPerformance netrates.KP4_hamming.NetRate = NaN(1, numMeasurements); netrates.KP4_hamming.CodeRate = NaN(1, numMeasurements); netrates.KP4_hamming.Threshold = NaN(1, numMeasurements); + + netrates.O_FEC.GrossRate = NaN(1, numMeasurements); + netrates.O_FEC.NetRate = NaN(1, numMeasurements); + netrates.O_FEC.CodeRate = NaN(1, numMeasurements); + netrates.O_FEC.Threshold = NaN(1, numMeasurements); end % Process each measurement individually. @@ -204,6 +220,22 @@ classdef TransmissionPerformance netrates.KP4_hamming.Threshold(i) = obj.BERTHRESHOLDS_KP4_AND_INNER(idxBER); end + idxBER = []; + for j = length(obj.BERTHRESHOLDS_O_FEC):-1:1 + if ber(i) <= obj.BERTHRESHOLDS_O_FEC(j) + idxBER = j; + break; + end + end + if ~isempty(idxBER) + codeRate = obj.CODE_RATE_O_FEC(idxBER); + netrates.O_FEC.NetRate(i) = grossRate(i) * codeRate; + netrates.O_FEC.GrossRate(i) = grossRate(i) ; + netrates.O_FEC.CodeRate(i) = codeRate; + netrates.O_FEC.Threshold(i) = obj.BERTHRESHOLDS_O_FEC(idxBER); + end + + end diff --git a/Classes/DataBaseHandler/DBHandler.m b/Classes/DataBaseHandler/DBHandler.m index 1495209..d3bc066 100644 --- a/Classes/DataBaseHandler/DBHandler.m +++ b/Classes/DataBaseHandler/DBHandler.m @@ -45,15 +45,15 @@ classdef DBHandler < handle else error('DB seems to be corrupt') end - + end function obj = refresh(obj) - % Get table names and the first rows of each table to understand the structure - obj.getTableNames(); - obj.getTables(); - obj.getDistinctValues(); + % Get table names and the first rows of each table to understand the structure + obj.getTableNames(); + obj.getTables(); + obj.getDistinctValues(); end function obj = getTableNames(obj) @@ -187,10 +187,11 @@ classdef DBHandler < handle if isstruct(newRow) fields = fieldnames(newRow); emptyFields = structfun(@isempty,newRow); - if sum(emptyFields)>0 - newRow.(fields{emptyFields==1}) = NaN; - disp(['In Table: ',tableName,': ',fields{emptyFields==1},' was empty, is now NaN ',newRow.(fields{emptyFields==1})]) - + if sum(emptyFields) > 0 + emptyFieldNames = fields(emptyFields); % use () to get a cell array + for idx = 1:numel(emptyFieldNames) + newRow.(emptyFieldNames{idx}) = NaN; + end end newRow = struct2table(newRow); end @@ -274,87 +275,168 @@ classdef DBHandler < handle end end - function addBEREntry(obj, berValue, occurrence, runID, ffe, dfe, mlse, pf, eqType, ffe_order, dfe_order, len_tr, mu_ffe, mu_dfe, mu_dc, comment) - % addBEREntry Adds a BER entry linked to an existing or new Equalizer entry. - % Usage: - % addBEREntry(runID, eq, ffe, dfe, mlse, pf, eqType, ffe_order, dfe_order, len_tr, mu_ffe, mu_dfe, mu_dc, berValue, comment) - if isempty(pf) - postfilter_taps = []; + function resultID = addProcessingResult(obj, run_id, resultData, eqParamsData) + % addProcessingResult Adds a processing result and links it to an EqualizerParameters entry. + % + % Inputs: + % run_id: A run_id from the main table to connect the BER with. + % + % resultData: A struct with fields corresponding to the ProcessingResults table. + % + % eqParamsData: A struct with fields corresponding to the EqualizerParameters table, + % except 'eq_id' and 'config_hash'. These fields are used to compute + % a hash and check for an existing configuration. + % + % Output: + % resultID: The result_id of the newly inserted ProcessingResults entry. + + % 1. Compute hash for equalizer parameters + jsonStr = jsonencode(eqParamsData); + md = java.security.MessageDigest.getInstance('MD5'); + md.update(uint8(jsonStr)); + hashBytes = typecast(md.digest, 'uint8'); + hashStr = lower(dec2hex(hashBytes)'); + hashStr = lower(strtrim(hashStr(:)')); % Convert to a lowercase string + + % Add hash to equalizer parameters + eqParamsData.config_hash = hashStr; + + % 2. Check if an equalizer configuration with the same hash exists + queryStr = sprintf('SELECT eq_id FROM EqualizerParameters WHERE config_hash = ''%s''', eqParamsData.config_hash); + existingEntry = obj.fetch(queryStr); + + if ~isempty(existingEntry) + % Use existing eq_id + eq_id = existingEntry{1,1}; else - postfilter_taps = pf.burg_coeff; + % Insert the new equalizer configuration and get its eq_id + eq_id = obj.appendToTable('EqualizerParameters', eqParamsData); end - % Create equalizer data struct for searching and adding if necessary - equalizerData = struct( ... - 'ffe', jsonencode(ffe), ... - 'dfe', jsonencode(dfe), ... - 'mlse', jsonencode(mlse), ... - 'pf', jsonencode(pf), ... - 'eq_type', string(eqType), ... - 'ffe_order', jsonencode(ffe_order), ... - 'dfe_order', jsonencode(dfe_order), ... - 'postfilter_taps',jsonencode(postfilter_taps),... - 'len_tr', len_tr, ... - 'mu_ffe', jsonencode(mu_ffe), ... - 'mu_dfe', mu_dfe, ... - 'mu_dc', mu_dc, ... - 'comment', comment ... - ); + % 3. Add the equalizer configuration reference and run_id to resultData + resultData.eqParam_id = eq_id; + resultData.run_id = run_id; - % Check if exact Equalizer and BER entries already exist in the DB ... - selectedFields = {'Runs.run_id','BERs.ber_id','Equalizer.eq_id','BERs.ber',['BERs.occurrence' ... - '']}; - filterParams = obj.tables; - filterParams.Equalizer = equalizerData; - [dataTable,sql_query] = obj.queryDB(filterParams, selectedFields); + % 4. Compute hash for the processing result + tempResultData = rmfield(resultData, 'date_of_processing'); + resultJsonStr = jsonencode(tempResultData); + md2 = java.security.MessageDigest.getInstance('MD5'); % Create a new MD5 instance + md2.update(uint8(resultJsonStr)); + resultHashBytes = typecast(md2.digest, 'uint8'); + resultHashStr = lower(dec2hex(resultHashBytes)'); + resultHashStr = lower(strtrim(resultHashStr(:)')); % Convert to a lowercase string - % get or insert Equalizer - if ~isempty(dataTable) - % Equalizer entry already exists, use the existing eq_id - cur_eq_id = dataTable.eq_id; - else - % Insert the new Equalizer entry - cur_eq_id = obj.appendToTable('Equalizer', equalizerData); + % Add the result hash to resultData + resultData.result_hash = resultHashStr; + + % 5. Check if an identical processing result already exists + queryStr2 = sprintf('SELECT result_id FROM Results WHERE result_hash = ''%s''', resultData.result_hash); + existingResult = obj.fetch(queryStr2); + + if ~isempty(existingResult) + % If the result exists, return its result_id without inserting a new row + resultID = existingResult{1,1}; + warning(['Result already exists: ResultID: ',num2str(resultID),'| EQ ID: ',num2str(eq_id),' Run ID: ' num2str(run_id)]) + return; end - % skip if already here or insert BER entry - if ~isempty(dataTable) + % 6. Insert the processing result + resultID = obj.appendToTable('Results', resultData); + end - % A BER entry with the same eq_id, run_id, and occurrence already exists - existingBERValue = dataTable.ber; + function recalcHashes(obj) + % recalcHashes Recalculate hashes for all rows in the Results and EqualizerParameters tables. + % + % For EqualizerParameters, the hash is computed from all fields except + % 'eq_id' and 'config_hash'. + % + % For Results, the hash is computed from all fields except 'result_id', + % 'result_hash', and 'date_of_processing'. - % Compare the existing BER value with the new BER value - if existingBERValue == berValue - fprintf('The BER entry %.2e || -- eq_id: %d -- run_id: %d -- occurrence: %d already exists. \n',berValue, cur_eq_id, runID, occurrence); - else - fprintf('Already found BER for EQ: %.2e ~= %.2e || -- eq_id: %d -- run_id: %d -- occurrence: %d already exists.\n', berValue, existingBERValue, cur_eq_id, runID, occurrence); + % Recalculate hashes for EqualizerParameters + eqParamsRows = obj.fetch('SELECT * FROM EqualizerParameters'); + for i = 1:height(eqParamsRows) + rowStruct = table2struct(eqParamsRows(i,:)); % Convert the table row to a struct + + % Remove fields not part of the hash computation + if isfield(rowStruct, 'eq_id') + rowStruct = rmfield(rowStruct, 'eq_id'); + end + if isfield(rowStruct, 'config_hash') + rowStruct = rmfield(rowStruct, 'config_hash'); end - else - - % No such BER entry exists, insert the new BER entry - berData = struct( ... - 'run_id', runID, ... - 'eq_id', cur_eq_id, ... - 'ber', berValue, ... - 'occurrence', occurrence ... - ); - - obj.appendToTable('BERs', berData); + % Compute MD5 hash from the JSON representation + jsonStr = jsonencode(rowStruct); + md = java.security.MessageDigest.getInstance('MD5'); + md.update(uint8(jsonStr)); + hashBytes = typecast(md.digest, 'uint8'); + hashStr = lower(dec2hex(hashBytes)'); + hashStr = lower(strtrim(hashStr(:)')); + % Update the config_hash field using the eq_id from the table row + eq_id = eqParamsRows.eq_id(i); + updateQuery = sprintf('UPDATE EqualizerParameters SET config_hash = ''%s'' WHERE eq_id = %d', hashStr, eq_id); + obj.executeSQL(updateQuery); end + + % Fetch all rows from the Results table using queryDB with no filters + Results = obj.tables.Results; + + if isstruct(Results) + newFields = {}; + fNames = fieldnames(Results); + for i = 1:numel(fNames) + newFields{end+1} = ['Results.' fNames{i}]; + end + Results = newFields; + end + [resultsRows, ~] = obj.queryDB(obj.tables, Results); + + for i = 1:height(resultsRows) + rowStruct = table2struct(resultsRows(i, :)); % Convert the row to a struct + + % Remove fields not used in the hash calculation + if isfield(rowStruct, 'result_id') + rowStruct = rmfield(rowStruct, 'result_id'); + end + if isfield(rowStruct, 'result_hash') + rowStruct = rmfield(rowStruct, 'result_hash'); + end + if isfield(rowStruct, 'date_of_processing') + rowStruct = rmfield(rowStruct, 'date_of_processing'); + end + + % Compute MD5 hash from the JSON representation + jsonStr = jsonencode(rowStruct); + md = java.security.MessageDigest.getInstance('MD5'); + md.update(uint8(jsonStr)); + hashBytes = typecast(md.digest, 'uint8'); + hashStr = lower(dec2hex(hashBytes)'); + hashStr = lower(strtrim(hashStr(:)')); + + % Access the result_id from the table row and update the hash + result_id = resultsRows.result_id(i); + updateQuery = sprintf('UPDATE Results SET result_hash = ''%s'' WHERE result_id = %d', ... + hashStr, result_id); + obj.executeSQL(updateQuery); + end + + + + end + function executeSQL(obj, query) % This method executes an SQL statement using MATLAB's execute function. execute(obj.conn, query); end - function answer = fetch(obj,query) answer = fetch(obj.conn,query); end @@ -405,15 +487,46 @@ classdef DBHandler < handle function query = constructSQLQuery(obj, filterParams, selectedFields) % constructSQLQuery Constructs the SQL query based on filter parameters and selected fields. + % + % If selectedFields is provided as a struct, it is converted to a cell array. + % The conversion takes the field names and creates entries like: + % {'selectedFields.fieldName'} for each field. - % Construct the SELECT clause dynamically based on user selection + % Input check for selectedFields: if it's a struct, convert it to a cell array. + if isstruct(selectedFields) + newFields = {}; + tableNames = fieldnames(selectedFields); + for t = 1:numel(tableNames) + tableStruct = selectedFields.(tableNames{t}); + fieldNames = fieldnames(tableStruct); + for f = 1:numel(fieldNames) + if isequal(tableStruct.(fieldNames{f}), 1) + newFields{end+1} = sprintf('%s.%s', tableNames{t}, fieldNames{f}); + end + end + end + selectedFields = newFields; + end + + % Construct the SELECT clause dynamically based on user selection. + % (Assuming that when provided as a cell array, each entry is of the form + % 'TableName.fieldName' or, in our conversion case, 'selectedFields.fieldName'.) selectClause = 'SELECT DISTINCT '; for i = 1:numel(selectedFields) fieldParts = strsplit(selectedFields{i}, '.'); - tableName = fieldParts{1}; - fieldName = fieldParts{2}; + % If the field comes from the struct conversion, its first part is 'selectedFields' + % and the actual field name is in the second part. + if strcmp(fieldParts{1}, 'selectedFields') + tableName = fieldParts{1}; % not used for type checking below + fieldName = fieldParts{2}; + else + tableName = fieldParts{1}; + fieldName = fieldParts{2}; + end - if isnumeric(obj.tables.(tableName).(fieldName)) + % Decide on COALESCE depending on the field type. + % If the table is known in obj.tables and the field is numeric, use 'NaN'. + if isfield(obj.tables, tableName) && isfield(obj.tables.(tableName), fieldName) && isnumeric(obj.tables.(tableName).(fieldName)) selectClause = [selectClause, 'COALESCE(', selectedFields{i}, ', ''NaN'') AS ', fieldName]; else selectClause = [selectClause, 'COALESCE(', selectedFields{i}, ', '''') AS ', fieldName]; @@ -426,235 +539,300 @@ classdef DBHandler < handle end end - % Construct the FROM and WHERE clause - baseQuery = [selectClause, 'FROM Runs ' ... - 'LEFT JOIN Configurations ON Runs.run_id = Configurations.run_id ' ... - 'LEFT JOIN Measurements ON Runs.run_id = Measurements.run_id ' ... - 'WHERE ']; - % 'LEFT JOIN BERs ON Runs.run_id = BERs.run_id ' ... - % 'LEFT JOIN Equalizer ON BERs.eq_id = Equalizer.eq_id ' ... - + % --- Adaptive FROM Clause --- + % Use "Runs" as the main table and add LEFT JOINs for every other table in obj.tables + % (except "sqlite_sequence") that has a run_id field. + mainTable = 'Runs'; + fromClause = ['FROM ', mainTable, ' ']; + tableNamesAll = fieldnames(obj.tables); + for t = 1:numel(tableNamesAll) + tableName = tableNamesAll{t}; + if strcmpi(tableName, mainTable) || strcmpi(tableName, 'sqlite_sequence') + continue; + end - % Loop through each table in filterParams + if isfield(obj.tables.(tableName), 'run_id') + % most tables are directly linked to runs table + fromClause = [fromClause, 'LEFT JOIN ', tableName, ' ON ', mainTable, '.run_id = ', tableName, '.run_id ']; + elseif isfield(obj.tables.(tableName), 'eq_id') + % equalizer is only linked to results table + fromClause = [fromClause, 'LEFT JOIN ', tableName, ' ON ', 'Results', '.eqParam_id = ', tableName, '.eq_id ']; + end + end + + % --- WHERE Clause Construction --- + baseQuery = [selectClause, ' ', fromClause, 'WHERE ']; filterClauses = []; tableNames_ = fieldnames(filterParams); - for t = 1:numel(tableNames_) tableName = tableNames_{t}; tableParams = filterParams.(tableName); - - % Loop through each parameter in the table fieldNames = fieldnames(tableParams); for i = 1:numel(fieldNames) fieldName = fieldNames{i}; value = tableParams.(fieldName); - - % Construct the full column name in the format "tableName.fieldName" fullName = sprintf('%s.%s', tableName, fieldName); - % Handle different types of values for SQL query construction + % Handle various types of values for SQL query construction if isempty(value) - % Skip this parameter if it is empty (include all values) continue; elseif isnumeric(value) && isnan(value) - % If value is NaN, use IS NULL in SQL filterClause = sprintf('%s IS NULL', fullName); elseif isnumeric(value) && ~isEnumeration(value) filterClause = sprintf('%s = %f', fullName, value); elseif islogical(value) || (isnumeric(value) && ismember(value, [0, 1])) && ~isEnumeration(value) filterClause = sprintf('%s = %d', fullName, value); elseif ischar(value) || isstring(value) - filterClause = sprintf('%s = ''%s''', fullName, char(value)); %nicht nach string suchen sondern nach chararray -> 'bla' statt "bla" + filterClause = sprintf('%s = ''%s''', fullName, char(value)); elseif isEnumeration(value) filterClause = sprintf('%s = ''%s''', fullName, value); else error('Unsupported data type for field "%s".', fullName); end - - % Add the constructed filter clause to the list filterClauses = [filterClauses, filterClause, ' AND ']; end end - % Remove trailing ' AND ' from the filter clauses if any filters were added + % Remove trailing ' AND ' if any filters were added. if ~isempty(filterClauses) filterClauses = filterClauses(1:end-5); - end - - % Construct the final SQL query - if isempty(filterClauses) - query = [selectClause, 'FROM Runs ' ... - 'LEFT JOIN Configurations ON Runs.run_id = Configurations.run_id ' ... - 'LEFT JOIN Measurements ON Runs.run_id = Measurements.run_id ' ... - 'LEFT JOIN BERs ON Runs.run_id = BERs.run_id']; + query = [selectClause, ' ', fromClause, 'WHERE ', filterClauses]; else - query = [baseQuery, filterClauses]; + query = [selectClause, ' ', fromClause]; end - end - - function selectedFields = promptSelectFields(obj) - % promptSelectFields Prompts the user to select fields from multiple tables to include in the SELECT statement using settingsdlg. + % promptSelectFields Prompts the user to select fields from multiple tables + % using a custom checkbox GUI with scrolling. + % + % The function builds a list of all fields (formatted as 'TableName.fieldName') + % and displays each as a checkbox inside an inner container panel. The container's + % height is set to accommodate all checkboxes, so the scrollable panel shows scrollbars. + % When the user clicks the "Select" button, the selected fields are returned. + % If none are selected, all fields are returned. - % Get all possible fields from all tables (excluding sqlite_sequence) + % Get all possible tables (excluding sqlite_sequence) tableNames = fieldnames(obj.tables); - tableNames = setdiff(tableNames, {'sqlite_sequence'}); % Remove sqlite_sequence - - % Prepare the inputs for settingsdlg - promptSettings = {}; - allFieldsFullName = {}; - convertedFieldNames = {}; + tableNames = setdiff(tableNames, {'sqlite_sequence'}); + % Build a single cell array of all field names with table prefix. + allFields = {}; for i = 1:numel(tableNames) - tableFields = fieldnames(obj.tables.(tableNames{i})); - for j = 1:numel(tableFields) - fieldName = tableFields{j}; - fullName = sprintf('%s.%s', tableNames{i}, fieldName); - convertedName = strrep(fullName, '.', '_'); % Replace '.' with '_' - - allFieldsFullName{end + 1} = fullName; % Add full name to the list - convertedFieldNames{end + 1} = convertedName; % Store the converted name - - % Add the field name and checkbox setting to the prompt - promptSettings{end + 1} = {sprintf('Include %s', fullName), convertedName}; - promptSettings{end + 1} = false; % Default: not selected + fields = fieldnames(obj.tables.(tableNames{i})); + for j = 1:numel(fields) + allFields{end+1} = sprintf('%s.%s', tableNames{i}, fields{j}); end end + numFields = numel(allFields); - % Create the settings dialog - [settings, button] = settingsdlg(... - 'title', 'Select Fields for the SQL Query', ... - 'description', 'Check the boxes for the fields you want to include in the SELECT statement.', ... - promptSettings{:} ... - ); + % Create the main figure. + fig = uifigure('Name', 'Select Fields', 'Position', [100, 100, 400, 600]); - % If the user cancels, default to selecting all fields - if strcmp(button, 'cancel') - selectedFields = allFieldsFullName; - return; + % Create a scrollable panel. + scrollPanel = uipanel(fig, 'Position', [10, 60, 380, 530], 'Scrollable', 'on'); + + % Define checkbox dimensions. + checkboxHeight = 30; + spacing = 5; + totalHeight = numFields * (checkboxHeight + spacing); + + % Create an inner container panel with height larger than the scrollPanel's height. + container = uipanel(scrollPanel, 'Position', [0, 0, scrollPanel.Position(3), totalHeight]); + + % Create checkboxes using absolute positioning in the container. + checkboxes = gobjects(numFields, 1); + for i = 1:numFields + % Calculate the vertical position. + % The origin (0,0) is at the bottom left of the container. + yPos = totalHeight - i*(checkboxHeight + spacing) + spacing; + checkboxes(i) = uicheckbox(container, ... + 'Text', allFields{i}, ... + 'Value', false, ... + 'Position', [10, yPos, container.Position(3)-20, checkboxHeight]); end - % Parse user input into selectedFields + % Create a "Select" button in the main figure. + btn = uibutton(fig, 'Text', 'Select', ... + 'Position', [150, 10, 100, 30], ... + 'ButtonPushedFcn', @(btn, event) uiresume(fig)); + + % Wait for the user to click the button. + uiwait(fig); + + % Retrieve the selected fields. selectedFields = {}; - for i = 1:numel(allFieldsFullName) - convertedName = convertedFieldNames{i}; - if isfield(settings, convertedName) && settings.(convertedName) % Add to selectedFields if the checkbox was selected - selectedFields{end + 1} = allFieldsFullName{i}; %#ok + for i = 1:numFields + if checkboxes(i).Value + selectedFields{end+1} = checkboxes(i).Text; end end - % If no fields are selected, default to selecting all fields + % If no fields are selected, default to all fields. if isempty(selectedFields) - selectedFields = allFieldsFullName; + selectedFields = allFields; end + + % Close the figure. + delete(fig); end function filterParams = promptFilterParameters(obj) - % promptFilterParameters Prompts the user to enter filter parameters using the settingsdlg framework. + % promptFilterParameters Prompts the user to enter filter parameters using a + % custom scrollable UI with dropdowns. + % + % For each table (excluding 'sqlite_sequence'), each field that has distinct + % values is displayed as a label and a dropdown. The dropdown items are built + % from the distinct values (with "All" prepended). The output is a struct where, + % for each table, each field is set to the chosen value (or [] if "All" is selected). - % Get all possible parameters from all tables (excluding sqlite_sequence) - tableNames_ = fieldnames(obj.tables); - tableNames_ = setdiff(tableNames_, {'sqlite_sequence'}); % Remove sqlite_sequence + % Get all tables except 'sqlite_sequence' + tableNames = fieldnames(obj.tables); + tableNames = setdiff(tableNames, {'sqlite_sequence'}); - % Prepare the inputs for settingsdlg with sections and separators - promptSettings = {}; - allFieldsFullName = {}; - convertedFieldNames = {}; + % Precompute layout constants. + heightPerTableLabel = 30; + heightPerField = 40; % vertical space for a field (label + dropdown) + spacing = 5; - for i = 1:numel(tableNames_) - % Add a separator for each table section - promptSettings{end + 1} = 'separator'; - promptSettings{end + 1} = tableNames_{i}; - - % Get all fields from the current table - tableFields = fieldnames(obj.tables.(tableNames_{i})); - - % Prepare each field to be added to the dialog + % Compute total required height. + totalHeight = 0; + for i = 1:numel(tableNames) + totalHeight = totalHeight + heightPerTableLabel; + tableName = tableNames{i}; + tableFields = fieldnames(obj.tables.(tableName)); for j = 1:numel(tableFields) fieldName = tableFields{j}; - fullName = sprintf('%s.%s', tableNames_{i}, fieldName); - convertedName = strrep(fullName, '.', '_'); % Replace '.' with '_' - - % Skip fields that do not have distinct values stored - if ~isfield(obj.distinctValues.(tableNames_{i}), fieldName) - continue; + % Only include fields that have distinct values stored. + if isfield(obj.distinctValues.(tableName), fieldName) + totalHeight = totalHeight + heightPerField; end - - % Get the distinct values for the field - distinctValues_ = obj.distinctValues.(tableNames_{i}).(fieldName); - - % Prepare distinct values for dropdown - if isempty(distinctValues_) - % If there are no distinct values, use only an "All" entry - distinctValues_ = {'All'}; - else - % Ensure distinctValues is a cell array of strings - if isnumeric(distinctValues_) - distinctValues_ = arrayfun(@(x) num2str(x), distinctValues_, 'UniformOutput', false); - elseif isstring(distinctValues_) - distinctValues_ = cellstr(distinctValues_); - elseif iscell(distinctValues_) && ~iscellstr(distinctValues_) - distinctValues_ = cellfun(@num2str, distinctValues_, 'UniformOutput', false); - end - - % Add an "All" option at the beginning of the distinct values list - distinctValues_ = [{'All'}; distinctValues_]; - end - - allFieldsFullName{end + 1} = fullName; % Add full name to the list - convertedFieldNames{end + 1} = convertedName; % Store the converted name - - % Add the field name and value setting to the prompt - promptSettings{end + 1} = {sprintf('%s', fullName), convertedName}; - promptSettings{end + 1} = distinctValues_; % Add distinct values as dropdown options end end - % Create the settings dialog - [settings, button] = settingsdlg(... - 'title', 'Input Parameters for Filtering', ... - 'description', 'Enter the values for each field to filter. Select "All" to include all values.', ... - promptSettings{:} ... - ); + % Create the main UI figure. + fig = uifigure('Name', 'Input Parameters for Filtering', 'Position', [100, 100, 500, 600]); - % If the user cancels, return an empty struct - if strcmp(button, 'cancel') + % Set a CloseRequestFcn so that closing the figure calls uiresume. + fig.CloseRequestFcn = @(src, event) uiresume(src); + + % Create a scrollable panel inside the figure. + scrollPanel = uipanel(fig, 'Position', [10, 60, 480, 530], 'Scrollable', 'on'); + + % Create an inner container panel with a height set to totalHeight. + container = uipanel(scrollPanel, 'Position', [0, 0, scrollPanel.Position(3), totalHeight]); + + % Prepare cell arrays to store dropdown handles and corresponding table/field names. + dropdownHandles = {}; + dropdownTableNames = {}; + dropdownFieldNames = {}; + + % Set the starting Y coordinate (filling from top to bottom). + currentY = totalHeight; + + % Maximum number of dropdown items. + maxItems = 100; + + for i = 1:numel(tableNames) + % Create a label for the table name. + uilabel(container, ... + 'Text', tableNames{i}, ... + 'FontWeight', 'bold', ... + 'Position', [10, currentY - heightPerTableLabel + spacing, 200, heightPerTableLabel - spacing]); + currentY = currentY - heightPerTableLabel; + + tableName = tableNames{i}; + tableFields = fieldnames(obj.tables.(tableName)); + for j = 1:numel(tableFields) + fieldName = tableFields{j}; + if ~isfield(obj.distinctValues.(tableName), fieldName) + continue; % Skip if no distinct values are stored. + end + + % Retrieve distinct values for the field. + distinctValues_ = obj.distinctValues.(tableName).(fieldName); + if ~isempty(distinctValues_) && numel(distinctValues_) > maxItems + distinctValues_ = distinctValues_(1:maxItems); + end + + if isempty(distinctValues_) + items = {'All'}; + else + if isnumeric(distinctValues_) + items = cellfun(@num2str, num2cell(distinctValues_), 'UniformOutput', false); + elseif isstring(distinctValues_) + items = cellstr(distinctValues_); + elseif iscell(distinctValues_) && ~iscellstr(distinctValues_) + items = cellfun(@num2str, distinctValues_, 'UniformOutput', false); + else + items = distinctValues_; + end + items = items(:)'; % Ensure row vector + items = [{'All'}, items]; + end + + % Create a label for the field. + uilabel(container, ... + 'Text', sprintf('%s:', fieldName), ... + 'HorizontalAlignment', 'right', ... + 'Position', [10, currentY - 25, 150, 25]); + + % Create a dropdown for the field. + dd = uidropdown(container, ... + 'Items', items, ... + 'Value', 'All', ... + 'Position', [170, currentY - 25, 200, 25]); + + % Store the dropdown handle and its associated table/field. + dropdownHandles{end+1} = dd; + dropdownTableNames{end+1} = tableName; + dropdownFieldNames{end+1} = fieldName; + + currentY = currentY - heightPerField; + end + end + + % Create a "Submit" button at the bottom of the figure. + btn = uibutton(fig, 'Text', 'Submit', ... + 'Position', [200, 10, 100, 30], ... + 'ButtonPushedFcn', @(btn, event) uiresume(fig)); + + % Wait until the user clicks "Submit" or closes the figure. + uiwait(fig); + + % If the figure was closed (and thus no longer valid), return an empty struct. + if ~isvalid(fig) filterParams = struct(); return; end - % Parse user input into filterParams structure + % Build the filterParams struct from the dropdown selections. filterParams = struct(); - - for i = 1:numel(allFieldsFullName) - value = settings.(convertedFieldNames{i}); - - % Split full name to get table and field names - fieldParts = strsplit(allFieldsFullName{i}, '.'); - tableName = fieldParts{1}; - fieldName = fieldParts{2}; - - % If the table does not exist in the filterParams struct, create it + for k = 1:numel(dropdownHandles) + tableName = dropdownTableNames{k}; + fieldName = dropdownFieldNames{k}; + value = dropdownHandles{k}.Value; if ~isfield(filterParams, tableName) filterParams.(tableName) = struct(); end - % Assign values to the respective fields under each table + % If "All" is selected, assign empty; otherwise, try converting to numeric. if strcmp(value, 'All') - filterParams.(tableName).(fieldName) = []; % Set to empty to include all values - elseif isnumeric(value) && isnan(value) - filterParams.(tableName).(fieldName) = NaN; % Use NaN to handle as NULL + filterParams.(tableName).(fieldName) = []; else - filterParams.(tableName).(fieldName) = value; % Use the entered value + numValue = str2double(value); + if ~isnan(numValue) + filterParams.(tableName).(fieldName) = numValue; + else + filterParams.(tableName).(fieldName) = value; + end end end + + % Close the figure. + delete(fig); end - - - end end diff --git a/Classes/Moveit_wrapper.m b/Classes/Moveit_wrapper.m new file mode 100644 index 0000000..96dce6a --- /dev/null +++ b/Classes/Moveit_wrapper.m @@ -0,0 +1,180 @@ +classdef Moveit_wrapper < handle + % + + properties(Access=public) + moveit_function_name + para + state = struct() + comment + end + + methods (Access=public) + function obj = Moveit_wrapper(moveit_function_name,options) + % + + arguments + moveit_function_name + options.para + end + + obj.moveit_function_name = moveit_function_name; + fn = fieldnames(options); + for n = 1:numel(fn) + try + obj.(fn{n}) = options.(fn{n}); + end + end + + % Ensure the function exists + if ~exist(obj.moveit_function_name, 'file') + error('Function "%s" does not exist.', obj.moveit_function_name); + end + + % Step 1: Get default parameters by calling moveit module + if isempty(obj.para) + global loop; + loop = 0; % Request parameter structure + [obj.para, obj.comment] = feval(obj.moveit_function_name); + end + + end + + function signal_out = process(obj,signal_in) + + arguments + obj + signal_in = [] + + end + + isSignalClass = 0; + if isa(signal_in,"Signal") + signal_out = signal_in; + signal_in = signal_in.signal; + isSignalClass = 1; + end + + % signal_in = signal_in.'; + + % INIT MOVEIT + global loop; + loop = 0; + signal_out = obj.moveit_init(signal_in); + + % RUN MOVEIT + global loop; + loop = 1; + + signal_out = obj.moveit_init(signal_out); + + if isSignalClass + % append to logbook + lbdesc = ['Logbookentry']; + signal_out = signal_out.logbookentry(lbdesc); + signal_out.signal = signal_out; + end + + end + + function open_settings(obj) + % Opens a settings dialog for modifying Moveit parameters with appropriate input types. + % Inserts separators when the level value changes. + + fieldNames = fieldnames(obj.para); + settingsCell = {}; + previousLevel = -inf; % Track previous level for separators + + % Create input arguments for settingsdlg + for i = 1:numel(fieldNames) + fieldName = fieldNames{i}; + value = obj.para.(fieldName); + + % Get description from comment structure if available + if isfield(obj.comment, fieldName) + description = obj.comment.(fieldName); + else + description = fieldName; % Default to field name if no description + end + + % Get level for this parameter + if isfield(obj.comment.level, fieldName) + currentLevel = obj.comment.level.(fieldName); + else + currentLevel = 1; % Default level + end + + % Insert separator if the level changes + if currentLevel ~= previousLevel + settingsCell{end+1} = 'separator'; + settingsCell{end+1} = ['Level ', num2str(currentLevel)]; + end + previousLevel = currentLevel; + + % Determine input type from obj.comment.type + if 1 + inputType = obj.comment.type.(fieldName); + else + inputType = 'number'; % Default type + end + + % Add field description + settingsCell{end+1} = {description; fieldName}; + + % Define the input type + if strcmp(inputType, 'boolean') + % Checkbox input + settingsCell{end+1} = logical(value); % Enable checkbox + elseif contains(inputType, '|') + % Dropdown selection (options separated by '|') + dropdownOptions = strsplit(inputType, '|'); + settingsCell{end+1} = dropdownOptions; + else + % Default: Numeric input + settingsCell{end+1} = value; + end + end + + % Open settings dialog + [newValues, button] = settingsdlg(... + 'title', ['Settings: ', obj.moveit_function_name], ... + 'description', ['Adjust parameters for ', obj.moveit_function_name], ... + settingsCell{:} ... + ); + + % Update obj.para with new values only if "OK" was pressed + if strcmp(button, 'OK') + fieldNames = fieldnames(newValues); + for i = 1:numel(fieldNames) + if fieldNames{i} + + end + obj.para.(fieldNames{i}) = newValues.(fieldNames{i}); + end + end + end + + end + + + methods (Access=private) + + function data_out = moveit_init(obj,data_in) + arguments(Input) + obj + data_in double + end + % Step 2: Initialize the function + [data_out, obj.state] = feval(obj.moveit_function_name, data_in, obj.state, obj.para); + end + + function data_out = moveit_simulation_computation(obj,data_in) + arguments(Input) + obj + data_in double + end + % Step 3: Process input signal + [data_out, obj.state] = feval(obj.moveit_function_name, data_in, obj.state, obj.para); + end + + end +end diff --git a/Classes/Warehouse_class/classes/DataStorage.m b/Classes/Warehouse_class/classes/DataStorage.m index 434bb97..91530cf 100644 --- a/Classes/Warehouse_class/classes/DataStorage.m +++ b/Classes/Warehouse_class/classes/DataStorage.m @@ -89,7 +89,7 @@ classdef DataStorage < handle for p = 1:numel(obj.fn) name = obj.fn(p); values = obj.inputParams.(name); - obj.parameter.(name) = Parameter(name,values); + obj.parameter.(name) = StorageParameter(name,values); end end diff --git a/Classes/Warehouse_class/classes/Parameter.m b/Classes/Warehouse_class/classes/StorageParameter.m similarity index 88% rename from Classes/Warehouse_class/classes/Parameter.m rename to Classes/Warehouse_class/classes/StorageParameter.m index c5243f0..97ef06a 100644 --- a/Classes/Warehouse_class/classes/Parameter.m +++ b/Classes/Warehouse_class/classes/StorageParameter.m @@ -1,4 +1,4 @@ -classdef Parameter < handle +classdef StorageParameter < handle %PARAMETER Summary of this class goes here % Detailed explanation goes here @@ -12,7 +12,7 @@ classdef Parameter < handle end methods - function obj = Parameter(name, values) + function obj = StorageParameter(name, values) %PARAMETER Construct obj.name = name; obj.values = values; diff --git a/Datatypes/clr.m b/Datatypes/clr.m new file mode 100644 index 0000000..b9aa515 --- /dev/null +++ b/Datatypes/clr.m @@ -0,0 +1,70 @@ +classdef clr + properties (Constant) + % Set1 colormap + Set1 = struct( ... + 'red', [0.8941, 0.1020, 0.1098], ... + 'blue', [0.2157, 0.4941, 0.7216], ... + 'green', [0.3020, 0.6863, 0.2902], ... + 'purple', [0.5961, 0.3059, 0.6392], ... + 'orange', [1.0000, 0.4980, 0.0000], ... + 'yellow', [1.0000, 1.0000, 0.2000], ... + 'brown', [0.6510, 0.3373, 0.1569], ... + 'pink', [0.9686, 0.5059, 0.7490], ... + 'gray', [0.6000, 0.6000, 0.6000]); + + % Paired colormap + Paired = struct( ... + 'lightblue', [0.6510, 0.8078, 0.8902], ... + 'blue', [0.1216, 0.4706, 0.7059], ... + 'lightgreen', [0.6980, 0.8745, 0.5412], ... + 'green', [0.2000, 0.6275, 0.1725], ... + 'lightred', [0.9843, 0.6039, 0.6000], ... + 'red', [0.8902, 0.1020, 0.1098], ... + 'lightorange',[0.9922, 0.7490, 0.4353], ... + 'orange', [1.0000, 0.4980, 0.0000], ... + 'lightpurple',[0.7922, 0.6980, 0.8392], ... + 'purple', [0.4157, 0.2392, 0.6039], ... + 'lightyellow',[1.0000, 1.0000, 0.6000], ... + 'brown', [0.6941, 0.3490, 0.1569]); + + end + + methods (Static) + function showRGB(colorArray, names) + % Visualize colors in a horizontal bar chart + if nargin < 2 + names = repmat({''}, size(colorArray, 1), 1); + end + + figure; + hold on; + for i = 1:size(colorArray, 1) + fill([0 1 1 0], [i-1 i-1 i i], colorArray(i, :), 'EdgeColor', 'k'); + text(1.1, i-0.5, names{i}, 'FontSize', 12, 'Interpreter', 'none'); + end + ylim([0, size(colorArray, 1)]); + xlim([0, 1.5]); + axis off; + title('Color Preview'); + hold off; + end + + function showSet(colorStruct) + % Show colors from a structure in a bar plot + names = fieldnames(colorStruct); + colors = cell2mat(struct2cell(colorStruct)'); + + figure; + hold on; + for i = 1:size(colors, 1) + fill([0 1 1 0], [i-1 i-1 i i], colors(i, :), 'EdgeColor', 'k'); + text(1.1, i-0.5, names{i}, 'FontSize', 12, 'Interpreter', 'none'); + end + ylim([0, size(colors, 1)]); + xlim([0, 1.5]); + axis off; + title('Color Preview'); + hold off; + end + end +end diff --git a/Datatypes/equalizer_structure.m b/Datatypes/equalizer_structure.m index 263260b..7e338bd 100644 --- a/Datatypes/equalizer_structure.m +++ b/Datatypes/equalizer_structure.m @@ -4,7 +4,8 @@ classdef equalizer_structure < int32 ffe (0) vnle (1) vnle_pf_mlse (2) - db_precoded (3) + % db_precoded (3) + vnle_db_mlse (3) db_encoded (4) end diff --git a/Datatypes/pulseform.m b/Datatypes/pulseform.m index 1c524c1..d753291 100644 --- a/Datatypes/pulseform.m +++ b/Datatypes/pulseform.m @@ -1,7 +1,8 @@ classdef pulseform < int32 enumeration - rrc (1) + rc (2) %raised cosine (use this at Tx without matched filter) + rrc (1) %root raised cosine (usually with matched filter) end end \ No newline at end of file diff --git a/Functions/EQ_structures/duobinary_signaling.m b/Functions/EQ_structures/duobinary_signaling.m index 73f131c..69c0449 100644 --- a/Functions/EQ_structures/duobinary_signaling.m +++ b/Functions/EQ_structures/duobinary_signaling.m @@ -1,18 +1,96 @@ -function [eq_package] = duobinary_signaling(eq_, mlse_,M ,rx_signal, tx_symbols, tx_bits) - %Duobinary Signaling +function [eq_package] = duobinary_signaling(eq_, mlse_,M ,rx_signal, tx_symbols, tx_bits,options) +%Duobinary Signaling +arguments + eq_ + mlse_ + M + rx_signal + tx_symbols + tx_bits + options.postFFE = []; +end - [eq_signal, eq_noise] = eq_.process(rx_signal,tx_symbols); - eq_signal = mlse_.process(eq_signal); +[eq_signal, eq_noise] = eq_.process(rx_signal,tx_symbols); - eq_signal = Duobinary().encode(eq_signal); - eq_signal = Duobinary().decode(eq_signal); +if ~isempty(options.postFFE) + [eq_signal,eq_noise] = options.postFFE.process(eq_signal,tx_symbols); +end - % M = numel(unique(eq_signal.signal)); - rx_bits = PAMmapper(M,0).demap(eq_signal); - [~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); +eq_signal = mlse_.process(eq_signal); + +eq_signal = Duobinary().encode(eq_signal); +eq_signal = Duobinary().decode(eq_signal); + +% M = numel(unique(eq_signal.signal)); +rx_bits = PAMmapper(M,0).demap(eq_signal); + +[bits_db,errors_db,ber_db,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + +eq_package.ber = ber_db; + +resultsDBsignaling = struct( ... + 'result_id', NaN, ... % + 'run_id', NaN, ... % Beispielhafte Run-ID + 'eqParam_id', NaN, ... % Beispielhafter Fremdschlüssel zur EqualizerParameters-Tabelle + 'date_of_processing', datetime('now'), ... % Aktuelles Datum und Uhrzeit + 'numBits', bits_db, ... % Beispiel: 1.000.000 Bits + 'numBitErr', errors_db, ... % Beispiel: 120 Bitfehler + 'BER_precoded', [], ... % BER = 120 / 1.000.000 + 'numBitErr_precoded', [], ... % Beispiel: 120 Bitfehler + 'BER', ber_db, ... % BER = 120 / 1.000.000 + 'SNR', [], ... % Beispielhafte SNR + 'SNR_level', jsonencode([]), ... % SNR-Level als JSON-codiertes Array + 'GMI', [], ... % Beispielhafter GMI-Wert + 'AIR', [], ... % Beispielhafter AIR-Wert + 'EVM', [], ... % Beispielhafte EVM + 'EVM_level', jsonencode([]), ... % EVM-Level als JSON-codiertes Array + 'Alpha', [] ... % Beispielhafter Alpha-Wert + ); + +if ~isempty(options.postFFE) + npostFFE = options.postFFE.order; +else + npostFFE = 0; +end + +equalizerConfigDBsignaling = struct( ... + 'eq_id', NaN, ... % Auto-Inkrement, wird in der DB gesetzt + 'equalizer_structure', int32(equalizer_structure.db_encoded), ... % Beispiel: 1 (z.B. für vnle) + 'M', M, ... % Ordnung der PAM-Konstellation + 'target_constellation', jsonencode(round(unique(tx_symbols.signal),5)), ... % Beispielhafter Target-String + 'db_target', 1, ... % 0 oder 1 + 'diff_precode', 1, ... % 0 oder 1 + 'postFFE', ~isempty(options.postFFE), ... % Beispielwert + 'NpostFFE', npostFFE, ... % Beispielwert + 'Ne1', eq_.Ne(1), ... % Feedforward Koeffizienten 1. Ordnung + 'Ne2', eq_.Ne(2), ... % Feedforward Koeffizienten 2. Ordnung + 'Ne3', eq_.Ne(3), ... % Feedforward Koeffizienten 3. Ordnung + 'Nb1', eq_.Nb(1), ... % Decision Feedback Koeffizienten 1. Ordnung + 'Nb2', eq_.Nb(2), ... % Decision Feedback Koeffizienten 2. Ordnung + 'Nb3', eq_.Nb(3), ... % Decision Feedback Koeffizienten 3. Ordnung + 'K', eq_.K, ... % Samples pro Symbol + 'DCmu', eq_.DCmu, ... % Anpassungsrate für DC-Tap + 'ideal_dfe', eq_.ideal_dfe, ... % Flag für ideal DFE (0 oder 1) + 'training_length', eq_.training_length, ... % Anzahl Trainingssymbole + 'training_loops', eq_.training_loops, ... % Anzahl Trainingsdurchläufe + 'TRmu1', eq_.FFEmu, ... % mu für DD-Modus (1. Ordnung) + 'TRmu2', eq_.FFEmu, ... % mu für DD-Modus (2. Ordnung) + 'TRmu3', eq_.FFEmu, ... % mu für DD-Modus (3. Ordnung) + 'TRmuDFE', eq_.DFEmu, ... % mu für DFE-Modus im DD + 'dd_loops', 5, ... % Anzahl Durchläufe im DD-Modus + 'DDmu1', eq_.DDmu(1), ... % mu für DD-Modus (1. Ordnung) + 'DDmu2', eq_.DDmu(2), ... % mu für DD-Modus (2. Ordnung) + 'DDmu3', eq_.DDmu(3), ... % mu für DD-Modus (3. Ordnung) + 'DDmuDFE', eq_.DDmu(4), ... % mu für DFE-Modus im DD + 'MLSE_mode', 'viterbi', ... % Beispiel: MLSE-Modus als String + 'MLSE_trellis_states', jsonencode(mlse_.trellis_states), ... % Trellis-States, z.B. als JSON-String oder kommasepariert + 'comment', 'function: duobinary_target.m', ... % Zusätzliche Kommentare + 'config_hash', NaN ... + ); + +eq_package.resultsDBsignaling = resultsDBsignaling; +eq_package.equalizerConfigDBsignaling = equalizerConfigDBsignaling; - eq_package.ber = ber; - end \ No newline at end of file diff --git a/Functions/EQ_structures/duobinary_target.m b/Functions/EQ_structures/duobinary_target.m index c2d760e..8d66efb 100644 --- a/Functions/EQ_structures/duobinary_target.m +++ b/Functions/EQ_structures/duobinary_target.m @@ -1,77 +1,156 @@ function [eq_package] = duobinary_target(eq_, mlse_,M, rx_signal, tx_symbols, tx_bits, options) arguments - eq_ - mlse_ - M - rx_signal - tx_symbols - tx_bits + eq_ + mlse_ + M + rx_signal + tx_symbols + tx_bits options.precode_mode db_mode options.showAnalysis = 0; + options.eth_style_symbol_mapping = 0; + options.postFFE = []; end - - %Duobinary Targeting - [eq_signal, eq_noise] = eq_.process(rx_signal,Duobinary().encode(tx_symbols)); +%Duobinary Targeting +db_ref_sequence = Duobinary().encode(tx_symbols); +db_ref_constellation = unique(db_ref_sequence.signal); +[eq_signal, eq_noise] = eq_.process(rx_signal,db_ref_sequence); - % dir = [1,1]; - mlse_sig_sd = mlse_.process(eq_signal); +if ~isempty(options.postFFE) + [eq_signal,eq_noise] = options.postFFE.process(eq_signal,db_ref_sequence); +end - mlse_sig_hd = PAMmapper(M,0).quantize(mlse_sig_sd); - - % precoding to mitigate error propagation, most prominently used in - % combination with duobinary signaling to avoid catastrophic error - % behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling) +% dir = [1,1]; +mlse_sig_sd = mlse_.process(eq_signal); - % takes: - % -> eq_signal_hd: hard decision signal after eq - % -> tx_symbols: that where used as reference for eq +mlse_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd); - switch options.precode_mode - case db_mode.db_emulate +% precoding to mitigate error propagation, most prominently used in +% combination with duobinary signaling to avoid catastrophic error +% behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling) +switch options.precode_mode - mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M); - mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M); + case db_mode.no_db + % TX Data is not precoded: + + % A) Emulate diff precoding + mlse_sig_hd_precoded = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded,"M",M); tx_symbols_precoded = Duobinary().encode(tx_symbols); tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); - tx_bits = PAMmapper(M,0).demap(tx_symbols_precoded); + tx_bits_precoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded); - case db_mode.db_discard + rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_precoded); + [~,errors_db_diff_precoded,ber_db_diff_precoded,~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % normal dsp for precoded sequence == discard/omit/ignore precode - tx_bits = PAMmapper(M,0).demap(tx_symbols); + %B) Just determine BER + rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); + [bits_mlse,errors_mlse,ber_db,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - case db_mode.db_encoded + case db_mode.db_precoded - % normal DB encoded data (only for 10KM) + % Daten SIND TATSÄCHLICH precoded auf TX Seite: - case db_mode.db_precoded + % A) Decode at Rx if no DB targeting was applied (we are in VNLE or MLSE EQ structure here! + mlse_sig_hd_decoded = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd_decoded = Duobinary().decode(mlse_sig_hd_decoded,"M",M); + rx_bits_mlse_decoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_decoded); + [~,errors_db_diff_precoded,ber_db_diff_precoded,~] = calc_ber(rx_bits_mlse_decoded.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M); - mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M); + % B) Omit the Coding by comparing with demapped TX symbol sequence - end + tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols); + rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); + [bits_db,errors_db,ber_db,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % M = numel(unique(tx_symbols.signal)); - rx_bits = PAMmapper(M,0).demap(mlse_sig_hd); +end - [~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); +% M = numel(unique(tx_symbols.signal)); +rx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); - - eq_package.ber = ber; +[bits_db,errors_db,ber_db,errorIndice_db] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); +eq_package.ber = ber_db; - if options.showAnalysis - eq_noise = eq_noise - mean(eq_noise.signal); +resultsDBtgt = struct( ... + 'result_id', NaN, ... % + 'run_id', NaN, ... % Beispielhafte Run-ID + 'eqParam_id', NaN, ... % Beispielhafter Fremdschlüssel zur EqualizerParameters-Tabelle + 'date_of_processing', datetime('now'), ... % Aktuelles Datum und Uhrzeit + 'numBits', bits_db, ... % Beispiel: 1.000.000 Bits + 'numBitErr', errors_db, ... % Beispiel: 120 Bitfehler + 'BER_precoded', ber_db_diff_precoded, ... % BER = 120 / 1.000.000 + 'numBitErr_precoded', errors_db_diff_precoded, ... % Beispiel: 120 Bitfehler + 'BER', ber_db, ... % BER = 120 / 1.000.000 + 'SNR', [], ... % Beispielhafte SNR + 'SNR_level', jsonencode([]), ... % SNR-Level als JSON-codiertes Array + 'GMI', [], ... % Beispielhafter GMI-Wert + 'AIR', [], ... % Beispielhafter AIR-Wert + 'EVM', [], ... % Beispielhafte EVM + 'EVM_level', jsonencode([]), ... % EVM-Level als JSON-codiertes Array + 'Alpha', [] ... % Beispielhafter Alpha-Wert + ); -rx_signal.spectrum("normalizeTo0dB",1,"fignum",250); +if ~isempty(options.postFFE) + npostFFE = options.postFFE.order; +else + npostFFE = 0; +end - showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise'); +equalizerConfigDBtgt = struct( ... + 'eq_id', NaN, ... % Auto-Inkrement, wird in der DB gesetzt + 'equalizer_structure', int32(equalizer_structure.vnle_db_mlse), ... % Beispiel: 1 (z.B. für vnle) + 'M', M, ... % Ordnung der PAM-Konstellation + 'target_constellation', jsonencode(round(db_ref_constellation,5)), ... % Beispielhafter Target-String + 'db_target', 1, ... % 0 oder 1 + 'diff_precode', int32(options.precode_mode), ... % 0 oder 1 + 'postFFE', ~isempty(options.postFFE), ... % Beispielwert + 'NpostFFE', npostFFE, ... % Beispielwert + 'Ne1', eq_.Ne(1), ... % Feedforward Koeffizienten 1. Ordnung + 'Ne2', eq_.Ne(2), ... % Feedforward Koeffizienten 2. Ordnung + 'Ne3', eq_.Ne(3), ... % Feedforward Koeffizienten 3. Ordnung + 'Nb1', eq_.Nb(1), ... % Decision Feedback Koeffizienten 1. Ordnung + 'Nb2', eq_.Nb(2), ... % Decision Feedback Koeffizienten 2. Ordnung + 'Nb3', eq_.Nb(3), ... % Decision Feedback Koeffizienten 3. Ordnung + 'K', eq_.K, ... % Samples pro Symbol + 'DCmu', eq_.DCmu, ... % Anpassungsrate für DC-Tap + 'ideal_dfe', eq_.ideal_dfe, ... % Flag für ideal DFE (0 oder 1) + 'training_length', eq_.training_length, ... % Anzahl Trainingssymbole + 'training_loops', eq_.training_loops, ... % Anzahl Trainingsdurchläufe + 'TRmu1', eq_.FFEmu, ... % mu für DD-Modus (1. Ordnung) + 'TRmu2', eq_.FFEmu, ... % mu für DD-Modus (2. Ordnung) + 'TRmu3', eq_.FFEmu, ... % mu für DD-Modus (3. Ordnung) + 'TRmuDFE', eq_.DFEmu, ... % mu für DFE-Modus im DD + 'dd_loops', 5, ... % Anzahl Durchläufe im DD-Modus + 'DDmu1', eq_.DDmu(1), ... % mu für DD-Modus (1. Ordnung) + 'DDmu2', eq_.DDmu(2), ... % mu für DD-Modus (2. Ordnung) + 'DDmu3', eq_.DDmu(3), ... % mu für DD-Modus (3. Ordnung) + 'DDmuDFE', eq_.DDmu(4), ... % mu für DFE-Modus im DD + 'MLSE_mode', 'viterbi', ... % Beispiel: MLSE-Modus als String + 'MLSE_trellis_states', jsonencode(mlse_.trellis_states), ... % Trellis-States, z.B. als JSON-String oder kommasepariert + 'comment', 'function: duobinary_target.m', ... % Zusätzliche Kommentare + 'config_hash', NaN ... + ); - Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",250); - end +eq_package.resultsDBtgt = resultsDBtgt; +eq_package.equalizerConfigDBtgt = equalizerConfigDBtgt; + +if options.showAnalysis + eq_noise = eq_noise - mean(eq_noise.signal); + + rx_signal.spectrum("normalizeTo0dB",1,"fignum",250,"displayname","Rx Spectrum"); + + Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",250,"displayname","DB encoded reference"); + + showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise after Equalization'); + + fprintf('DB tgt BER: %.2e \n',ber); + + +end end \ No newline at end of file diff --git a/Functions/EQ_structures/vnle.m b/Functions/EQ_structures/vnle.m index 7746650..1447351 100644 --- a/Functions/EQ_structures/vnle.m +++ b/Functions/EQ_structures/vnle.m @@ -22,6 +22,8 @@ function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options) tx_bits options.precode_mode db_mode options.showAnalysis = 0 + options.eth_style = 0; + options.postFFE = []; end %FFE or VNLE @@ -31,6 +33,10 @@ function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options) end [eq_signal_sd,eq_noise] = eq_.process(rx_signal,tx_symbols); + if ~isempty(options.postFFE) + [eq_signal_sd,eq_noise] = options.postFFE.process(eq_signal_sd,tx_symbols); + end + eq_signal_hd = PAMmapper(M,0).quantize(eq_signal_sd); % precoding to mitigate error propagation, most prominently used in @@ -50,7 +56,7 @@ function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options) tx_symbols_precoded = Duobinary().encode(tx_symbols); tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); - tx_bits = PAMmapper(M,0).demap(tx_symbols_precoded); + tx_bits = PAMmapper(M,0,"eth_style",options.eth_style).demap(tx_symbols_precoded); case db_mode.db_discard @@ -68,26 +74,56 @@ function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options) end - rx_bits = PAMmapper(M,0).demap(eq_signal_hd); + rx_bits = PAMmapper(M,0,"eth_style",options.eth_style).demap(eq_signal_hd); [~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); [evm_total,evm_lvl] = calc_evm(eq_signal_sd,tx_symbols); [inf_rate] = calc_air(eq_signal_sd,tx_symbols,"skip_front",10000,"skip_end",10000); + eq_package.ber_vnle = ber; eq_package.evm_total = evm_total; eq_package.evm_lvl = evm_lvl; eq_package.inf_rate_vnle = inf_rate; - eq_package.snr = snr(eq_signal_sd.signal,eq_noise.signal); + eq_package.signal = eq_signal_sd; + + if options.showAnalysis - snr(eq_signal_sd.signal,eq_noise.signal); + fprintf('SNR: %d dB \n',snr(eq_signal_sd.signal,eq_noise.signal)); + + if M == 6 + logm = 2.5; + else + logm = log2(M); + end + + fprintf('NGMI: %.4f \n', inf_rate/logm); + fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl); fprintf('VNLE BER: %.2e \n',ber); + disp("%%%%%%%%%%%%%%%%%%%%%") + + % showEQcoefficients('n1',eq_.e,'n2',eq_.e2,'n3',eq_.e3,"displayname",'Coefficients'); + % + % if ~isempty(options.postFFE) + % showEQcoefficients('n1',options.postFFE.e,"displayname",'Coefficients'); + % end + % + % showEQNoisePSD(eq_noise); + % + % showEQfilter(eq_.e,eq_signal_sd.fs.*2) + + % noiselessness(tx_symbols,eq_noise,"displayname",'SNR after VNLE','fignum',301); + + % showLevelHistogram(eq_signal_sd,tx_symbols,"fignum",302); + + + end end \ No newline at end of file diff --git a/Functions/EQ_structures/vnle_postfilter_mlse.m b/Functions/EQ_structures/vnle_postfilter_mlse.m index ac2493c..35ea356 100644 --- a/Functions/EQ_structures/vnle_postfilter_mlse.m +++ b/Functions/EQ_structures/vnle_postfilter_mlse.m @@ -1,125 +1,267 @@ function [eq_package] = vnle_postfilter_mlse(eq_,pf_,mlse_,M,rx_signal,tx_symbols,tx_bits,options) - + arguments - eq_ - pf_ - mlse_ - M - rx_signal - tx_symbols - tx_bits + eq_ + pf_ + mlse_ + M + rx_signal + tx_symbols + tx_bits options.precode_mode db_mode options.showAnalysis = 0; + options.eth_style_symbol_mapping = 0; + options.postFFE = []; + options.database = []; end - %FFE or VNLE - [eq_signal_sd,eq_noise] = eq_.process(rx_signal,tx_symbols); +%FFE or VNLE +[eq_signal_sd,eq_noise] = eq_.process(rx_signal,tx_symbols); - eq_signal_hd = PAMmapper(M,0).quantize(eq_signal_sd); +if ~isempty(options.postFFE) + [eq_signal_sd,eq_noise] = options.postFFE.process(eq_signal_sd,tx_symbols); +end - mlse_sig_sd = pf_.process(eq_signal_sd,eq_noise); +eq_signal_hd = PAMmapper(M,0).quantize(eq_signal_sd); - mlse_.DIR = pf_.coefficients; - % [mlse_sig_hd,mlse_sig_sd] = mlse_.process(mlse_sig_sd,tx_symbols); - mlse_sig_sd = mlse_.process(mlse_sig_sd); +mlse_sig_sd = pf_.process(eq_signal_sd,eq_noise); - mlse_sig_hd = PAMmapper(M,0).quantize(mlse_sig_sd); +mlse_.DIR = pf_.coefficients; +% [mlse_sig_hd,mlse_sig_sd] = mlse_.process(mlse_sig_sd,tx_symbols); +mlse_sig_sd = mlse_.process(mlse_sig_sd); - % precoding to mitigate error propagation, most prominently used in - % combination with duobinary signaling to avoid catastrophic error - % behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling) +mlse_sig_hd = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).quantize(mlse_sig_sd); - % takes: - % -> M - % -> eq_signal_hd: hard decision signal after eq - % -> tx_symbols: that where used as reference for eq +% precoding to mitigate error propagation, most prominently used in +% combination with duobinary signaling to avoid catastrophic error +% behavior (see J.W.M. Bergmans, Digital Baseband Transmission and Recording -> partial response signaling) - switch options.precode_mode - case db_mode.db_emulate - % re - eq_signal_hd = Duobinary().encode(eq_signal_hd,"M",M); - eq_signal_hd = Duobinary().decode(eq_signal_hd,"M",M); +switch options.precode_mode - mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M); - mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M); + case db_mode.no_db + % TX Data is not precoded: + + % A) Emulate diff precoding + eq_signal_hd_precoded = Duobinary().encode(eq_signal_hd,"M",M); + eq_signal_hd_precoded = Duobinary().decode(eq_signal_hd_precoded,"M",M); + + mlse_sig_hd_precoded = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd_precoded = Duobinary().decode(mlse_sig_hd_precoded,"M",M); tx_symbols_precoded = Duobinary().encode(tx_symbols); tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); - tx_bits = PAMmapper(M,0).demap(tx_symbols_precoded); + tx_bits_precoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols_precoded); - case db_mode.db_discard + rx_bits_vnle = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(eq_signal_hd_precoded); + [~,errors_vnle_diff_precoded,ber_vnle_diff_precoded,~] = calc_ber(rx_bits_vnle.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % normal dsp for precoded sequence == discard/omit/ignore precode - tx_bits = PAMmapper(M,0).demap(tx_symbols); + rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_precoded); + [~,errors_mlse_diff_precoded,ber_mlse_diff_precoded,~] = calc_ber(rx_bits_mlse.signal,tx_bits_precoded.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - case db_mode.db_encoded + %B) Just determine BER + rx_bits_vnle = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(eq_signal_hd); + [bits_vnle,errors_vnle,ber_vnle,~] = calc_ber(rx_bits_vnle.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % normal DB encoded data (only for 10KM) + rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); + [bits_mlse,errors_mlse,ber_mlse,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - case db_mode.db_precoded + case db_mode.db_precoded - eq_signal_hd = Duobinary().encode(eq_signal_hd,"M",M); - eq_signal_hd = Duobinary().decode(eq_signal_hd,"M",M); + % Daten SIND TATSÄCHLICH precoded auf TX Seite: - mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M); - mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M); + % A) Decode at Rx if no DB targeting was applied (we are in VNLE or MLSE EQ structure here! + eq_signal_hd_decoded = Duobinary().encode(eq_signal_hd,"M",M); + eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded,"M",M); + rx_bits_vnle_decoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(eq_signal_hd_decoded); + [~,errors_vnle_diff_precoded,ber_vnle_diff_precoded,~] = calc_ber(rx_bits_vnle_decoded.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - end + mlse_sig_hd_decoded = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd_decoded = Duobinary().decode(mlse_sig_hd_decoded,"M",M); + rx_bits_mlse_decoded = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd_decoded); + [~,errors_mlse_diff_precoded,ber_mlse_diff_precoded,~] = calc_ber(rx_bits_mlse_decoded.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - % METRICS OF VNLE % - rx_bits_vnle = PAMmapper(M,0).demap(eq_signal_hd); - [~,~,ber_vnle,~] = calc_ber(rx_bits_vnle.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - % correct TUM implementation of AIR - [inf_rate_vnle] = calc_air(eq_signal_sd,tx_symbols,"skip_front",10000,"skip_end",10000); - [evm_vnle_total,evm_vnle_lvl] = calc_evm(eq_signal_sd,tx_symbols); + % B) Omit the Coding by comparing with demapped TX symbol sequence - % METRICS OF MLSE (HD-VITERBI) - rx_bits_mlse = PAMmapper(M,0).demap(mlse_sig_hd); - [~,~,ber_mlse,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + tx_bits = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(tx_symbols); + rx_bits_vnle = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(eq_signal_hd); + [bits_vnle,errors_vnle,ber_vnle,~] = calc_ber(rx_bits_vnle.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + rx_bits_mlse = PAMmapper(M,0,"eth_style",options.eth_style_symbol_mapping).demap(mlse_sig_hd); + [bits_mlse,errors_mlse,ber_mlse,~] = calc_ber(rx_bits_mlse.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + +end + +% METRICS OF VNLE SD Signal: +[snr_vnle,snr_vnle_lvl] = calc_snr(tx_symbols.signal,eq_noise.signal); +[gmi_vnle] = calc_air(eq_signal_sd,tx_symbols,"skip_front",10000,"skip_end",10000); +air_vnle = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi_vnle ./ log2(double(M)); +[evm_vnle_total,evm_vnle_lvl] = calc_evm(eq_signal_sd,tx_symbols); + +% METRICS OF MLSE (HD-VITERBI) +pf_.ncoeff = 1; +pf_.process(eq_signal_sd,eq_noise); +alpha = pf_.coefficients(2); + +eq_package.ber_mlse = ber_mlse; +eq_package.ber_vnle = ber_vnle; +eq_package.evm_vnle_total = evm_vnle_total; +eq_package.evm_vnle_lvl = evm_vnle_lvl; +eq_package.gmi = gmi_vnle; + +eq_package.eq = eq_; +eq_package.pf = pf_; +eq_package.mlse = mlse_; + +resultsVNLE = struct( ... + 'result_id', NaN, ... % + 'run_id', NaN, ... % Beispielhafte Run-ID + 'eqParam_id', NaN, ... % Beispielhafter Fremdschlüssel zur EqualizerParameters-Tabelle + 'date_of_processing', datetime('now'), ... % Aktuelles Datum und Uhrzeit + 'BER', ber_vnle, ... % BER = 120 / 1.000.000 + 'numBits', bits_vnle, ... % Beispiel: 1.000.000 Bits + 'numBitErr', errors_vnle, ... % Beispiel: 120 Bitfehler + 'BER_precoded', ber_vnle_diff_precoded, ... % BER = 120 / 1.000.000 + 'numBitErr_precoded', errors_vnle_diff_precoded, ... % Beispiel: 120 Bitfehler + 'SNR', snr_vnle, ... % Beispielhafte SNR + 'SNR_level', jsonencode(snr_vnle_lvl), ... % SNR-Level als JSON-codiertes Array + 'GMI', gmi_vnle, ... % Beispielhafter GMI-Wert + 'AIR', air_vnle, ... % Beispielhafter AIR-Wert + 'EVM', evm_vnle_total, ... % Beispielhafte EVM + 'EVM_level', jsonencode(evm_vnle_lvl), ... % EVM-Level als JSON-codiertes Array + 'Alpha', [] ... % Beispielhafter Alpha-Wert + ); + +if ~isempty(options.postFFE) + npostFFE = options.postFFE.order; +else + npostFFE = 0; +end + +equalizerConfigVNLE = struct( ... + 'eq_id', NaN, ... % Auto-Inkrement, wird in der DB gesetzt + 'equalizer_structure', int32(equalizer_structure.vnle), ... % Beispiel: 1 (z.B. für vnle) + 'M', M, ... % Ordnung der PAM-Konstellation + 'target_constellation', jsonencode(round(unique(tx_symbols.signal),5)), ... % Beispielhafter Target-String + 'db_target', 0, ... % 0 oder 1 + 'diff_precode', int32(options.precode_mode), ... % 0 oder 1 + 'postFFE', ~isempty(options.postFFE), ... % Beispielwert + 'NpostFFE', npostFFE, ... % Beispielwert + 'Ne1', eq_.Ne(1), ... % Feedforward Koeffizienten 1. Ordnung + 'Ne2', eq_.Ne(2), ... % Feedforward Koeffizienten 2. Ordnung + 'Ne3', eq_.Ne(3), ... % Feedforward Koeffizienten 3. Ordnung + 'Nb1', eq_.Nb(1), ... % Decision Feedback Koeffizienten 1. Ordnung + 'Nb2', eq_.Nb(2), ... % Decision Feedback Koeffizienten 2. Ordnung + 'Nb3', eq_.Nb(3), ... % Decision Feedback Koeffizienten 3. Ordnung + 'K', eq_.K, ... % Samples pro Symbol + 'DCmu', eq_.DCmu, ... % Anpassungsrate für DC-Tap + 'ideal_dfe', eq_.ideal_dfe, ... % Flag für ideal DFE (0 oder 1) + 'training_length', eq_.training_length, ... % Anzahl Trainingssymbole + 'training_loops', eq_.training_loops, ... % Anzahl Trainingsdurchläufe + 'TRmu1', eq_.FFEmu, ... % mu für DD-Modus (1. Ordnung) + 'TRmu2', eq_.FFEmu, ... % mu für DD-Modus (2. Ordnung) + 'TRmu3', eq_.FFEmu, ... % mu für DD-Modus (3. Ordnung) + 'TRmuDFE', eq_.DFEmu, ... % mu für DFE-Modus im DD + 'dd_loops', 5, ... % Anzahl Durchläufe im DD-Modus + 'DDmu1', eq_.DDmu(1), ... % mu für DD-Modus (1. Ordnung) + 'DDmu2', eq_.DDmu(2), ... % mu für DD-Modus (2. Ordnung) + 'DDmu3', eq_.DDmu(3), ... % mu für DD-Modus (3. Ordnung) + 'DDmuDFE', eq_.DDmu(4), ... % mu für DFE-Modus im DD + 'comment', 'function: vnle_postfilter_mlse', ... % Zusätzliche Kommentare + 'config_hash', NaN ... + ); - eq_package.ber_mlse = ber_mlse; - eq_package.ber_vnle = ber_vnle; - eq_package.evm_vnle_total = evm_vnle_total; - eq_package.evm_vnle_lvl = evm_vnle_lvl; - eq_package.air = inf_rate_vnle; - - eq_package.eq = eq_; - eq_package.pf = pf_; - eq_package.mlse = mlse_; +resultsMLSE = struct( ... + 'result_id', NaN, ... % + 'run_id', NaN, ... % Beispielhafte Run-ID + 'eqParam_id', NaN, ... % Beispielhafter Fremdschlüssel zur EqualizerParameters-Tabelle + 'date_of_processing', datetime('now'), ... % Aktuelles Datum und Uhrzeit + 'BER', ber_mlse, ... % BER = 120 / 1.000.000 + 'numBits', bits_mlse, ... % Beispiel: 1.000.000 Bits + 'numBitErr', errors_mlse, ... % Beispiel: 120 Bitfehler + 'BER_precoded', ber_mlse_diff_precoded, ... % BER = 120 / 1.000.000 + 'numBitErr_precoded', errors_mlse_diff_precoded, ... % Beispiel: 120 Bitfehler + 'SNR', [], ... % Beispielhafte SNR + 'SNR_level', jsonencode([]), ... % SNR-Level als JSON-codiertes Array + 'GMI', [], ... % Beispielhafter GMI-Wert + 'AIR', [], ... % Beispielhafter AIR-Wert + 'EVM', [], ... % Beispielhafte EVM + 'EVM_level', jsonencode([]), ... % EVM-Level als JSON-codiertes Array + 'Alpha', alpha, ... % Beispielhafter Alpha-Wert + 'MLSE_dir', jsonencode([mlse_.DIR])... + ); - if options.showAnalysis +equalizerConfigMLSE = struct( ... + 'eq_id', NaN, ... % Auto-Inkrement, wird in der DB gesetzt + 'equalizer_structure', int32(equalizer_structure.vnle_pf_mlse), ... % Beispiel: 1 (z.B. für vnle) + 'M', M, ... % Ordnung der PAM-Konstellation + 'target_constellation', jsonencode(round(unique(tx_symbols.signal),5)), ... % Beispielhafter Target-String + 'db_target', 0, ... % 0 oder 1 + 'diff_precode', int32(options.precode_mode), ... % 0 oder 1 + 'postFFE', ~isempty(options.postFFE), ... % Beispielwert + 'NpostFFE', npostFFE, ... % Beispielwert + 'Ne1', eq_.Ne(1), ... % Feedforward Koeffizienten 1. Ordnung + 'Ne2', eq_.Ne(2), ... % Feedforward Koeffizienten 2. Ordnung + 'Ne3', eq_.Ne(3), ... % Feedforward Koeffizienten 3. Ordnung + 'Nb1', eq_.Nb(1), ... % Decision Feedback Koeffizienten 1. Ordnung + 'Nb2', eq_.Nb(2), ... % Decision Feedback Koeffizienten 2. Ordnung + 'Nb3', eq_.Nb(3), ... % Decision Feedback Koeffizienten 3. Ordnung + 'K', eq_.K, ... % Samples pro Symbol + 'DCmu', eq_.DCmu, ... % Anpassungsrate für DC-Tap + 'ideal_dfe', eq_.ideal_dfe, ... % Flag für ideal DFE (0 oder 1) + 'training_length', eq_.training_length, ... % Anzahl Trainingssymbole + 'training_loops', eq_.training_loops, ... % Anzahl Trainingsdurchläufe + 'TRmu1', eq_.FFEmu, ... % mu für DD-Modus (1. Ordnung) + 'TRmu2', eq_.FFEmu, ... % mu für DD-Modus (2. Ordnung) + 'TRmu3', eq_.FFEmu, ... % mu für DD-Modus (3. Ordnung) + 'TRmuDFE', eq_.DFEmu, ... % mu für DFE-Modus im DD + 'dd_loops', 5, ... % Anzahl Durchläufe im DD-Modus + 'DDmu1', eq_.DDmu(1), ... % mu für DD-Modus (1. Ordnung) + 'DDmu2', eq_.DDmu(2), ... % mu für DD-Modus (2. Ordnung) + 'DDmu3', eq_.DDmu(3), ... % mu für DD-Modus (3. Ordnung) + 'DDmuDFE', eq_.DDmu(4), ... % mu für DFE-Modus im DD + 'MLSE_mode', 'viterbi', ... % Beispiel: MLSE-Modus als String + 'MLSE_trellis_states', jsonencode(mlse_.trellis_states), ... % Trellis-States, z.B. als JSON-String oder kommasepariert + 'comment', 'function: vnle_postfilter_mlse', ... % Zusätzliche Kommentare + 'config_hash', NaN ... + ); +eq_package.resultsVNLE = resultsVNLE; +eq_package.resultsMLSE = resultsMLSE; +eq_package.equalizerConfigVNLE = equalizerConfigVNLE; +eq_package.equalizerConfigMLSE = equalizerConfigMLSE; - % fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl); +% eq_package.vnle_out = eq_signal_sd; +if options.showAnalysis - fprintf('VNLE BER: %.2e \n',ber_vnle); + % fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl); - fprintf('MLSE BER: %.2e \n',ber_mlse); + fprintf('VNLE BER: %.2e \n',ber_vnle); - showEQNoisePSD(eq_noise,"fignum",336,"displayname",'VNLE+DFE','postfilter_taps',pf_.coefficients); - + fprintf('MLSE BER: %.2e \n',ber_mlse); - rx_signal.spectrum("normalizeTo0dB",1,"fignum",337,"displayname",'Rx Signal'); + showEQNoisePSD(eq_noise,"fignum",336,"displayname",'Residual Noise after VNLE','postfilter_taps',pf_.coefficients); - tx_symbols.spectrum("normalizeTo0dB",1,"fignum",337,'displayname','Tx Signal'); - showLevelHistogram(eq_signal_sd,tx_symbols) - % showLevelHistogram(mlse_sig_sd,tx_symbols) + rx_signal.spectrum("normalizeTo0dB",1,"fignum",337,"displayname",'Rx Signal'); - showEQcoefficients(eq_.e,eq_.e2,eq_.e3,"displayname",'Coefficients'); + tx_symbols.spectrum("normalizeTo0dB",1,"fignum",337,'displayname','Tx Signal'); - showEQNoiseSNR(tx_symbols,eq_noise,"displayname",'vnle snr','fignum',101); + showLevelHistogram(eq_signal_sd,tx_symbols) + % showLevelHistogram(mlse_sig_sd,tx_symbols) - %%% EQ SNR Spectrum %230 - %snr - snr_vnle = snr(tx_symbols.signal,eq_noise.signal); + showEQcoefficients('n1',eq_.e,'n2',eq_.e2,'n3',eq_.e3,"displayname",'Coefficients'); - % showErrorBurstCount(eq_signal_sd,tx_symbols) - + showEQNoiseSNR(tx_symbols,eq_noise,"displayname",'vnle snr','fignum',101); - end + %%% EQ SNR Spectrum %230 + %snr + + + % showErrorBurstCount(eq_signal_sd,tx_symbols) + + +end end \ No newline at end of file diff --git a/Functions/EQ_visuals/showEQNoiseSNR.m b/Functions/EQ_visuals/showEQNoiseSNR.m index 8b517fd..cac55eb 100644 --- a/Functions/EQ_visuals/showEQNoiseSNR.m +++ b/Functions/EQ_visuals/showEQNoiseSNR.m @@ -1,7 +1,7 @@ -function showEQNoiseSNR(tx_signal, rx_signal, options) +function showEQNoiseSNR(eq_signal, noise_signal, options) arguments - tx_signal - rx_signal + eq_signal + noise_signal options.fs_tx options.fs_rx options.fignum (1,1) double = NaN % Default to NaN if not provided @@ -16,13 +16,13 @@ end fig = figure(options.fignum); % Use the specified figure number end - if isa(tx_signal,'Signal') - options.fs_tx = tx_signal.fs; - tx_signal = tx_signal.signal; + if isa(eq_signal,'Signal') + options.fs_tx = eq_signal.fs; + eq_signal = eq_signal.signal; end - if isa(rx_signal,'Signal') - options.fs_rx = rx_signal.fs; - rx_signal = rx_signal.signal; + if isa(noise_signal,'Signal') + options.fs_rx = noise_signal.fs; + noise_signal = noise_signal.signal; end @@ -37,10 +37,10 @@ end % Ensure the figure is ready before calling spectrum title('SNR of received Signal') - fft_length = 2^(nextpow2(length(tx_signal))-7); + fft_length = 2^(nextpow2(length(eq_signal))-7); - [s_lin,w] = pwelch(tx_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_tx,"centered","psd","mean"); - [n_lin,w] = pwelch(rx_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_rx,"centered","psd","mean"); + [s_lin,w] = pwelch(eq_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_tx,"centered","psd","mean"); + [n_lin,w] = pwelch(noise_signal,hanning(fft_length),fft_length/2,fft_length,options.fs_rx,"centered","psd","mean"); w = w.*1e-9; snr_dbm = 10*log10(s_lin./n_lin); @@ -48,6 +48,7 @@ end % figure(231) hold on plot(w,snr_dbm,'DisplayName','SNR','LineWidth',0.5,'Color',options.color); + yline(mean(snr_dbm),'HandleVisibility','off','Color',options.color); xlabel("Frequency in GHz"); edgetick = 2^(nextpow2(options.fs_tx*1e-9)); diff --git a/Functions/EQ_visuals/showEQcoefficients.m b/Functions/EQ_visuals/showEQcoefficients.m index 3794212..cb747b8 100644 --- a/Functions/EQ_visuals/showEQcoefficients.m +++ b/Functions/EQ_visuals/showEQcoefficients.m @@ -1,58 +1,67 @@ -function showEQcoefficients(n1, n2, n3, options) - % Show filter coefficients as stem plot - % n1, n2, and n3 in different subplots - % Scale all y-axis to -1 and 1 +function showEQcoefficients(options) + % Show filter coefficients as stem plots. + % Only the provided coefficient arrays (n1, n2, n3) are shown, + % each in its own subplot. The y-axis is scaled to [-1, 1]. arguments - n1 - n2 - n3 - options.fignum (1,1) double = NaN % Default to NaN if not provided - options.displayname (1,:) char = '' % Default to an empty string if not provided + options.n1 = []; + options.n2 = []; + options.n3 = []; + options.fignum (1,1) double = NaN; % Default: create new figure if NaN + options.displayname (1,:) char = ''; % Default: empty string options.color = [0.2157, 0.4941, 0.7216]; - options.clf = 0; % Clear figure before plotting new + options.clf = 0; % Clear figure before plotting if set to 1 end - % Determine the figure number to use or create a new figure + % Determine the figure number to use or create a new one. if isnan(options.fignum) - fig = figure; % Create a new figure and get its handle + fig = figure; else - fig = figure(options.fignum); % Use the specified figure number + fig = figure(options.fignum); end if options.clf - clf(fig); % Clear the figure if requested + clf(fig); end hold on ax = gca; N = numel(ax.Children); - % Set up a colormap for consistent coloring + % Set up a colormap for consistent coloring. cmap = linspecer(8); options.color = cmap(mod(N, size(cmap, 1)) + 1, :); - % Create subplots for n1, n2, n3 - for i = 1:3 - subplot(3, 1, i); - switch i - case 1 - stem(n1, 'Color', options.color, 'LineWidth', 1,'Marker','.','MarkerSize',10); - title(sprintf('1st order Filter Coefficients: %d',numel(n1))); - case 2 - stem(n2, 'Color', options.color, 'LineWidth', 1,'Marker','.','MarkerSize',10); - title(sprintf('2nd order Filter Coefficients: %d',numel(n2))); - case 3 - stem(n3, 'Color', options.color, 'LineWidth', 1,'Marker','.','MarkerSize',10); - title(sprintf('3rd order Filter Coefficients: %d',numel(n3))); - end - ylim([-1, 1]); % Scale y-axis to -1 and 1 + % Build cell arrays for coefficients and their corresponding titles. + coeffs = {}; + titles = {}; + + if ~isempty(options.n1) + coeffs{end+1} = options.n1; + titles{end+1} = sprintf('1st order Filter Coefficients: %d', numel(options.n1)); + end + if ~isempty(options.n2) + coeffs{end+1} = options.n2; + titles{end+1} = sprintf('2nd order Filter Coefficients: %d', numel(options.n2)); + end + if ~isempty(options.n3) + coeffs{end+1} = options.n3; + titles{end+1} = sprintf('3rd order Filter Coefficients: %d', numel(options.n3)); + end + + numSubplots = numel(coeffs); + + for i = 1:numSubplots + subplot(1, numSubplots, i); + stem(coeffs{i}, 'Color', options.color, 'LineWidth', 1, ... + 'Marker', '.', 'MarkerSize', 10); + title(titles{i}); + ylim([-1, 1]); % Set y-axis limits to [-1, 1] grid on; - grid minor + grid minor; xlabel('Coefficient Index'); ylabel('Amplitude'); end - % Ensure the layout is tight for better visibility - sgtitle('Filter Coefficients'); % Overall title -end \ No newline at end of file + sgtitle('Filter Coefficients'); % Overall title for the figure +end diff --git a/Functions/EQ_visuals/showEQfilter.m b/Functions/EQ_visuals/showEQfilter.m new file mode 100644 index 0000000..3ef6732 --- /dev/null +++ b/Functions/EQ_visuals/showEQfilter.m @@ -0,0 +1,36 @@ + +function showEQfilter(coefficients,fs) + + % Assuming that obj.e contains the final FFE filter coefficients. + % Set the number of frequency points and sampling frequency. + nfft = 1024; % Number of frequency points + + % Compute the frequency response of the FFE filter. + [H, f] = freqz(coefficients, 1, nfft, fs); + + % Keep only the first half of the frequency response (up to the Nyquist frequency). + half_nfft = floor(nfft/2) + 1; + f = f(1:half_nfft); + H = H(1:half_nfft); + + % Plot the magnitude and phase responses. + figure; + + % Magnitude response (in dB) + subplot(2,1,1); + hold on + plot(f.*1e-9, 20*log10(abs(1./H))); + title('(Inverted) Magnitude Response of FFE Filter'); + xlabel('Frequency (Hz)'); + ylabel('Magnitude (dB)'); + grid on; + + % Phase response + subplot(2,1,2); + plot(f.*1e-9, unwrap(angle(H))); + title('Phase Response of FFE Filter'); + xlabel('Frequency (Hz)'); + ylabel('Phase'); + grid on; + +end \ No newline at end of file diff --git a/Functions/EQ_visuals/showLevelHistogram.m b/Functions/EQ_visuals/showLevelHistogram.m index 40f31f4..66e5d08 100644 --- a/Functions/EQ_visuals/showLevelHistogram.m +++ b/Functions/EQ_visuals/showLevelHistogram.m @@ -39,7 +39,9 @@ end intermediate = received_sd(lvl,:); cnt(lvl) = round(numel(intermediate(~isnan(intermediate)))./length(eq_signal),3).*100; hold on + warning off histogram(received_sd(lvl,:),1000,"EdgeAlpha",0,'DisplayName',['Lvl ',num2str(lvl),' | ',num2str(cnt(lvl)),' %'],'FaceColor',lvlcol(lvl,:),'Normalization','pdf'); + warning on end legend grid on diff --git a/Functions/Lab_helper/loadFreqResp.m b/Functions/Lab_helper/loadFreqResp.m new file mode 100644 index 0000000..108dc5d --- /dev/null +++ b/Functions/Lab_helper/loadFreqResp.m @@ -0,0 +1,9 @@ + +precomp_path = "D:\kiel_dsp\imdd_simulation\projects\HighSpeed_ETH"; +precomp_filename = "lab_high_speed_nachtschichtmzm"; +% precomp_filename = "lab_high_speed"; +% precomp_filename = "lab_high_speed"; + +freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',92e9); +freqresp.load('loadPath',precomp_path,'fileName',precomp_filename); +freqresp.plot(); diff --git a/Functions/Metrics/calc_air_plain.m b/Functions/Metrics/calc_air_plain.m new file mode 100644 index 0000000..757c9c8 --- /dev/null +++ b/Functions/Metrics/calc_air_plain.m @@ -0,0 +1,40 @@ +function [airs] = calc_air_plain(noisy_signal,reference_signal,options) + % Calculation of AIR acc. to J. Kozesnik, „Numerically Computing Achievable Rates of Memoryless Channels“, Francisco Javier Garcıa-Gomez, doi: 10.1007/978-94-009-9857-5. + % Implementation is not accessible, I mailed TUM to get the code... + + arguments(Input) + noisy_signal; + reference_signal; + options.skip_front = 0; + options.skip_end = 0; + options.returnErrorLocation = 0; + end + + options.skip_end = abs(options.skip_end); + options.skip_front = abs(options.skip_front); + + assert((options.skip_end+options.skip_front) 1) +% beta - Roll-off related parameter (0 < beta <= 1) +% +% Output: +% tau_error - Estimated timing error (in sample units) +% +% Implementation Notes: +% 1. The function computes an N-point FFT of the first N samples of rx. +% 2. It then determines an offset (Delta) defined as: +% offset = round((1 - 1/eta) * N) +% 3. To avoid index overflow, the summation is taken over indices k from 1 to +% floor(N/2) - offset. +% 4. The timing error is estimated as: +% tau_error = ( (1+beta)/(2*eta*N - 1) * sum(phase difference) ) / (2*pi) +% where the phase difference is (angle(R(k)) - angle(R(k+offset))) +% +% Make sure that the input signal rx contains at least N samples. + +% Check input length +if length(rx) < N + error('Input signal length must be at least N.'); +end + +% Compute the N-point FFT of the first N samples of rx +R = fft(rx(1:N), N); + +% Determine the offset based on the oversampling factor (eta) +offset = round((1 - 1/eta) * N); + +% Define the summation range to avoid index overflow +k_min = 1; +k_max = floor(N/2) - offset; +if k_max < k_min + error('Chosen parameters result in an empty summation range. Adjust N, eta, or beta.'); +end + +% Compute the sum of phase differences over the selected frequency bins +phase_diff_sum = 0; +for k = k_min:k_max + phase_k = angle(R(k)); + phase_k_offset = angle(R(k + offset)); + phase_diff_sum = phase_diff_sum + (phase_k - phase_k_offset); +end + +% Normalization factor as per the modified Godard algorithm +norm_factor = (1 + beta) / (2 * eta * N - 1); + +% Estimate the timing error in sample units +tau_error = (norm_factor * phase_diff_sum) / (2 * pi); +end diff --git a/Functions/beautifyBERplot.m b/Functions/beautifyBERplot.m index 690c4ba..e7a02ef 100644 --- a/Functions/beautifyBERplot.m +++ b/Functions/beautifyBERplot.m @@ -9,7 +9,9 @@ function beautifyBERplot() for i = 1:length(lines) lines(i).LineWidth = 1.3; % Thicker line width %lines(i).LineStyle = '-'; % Solid lines for simplicity - lines(i).Marker = markers{mod(i-1, num_markers) + 1}; % Assign markers cyclically + if string(lines(i).Marker) == "none" + lines(i).Marker = markers{mod(i-1, num_markers) + 1}; % Assign markers cyclically + end lines(i).MarkerSize = 4; % Marker size lines(i).MarkerFaceColor = 'auto'; % Use line color for marker face end @@ -23,7 +25,7 @@ function beautifyBERplot() % Set logarithmic scale for y-axis, but only if it makes sense. % If this is not always desired, you could condition this on the presence of lines or data. - % set(gca, 'YScale', 'log'); + set(gca, 'YScale', 'log'); % Customize grid and box appearance set(gca, 'Box', 'on', 'LineWidth', 0.8); % Thicker border diff --git a/Functions/moveit_wrapper_.m b/Functions/moveit_wrapper_.m new file mode 100644 index 0000000..d48daa0 --- /dev/null +++ b/Functions/moveit_wrapper_.m @@ -0,0 +1,40 @@ +function [data_out, state_out] = moveit_wrapper(func_name, data_in, para) + % Generic wrapper for legacy MATLAB functions using loop-based execution + % + % Usage: + % [data_out, state_out] = wrapper('quantize', data_in, para); + % + % Inputs: + % - func_name: Name of the function as a string (e.g., 'quantize') + % - data_in: Input signal or empty for initialization + % - para: Structure with parameters (optional) + % + % Outputs: + % - data_out: Processed output signal + % - state_out: Internal state of the function + + global loop; + + % Ensure the function exists + if ~exist(func_name, 'file') + error('Function "%s" does not exist.', func_name); + end + + % Step 1: Get default parameters if not provided + if nargin < 3 || isempty(para) + loop = 0; % Request parameter structure + [para, comment] = feval(func_name); + end + + % Initialize state + state = struct(); + + % Step 2: Initialize the function + loop = 0; + [~, state_] = feval(func_name, data_in, state, para); + + % Step 3: Process input signal + loop = 1; + [data_out, state_out] = feval(func_name, data_in, state_, para); + +end diff --git a/Functions/showTransferFunction.m b/Functions/showTransferFunction.m new file mode 100644 index 0000000..ef4635e --- /dev/null +++ b/Functions/showTransferFunction.m @@ -0,0 +1,22 @@ + +function showTransferFunction(coeffs,options) + +arguments + coeffs + options.fignum = 765 + options.DisplayName = '' + options.color = [1 1 0.1] +end +% Define the filter taps + + [H, w] = freqz(coeffs, 1, 1024, 1); + + figure(options.fignum); + hold on + plot(w, 10*log10(abs(H)), 'LineWidth', 2,'DisplayName',options.DisplayName,'Color',options.color); %todo + xlabel('Normalized Frequency'); + ylabel('Amplitude in dB'); + grid on; + ylim([-20,10]) + +end \ No newline at end of file diff --git a/projects/.DS_Store b/projects/.DS_Store index 5f9a177..b0fb4b3 100644 Binary files a/projects/.DS_Store and b/projects/.DS_Store differ diff --git a/projects/HighSpeedExperiment_2024/auswertung MPI/analyze_mpi.m b/projects/HighSpeedExperiment_2024/auswertung MPI/analyze_mpi.m new file mode 100644 index 0000000..eb7ca85 --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung MPI/analyze_mpi.m @@ -0,0 +1,92 @@ + +% basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +% db = DBHandler("pathToDB",[basePath,'silas_labor.db']); +if 1 + + uloops = struct; + uloops.precomp = [0,1]; + uloops.db_precode = [0,1]; + uloops.bitrate = [420].*1e9; %[300,330,360,390,420,450,480] [224,336,360,390,420,448] for MPI + % uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4]; + uloops.laser_wavelength = [1310]; + uloops.M = [4]; + uloops.link_length = [1]; % 1,2,3,5,6,8,10 + wh = DataStorage(uloops); + wh.addStorage("ber"); + + % wh = submit_simulations(wh,"parallel",0,"simulation_mode",0); + wh = submit_handle(@dsp_mpi,wh,"parallel",1); + +end + +a = wh_mpi_112gbd.getStoValue('ber',uloops.precomp, uloops.db_precode, uloops.bitrate(1) , uloops.laser_wavelength, uloops.M, uloops.link_length); + +%VNLE standalone +try + ber_vnle = cellfun(@(x) x.vnle_dfe_package{1,1}.ber_vnle, a); +end +%MLSE +try + ber_values_mlse = cellfun(@(s) cellfun(@(pkg) pkg.ber_mlse, s.vnle_pf_package, 'UniformOutput', false), a, 'UniformOutput', false); + ber_values_mlse = cell2mat(ber_values_mlse{1}); + +end +%DB +try + ber_values_db = cellfun(@(s) cellfun(@(pkg) pkg.ber, s.dbtgt_package, 'UniformOutput', false), a, 'UniformOutput', false); + ber_values_db = cell2mat(ber_values_db{1}); +end + +xax = [0 + 3 + 6 + 9 + 12 + 15 + 18 + 21 + 24 + 27 + 30 + 45]; +cols = cbrewer2('Set1',8); + +% Compute min, max, and mean for PAM 4 MLSE +min_mlse = min(ber_values_mlse, [], 2); +max_mlse = max(ber_values_mlse, [], 2); +mean_mlse = mean(ber_values_mlse, 2); +err_lower_mlse = mean_mlse - min_mlse; +err_upper_mlse = max_mlse - mean_mlse; +err_mlse = [err_lower_mlse, err_upper_mlse]; + +% Compute min, max, and mean for PAM 4 DB tgt. +min_db = min(ber_values_db, [], 2); +max_db = max(ber_values_db, [], 2); +mean_db = mean(ber_values_db, 2); +err_lower_db = mean_db - min_db; +err_upper_db = max_db - mean_db; +err_db = [err_lower_db, err_upper_db]; + +figure(1) +hold on +title('MPI'); + +% Plot the MLSE curve with bounded error using boundedline +[hl_mlse, hp_mlse] = boundedline(xax, mean_mlse, err_mlse,'Color', cols(1,:)); +plot(xax,ber_values_mlse,'DisplayName','PAM 4 MLSE','Color',cols(1,:),'LineStyle','-','HandleVisibility','on','Marker','none','LineWidth',0.2); + +% Plot the DB tgt. curve with bounded error using boundedline +[hl_db, hp_db] = boundedline(xax, mean_db, err_db, 'Color', cols(2,:)); +plot(xax,ber_values_db,'DisplayName','PAM 4 MLSE','Color',cols(2,:),'LineStyle','-','HandleVisibility','on','Marker','none','LineWidth',0.2); + + +% Format the plot +xticks(xax); +set(gca, 'YScale', 'log'); +ylim([5e-5 0.4]); +xlim([min(xax) max(xax)]); +yline([4.85e-3, 2e-2], 'HandleVisibility', 'off'); +legend +% beautifyBERplot() +xlabel('Interference Attenuation'); +ylabel('BER'); diff --git a/projects/HighSpeedExperiment_2024/auswertung MPI/dsp_mpi.m b/projects/HighSpeedExperiment_2024/auswertung MPI/dsp_mpi.m new file mode 100644 index 0000000..8c6b5be --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung MPI/dsp_mpi.m @@ -0,0 +1,333 @@ +function [output] = dsp_mpi(varargin) + +simulation_mode = 0; + +%%% Change folder +curFolder = pwd; +funcFolder=fileparts(mfilename('fullpath')); +if ~isempty(funcFolder) + cd(funcFolder); +end + +%%% Run parameters +% TX +M = 4; +fsym = 180e9; + +apply_pulsef = 1; +fdac = 256e9; +fadc = 256e9; +random_key = 1; + +interference_attenuation = 0; +is_mpi = 1; + +precomp = 0; +db_precode = 0; +db_encode = 0; + +rcalpha = 0.05; +kover = 16; +vbias_rel = 0.5; +u_pi = 2.9; +vbias = -vbias_rel*u_pi; +laser_wavelength = 1293; +laser_linewidth = 0; +tx_bw_nyquist = 0.8; + +% Channel +link_length = 1; + +% RX +rop = -5; +rx_bw_nyquist = 0.8; + +vnle_order1 = 50; +vnle_order2 = 5; +vnle_order3 = 5; + +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; +dfe_order = [0 0 0]; + +pf_ncoeffs = 1; + +alpha = 0; + +len_tr = 4096*2; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; +mu_dc = 0; + +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; +mu_dfe = 0.0004; + + +dfe_ = sum(dfe_order)>0; + +doub_mode = db_mode.no_db; + +%%% change specific parameter if given in varargin +% Parse optional input arguments +if ~isempty(varargin) + var_s = varargin{1}; + if isstruct(var_s) + fields = fieldnames(var_s); + for i = 1:numel(fields) + if isnumeric(fields{i}) + eval([fields{i}, ' = ', num2str( var_s.(fields{i}) ), ';']); + fprintf("%s <-- %.2f \n", fields{i}, var_s.(fields{i})); + else + eval([fields{i}, ' = ', 'var_s.(fields{',num2str(i),'})' , ';']); + end + + end + else + error('Optional variables should be passed as a struct.'); + end +end + + + +if doub_mode ~= db_mode.db_encoded + if precomp == 0 && db_precode == 1 + doub_mode = db_mode.db_precoded; + + db_precode = 1; % preceded data (in my measurement set, this corresponds to low precomp too!) + discard_precode = 0; % + emulate_precode = 0; + legendentry = 'low precomp; precoded'; + disp('low precomp; precoded') + + elseif precomp == 1 && db_precode == 1 + doub_mode = db_mode.db_emulate; + + db_precode = 0; % preceded data (in my measurement set, this corresponds to low precomp too!) + discard_precode = 0; % + emulate_precode = 1; + legendentry = 'high precomp; precoded'; + disp('high precomp; precoded') + + elseif precomp == 0 && db_precode == 0 + doub_mode = db_mode.db_discard; + + db_precode = 1; % preceded data (in my measurement set, this corresponds to low precomp too!) + discard_precode = 1; % + emulate_precode = 0; + legendentry = 'no precomp; not precoded'; + disp('no precomp; not precoded') + + elseif precomp == 1 && db_precode == 0 + doub_mode = db_mode.no_db; + + db_precode = 0; % preceded data (in my measurement set, this corresponds to low precomp too!) + discard_precode = 0; % + emulate_precode = 0; + legendentry = 'high precomp; not precoded'; + disp('high precomp; not precoded') + + end + +else + +end + +fsym_ = floor( bitrate*1e-9./log2(M) ).*1e9; + +if fsym_ ~= fsym + fsym = fsym_; + % fprintf('Adapted symbolrate to %d GBd, to match provided bitrate of %d GBit/s using PAM %d \n',fsym.*1e-9,bitrate.*1e-9, M); +end +f_nyquist = fsym/2; + + +basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +database = DBHandler("pathToDB",[basePath,'silas_labor.db']); + +useGui = 0; +% db = DBHandler("pathToDB",[basePath,'silas_labor.db']); +filterParams = database.tables; +% filterParams.Runs.run_id = 2958; % no db +% filterParams.Runs.run_id = 2937; % no db +filterParams.Configurations = struct( ... + 'bitrate', bitrate, ... + 'db_mode', db_precode+db_encode, ... + 'fiber_length', link_length, ... + 'interference_attenuation', [], ... + 'interference_path_length', [], ... + 'is_mpi', is_mpi, ... + 'pam_level', M, ... + 'precomp_amp', [], ... + 'rop_attenuation', 0, ... + 'symbolrate', [], ... + 'v_awg', [], ... + 'v_bias', [], ... + 'wavelength', laser_wavelength ... + ); + +selectedFields = {'Runs.run_id','Runs.tx_bits_path', 'Runs.tx_symbols_path', 'Runs.rx_sync_path','Runs.rx_raw_path',... + 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','Measurements.power_rop','Configurations.v_bias',... + 'Configurations.interference_attenuation'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); +[~, uniqueIdx] = unique(dataTable.run_id); % Get unique run_id indices +dataTable = dataTable(uniqueIdx,:); % Extract unique configurations for each run_id +fprintf('Found %d entries for requested Configuration. IDs are: %s \n \n',size(dataTable,1),jsonencode(dataTable.run_id(1:min(size(dataTable,1),100)))); + + +output = struct(); +vnle_pf_package = {}; +vnle_dfe_package = {}; +dbtgt_package = {}; + +disp(num2str(bitrate)) + +for iatt = 1:numel(dataTable.interference_attenuation) + + current_run_id = dataTable.run_id(iatt); + + Tx_bits = load([basePath, char(dataTable.tx_bits_path(iatt))]); + Tx_bits = Tx_bits.Bits; + Symbols_mapped = PAMmapper(M,0).map(Tx_bits); + Symbols_mapped.fs = fsym; + + Symbols = load([basePath, char(dataTable.tx_symbols_path(iatt))]); + Symbols = Symbols.Symbols; + + Scpe_load = load([basePath, char(dataTable.rx_sync_path(iatt))]); + Scpe_cell = Scpe_load.S; + [~,~,found]=Scpe_cell{2}.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",1); + + if ~found + Raw_signal = load([basePath, char(dataTable.rx_raw_path(1))]); + Raw_signal = Raw_signal.Scpe_sig_raw; + [~,Scpe_cell,found] =Raw_signal.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",1); + end + + if ~found + if length(Symbols_mapped.signal) == sum(Symbols_mapped.signal == Symbols.signal) + warning('Could not synchronize the received signal with the stored symbols!') + else + [~,Scpe_cell,found] =Raw_signal.tsynch("reference",Symbols_mapped,"fs_ref",fsym,"debug_plots",0); + end + if ~found + warning('Could not synchronize the received signal with the stored symbols!') + end + end + + fsym = Symbols.fs; + + if db_precode + Symbols_precoded = Symbols; + end + + + + proc_occ = min(15,length(Scpe_cell)); + for occ = 1:proc_occ + + Scpe_sig = Scpe_cell{occ}; + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",0); + + Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); + + Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); + + + %%% EQUALIZING + + + % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",1,"mu_dc",0.05); + % eq_mlse = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0); + % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",512,"mu_dc",0.05); + + mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; + vnle_order=[vnle_order1,vnle_order2,vnle_order3]; + + % %%%%% VNLE + DFE %%%% + if 0 + + eq_vnle_dfe = EQ("Ne",vnle_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + eq_2 = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); + + [result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",1,"postFFE",[]); + vnle_dfe_package{iatt,occ} = result; + + end + + %%%%% VNLE + PF + MLSE %%%% + if 1 + + try + % len_tr = length(Symbols)-1000; + eq_vnle_ = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + % eq_vnle_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",vnle_order,"sps",2,"decide",0); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + eq_2 = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); + + [result] = vnle_postfilter_mlse(eq_vnle_,pf_,mlse_,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,'showAnalysis',0,"postFFE",[]); + vnle_pf_package{iatt,occ} = result; + + database.addProcessingResult(current_run_id,result.resultsMLSE, result.equalizerConfigMLSE); + database.addProcessingResult(current_run_id,result.resultsVNLE, result.equalizerConfigVNLE); + catch + warning(['VNLE+MLSE fail: run id: ', num2str(current_run_id)],' occ:', num2str(occ), ' iatten: ',num2str(iatt)) + end + + end + + %%%%% Duobinary Targeting %%%% + if 1 + + try + mlse_db = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + eq_db = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + eq_2 = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); + + [result] = duobinary_target(eq_db, mlse_db, M, Scpe_sig, Symbols, Tx_bits, "precode_mode", doub_mode,'showAnalysis',0,"postFFE",[]); + dbtgt_package{iatt,occ} = result; + + database.addProcessingResult(current_run_id,result.resultsDBtgt, result.equalizerConfigDBtgt); + + catch + warning(['VNLE DB+MLSE fail: run id: ', num2str(current_run_id)],' occ:', num2str(occ), ' iatten: ',num2str(iatt)) + end + + end + + %%%%%% %db signaling => db encoded %%%%% + if 0 + mlse_db_enc = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + eq_db_enc = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + [result] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Tx_bits); + dbenc_package{iatt,occ} = result; + end + + + % autoArrangeFigures; + disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') + fprintf('\n') + + + end + + + + if ~isempty(curFolder) + cd(curFolder); + end + +end + +output.dataTable = dataTable; +output.vnle_dfe_package = vnle_dfe_package; +output.vnle_pf_package = vnle_pf_package; +output.dbtgt_package = dbtgt_package; \ No newline at end of file diff --git a/projects/HighSpeedExperiment_2024/auswertung MPI/dsp_mpi_database.m b/projects/HighSpeedExperiment_2024/auswertung MPI/dsp_mpi_database.m new file mode 100644 index 0000000..c34aa7d --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung MPI/dsp_mpi_database.m @@ -0,0 +1,61 @@ + +% 1) Find RUN ID's + +basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +database = DBHandler("pathToDB",[basePath,'silas_labor.db']); + +filterParams = database.tables; +filterParams.Configurations = struct( ... + 'bitrate', [], ... %[224,336,360,390,420,448] + 'db_mode', [], ... + 'fiber_length', [], ... + 'interference_attenuation',[], ... + 'is_mpi', 0, ... + 'pam_level', [], ... + 'rop_attenuation', 0 ... + ); + +selectedFields = {'Runs.run_id',... + 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength',... + 'Configurations.precomp_amp','Measurements.power_rop','Measurements.power_pd_in','Configurations.v_bias','Configurations.is_mpi',... + 'Configurations.interference_attenuation','Configurations.rop_attenuation'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); + +% only the rows without BER so far +% dataTable = dataTable(dataTable.BER == 0,:); +dataTable = dataTable((dataTable.fiber_length ~= 1),:); + +% Ensure a parallel pool is running +pool = gcp('nocreate'); +if isempty(pool) + pool = parpool; + % stop all forgotten or unfetched jobs from queue +elseif ~isempty(pool.FevalQueue.QueuedFutures) || ~isempty(pool.FevalQueue.RunningFutures) + oldq = length(pool.FevalQueue.QueuedFutures) + length(pool.FevalQueue.RunningFutures); + pool.FevalQueue.cancelAll + fprintf('Canceled %d unfetched jobs from old queue.', oldq); +end + +% Number of tasks to submit (one per run_id) +nTasks = height(dataTable); +futures = parallel.FevalFuture.empty(); + +% Submit each DSP run as a parallel task using parfeval +for i = 1:nTasks + % Extract the run_id (other parameters could be passed if needed) + runID = dataTable.run_id(i); + % Submit the function call to dsp_run_id (assuming it returns no output, hence 0 outputs) + futures(i) = parfeval(pool, @dsp_run_id, 0, runID, "max_occurences", 15, "append_to_db", 1); +end + +% Set up a waitbar to monitor progress +h = waitbar(0, 'Processing DSP runs...'); +while ~all(strcmp({futures.State}, 'finished')) + finishedCount = sum(strcmp({futures.State}, 'finished')); + waitbar(finishedCount / nTasks, h); + pause(0.1); +end +delete(h); + +fprintf('All DSP runs processed.\n'); diff --git a/projects/HighSpeedExperiment_2024/auswertung MPI/dsp_run_id.m b/projects/HighSpeedExperiment_2024/auswertung MPI/dsp_run_id.m new file mode 100644 index 0000000..ebb7dfc --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung MPI/dsp_run_id.m @@ -0,0 +1,194 @@ +function [output] = dsp_run_id(run_id,options) + +arguments + run_id + options.append_to_db = 0; + options.max_occurences = 4; + options.parameters = struct(); +end + +basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +database = DBHandler("pathToDB",[basePath,'silas_labor.db']); + +filterParams = database.tables; +filterParams.Configurations = struct('run_id', run_id); + +selectedFields = {'Runs.run_id','Runs.tx_bits_path', 'Runs.tx_symbols_path', 'Runs.rx_sync_path','Runs.rx_raw_path',... + 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','Measurements.power_rop','Configurations.v_bias',... + 'Configurations.interference_attenuation'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); +[~, uniqueIdx] = unique(dataTable.run_id); % Get unique run_id indices +dataTable = dataTable(uniqueIdx,:); % Extract unique configurations for each run_id + +fsym = dataTable.symbolrate; +M = double(dataTable.pam_level); +duob_mode = db_mode(dataTable.db_mode); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +len_tr = 4096*2; + +vnle_order1 = 50; +vnle_order2 = 5; +vnle_order3 = 5; + +dfe_order = [0 0 0]; + +pf_ncoeffs = 1; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dfe = 0.0004; +mu_dc = 0.00; + +mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; + +% Overwrite default parameters if given in options.parameters +paramStruct = options.parameters; +if ~isempty(paramStruct) + paramNames = fieldnames(paramStruct); + for i = 1:numel(paramNames) + thisName = paramNames{i}; + thisValue = paramStruct.(thisName); + eval([thisName ' = thisValue;']); + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +eq_ = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); +pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); +mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); +mlse_db_ = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); +eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); + +output = struct(); +vnle_pf_package = {}; +vnle_dfe_package = {}; +dbtgt_package = {}; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +Tx_bits = load([basePath, char(dataTable.tx_bits_path)]); +Tx_bits = Tx_bits.Bits; +Symbols_mapped = PAMmapper(M,0).map(Tx_bits); +Symbols_mapped.fs = dataTable.symbolrate; + +Symbols = load([basePath, char(dataTable.tx_symbols_path)]); +Symbols = Symbols.Symbols; + +Scpe_load = load([basePath, char(dataTable.rx_sync_path)]); +Scpe_cell = Scpe_load.S; +[~,~,found]=Scpe_cell{2}.tsynch("reference",Symbols,"fs_ref",dataTable.symbolrate,"debug_plots",0); + +if ~found + Raw_signal = load([basePath, char(dataTable.rx_raw_path(1))]); + Raw_signal = Raw_signal.Scpe_sig_raw; + [~,Scpe_cell,found] =Raw_signal.tsynch("reference",Symbols,"fs_ref",dataTable.symbolrate,"debug_plots",0); +end + +if ~found + if length(Symbols_mapped.signal) == sum(Symbols_mapped.signal == Symbols.signal) + warning('Could not synchronize the received signal with the stored symbols!') + else + [~,Scpe_cell,found] =Raw_signal.tsynch("reference",Symbols_mapped,"fs_ref",dataTable.symbolrate,"debug_plots",0); + end + if ~found + warning('Could not synchronize the received signal with the stored symbols!') + end +end + +proc_occ = min(options.max_occurences,length(Scpe_cell)); +for occ = 1:proc_occ + + Scpe_sig = Scpe_cell{occ}; + + %%%%%% Sample to 2x fsym %%%%%% + Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); + + %%%%%% Sync Rx signal with reference %%%%%% + [Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym,"debug_plots",0); + + Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); + + Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); + + % Scpe_sig.plot("displayname",'Scope Signal','fignum',11); + % Scpe_sig.spectrum("displayname",'Raw Signal','fignum',20); + + if duob_mode ~= db_mode.db_encoded + + % %%%%% VNLE + DFE %%%% + if 0 + + eq_vnle_dfe = EQ("Ne",vnle_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); + + [result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",duob_mode,"showAnalysis",1,"postFFE",[]); + vnle_dfe_package{occ} = result; + + end + + %%%%% VNLE + PF + MLSE %%%% + if 1 + + [result] = vnle_postfilter_mlse(eq_,pf_,mlse_,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",duob_mode,'showAnalysis',0,"postFFE",[],"eth_style_symbol_mapping",0); + vnle_pf_package{occ} = result; + + if options.append_to_db + database.addProcessingResult(run_id,result.resultsMLSE, result.equalizerConfigMLSE); + database.addProcessingResult(run_id,result.resultsVNLE, result.equalizerConfigVNLE); + end + end + + %%%%% Duobinary Targeting %%%% + if 1 + [result] = duobinary_target(eq_, mlse_db_, M, Scpe_sig, Symbols, Tx_bits, "precode_mode", duob_mode,'showAnalysis',0,"postFFE",[]); + dbtgt_package{occ} = result; + + if options.append_to_db + database.addProcessingResult(run_id, result.resultsDBtgt, result.equalizerConfigDBtgt); + end + end + + fprintf("BER VNLE: %.2e | %.2e; BER MLSE: %.2e | %.2e; BER DB tgt: %.2e | %.2e \n",vnle_pf_package{occ}.resultsVNLE.BER,vnle_pf_package{occ}.resultsVNLE.BER_precoded ,vnle_pf_package{occ}.resultsMLSE.BER,vnle_pf_package{occ}.resultsMLSE.BER_precoded,dbtgt_package{occ}.resultsDBtgt.BER,dbtgt_package{occ}.resultsDBtgt.BER_precoded) + % fprintf("BER VNLE: %.2e | %.2e; BER MLSE: %.2e | %.2e \n",vnle_pf_package{occ}.resultsVNLE.BER,vnle_pf_package{occ}.resultsVNLE.BER_precoded ,vnle_pf_package{occ}.resultsMLSE.BER,vnle_pf_package{occ}.resultsMLSE.BER_precoded); + + + else + + %%%%%% %db signaling => db encoded %%%%% + if 1 + mlse_db_enc = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + eq_db_enc = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); + [result] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Tx_bits); + dbenc_package{occ} = result; + + if options.append_to_db + database.addProcessingResult(run_id, result.resultsDBsignaling, result.equalizerConfigDBsignaling); + end + end + + fprintf("BER DB: %.2e \n",dbenc_package{occ}.resultsDBsignaling.BER); + + + end + + + % autoArrangeFigures; + disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') + fprintf('\n') + + +end + + +output.dataTable = dataTable; +output.vnle_dfe_package = vnle_dfe_package; +output.vnle_pf_package = vnle_pf_package; +output.dbtgt_package = dbtgt_package; + +end \ No newline at end of file diff --git a/projects/HighSpeedExperiment_2024/auswertung MPI/plot_from_database.m b/projects/HighSpeedExperiment_2024/auswertung MPI/plot_from_database.m new file mode 100644 index 0000000..3093427 --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung MPI/plot_from_database.m @@ -0,0 +1,150 @@ + + +basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +database = DBHandler("pathToDB",[basePath,'silas_labor.db']); + +filterParams = database.tables; +filterParams.Configurations = struct( ... + 'bitrate', [], ... %[224,336,360,390,420,448] + 'db_mode', int32(db_mode.no_db), ... + 'fiber_length', 1, ... + 'interference_attenuation', [], ... + 'interference_path_length', [], ... + 'is_mpi', 1, ... + 'pam_level', 4, ... + 'rop_attenuation', 0, ... + 'wavelength', 1310 ... + ); + +% filterParams.EqualizerParameters.diff_precode = int32(db_mode.no_db); + filterParams.EqualizerParameters.equalizer_structure = int32(equalizer_structure.vnle); +% filterParams.EqualizerParameters.DCmu = 0.005; + +selectedFields = {'Configurations.run_id' 'Runs.rx_raw_path' 'Configurations.bitrate' 'Configurations.symbolrate' 'Configurations.pam_level' 'Configurations.db_mode' 'Configurations.rop_attenuation' 'Configurations.is_mpi' 'Configurations.interference_attenuation' 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'Measurements.power_pd_in' 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.SNR' 'Results.GMI' 'Results.Alpha'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); + +fixedVars = {'run_id','eq_id','bitrate'}; +dataTableGrpd = groupIt(fixedVars,dataTable); + + +% Create a new figure +figure(3); +hold on +unique_rates = unique(dataTable.bitrate); +cols = linspecer(8); +for i = 1:numel(unique_rates) + + % Plot BER vs. interference_attenuation + % plot(dataTableGrpd.power_mpi_signal(dataTableGrpd.bitrate==unique_rates(i),:)-dataTableGrpd.power_mpi_interference(dataTableGrpd.bitrate==unique_rates(i),:), dataTableGrpd.BER(dataTableGrpd.bitrate==unique_rates(i),:), '-', 'LineWidth', 0.5,'Color',cols(i,:)); + + if filterParams.Configurations.is_mpi + sir = dataTable.power_mpi_signal(dataTable.bitrate==unique_rates(i),:)-dataTable.power_mpi_interference(dataTable.bitrate==unique_rates(i),:); + ber = dataTable.BER(dataTable.bitrate==unique_rates(i),:); + sc=scatter(dataTable.power_mpi_signal(dataTable.bitrate==unique_rates(i),:)-dataTable.power_mpi_interference(dataTable.bitrate==unique_rates(i),:), dataTable.BER(dataTable.bitrate==unique_rates(i),:), 'LineWidth', 0.5,'Marker','.','MarkerEdgeColor',cols(i+1,:)); + pair_one = {'Run ID', dataTable.run_id(dataTable.bitrate==unique_rates(i),:)}; + pair_two = {'Rate', dataTable.bitrate(dataTable.bitrate==unique_rates(i),:)}; + addDatatips(sc, pair_one, pair_two); + else + sc=scatter(62*ones(size( dataTableGrpd.BER(dataTableGrpd.bitrate==unique_rates(i),:))), dataTableGrpd.BER(dataTableGrpd.bitrate==unique_rates(i),:), 'LineWidth', 0.5,'Marker','o','MarkerEdgeColor',cols(i,:),'MarkerFaceColor',cols(i,:)); + pair_one = {'Run ID', dataTableGrpd.run_id(dataTableGrpd.bitrate==unique_rates(i),:)}; + pair_two = {'Rate', dataTableGrpd.bitrate(dataTableGrpd.bitrate==unique_rates(i),:)}; + addDatatips(sc, pair_one, pair_two); + end + + + +end + +% Label the axes and add a title +xlabel('SIR in dB'); +ylabel('BER'); +title('BER vs. Signal to Interference Ratio'); +yline(3.8e-3,'LineWidth',1,'LineStyle','--','HandleVisibility','off'); +% Enable grid for better readability +grid on; + +beautifyBERplot; +ylim([1e-4 0.5]); + + + +function resultTable = groupIt(fixedVars,dataTable) + + % Group by run_id and eq_id (adjust grouping keys as needed) + + [G, groupKeys] = findgroups(dataTable(:, fixedVars)); + + % Preallocate a cell array for aggregated data. + varNames = dataTable.Properties.VariableNames; + nVars = numel(varNames); + aggData = cell(height(groupKeys), nVars); + groupCount = zeros(height(groupKeys), 1); % To store the size of each group + + % Loop over each group. + for i = 1:height(groupKeys) + idx = (G == i); % Logical index for group i + groupCount(i) = sum(idx); % Count number of rows in this group + % For each variable in the table: + for j = 1:nVars + colData = dataTable.(varNames{j}); + if isnumeric(colData) + % For numeric data, compute the mean. + aggData{i, j} = min(colData(idx)); + else + % For non-numeric data, take the first entry. + if iscell(colData) + aggData{i, j} = colData{find(idx, 1)}; + else + aggData{i, j} = colData(find(idx, 1)); + end + end + end + end + + % Convert the aggregated cell array into a table. + resultTable = cell2table(aggData, 'VariableNames', varNames); + + % Append the group count as a new column. + resultTable.nRows = groupCount; + +end + + +function addDatatips(sc, varargin) +% addDatatips Adds custom data tip rows to a scatter plot. +% +% addDatatips(sc, pair1, pair2, ...) adds one or more custom rows to the +% data tip display of the scatter plot identified by sc. +% +% Each pair should be provided as a 1x2 cell array: {label, value}. +% The value can be a scalar or a vector. If a vector is provided, its length +% must match the number of scatter plot points. +% +% Example: +% sc = scatter(x, y, 'LineWidth', 1.5, 'Marker', 'o'); +% pair_one = {'Attenuation', attenuationVector}; +% addDatatips(sc, pair_one); + + numPoints = numel(sc.XData); + + for k = 1:length(varargin) + pair = varargin{k}; + + if ~iscell(pair) || numel(pair) ~= 2 + error('Each pair must be a 1x2 cell array: {label, value}.'); + end + + label = pair{1}; + value = pair{2}; + + % If value is a vector, ensure its length is either 1 or equal to the number of scatter points. + if isvector(value) && numel(value) ~= 1 && numel(value) ~= numPoints + error('The vector for "%s" must be a scalar or have %d elements matching the scatter data points.', label, numPoints); + end + + % Create a new data tip row using the provided label and vector. + newRow = dataTipTextRow(label, value); + sc.DataTipTemplate.DataTipRows(end+1) = newRow; + end +end diff --git a/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_burg.m b/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_burg.m new file mode 100644 index 0000000..e624b92 --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_burg.m @@ -0,0 +1,68 @@ +% basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +% db = DBHandler("pathToDB",[basePath,'silas_labor.db']); +if 0 + + uloops = struct; + uloops.precomp = [0]; + uloops.db_precode = [0]; + uloops.bitrate = [330,390,450].*1e9; %[300,330,360,390,420,450,480] + uloops.laser_wavelength = [1310]; + uloops.M = [4]; + uloops.link_length = [2]; % 1,2,3,5,6,8,10 + %uloops.alpha = [0:0.1:1]; + + wh = DataStorage(uloops); + wh.addStorage("ber"); + + wh = submit_simulations(wh,"parallel",1,"simulation_mode",0); + +end + +col = cbrewer2('spectral',6);% +col=linspecer(5); + +cnt = 1; +for br = uloops.bitrate + a = wh_burg.getStoValue('ber',uloops.precomp, uloops.db_precode, br , uloops.laser_wavelength, uloops.M, uloops.link_length); + ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); + alpha = cellfun(@(x) x.vnle_pf_package{1,1}.pf.coefficients, a,'UniformOutput', false); + figure(23) + hold on + scatter(alpha{1}(2),ber_mlse,100,'MarkerEdgeColor',col(cnt,:),'Marker','x','LineWidth',2,'HandleVisibility','off'); + cnt = cnt+1; +end + + + +cnt = 1; +alpha = []; + +for br = uloops.bitrate + + a = wh_alphas.getStoValue('ber',uloops.precomp, uloops.db_precode, br , uloops.laser_wavelength, uloops.M, uloops.link_length, [0:0.1:1]); + ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); + + x_ax = [0:0.1:1]; + + figure(23) + hold on + % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,uloops.M)); + plot(x_ax,ber_mlse,'DisplayName',sprintf(' %d GBps PAM 4',br.*1e-9),'LineStyle','-','HandleVisibility','on','Color',col(cnt,:)); + + xticks(x_ax); + set(gca, 'YScale', 'log'); + ylim([1e-5 0.4]); + xlim([min(x_ax), max(x_ax) ]); + yline([3.8e-3, 2e-2],'HandleVisibility','off'); + legend + beautifyBERplot(); + xlabel('Channel $\alpha$'); + ylabel('BER'); + + cnt = cnt+1; + +end + + + + diff --git a/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_equalizer_taps.m b/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_equalizer_taps.m new file mode 100644 index 0000000..a393e25 --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_equalizer_taps.m @@ -0,0 +1,57 @@ +% basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +% db = DBHandler("pathToDB",[basePath,'silas_labor.db']); +if 0 + + uloops = struct; + uloops.precomp = [0]; + uloops.db_precode = [0]; + uloops.bitrate = [390].*1e9; %[300,330,360,390,420,450,480] + uloops.laser_wavelength = [1310]; + uloops.M = [4]; + uloops.link_length = [2]; % 1,2,3,5,6,8,10 + uloops.vnle_order1 = [5,10:10:100]; + uloops.vnle_order2 = [0]; + uloops.vnle_order3 = [0]; + + wh = DataStorage(uloops); + wh.addStorage("ber"); + + wh = submit_simulations(wh,"parallel",1,"simulation_mode",0); + +end + + +col = cbrewer2('spectral',6);% +col=linspecer(5); +cnt = 1; +for n3 = uloops.vnle_order3 + + a = wh.getStoValue('ber',uloops.precomp, uloops.db_precode, uloops.bitrate , uloops.laser_wavelength, uloops.M, uloops.link_length,... + uloops.vnle_order1,... + uloops.vnle_order2,... + n3); + + ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); + ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); + ber_db = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); + log_bers = log10(ber_vnle + 1e-12); + + figure(20) + hold on + title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,uloops.M)); + % plot(uloops.vnle_order1,ber_vnle,'DisplayName',sprintf('Tx precomp. + VNLE'),'LineStyle','-','HandleVisibility','on','Color',col(cnt,:)); + plot(uloops.vnle_order1,ber_mlse,'DisplayName',sprintf('VNLE + Postfilter + ;MLSE'),'LineStyle','-','HandleVisibility','on','Color',col(cnt,:)); + plot(uloops.vnle_order1,ber_db,'DisplayName',sprintf('DB tgt. VNLE + MLSE'),'LineStyle','-','HandleVisibility','on','Color',col(cnt+2,:)); + xticks(uloops.vnle_order1([1:1:end])); + set(gca, 'YScale', 'log'); + ylim([5e-5 0.4]); + xlim([min(uloops.vnle_order1), max(uloops.vnle_order1) ]); + yline([3.8e-3, 2e-2],'HandleVisibility','off'); + legend + beautifyBERplot(); + xlabel('Number of 1st order coeff.'); + ylabel('BER'); + + cnt = cnt+1; +end + diff --git a/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_lambda.m b/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_lambda.m new file mode 100644 index 0000000..6146a7d --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_lambda.m @@ -0,0 +1,50 @@ +if 0 + uloops = struct; + uloops.precomp = [0]; + uloops.db_precode = [0,1]; + uloops.bitrate = [420].*1e9; %[300,330,360,390,420,450,480] + % uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4]; + uloops.laser_wavelength = [1293, 1302,1310,1318,1327.4]; + uloops.M = [4,6,8]; + uloops.link_length = [5]; % 1,2,3,5,6,8,10 + wh = DataStorage(uloops); + wh.addStorage("ber"); + + wh = submit_simulations(wh,"parallel",1,"simulation_mode",0); +end + +col = cbrewer2('Set2',6);% +col=linspecer(5); +cnt = 1; + +m = 8; + +% a = wh.getStoValue('ber',1, 0, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); +% ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); + +a = wh.getStoValue('ber',0, 0, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); +ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); + +a = wh.getStoValue('ber',0, 1, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); +ber_db = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); + +x_ax = uloops.laser_wavelength; + +figure(21) +hold on +% title(sprintf('%d km | %d GBd | PAM %d',uloops.link_length,uloops.bitrate/log2(m).*1e-9,m)); +% plot(x_ax,ber_vnle,'DisplayName',sprintf('Tx precomp. + VNLE'),'LineStyle','-','HandleVisibility','on','Color',col(cnt+1,:)); +plot(x_ax,ber_mlse,'DisplayName',sprintf('VNLE + Postfilter + MLSE'),'LineStyle','-','HandleVisibility','on','Color',colorsets.DeepRed.RGB); +% plot(x_ax,ber_db,'DisplayName',sprintf('DB tgt. VNLE + MLSE'),'LineStyle','-','HandleVisibility','on','Color',col(cnt,:)); +xticks(x_ax([1:1:end])); +set(gca, 'YScale', 'log'); +ylim([5e-5 0.4]); +xlim([min(x_ax)-3, max(x_ax)+3 ]); +yline([3.8e-3, 2e-2],'HandleVisibility','off'); +legend +beautifyBERplot(); +xlabel('Number of 1st order coeff.'); +ylabel('BER'); + +cnt = cnt+1; + diff --git a/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_postfilter_taps.m b/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_postfilter_taps.m new file mode 100644 index 0000000..834b35a --- /dev/null +++ b/projects/HighSpeedExperiment_2024/auswertung deliv 06 feb 2025/analyze_postfilter_taps.m @@ -0,0 +1,70 @@ +% basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +% db = DBHandler("pathToDB",[basePath,'silas_labor.db']); +if 0 + + uloops = struct; + uloops.precomp = [0]; + uloops.db_precode = [0]; + uloops.bitrate = [300,330,360,390,420,450,480].*1e9; %[300,330,360,390,420,450,480] + uloops.laser_wavelength = [1310]; + uloops.M = [4]; + uloops.link_length = [2]; % 1,2,3,5,6,8,10 + uloops.vnle_order1 = [50]; + uloops.vnle_order2 = [7]; + uloops.vnle_order3 = [7]; + uloops.pf_ncoeffs = [1,2,3]; + + wh = DataStorage(uloops); + wh.addStorage("ber"); + + wh = submit_simulations(wh,"parallel",1,"simulation_mode",0); + +end + +col = cbrewer2('spectral',6);% +col=linspecer(5); +cnt = 1; +alpha = []; +for n = uloops.pf_ncoeffs + + a = wh.getStoValue('ber',uloops.precomp, uloops.db_precode, uloops.bitrate , uloops.laser_wavelength, uloops.M, uloops.link_length,... + uloops.vnle_order1,... + uloops.vnle_order2,... + uloops.vnle_order3,... + n); + + ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); + ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); + + % PF = cellfun(@(x) x.vnle_pf_package{1,1}.pf.coefficients, a,'UniformOutput',false); + % + % showTransferFunction(PF{3}.coefficients,"fignum",12,"color",clr.Set1.red,"DisplayName",['360 GBd']); + % + % showTransferFunction(PF{4}.coefficients,"fignum",12,"color",clr.Set1.blue,"DisplayName",['390 GBd']); + % + % showTransferFunction(PF{5}.coefficients,"fignum",12,'color',clr.Set1.green,"DisplayName",['420 GBd']); + % ber_db = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); + + x_ax = uloops.bitrate.*1e-9; + + figure(23) + hold on + title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,uloops.M)); + if n==1 + plot(x_ax,ber_vnle,'DisplayName',sprintf('Tx precomp. + VNLE'),'LineStyle','-','HandleVisibility','on','Color',col(cnt+4,:)); + end + plot(x_ax,ber_mlse,'DisplayName',sprintf('VNLE + Postfilter + ;MLSE'),'LineStyle','-','HandleVisibility','on','Color',col(cnt,:)); + % plot(x_ax,ber_db,'DisplayName',sprintf('DB tgt. VNLE + MLSE'),'LineStyle','-','HandleVisibility','on','Color',col(cnt+2,:)); + xticks(x_ax([1:1:end])); + set(gca, 'YScale', 'log'); + ylim([1e-5 0.4]); + xlim([min(x_ax(2:end)), max(x_ax) ]); + yline([3.8e-3, 2e-2],'HandleVisibility','off'); + legend + beautifyBERplot(); + xlabel('Gross Bitrate in Gbps'); + ylabel('BER'); + + cnt = cnt+1; +end + diff --git a/projects/HighSpeedExperiment_2024/auswertung/App/auswertungs_app.mlapp b/projects/HighSpeedExperiment_2024/auswertung/App/auswertungs_app.mlapp index 375605d..9fb4d4a 100644 Binary files a/projects/HighSpeedExperiment_2024/auswertung/App/auswertungs_app.mlapp and b/projects/HighSpeedExperiment_2024/auswertung/App/auswertungs_app.mlapp differ diff --git a/projects/HighSpeedExperiment_2024/db_auswertung/rate_vs_ber.m b/projects/HighSpeedExperiment_2024/db_auswertung/rate_vs_ber.m new file mode 100644 index 0000000..e53d7c6 --- /dev/null +++ b/projects/HighSpeedExperiment_2024/db_auswertung/rate_vs_ber.m @@ -0,0 +1,186 @@ + + +basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +database = DBHandler("pathToDB",[basePath,'silas_labor.db']); + +filterParams = database.tables; +filterParams.Configurations = struct( ... + 'bitrate', [], ... %[224,336,360,390,420,448] + 'db_mode', int32(db_mode.db_encoded), ... + 'fiber_length', 10, ... + 'interference_attenuation', [], ... + 'interference_path_length', [], ... + 'is_mpi', 0, ... + 'pam_level', 4, ... + 'rop_attenuation', 0, ... + 'wavelength', 1310 ... + ); + +% filterParams.EqualizerParameters.diff_precode = int32(db_mode.db_encoded); +% filterParams.EqualizerParameters.equalizer_structure = int32(equalizer_structure.vnle); +filterParams.EqualizerParameters.DCmu = 0.00; + +selectedFields = {'Configurations.run_id' 'Runs.rx_raw_path' 'Configurations.bitrate' 'Configurations.symbolrate' 'Configurations.pam_level'... + 'Configurations.db_mode' 'Configurations.rop_attenuation' 'Configurations.is_mpi' 'Configurations.interference_attenuation' ... + 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'Measurements.power_pd_in' ... + 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); + +fixedVars = {'eq_id','bitrate'}; +dataTableGrpd = groupIt(fixedVars,dataTable); + + +plotRealizations = 0; + +% Create a new figure +figure(); +hold on +unique_eq = unique(dataTable.eq_id); +cols = linspecer(8); +for i = 1:numel(unique_eq) + + idx = find(dataTableGrpd.eq_id == unique_eq(i), 1, 'first'); + equalizer_ = equalizer_structure(dataTableGrpd.equalizer_structure(idx)); + + + loop_filt = dataTableGrpd.eq_id==unique_eq(i); + + % Plot LINE: BER vs. interference_attenuation + switch equalizer_ + case equalizer_structure.vnle + bers = dataTableGrpd.BER(loop_filt,:); + case equalizer_structure.vnle_pf_mlse + bers = dataTableGrpd.BER(loop_filt,:); + case equalizer_structure.vnle_db_mlse + bers = dataTableGrpd.BER_precoded(loop_filt,:); + case equalizer_structure.db_encoded + bers = dataTableGrpd.BER(loop_filt,:); + end + + + name = sprintf('%s',equalizer_); + + p = plot(dataTableGrpd.bitrate(loop_filt,:).*1e-9, bers, '-', 'LineWidth', 0.5,'Color',cols(i,:),'DisplayName',name); + pair_one = {'Run ID', dataTableGrpd.run_id(loop_filt,:)}; + pair_two = {'Rate', dataTableGrpd.bitrate(loop_filt,:)}; + addDatatips(p, pair_one, pair_two); + xticks(unique(dataTableGrpd.bitrate(loop_filt,:).*1e-9)); + + % Plot SCATTERS: BER vs. interference_attenuation + + loop_filt = dataTable.eq_id==unique_eq(i); + if plotRealizations + switch equalizer_ + case equalizer_structure.vnle + bers = dataTable.BER(loop_filt,:); + case equalizer_structure.vnle_pf_mlse + bers = dataTable.BER(loop_filt,:); + case equalizer_structure.vnle_db_mlse + bers = dataTable.BER_precoded(loop_filt,:); + case equalizer_structure.db_encoded + bers = dataTable.BER(loop_filt,:); + end + + sc = scatter(dataTable.bitrate(loop_filt,:).*1e-9, bers, 'LineWidth', 0.5,'Marker','*','MarkerEdgeColor',cols(i,:),'HandleVisibility','off'); + pair_one = {'Run ID', dataTable.run_id(loop_filt,:)}; + pair_two = {'Rate', dataTable.bitrate(loop_filt,:)}; + addDatatips(sc, pair_one, pair_two); + xticks(unique(dataTable.bitrate(loop_filt,:).*1e-9)); + end + +end + + +% Label the axes and add a title +xlabel('Bitrate in Gbps'); +ylabel('BER'); +title('Line Rate vs. BER'); +yline(3.8e-3,'LineWidth',1,'LineStyle','--','HandleVisibility','off'); +% Enable grid for better readability +grid on; + +beautifyBERplot; +ylim([1e-4 0.5]); + + + +function resultTable = groupIt(fixedVars,dataTable) + + % Group by run_id and eq_id (adjust grouping keys as needed) + + [G, groupKeys] = findgroups(dataTable(:, fixedVars)); + + % Preallocate a cell array for aggregated data. + varNames = dataTable.Properties.VariableNames; + nVars = numel(varNames); + aggData = cell(height(groupKeys), nVars); + groupCount = zeros(height(groupKeys), 1); % To store the size of each group + + % Loop over each group. + for i = 1:height(groupKeys) + idx = (G == i); % Logical index for group i + groupCount(i) = sum(idx); % Count number of rows in this group + % For each variable in the table: + for j = 1:nVars + colData = dataTable.(varNames{j}); + if isnumeric(colData) + % For numeric data, compute the mean. + aggData{i, j} = min(colData(idx)); + else + % For non-numeric data, take the first entry. + if iscell(colData) + aggData{i, j} = colData{find(idx, 1)}; + else + aggData{i, j} = colData(find(idx, 1)); + end + end + end + end + + % Convert the aggregated cell array into a table. + resultTable = cell2table(aggData, 'VariableNames', varNames); + + % Append the group count as a new column. + resultTable.nRows = groupCount; + +end + + +function addDatatips(sc, varargin) +% addDatatips Adds custom data tip rows to a scatter plot. +% +% addDatatips(sc, pair1, pair2, ...) adds one or more custom rows to the +% data tip display of the scatter plot identified by sc. +% +% Each pair should be provided as a 1x2 cell array: {label, value}. +% The value can be a scalar or a vector. If a vector is provided, its length +% must match the number of scatter plot points. +% +% Example: +% sc = scatter(x, y, 'LineWidth', 1.5, 'Marker', 'o'); +% pair_one = {'Attenuation', attenuationVector}; +% addDatatips(sc, pair_one); + + numPoints = numel(sc.XData); + + for k = 1:length(varargin) + pair = varargin{k}; + + if ~iscell(pair) || numel(pair) ~= 2 + error('Each pair must be a 1x2 cell array: {label, value}.'); + end + + label = pair{1}; + value = pair{2}; + + % If value is a vector, ensure its length is either 1 or equal to the number of scatter points. + if isvector(value) && numel(value) ~= 1 && numel(value) ~= numPoints + error('The vector for "%s" must be a scalar or have %d elements matching the scatter data points.', label, numPoints); + end + + % Create a new data tip row using the provided label and vector. + newRow = dataTipTextRow(label, value); + sc.DataTipTemplate.DataTipRows(end+1) = newRow; + end +end diff --git a/projects/HighSpeedExperiment_2024/db_auswertung/realization_vs_ber.m b/projects/HighSpeedExperiment_2024/db_auswertung/realization_vs_ber.m new file mode 100644 index 0000000..b941680 --- /dev/null +++ b/projects/HighSpeedExperiment_2024/db_auswertung/realization_vs_ber.m @@ -0,0 +1,164 @@ + + +basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; +database = DBHandler("pathToDB",[basePath,'silas_labor.db']); + +filterParams = database.tables; +filterParams.Configurations = struct( ... + 'bitrate', 450e9, ... %[224,336,360,390,420,448] + 'db_mode', int32(db_mode.no_db), ... + 'fiber_length', 10, ... + 'interference_attenuation', [], ... + 'interference_path_length', [], ... + 'is_mpi', 0, ... + 'pam_level', 4, ... + 'rop_attenuation', 0, ... + 'wavelength', 1310 ... + ); + +% filterParams.EqualizerParameters.diff_precode = int32(db_mode.db_encoded); +% filterParams.EqualizerParameters.equalizer_structure = int32(equalizer_structure.vnle); +filterParams.EqualizerParameters.DCmu = 0.00; + +selectedFields = {'Configurations.run_id' 'Runs.rx_raw_path' 'Configurations.bitrate' 'Configurations.symbolrate' 'Configurations.pam_level'... + 'Configurations.db_mode' 'Configurations.rop_attenuation' 'Configurations.is_mpi' 'Configurations.interference_attenuation' ... + 'EqualizerParameters.equalizer_structure' 'EqualizerParameters.diff_precode' 'EqualizerParameters.eq_id' 'Measurements.power_pd_in' ... + 'Measurements.power_mpi_interference' 'Measurements.power_mpi_signal' 'Results.BER' 'Results.BER_precoded' 'Results.SNR' 'Results.GMI' 'Results.Alpha' 'Results.date_of_processing'}; + +[dataTable,sql_query] = database.queryDB(filterParams, selectedFields); + +fixedVars = {'eq_id','bitrate'}; +dataTableGrpd = groupIt(fixedVars,dataTable); + + +plotRealizations = 1; + +% Create a new figure +figure(); +hold on +unique_eq = unique(dataTable.eq_id); +cols = linspecer(8); + +for i = 1:numel(unique_eq) + + idx = find(dataTableGrpd.eq_id == unique_eq(i), 1, 'first'); + equalizer_ = equalizer_structure(dataTableGrpd.equalizer_structure(idx)); + + % Plot SCATTERS: timestamp vs. interference_attenuation + + loop_filt = dataTable.eq_id==unique_eq(i); + if plotRealizations + switch equalizer_ + case equalizer_structure.vnle + bers = dataTable.BER(loop_filt,:); + case equalizer_structure.vnle_pf_mlse + bers = dataTable.BER(loop_filt,:); + case equalizer_structure.vnle_db_mlse + bers = dataTable.BER_precoded(loop_filt,:); + case equalizer_structure.db_encoded + bers = dataTable.BER(loop_filt,:); + end + date_of_proc = datetime(dataTable.date_of_processing(loop_filt,:)); + + sc = scatter(date_of_proc, bers, 'LineWidth', 0.5,'Marker','*','MarkerEdgeColor',cols(i,:),'HandleVisibility','off'); + pair_one = {'Run ID', dataTable.run_id(loop_filt,:)}; + pair_two = {'Rate', dataTable.bitrate(loop_filt,:)}; + addDatatips(sc, pair_one, pair_two); + xticks(date_of_proc); + end + +end + + +% Label the axes and add a title +xlabel('Time of Processing'); +ylabel('BER'); +title('Line Rate vs. BER'); +yline(3.8e-3,'LineWidth',1,'LineStyle','--','HandleVisibility','off'); +% Enable grid for better readability +grid on; + +beautifyBERplot; +ylim([1e-4 0.5]); + + + +function resultTable = groupIt(fixedVars,dataTable) + + % Group by run_id and eq_id (adjust grouping keys as needed) + + [G, groupKeys] = findgroups(dataTable(:, fixedVars)); + + % Preallocate a cell array for aggregated data. + varNames = dataTable.Properties.VariableNames; + nVars = numel(varNames); + aggData = cell(height(groupKeys), nVars); + groupCount = zeros(height(groupKeys), 1); % To store the size of each group + + % Loop over each group. + for i = 1:height(groupKeys) + idx = (G == i); % Logical index for group i + groupCount(i) = sum(idx); % Count number of rows in this group + % For each variable in the table: + for j = 1:nVars + colData = dataTable.(varNames{j}); + if isnumeric(colData) + % For numeric data, compute the mean. + aggData{i, j} = min(colData(idx)); + else + % For non-numeric data, take the first entry. + if iscell(colData) + aggData{i, j} = colData{find(idx, 1)}; + else + aggData{i, j} = colData(find(idx, 1)); + end + end + end + end + + % Convert the aggregated cell array into a table. + resultTable = cell2table(aggData, 'VariableNames', varNames); + + % Append the group count as a new column. + resultTable.nRows = groupCount; + +end + + +function addDatatips(sc, varargin) +% addDatatips Adds custom data tip rows to a scatter plot. +% +% addDatatips(sc, pair1, pair2, ...) adds one or more custom rows to the +% data tip display of the scatter plot identified by sc. +% +% Each pair should be provided as a 1x2 cell array: {label, value}. +% The value can be a scalar or a vector. If a vector is provided, its length +% must match the number of scatter plot points. +% +% Example: +% sc = scatter(x, y, 'LineWidth', 1.5, 'Marker', 'o'); +% pair_one = {'Attenuation', attenuationVector}; +% addDatatips(sc, pair_one); + + numPoints = numel(sc.XData); + + for k = 1:length(varargin) + pair = varargin{k}; + + if ~iscell(pair) || numel(pair) ~= 2 + error('Each pair must be a 1x2 cell array: {label, value}.'); + end + + label = pair{1}; + value = pair{2}; + + % If value is a vector, ensure its length is either 1 or equal to the number of scatter points. + if isvector(value) && numel(value) ~= 1 && numel(value) ~= numPoints + error('The vector for "%s" must be a scalar or have %d elements matching the scatter data points.', label, numPoints); + end + + % Create a new data tip row using the provided label and vector. + newRow = dataTipTextRow(label, value); + sc.DataTipTemplate.DataTipRows(end+1) = newRow; + end +end diff --git a/projects/IMDD_base_system/imdd_it.m b/projects/IMDD_base_system/imdd_it.m index 01a40cd..d2fdf3c 100644 --- a/projects/IMDD_base_system/imdd_it.m +++ b/projects/IMDD_base_system/imdd_it.m @@ -5,141 +5,307 @@ if 1 uloops = struct; uloops.precomp = [1]; uloops.db_precode = [0]; - uloops.bitrate = [300,330,360,390,420,450,480].*1e9; %[300,330,360,390,420,450,480] - % uloops.bitrate = 390e9; + uloops.bitrate = [224].*1e9; %[300,330,360,390,420,450,480] [224,336,360,390,420,448] for MPI % uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4]; uloops.laser_wavelength = [1310]; - uloops.M = [6]; - uloops.link_length = [2]; % 1,2,3,5,6,8,10 + uloops.M = [4]; + uloops.link_length = [1]; % 1,2,3,5,6,8,10 + uloops.interference_attenuation = [0,3,6,9,12,15,18,21,24,27,30,45]; wh = DataStorage(uloops); wh.addStorage("ber"); - - wh = submit_simulations(wh,"parallel",0,"simulation_mode",0); + + % wh = submit_simulations(wh,"parallel",0,"simulation_mode",0); + wh = submit_handle(@imdd_model,wh,"parallel",1); + end -wh_ana = wh; +wh_ana = wh_master; -ber_mlse = {}; -ber_vnle = {}; -inf_rate_vnle ={}; +cols = cbrewer2('Paired',8); -ber_dbtgt ={}; +figure() -ber_dbenc ={}; -alpha = {}; -ber_dfe = {}; -ngmi = []; +for precomp = [0,1] + wavelength=uloops.laser_wavelength; + for m = [6] -wavelength=uloops.laser_wavelength; -for m = [4,6,8] - %1302 - %VNLE - precomp = 1; - precode = 0; - a = wh_ana.getStoValue('ber',precomp, precode, uloops.bitrate , wavelength, m, uloops.link_length); - ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); - - %MLSE - % precomp = 0; - % precode = 1; - a = wh_ana.getStoValue('ber',precomp, precode, uloops.bitrate , wavelength, m, uloops.link_length); - ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); - - %DB - % precomp = 0; - % precode = 1; - a = wh_ana.getStoValue('ber',precomp, precode, uloops.bitrate , wavelength, m, uloops.link_length); - ber_dbtgt = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); - - figure(m+20) - hold on - title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); - - plot(uloops.bitrate,ber_vnle,'DisplayName',sprintf('VNLE'),'Color',cols(3,:),'LineStyle','-','HandleVisibility','on'); - - plot(uloops.bitrate,ber_dbtgt,'DisplayName',sprintf('DB tgt. + MLSE',uloops.link_length,uloops.M),'Color',cols(1,:),'LineStyle','-','HandleVisibility','on'); - - plot(uloops.bitrate,ber_mlse,'DisplayName',sprintf('VNLE + 1 tap post-filter + MLSE',uloops.link_length,uloops.M),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); - - set(gca, 'YScale', 'log'); - ylim([5e-5 0.3]); - % xlim([min(uloops.bitrate.*1e-9), max(uloops.bitrate.*1e-9) ]); - yline([3.8e-3, 2e-2],'HandleVisibility','off'); - legend - beautifyBERplot() - xlabel('Bit Rate in Gbps'); - ylabel('BER'); -end + baudrate = wh_ana.parameter.bitrate.values; -cols = linspecer(7);%cbrewer2('Set2',10); + %VNLE + precode = 1; + a = wh_ana.getStoValue('ber',precomp, precode, baudrate , wavelength, m, uloops.link_length); + ber_vnle_pc = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); + %MLSE + ber_mlse_pc = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); + %DB + ber_dbtgt_pc = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); + precode = 0; + a = wh_ana.getStoValue('ber',precomp, precode, baudrate , wavelength, m, uloops.link_length); + ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); + %MLSE + ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); + %DB + ber_dbtgt = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); -for w = uloops.laser_wavelength - - figure(w) - figcnt = 0; - - for precode = uloops.db_precode - - for precomp = uloops.precomp - - for m = uloops.M - - a = wh_ana.getStoValue('ber',precomp, precode, uloops.bitrate , w, m, uloops.link_length); - ber_dbtgt = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); - - ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); - - ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); - - figcnt = figcnt+1; - subplot(4,3,figcnt); - hold on - title(sprintf('precomp = %d | precode = %d | %d km | %d nm | PAM %d',precomp,precode,uloops.link_length,w,m)); - - plot(uloops.bitrate,ber_dbtgt,'DisplayName',sprintf('DB tgt. + MLSE',uloops.link_length,uloops.M),'Color',cols(1,:),'LineStyle','-','HandleVisibility','on'); - - plot(uloops.bitrate,ber_vnle,'DisplayName',sprintf('VNLE'),'Color',cols(3,:),'LineStyle','-','HandleVisibility','on'); - - plot(uloops.bitrate,ber_mlse,'DisplayName',sprintf('VNLE + 1 tap post-filter + MLSE',uloops.link_length,uloops.M),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); - - % plot(uloops.bitrate,cellfun(@min, ber_dfe),'DisplayName',sprintf('VNLE + DFE',uloops.link_length,uloops.M),'Color',cols(2,:),'LineStyle','--'); - - % plot(uloops.bitrate,cellfun(@min, ber_dbenc),'DisplayName',sprintf('DB Encoded',uloops.link_length,uloops.M),'Color',cols(5,:),'LineStyle','-'); - - set(gca, 'YScale', 'log'); - ylim([5e-5 0.5]); - % xlim([min(uloops.bitrate.*1e-9), max(uloops.bitrate.*1e-9) ]); - yline([3.8e-3, 2e-2],'HandleVisibility','off'); - legend - beautifyBERplot() - xlabel('Bit Rate in Gbps'); - ylabel('Channel Wavelength (nm)'); - - end + if precomp + legndname1 = ['Pre-Emphasis']; + else + legndname1 = ['No Pre-Emphasis']; end + + + baudrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* 2.5 .* 1e9; + + subplot(1,3,1) + hold on + title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); + title(sprintf('PAM %d',m)); + plot(baudrate.*1e-9,ber_vnle,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(1+precomp,:),'LineStyle','-','HandleVisibility','on'); + plot(baudrate.*1e-9,ber_vnle_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(1+precomp,:),'LineStyle','-.','HandleVisibility','on'); + xticks(baudrate.*1e-9); + set(gca, 'YScale', 'log'); + ylim([5e-5 0.4]); + xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); + yline([3.8e-3, 2e-2],'HandleVisibility','off'); + legend + beautifyBERplot() + xlabel('Bit Rate in Gbps'); + ylabel('BER'); + + + + subplot(1,3,2) + hold on + % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); + title(sprintf('PAM %d',m)); + plot(baudrate.*1e-9,ber_dbtgt,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(5+precomp,:),'LineStyle','-','HandleVisibility','on'); + plot(baudrate.*1e-9,ber_dbtgt_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(5+precomp,:),'LineStyle','-.','HandleVisibility','on'); + xticks(baudrate.*1e-9); + set(gca, 'YScale', 'log'); + ylim([5e-5 0.4]); + xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); + yline([3.8e-3, 2e-2],'HandleVisibility','off'); + legend + beautifyBERplot() + xlabel('Bit Rate in Gbps'); + ylabel('BER'); + + + subplot(1,3,3) + hold on + % title(sprintf('%d km | %d nm | PAM %d',uloops.link_length,wavelength,m)); + title(sprintf('PAM %d',m)); + plot(baudrate.*1e-9,ber_mlse,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 0'],'Color',cols(3+precomp,:),'LineStyle','-','HandleVisibility','on'); + plot(baudrate.*1e-9,ber_mlse_pc,'DisplayName',['Pre-Emphasis: ', num2str(precomp), '| Diff.-Code: 1'],'Color',cols(3+precomp,:),'LineStyle','-.','HandleVisibility','on'); + xticks(baudrate.*1e-9); + set(gca, 'YScale', 'log'); + ylim([5e-5 0.4]); + xlim([min(baudrate(2).*1e-9), max(baudrate.*1e-9) ]); + yline([3.8e-3, 2e-2],'HandleVisibility','off'); + legend + beautifyBERplot() + xlabel('Bit Rate in Gbps'); + ylabel('BER'); end end +% +% +% cols = linspecer(7);%cbrewer2('Set2',10); +% +% +% for w = uloops.laser_wavelength +% +% figure(w) +% figcnt = 0; +% +% for precode = uloops.db_precode +% +% for precomp = uloops.precomp +% +% for m = uloops.M +% +% a = wh_ana.getStoValue('ber',precomp, precode, uloops.bitrate , w, m, uloops.link_length); +% ber_dbtgt = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); +% +% ber_vnle = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); +% +% ber_mlse = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); +% +% figcnt = figcnt+1; +% subplot(4,3,figcnt); +% hold on +% title(sprintf('precomp = %d | precode = %d | %d km | %d nm | PAM %d',precomp,precode,uloops.link_length,w,m)); +% +% plot(uloops.bitrate,ber_dbtgt,'DisplayName',sprintf('DB tgt. + MLSE',uloops.link_length,uloops.M),'Color',cols(1,:),'LineStyle','-','HandleVisibility','on'); +% +% plot(uloops.bitrate,ber_vnle,'DisplayName',sprintf('VNLE'),'Color',cols(3,:),'LineStyle','-','HandleVisibility','on'); +% +% plot(uloops.bitrate,ber_mlse,'DisplayName',sprintf('VNLE + 1 tap post-filter + MLSE',uloops.link_length,uloops.M),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); +% +% % plot(uloops.bitrate,cellfun(@min, ber_dfe),'DisplayName',sprintf('VNLE + DFE',uloops.link_length,uloops.M),'Color',cols(2,:),'LineStyle','--'); +% +% % plot(uloops.bitrate,cellfun(@min, ber_dbenc),'DisplayName',sprintf('DB Encoded',uloops.link_length,uloops.M),'Color',cols(5,:),'LineStyle','-'); +% +% set(gca, 'YScale', 'log'); +% ylim([5e-5 0.5]); +% % xlim([min(uloops.bitrate.*1e-9), max(uloops.bitrate.*1e-9) ]); +% yline([3.8e-3, 2e-2],'HandleVisibility','off'); +% legend +% beautifyBERplot() +% xlabel('Bit Rate in Gbps'); +% ylabel('Channel Wavelength (nm)'); +% +% end +% end +% end +% end -tp = TransmissionPerformance; -netRatesVNLE = tp.calculateNetRate(uloops.bitrate, 'NGMI', cellfun(@min, inf_rate_vnle)./log2(uloops.M), 'BER', cellfun(@min, ber_vnle)); -pam8= [2.9515 2.9284 2.9203 2.9311 2.8473 2.7740 2.6253]; -pam6 = [2.5280 2.5452 2.5579 2.5549 2.5272 2.4243 2.2617]; -pam4 = [ 1.9982 1.9972 1.9690 1.7909 1.2493 0.8014 0.6385]; -figure(6) -hold on -title(sprintf('Performance at 1310 for all lengths')); -% plot(uloops.bitrate.*1e-9,cellfun(@min, inf_rate_vnle),'DisplayName',sprintf('NGMI VNLE; %d km',uloops.link_length),'Color',cols(3,:),'LineStyle',':'); -plot(uloops.bitrate.*1e-9,pam4/2,'DisplayName',sprintf('GMI VNLE; PAM 4'),'Color',cols(1,:),'LineStyle',':'); -plot(uloops.bitrate.*1e-9,pam6/log2(6),'DisplayName',sprintf('GMI VNLE; PAM 6'),'Color',cols(2,:),'LineStyle',':'); -plot(uloops.bitrate.*1e-9,pam8/3,'DisplayName',sprintf('GMI VNLE; PAM 8'),'Color',cols(3,:),'LineStyle',':'); -xlabel('Gross Bitrate in Gbps'); -ylabel('NGMI') -beautifyBERplot() -ylim([0,1]); +m = 6; +ir = [2,2.5,3]; +cols = linspecer(6); +baudrate_gather = []; +for i = 1:3 + m = uloops.M(i); + + %%% GET VNLE VALS + precode = 0; + precomp = 1; + a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); + ber_vnle(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.ber_vnle, a); + + %%% GET DB VALS + precode = 1; + precomp = 0; + a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); + ber_db(i,:) = cellfun(@(x) x.dbtgt_package{1,1}.ber, a); + + %%% GET MLSE VALS + precode = 0; + precomp = 0; + a = wh_master.getStoValue('ber',precomp, precode, uloops.bitrate , uloops.laser_wavelength, m, uloops.link_length); + ber_mlse(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.ber_mlse, a); + + inf_rate_pam(i,:) = cellfun(@(x) x.vnle_pf_package{1,1}.air, a); + inf_rate_pam(i,:) = inf_rate_pam(i,:)./log2(m); + + bitrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* ir(i) .* 1e9; + baudrate = floor( uloops.bitrate.*1e-9 ./log2(m) ) .* 1e9; + baudrate_gather = union(baudrate_gather,baudrate); + baudrate_ticks = 100:20:240; + bitrate_ticks = 300:30:480; + + tp = TransmissionPerformance; + netRatesVNLE = tp.calculateNetRate(bitrate, 'NGMI', inf_rate_pam(i,:), 'BER', ber_vnle(i,:)); + + %%% NGMI + figure(12) + hold on + title(sprintf('Performance at 1310 nm')); + plot(baudrate.*1e-9,inf_rate_pam(i,:),'DisplayName',sprintf('NGMI; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); + xlabel('Baud rate in GBd'); + ylabel('NGMI') + beautifyBERplot() + xticks(baudrate_ticks); + xlim([min(baudrate_ticks) max(baudrate_ticks)]); + + + %%% AIR + figure(14) + hold on + title(sprintf('Performance at 1310 nm')); + plot(baudrate.*1e-9,inf_rate_pam(i,:).*bitrate.*1e-9,'DisplayName',sprintf('AIR; PAM %d',m),'Color',cols(i,:),'LineStyle','-'); + xlabel('Baud rate in GBd'); + ylabel('AIR'); + beautifyBERplot() + xticks(baudrate_ticks); + xlim([min(baudrate_ticks) max(baudrate_ticks)]); + + %%% RATES + figure(16) + hold on + title(sprintf('Performance at 1310 nm')); + if i == 1 + hv = 'on'; + else + hv = 'off'; + end + plot(baudrate*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); + plot(baudrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); + plot(baudrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate*1e-9,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); + xticks(baudrate_ticks); + xlim([min(baudrate_ticks) max(baudrate_ticks)]); + xlabel('Baud rate in GBd'); + ylabel('Net Bitrate in Gbps') + beautifyBERplot() + ylim([250 410]) + + %%% CODE OVERHEAD IN % + figure(18) + hold on + title(sprintf('Performance at 1310 nm')); + plot(baudrate*1e-9,100.*(1-netRatesVNLE.SDHD.CodeRate)./netRatesVNLE.SDHD.CodeRate,'DisplayName',sprintf('SD+HD FEC',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility',hv,'Marker','o'); + plot(baudrate.*1e-9,100.*(1-netRatesVNLE.HD.CodeRate)./netRatesVNLE.HD.CodeRate,'DisplayName',sprintf('HD FEC',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','diamond'); + plot(baudrate.*1e-9,100.*(1-netRatesVNLE.KP4_hamming.CodeRate)./netRatesVNLE.KP4_hamming.CodeRate,'DisplayName',sprintf('KP4+Hamming',m),'Color',cols(i,:),'LineStyle','-.','HandleVisibility',hv,'Marker','square'); + xticks(baudrate_ticks); + xlim([min(baudrate_ticks) max(baudrate_ticks)]); + xlabel('Baud rate in GBd'); + ylabel('FEC Overhead in %') + beautifyBERplot() + + %%% CLASSIC BER + figure(22) + subplot(1,4,i) + hold on + plot(baudrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('Tx precomp + VNLE',m),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); + plot(baudrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('VNLE + PF + MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility','on','Marker','square'); + plot(baudrate*1e-9,ber_db(i,:),'DisplayName',sprintf('Diff. Code + DB tgt.',m),'Color',cols(i,:),'LineStyle','--','HandleVisibility','on','Marker','diamond'); + yline(4.85e-3,'HandleVisibility','off'); + yline(2e-2,'HandleVisibility','off'); + xticks(baudrate*1e-9); + xlim([min(baudrate*1e-9) max(baudrate*1e-9)]); + ylim([1e-4 0.3]); + xlabel('Baudrate in GBd'); + if i == 1 + ylabel('BER') + end + beautifyBERplot() + set(gca, 'YScale', 'log'); +legend + subplot(1,4,4) + hold on + if m == 4 + + plot(bitrate*1e-9,ber_db(i,:),'DisplayName',sprintf('Diff. Code + DB tgt.'),'Color',cols(i,:),'LineStyle','--','HandleVisibility','on','Marker','diamond'); + + elseif m == 6 + + plot(bitrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('VNLE + PF + MLSE'),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); + % plot(bitrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','square'); + + elseif m ==8 + + plot(bitrate*1e-9,ber_vnle(i,:),'DisplayName',sprintf('Tx precomp + VNLE'),'Color',cols(i,:),'LineStyle','-','HandleVisibility','on','Marker','o'); + % plot(bitrate*1e-9,ber_mlse(i,:),'DisplayName',sprintf('MLSE',m),'Color',cols(i,:),'LineStyle',':','HandleVisibility',hv,'Marker','square'); + + end + yline(4.85e-3,'HandleVisibility','off'); + yline(2e-2,'HandleVisibility','off'); + xticks(bitrate_ticks); + xlim([min(bitrate_ticks) max(bitrate_ticks)]); + ylim([1e-4 0.3]); + xlabel('Gross Bitrate in Gbps'); + % ylabel('BER') + beautifyBERplot() + set(gca, 'YScale', 'log'); + + +end + + figure() title(sprintf('%d km | 1310 nm | PAM %d | VNLE',uloops.link_length,uloops.M)); @@ -160,7 +326,7 @@ legend % plot(uloops.bitrate.*1e-9,ber_dfe,'DisplayName',sprintf('NGMI MLSE; %d km',len),'Color',cols(1,:),'LineStyle','-'); % plot(uloops.bitrate.*1e-9,ber_mlse,'DisplayName',sprintf('NGMI MLSE; %d km',len),'Color',cols(1,:),'LineStyle','-'); % plot(uloops.bitrate.*1e-9,ber_db,'DisplayName',sprintf('NGMI VNLE; %d km',len),'Color',cols(2,:),'LineStyle','-'); -% +% % set(gca, 'YScale', 'log'); % ylim([1e-5, 0.1]); % yline([3.8e-3;1e-2],'LineWidth',1,'HandleVisibility','off'); @@ -172,13 +338,13 @@ legend % for alpha = uloops.alpha % for precode = uloops.db_precode % for precomp = uloops.precomp -% -% +% +% % cols = linspecer(6);%cbrewer2('Set2',10); % cnt = 1; % a = wh.getStoValue('ber',alpha,uloops.vnle_order2,uloops.vnle_order3,precomp, precode, uloops.bitrate,uloops.M); -% -% +% +% % for i = 1:numel(a) % ber_db(i) = mean((a{i}.ber_db)); % ber_vnle(i) = mean((a{i}.ber_vnle)); @@ -186,98 +352,98 @@ legend % ber_mlse(i) = mean((a{i}.ber_mlse)); % pf_taps(i) = a{i}.pf_taps{1}(2); % vnle_taps(i,:) = a{i}.eq_vnle{1}.e; -% +% % figure(222) -% -% +% +% % subplot(7,3,3*(i-1)+1) % ylim([-0.5 1]); % title(sprintf('%d GBps',uloops.bitrate(i).*1e-9)); % hold on % stem(a{i}.eq_vnle{1}.e) -% +% % subplot(7,3,3*(i-1)+2) % ylim([-0.5 1]); % title(sprintf('%d GBps',uloops.bitrate(i).*1e-9)); % hold on % stem(a{i}.eq_vnle{1}.e2) -% +% % subplot(7,3,3*(i-1)+3) % ylim([-0.5 1]); % title(sprintf('%d GBps',uloops.bitrate(i).*1e-9)); % hold on % stem(a{i}.eq_vnle{1}.e3) -% +% % showEQNoisePSD(a{i}.noise_vnle{1},"fignum",220,"displayname",sprintf('%d 2nd order taps; %d GBps',vnle_order2,uloops.bitrate(i).*1e-9),"postfilter_taps",a{i}.pf_taps{1}); -% +% % eq_sig = a{i}.signal_vnle{1}; -% -% +% +% % end -% +% % figure() % hold on -% +% % if precomp % lsty = '-'; % else % lsty = '-'; % end -% +% % if precode % coloffset = 1; % else % coloffset=1; % end -% +% % % plot(uloops.bitrate,ber_mlse,'DisplayName',sprintf('Precomp: %d; Precode %d',precomp,precode),'Color',cols(cnt,:)); % title(sprintf('Precomp: %d; Precode %d',precomp,precode)); % plot(uloops.bitrate.*1e-9,ber_vnle,'DisplayName','VNLE','Color',cols(4,:),'LineStyle',lsty); % plot(uloops.bitrate.*1e-9,ber_vnle_dfe,'DisplayName','VNLE+DFE','Color',cols(2,:),'LineStyle',lsty); % plot(uloops.bitrate.*1e-9,ber_mlse,'DisplayName','VNLE+PF+MLSE','Color',cols(3,:),'LineStyle',lsty); % plot(uloops.bitrate.*1e-9,ber_db,'DisplayName','DB tgt.','Color',cols(1,:),'LineStyle',lsty); -% +% % cnt = cnt+1; % yline([3.8e-3;1e-2],'LineWidth',1,'HandleVisibility','off'); % beautifyBERplot() % end % end % end -% -% -% +% +% +% % for vnle_order3 = uloops.vnle_order3([1,3,5,7,9,10]) % cnt = 0; % for vnle_order2 = 3%= uloops.vnle_order2 % for precode = uloops.db_precode % for precomp = uloops.precomp -% -% +% +% % cols = linspecer(10);%cbrewer2('Set2',10); % cnt = cnt+1; % a = wh.getStoValue('ber',vnle_order2,vnle_order3,precomp, precode, uloops.bitrate,uloops.M); -% -% +% +% % ber_db(cnt) = mean((a{1}.ber_db)); % ber_vnle(cnt) = mean((a{1}.ber_vnle)); % ber_vnle_dfe(cnt) = mean((a{1}.ber_vnle_dfe)); % ber_mlse(cnt) = mean((a{1}.ber_mlse)); % pf_taps(cnt) = a{1}.pf_taps{1}(2); -% +% % showEQNoisePSD(a{1}.noise_vnle{1},"fignum",220,"displayname",sprintf('%d 3rd order taps; %d GBps',vnle_order3,uloops.bitrate(i).*1e-9)); -% +% % eq_sig = a{1}.signal_vnle{1}; -% -% +% +% % end % end % end -% +% % figure(180) % hold on % scatter(uloops.vnle_order2,ber_mlse,'MarkerEdgeColor',cols(vnle_order3,:),'LineWidth',1,'DisplayName',sprintf('%d 3rd order taps',vnle_order3)); % ylim([1e-3 0.5]); -% +% % end % title(sprintf('%d GBps',uloops.bitrate(i).*1e-9)); % ylabel('BER'); @@ -287,7 +453,7 @@ legend % r = 1; % for rate = uloops.bitrate -% +% % i = 1; % for alpha = uloops.alpha % cols = linspecer(7);%cbrewer2('Set2',10); @@ -297,24 +463,24 @@ legend % pf_taps(i) = a{1}.pf_taps{1}(2); % i = i+1; % end -% +% % figure(150) % hold on % scatter(pf_taps(1),ber_mlse(1),200,'DisplayName','Burg','MarkerEdgeColor',cols(r,:),'Marker','x','LineWidth',2,'HandleVisibility','off'); % plot(uloops.alpha(2:end),ber_mlse(2:end),'DisplayName',sprintf('%d GBps',rate.*1e-9),'Color',cols(r,:),'LineStyle','-'); % r=r+1; -% +% % end -% -% +% +% % % title(sprintf('%d GBps',uloops.bitrate.*1e-9)); % ylabel('BER'); % xlabel('Alpha'); % yline([3.8e-3;1e-2],'LineWidth',1,'HandleVisibility','off'); % beautifyBERplot() -% -% -% +% +% +% % figure(2024) % hold on % for j = 1:numel(uloops.vnle_order3) @@ -323,23 +489,23 @@ legend % ber_mlse(i) = all2nd{i}.ber_mlse; % ber_vnle(i) = all2nd{i}.ber_vnle; % pf_taps(:,i) = all2nd{i}.pf_taps; -% +% % end % plot(uloops.vnle_order2,ber_mlse,'DisplayName',sprintf('%d 3rd order',uloops.vnle_order3(j))); % end -% +% % yline(3.8e-3,'LineWidth',2,'DisplayName','3.8e-3'); % yline(2e-2,'LineWidth',2,'LineStyle','--','DisplayName','2e-2'); % beautifyBERplot() % legend -% +% % figure(2025) % clf; % Clear figure so we start fresh -% +% % numJ = numel(uloops.vnle_order3); % numI = numel(uloops.vnle_order2); % ber_mlse_mat = zeros(numJ, numI); -% +% % % Gather data into a 2D matrix % for j = 1:numJ % all2nd = wh.getStoValue('ber', uloops.vnle_order3(j), uloops.vnle_order2, uloops.bitrate, uloops.M); @@ -347,21 +513,21 @@ legend % ber_mlse_mat(j,i) = all2nd{i}.ber_mlse; % end % end -% +% % % Create a 2D plot % % 'imagesc' displays the matrix as an image with a colorbar. % imagesc(uloops.vnle_order2, uloops.vnle_order3, ber_mlse_mat); % set(gca,'YDir','normal'); % Ensure that lower vnle_order3 is at the bottom % colorbar; % Add a colorbar to show BER scale -% +% % xlabel('VNLE Order 2'); % ylabel('VNLE Order 3'); % title('BER MLSE as a function of VNLE Orders'); -% +% % % If you want to highlight certain BER levels, you can add contour lines: % hold on; % [C,h] = contour(uloops.vnle_order2, uloops.vnle_order3, ber_mlse_mat, [3.8e-3, 2e-2], 'LineWidth',2,'LineColor','k'); % clabel(C,h,'Color','k','FontWeight','bold'); -% +% % beautifyBERplot(); % legend('BER contour lines'); \ No newline at end of file diff --git a/projects/IMDD_base_system/imdd_model.m b/projects/IMDD_base_system/imdd_model.m index 800e946..c15a919 100644 --- a/projects/IMDD_base_system/imdd_model.m +++ b/projects/IMDD_base_system/imdd_model.m @@ -1,4 +1,6 @@ -function [output] = imdd_model(simulation_mode,varargin) +function [output] = imdd_model(varargin) + +simulation_mode = 0; %%% Change folder curFolder = pwd; @@ -17,6 +19,9 @@ fdac = 256e9; fadc = 256e9; random_key = 1; +interference_attenuation = 0; +is_mpi = 1; + precomp = 0; db_precode = 0; @@ -35,7 +40,7 @@ tx_bw_nyquist = 0.8; link_length = 1; % RX -rop = -8; +rop = -5; rx_bw_nyquist = 0.8; vnle_order1 = 50; @@ -45,6 +50,7 @@ vnle_order3 = 7; vnle_order=[vnle_order1,vnle_order2,vnle_order3]; dfe_order = [0 0 0]; +pf_ncoeffs = 1; alpha = 0; @@ -54,6 +60,7 @@ mu_ffe1 = 0.0001; mu_ffe2 = 0.0008; mu_ffe3 = 0.001; mu_dc = 0.005; +% mu_dc = 0; mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; mu_dfe = 0.0004; @@ -61,7 +68,7 @@ mu_dfe = 0.0004; dfe_ = sum(dfe_order)>0; -doub_mode = db_mode.db_precoded; +doub_mode = db_mode.no_db; %%% change specific parameter if given in varargin % Parse optional input arguments @@ -132,10 +139,15 @@ f_nyquist = fsym/2; %%% run the simulation or measurement or ... if simulation_mode - Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rcalpha); + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); + rcalpha = 1; + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha); + db_precode = 0; + db_encode = 0; + apply_pulsef = 1; [Digi_sig,Symbols,Tx_bits] = PAMsource(... - "fsym",fsym,"M",M,"order",17,"useprbs",1,... + "fsym",fsym,"M",M,"order",18,"useprbs",0,... "fs_out",fdac,... "applyclipping",0,"clipfactor",1.5,... "applypulseform",apply_pulsef,"pulseformer",Pform,... @@ -143,7 +155,9 @@ if simulation_mode "db_precode",db_precode,"db_encode",db_encode,... "mrds_code",0,"mrds_blocklength",512).process(); - % Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); + Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); + + %%%%% AWG % El_sig = M8199A("kover",kover).process(Digi_sig); @@ -207,7 +221,7 @@ else 'fiber_length', link_length, ... 'interference_attenuation', [], ... 'interference_path_length', [], ... - 'is_mpi', 0, ... + 'is_mpi', is_mpi, ... 'pam_level', M, ... 'precomp_amp', [], ... 'rop_attenuation', 0, ... @@ -218,7 +232,8 @@ else ); selectedFields = {'Runs.run_id','Runs.tx_bits_path', 'Runs.tx_symbols_path', 'Runs.rx_sync_path','Runs.rx_raw_path',... - 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','Measurements.power_rop','Configurations.v_bias'}; + 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','Measurements.power_rop','Configurations.v_bias',... + 'Configurations.interference_attenuation'}; [dataTable,sql_query] = database.queryDB(filterParams, selectedFields); [~, uniqueIdx] = unique(dataTable.run_id); % Get unique run_id indices @@ -236,13 +251,14 @@ else % Raw_signal = load([basePath, char(dataTable.rx_raw_path(1))]); + % Raw_signal.Scpe_sig_raw.plot("displayname",'0db atten','fignum',10101) % Raw_signal = Raw_signal.Scpe_sig_raw; % % Raw_signal = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.55,"fs",Raw_signal.fs,"filterType",filtertypes.gaussian,"active",true).process(Raw_signal); % - Scpe_cell{1}.eye(fsym,M,"displayname",'eye','fignum',227); + % Scpe_cell{1}.eye(fsym,M,"displayname",'eye','fignum',227); % - % Raw_signal.spectrum("normalizeTo0dB",0,"fignum",336,"fft_length",2^12); + % Raw_signal.spectrum("normalizeTo0dB",0,"fignum",11,"fft_length",2^12); % Raw_signal.move_it_spectrum("fignum",334); % Raw_signal.move_it_spectrum("fignum",334); @@ -261,7 +277,7 @@ dbtgt_package = {}; proc_occ = min(1,length(Scpe_cell)); -for occ = 1:proc_occ +for occ = 1%:proc_occ Scpe_sig = Scpe_cell{occ}; @@ -269,11 +285,18 @@ for occ = 1:proc_occ Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); %%%%%% Sync Rx signal with reference %%%%%% - % [Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); + [Scpe_sig,~] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); + % + % Pform = Pulseformer("fsym",Scpe_sig.fs,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha,"matched",0); + % + % Scpe_sig_matched = Pform.process(Scpe_sig); + % + % Scpe_sig.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","scope "); + % Scpe_sig_matched.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","matched"); %%% EQUALIZING @@ -283,15 +306,17 @@ for occ = 1:proc_occ % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",512,"mu_dc",0.05); mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; + vnle_order=[vnle_order1,vnle_order2,vnle_order3]; % %%%%% VNLE + DFE %%%% - if 1 - eq_vnle_dfe = EQ("Ne",vnle_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",1,"ideal_dfe",0); - - [result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",1); + if 0 + + eq_vnle_dfe = EQ("Ne",vnle_order,"Nb",[0,0,0],"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",0); + eq_2 = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); + + [result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",1,"postFFE",[]); vnle_dfe_package{occ} = result; - end %%%%% VNLE + PF + MLSE %%%% if 1 @@ -299,7 +324,7 @@ for occ = 1:proc_occ % len_tr = length(Symbols)-1000; eq_vnle_ = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); % eq_vnle_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",vnle_order,"sps",2,"decide",0); - pf_ = Postfilter("ncoeff",1,"useBurg",1); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); [result] = vnle_postfilter_mlse(eq_vnle_,pf_,mlse_,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,'showAnalysis',1); @@ -314,7 +339,7 @@ for occ = 1:proc_occ mlse_db = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); eq_db = EQ("Ne",vnle_order,"Nb",dfe_order,"training_length",len_tr,"training_loops",5,"dd_loops",5,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005,"FFEmu",0,"plotfinal",0,"ideal_dfe",1); - [result] = duobinary_target(eq_db, mlse_db, M, Scpe_sig, Symbols, Tx_bits, "precode_mode", doub_mode,'showAnalysis',1); + [result] = duobinary_target(eq_db, mlse_db, M, Scpe_sig, Symbols, Tx_bits, "precode_mode", doub_mode,'showAnalysis',0); dbtgt_package{occ} = result; diff --git a/projects/Lab_2024/offline_dsp_analysis/load_n_plot_transfer_characteristic.m b/projects/Lab_2024/offline_dsp_analysis/load_n_plot_transfer_characteristic.m index a23d25a..17d04f5 100644 --- a/projects/Lab_2024/offline_dsp_analysis/load_n_plot_transfer_characteristic.m +++ b/projects/Lab_2024/offline_dsp_analysis/load_n_plot_transfer_characteristic.m @@ -1,7 +1,7 @@ -precomp_path = "C:\Users\Silas\Documents\MATLAB\imdd_simulation\projects\standard_system"; -precomp_filename = "lab_mpi_setup_2"; +precomp_path = "C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\precomp"; +precomp_filename = "lab_high_speed"; freqresp = ChannelFreqResp("Nacq",1024,"Navg",64,"Ncp",63,'f_ref',92e9); freqresp.load('loadPath',precomp_path,'fileName',precomp_filename); diff --git a/projects/Messung_Zürich/AIR_simp_new.fig b/projects/Messung_Zürich/AIR_simp_new.fig new file mode 100644 index 0000000..dbc22cc Binary files /dev/null and b/projects/Messung_Zürich/AIR_simp_new.fig differ diff --git a/projects/Messung_Zürich/BER_ief.fig b/projects/Messung_Zürich/BER_ief.fig new file mode 100644 index 0000000..51692ef Binary files /dev/null and b/projects/Messung_Zürich/BER_ief.fig differ diff --git a/projects/Messung_Zürich/BER_simp_new.fig b/projects/Messung_Zürich/BER_simp_new.fig new file mode 100644 index 0000000..5f7044f Binary files /dev/null and b/projects/Messung_Zürich/BER_simp_new.fig differ diff --git a/projects/Messung_Zürich/NDR_simp_new.fig b/projects/Messung_Zürich/NDR_simp_new.fig new file mode 100644 index 0000000..9ba7d17 Binary files /dev/null and b/projects/Messung_Zürich/NDR_simp_new.fig differ diff --git a/projects/Messung_Zürich/NGMI_ief.fig b/projects/Messung_Zürich/NGMI_ief.fig new file mode 100644 index 0000000..eab9017 Binary files /dev/null and b/projects/Messung_Zürich/NGMI_ief.fig differ diff --git a/projects/Messung_Zürich/NGMI_simp_new.fig b/projects/Messung_Zürich/NGMI_simp_new.fig new file mode 100644 index 0000000..8dff0cc Binary files /dev/null and b/projects/Messung_Zürich/NGMI_simp_new.fig differ diff --git a/projects/Messung_Zürich/change_filenames.m b/projects/Messung_Zürich/change_filenames.m new file mode 100644 index 0000000..0d3afd9 --- /dev/null +++ b/projects/Messung_Zürich/change_filenames.m @@ -0,0 +1,85 @@ +%% 1) Read all files in the TR folder +pathToTimingRecov = "C:\Users\Silas\Documents\MATLAB\Datensätze\IEF_Polariton_2025\36_IMDD_Kiel\Data\TR_ZIP"; +trFiles = dir(fullfile(pathToTimingRecov, 'TR_SILAS_*.mat')); + +%% Directory for measurement files (used to extract metadata) +pathToMeasurement = "C:\Users\Silas\Documents\MATLAB\Datensätze\IEF_Polariton_2025\36_IMDD_Kiel\Data\20250221"; + +%% Initialize lists for different PAM types +listPAM2 = {}; +listPAM4 = {}; +listPAM6 = {}; +listPAM8 = {}; + +%% Loop over each TR file +for k = 1:length(trFiles) + % Get current TR file name + trFileName = trFiles(k).name; + + % 2) Extract file code from TR file name. + % For a filename like "TR_SILAS_20250221T001800.mat" the code is "20250221T001800". + filecode = extractBetween(trFileName, 'TR_SILAS_', '.mat'); + + % Find corresponding measurement file by code (custom function) + measurementFilename = findFileByCode(pathToMeasurement, filecode{1}); + + % 3) Extract parameters from the measurement filename using regex. + % Expected measurement filename format (example): + % "Pmod_-10p000dBm_P_PD_-20p000dBm_..._32GBd_4PAM__1234T5678" + tokens = regexp(measurementFilename, ... + 'Pmod_([-0-9p]+)dBm_P_PD_([-0-9p]+)dBm_.*?_(\d+)GBd_(\d+)PAM__\d+T\d+', ... + 'tokens'); + if isempty(tokens) + error('Filename format not recognized for measurement file: %s', measurementFilename); + end + tokens = tokens{1}; + + % Convert token strings to numbers + config.P_laser = str2double(strrep(tokens{1}, 'p', '.')); + config.P_pd = str2double(strrep(tokens{2}, 'p', '.')); + config.fsym = str2double(tokens{3}) * 1e9; % Convert from GBd to Hz + config.M = str2double(tokens{4}); + + % Display loaded metadata + fprintf('Loaded measurement file: %s\n', measurementFilename); + fprintf('P_laser: %.3f dBm\n', config.P_laser); + fprintf('P_pd: %.3f dBm\n', config.P_pd); + fprintf('fsym: %.1f GBd\n', config.fsym * 1e-9); + fprintf('M: %d\n', config.M); + + % 4) Rename the TR file to include the metadata. + % New filename format: TR_SILAS__Pmod_dBm_P_PD_dBm_GBd_PAM.mat + newTRname = sprintf('TR_SILAS_%s_Pmod_%.3fdBm_P_PD_%.3fdBm_%dGBd_%dPAM', ... + filecode{1}, config.P_laser, config.P_pd, config.fsym/1e9, config.M); + + newTRname = strrep(newTRname,'.','p'); + newTRname = [newTRname, '.mat']; + movefile(fullfile(pathToTimingRecov, trFileName), fullfile(pathToTimingRecov, newTRname)); + + % Append the file code to the corresponding PAM list based on config.M + switch config.M + case 2 + listPAM2{end+1} = filecode{1}; + case 4 + listPAM4{end+1} = filecode{1}; + case 6 + listPAM6{end+1} = filecode{1}; + case 8 + listPAM8{end+1} = filecode{1}; + otherwise + warning('Unexpected PAM value %d in file %s', config.M, measurementFilename); + end +end + +%% Display the lists of file codes for each PAM type +disp('List of file codes for PAM2:'); +disp(listPAM2); + +disp('List of file codes for PAM4:'); +disp(listPAM4); + +disp('List of file codes for PAM6:'); +disp(listPAM6); + +disp('List of file codes for PAM8:'); +disp(listPAM8); diff --git a/projects/Messung_Zürich/dsp_all_files.m b/projects/Messung_Zürich/dsp_all_files.m new file mode 100644 index 0000000..9ec98ab --- /dev/null +++ b/projects/Messung_Zürich/dsp_all_files.m @@ -0,0 +1,366 @@ + +M = 6; +if M == 2 + file_codes = {'20250221T035221' '20250221T035354' '20250221T035824' '20250221T035931' '20250221T040035' '20250221T040132' '20250221T040226' '20250221T040523' '20250221T040646' '20250221T040723' '20250221T040843' '20250221T041011' '20250221T041101' '20250221T041244'}; +elseif M == 4 + file_codes = {'20250221T030844' '20250221T032043' '20250221T032312' '20250221T032424' '20250221T032529' '20250221T032632' '20250221T032800' '20250221T033035' '20250221T033138' '20250221T033246' '20250221T033425' '20250221T033527' '20250221T033642' '20250221T033743' '20250221T033851' '20250221T034314' '20250221T034529' '20250221T034647' '20250221T034756' '20250221T034915'}; +elseif M == 6 + file_codes = {'20250221T041445' '20250221T041512' '20250221T041539' '20250221T041607' '20250221T041633' '20250221T041702' '20250221T041729' '20250221T041758' '20250221T041825' '20250221T041854' '20250221T041922' '20250221T041951' '20250221T042019' '20250221T042048' '20250221T042117' '20250221T042147' '20250221T042215'}; +elseif M ==8 + file_codes = {'20250221T004926' '20250221T023534' '20250221T024256' '20250221T024629' '20250221T024929' '20250221T025305' '20250221T025505' '20250221T025856' '20250221T030122' '20250221T030311' '20250221T030513'}; +end + +if 0 + uloops = struct; + uloops.filecode = file_codes; + uloops.vnle_order1 = [50]; + uloops.vnle_order2 = [5]; + uloops.vnle_order3 = [5]; + wh = DataStorage(uloops); + wh.addStorage("ber"); + + wh = submit_handle(@dsp_ief_file,wh,"parallel",1); +end + +if 0 + % Bring figure(5) to focus + fig = figure(2); + + % Get handles to all line objects in the figure + lines = findobj(fig, 'Type', 'line'); + + % Preallocate cell arrays to store the data for each line + xData = cell(numel(lines),1); + yData = cell(numel(lines),1); + + % Loop through each line and extract its data + + for k = 1:numel(lines) + xData{k} = get(lines(k), 'XData'); + yData{k} = get(lines(k), 'YData'); + end + + ief.M2.baudr_new = xData{4}; + ief.M2.ngmi_new = yData{4}; + + ief.M4.baudr_new = xData{3}; + ief.M4.ngmi_new = yData{3}; + + ief.M6.baudr_new = xData{2}; + ief.M6.ngmi_new = yData{2}; + + ief.M8.baudr_new = xData{1}; + ief.M8.ngmi_new = yData{1}; + + ief.M2.ber_new = yData{4}; + ief.M4.ber_new = yData{3}; + ief.M6.ber_new = yData{2}; + ief.M8.ber_new = yData{1}; + + pam2_baudr = xData{4}; + pam2_ber = yData{4}; + pam2_ngmi = yData{4}; + + pam4_baudr = xData{3}; + pam4_ber = yData{3}; + pam4_ngmi = yData{3}; + + pam6_baudr = xData{2}; + pam6_ber = yData{2}; + pam6_ngmi = yData{2}; + + pam8_baudr = xData{1}; + pam8_ber = yData{1}; + pam8_ngmi = yData{1}; + + pam2_baudr_lowdsp = xData{4}; + pam2_ber_lowdsp = yData{4}; + pam2_ngmi_lowdsp = yData{4}; + + pam4_baudr_lowdsp = xData{3}; + pam4_ber_lowdsp = yData{3}; + pam4_ngmi_lowdsp = yData{3}; + + pam6_baudr_lowdsp = xData{2}; + pam6_ber_lowdsp = yData{2}; + pam6_ngmi_lowdsp = yData{2}; + + pam8_baudr_lowdsp = xData{1}; + pam8_ber_lowdsp = yData{1}; + pam8_ngmi_lowdsp = yData{1}; + + + +end + + +close all +figure(1) +cnt = 1; +for M = [2,4,6,8] + % for i = 1:numel(uloops.vnle_order2) + % for j = 1:numel(uloops.vnle_order3) + wh = eval(sprintf('wh_pam%d',M)); + a =wh.getStoValue('ber',wh.parameter.filecode.values, wh.parameter.vnle_order1.values, wh.parameter.vnle_order2.values,wh.parameter.vnle_order3.values); + % a =wh_mit_1001.getStoValue('ber',uloops.filecode, wh_mit_1001.parameter.vnle_order1.values, wh_mit_1001.parameter.vnle_order2.values,wh_mit_1001.parameter.vnle_order3.values); + baudrate = cell2mat(cellfun(@(a) a.config.fsym, a, 'UniformOutput', false)); + [baudrate, idx] = sort(baudrate); + + % get best results per baudrate + vnle = 0; + try + ber_vnle_values = cellfun(@(a) cellfun(@(y) y.ber_vnle, a.vnle_package, 'UniformOutput', false), a, 'UniformOutput', false); + ber_vnle_infrate = cellfun(@(a) cellfun(@(y) y.inf_rate_vnle, a.vnle_package, 'UniformOutput', false), a, 'UniformOutput', false); + + best_vnle = cellfun(@(x) min(cell2mat(x)), ber_vnle_values); + best_vnle = best_vnle(idx); + best_gmi_vnle = cellfun(@(x) min(cell2mat(x)), ber_vnle_infrate); + best_gmi_vnle = best_gmi_vnle(idx); + + % netRateVNLE = TransmissionPerformance().calculateNetRate(log2(M)*baudrate,'NGMI',best_gmi_vnle./log2(M), 'BER',best_vnle); + + vnle = 1; + end + + mlse = 0; + try + ber_vnle_values = cellfun(@(a) cellfun(@(y) y.ber_vnle, a.vnle_pf_package, 'UniformOutput', false), a, 'UniformOutput', false); + best_vnle_2 = cellfun(@(x) min(cell2mat(x)), ber_vnle_values); + best_vnle_2 = best_vnle_2(idx); + + gmi_vnle_values = cellfun(@(a) cellfun(@(y) y.gmi, a.vnle_pf_package, 'UniformOutput', false), a, 'UniformOutput', false); + best_gmi_vnle = cellfun(@(x) min(cell2mat(x)), gmi_vnle_values); + best_gmi_vnle = best_gmi_vnle(idx); + + ber_mlse_values = cellfun(@(a) cellfun(@(y) y.ber_mlse, a.vnle_pf_package, 'UniformOutput', false), a, 'UniformOutput', false); + best_mlse = cellfun(@(x) min(cell2mat(x)), ber_mlse_values); + best_mlse = best_mlse(idx); + + if M == 8 + best_vnle_2(2) = []; + best_mlse(2) = []; + best_gmi_vnle(2) = []; + best_db(2) = []; + + baudrate(2) = []; + elseif M == 4 + + end + + mlse = 1; + end + + db = 0; + try + ber_db_values = cellfun(@(a) cellfun(@(y) y.ber, a.dbtgt_package, 'UniformOutput', false), a, 'UniformOutput', false); + best_db = cellfun(@(x) min(cell2mat(x)), ber_db_values); + best_db = best_db(idx); + if M == 4 + best_db(end-1:end) = []; + baudrate(end-1:end) = []; + end + + db = 1; + end + + %% BER PLOT + + % figure(1) + subplot(1,2,1) + cols = cbrewer2('Set1',6); + hold on + title(sprintf('BER')); + % title(sprintf('%d 1st order',uloops.vnle_order1(i))); + xax = baudrate.*1e-9; + + if M == 4 || M == 2 + plot(xax,best_db,'DisplayName',sprintf('PAM %d | DB+MLSE',M),'Color',cols(M/2,:),'LineStyle',':','HandleVisibility','on','Marker','o','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + plot(ief.(sprintf('M%d', M)).baudr_new,ief.(sprintf('M%d', M)).ber_new,'DisplayName',sprintf('PAM %d |VNLE',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + elseif M == 6 + plot(xax,best_vnle_2,'DisplayName',sprintf('PAM %d |VNLE',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + % plot(ief.(sprintf('M%d', M)).baudr_new,ief.(sprintf('M%d', M)).ber_new,'DisplayName',sprintf('PAM %d |VNLE',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + else + % plot(xax,best_vnle_2,'DisplayName',sprintf('PAM %d |VNLE',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + plot(ief.(sprintf('M%d', M)).baudr_new,ief.(sprintf('M%d', M)).ber_new,'DisplayName',sprintf('PAM %d |VNLE',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + end + + if 0 + if vnle + plot(xax,best_vnle,'DisplayName',sprintf('VNLE'),'Color',cols(2,:),'LineStyle','-','HandleVisibility','on'); + end + if mlse + plot(xax,best_vnle_2,'DisplayName',sprintf('VNLE'),'Color',cols(2,:),'LineStyle','-','HandleVisibility','on'); + plot(xax,best_mlse,'DisplayName',sprintf('VNLE + 1 tap post-filter + MLSE'),'Color',cols(3,:),'LineStyle','-','HandleVisibility','on'); + end + if db + plot(xax,best_db,'DisplayName',sprintf('DB tgt. + MLSE'),'Color',cols(1,:),'LineStyle','-','HandleVisibility','on'); + end + + plot(ief.(sprintf('M%d', M)).baudr,ief.(sprintf('M%d', M)).ber,'DisplayName',sprintf('VNLE [100,15,15] (ETH)'),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); + plot(ief.(sprintf('M%d', M)).baudr_lowdsp,ief.(sprintf('M%d', M)).ber_lowdsp,'DisplayName',sprintf('VNLE [100,15,15] (ETH)'),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); + end + + set(gca, 'YScale', 'log'); + ylim([1e-6 0.3]); + xlim([92, 260 ]); + xticks([0:16:280]); + yline([4.85e-3, 2e-2],'HandleVisibility','off'); + % legend + % beautifyBERplot() + xlabel('baudrate in GBd'); + ylabel('BER'); + set(findall(gca, '-property', 'Interpreter'), 'Interpreter', 'latex'); + + set(gcf, 'Color', 'w'); + set(gca, 'Box', 'on', 'LineWidth', 0.8); % Thicker border + grid on; + set(gca, 'FontSize', 10, 'FontName', 'Times New Roman'); + % + + %% NGMI PLOT + + % figure(2) + % cols = cbrewer2('Set1',6); + % hold on + % title(sprintf('%d km ; %.1f nm ; PAM %d',1,1313,M)); + % xax = baudrate.*1e-9; + % + % plot(xax,best_gmi_vnle./log2(M),'DisplayName',sprintf('NGMI VNLE'),'Color',cols(2,:),'LineStyle','-','HandleVisibility','on'); + % + % plot(ief.(sprintf('M%d', M)).baudr,ief.(sprintf('M%d', M)).ngmi,'DisplayName',sprintf('VNLE [100,15,15] (ETH)'),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); + % plot(ief.(sprintf('M%d', M)).baudr_lowdsp,ief.(sprintf('M%d', M)).ngmi_lowdsp,'DisplayName',sprintf('VNLE [100,15,15] (ETH)'),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); + % + % set(gca, 'YScale', 'log'); + % ylim([0.7 1]); + % xlim([min(xax), max(xax) ]); + % xticks(xax); + % yline([0.8],'HandleVisibility','off'); + % legend + % beautifyBERplot() + % xlabel('Baudrate in GBd'); + % ylabel('NGMI'); + + %% AIR Plot + % subplot(1,3,2) + % + % netRateMLSE = TransmissionPerformance().calculateNetRate(log2(M)*baudrate,'BER',best_mlse); + % netRateVNLE = TransmissionPerformance().calculateNetRate(log2(M)*baudrate,'NGMI',best_gmi_vnle./log2(M), 'BER',best_vnle_2); + % netRateDB = TransmissionPerformance().calculateNetRate(log2(M)*baudrate,'BER',best_db); + % + % netRateIEF = TransmissionPerformance().calculateNetRate(log2(M)*ief.(sprintf('M%d', M)).baudr,'BER',ief.(sprintf('M%d', M)).ber,'NGMI',ief.(sprintf('M%d', M)).ngmi); + % netRateIEF_lowdsp = TransmissionPerformance().calculateNetRate(log2(M)*ief.(sprintf('M%d', M)).baudr_lowdsp,'BER',ief.(sprintf('M%d', M)).ber_lowdsp,'NGMI',ief.(sprintf('M%d', M)).ngmi_lowdsp); + % + % cols = cbrewer2('Set1',8); + % hold on + % title(sprintf('AIR')); + % xax = baudrate.*1e-9; + % + % if M == 2 || M == 4 + % + % % plot(xax,netRateMLSE.HD.NetRate.*1e-9,'DisplayName',sprintf('MLSE - HD'),'Color',cols(1,:),'LineStyle',':','HandleVisibility','on'); + % % plot(xax,netRateVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('VNLE SD+HD'),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); + % plot(xax,netRateDB.HD.NetRate.*1e-9,'DisplayName',sprintf('PAM %d | DB+MLSE HD',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on'); + % plot(xax,netRateDB.O_FEC.NetRate.*1e-9,'DisplayName',sprintf('PAM %d | DB+MLSE O-FEC',M),'Color',cols(M/2,:),'LineStyle',':','HandleVisibility','on'); + % % plot(xax,netRateDB.KP4_hamming.NetRate.*1e-9,'DisplayName',sprintf('PAM %d | DB+MLSE KP4-FEC',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on'); + % + % else + % + % % plot(xax,netRateMLSE.HD.NetRate.*1e-9,'DisplayName',sprintf('MLSE - HD'),'Color',cols(1,:),'LineStyle',':','HandleVisibility','on'); + % % plot(xax,netRateVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('VNLE SD+HD'),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); + % % plot(ief.(sprintf('M%d', M)).baudr,netRateIEF.HD.NetRate,'DisplayName',sprintf('PAM %d | HD IEF',M),'Color',cols(M/2,:),'LineStyle',':','HandleVisibility','on'); + % plot(ief.(sprintf('M%d', M)).baudr,netRateIEF.SDHD.NetRate,'DisplayName',sprintf('PAM %d | SD+HD IEF',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on'); + % + % end + % set(gca, 'YScale', 'log'); + % ylim([92 450]); + % xlim([min(xax), max(xax) ]); + % xlim([92, 260 ]); + % xticks([0:16:280]); + % yline([0.8],'HandleVisibility','off'); + % legend + % beautifyBERplot() + % xlabel('Baudrate in GBd'); + % ylabel('Net Rate in Gbps'); + % + % + + + %% NDR PLOT + + subplot(1,2,2) + netRateIEF = TransmissionPerformance().calculateNetRate(log2(M)*ief.(sprintf('M%d', M)).baudr,'BER',ief.(sprintf('M%d', M)).ber,'NGMI',ief.(sprintf('M%d', M)).ngmi); + netRateIEF_lowdsp = TransmissionPerformance().calculateNetRate(log2(M)*ief.(sprintf('M%d', M)).baudr_lowdsp,'BER',ief.(sprintf('M%d', M)).ber_lowdsp,'NGMI',ief.(sprintf('M%d', M)).ngmi_lowdsp); + netRateIEF_new = TransmissionPerformance().calculateNetRate(log2(M)*ief.(sprintf('M%d', M)).baudr_new,'BER',ief.(sprintf('M%d', M)).ber_new,'NGMI',ief.(sprintf('M%d', M)).ngmi_new); + + + cols = cbrewer2('Set1',8); + hold on + title(sprintf('Net Bitrate')); + xax = baudrate.*1e-9; + + if M == 2 || M == 4 + + log2M = log2(M); + + netRateDB = TransmissionPerformance().calculateNetRate(log2(M)*baudrate,'BER',best_db); + + % plot(xax,netRateMLSE.HD.NetRate.*1e-9,'DisplayName',sprintf('MLSE - HD'),'Color',cols(1,:),'LineStyle',':','HandleVisibility','on'); + % plot(xax,netRateVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('VNLE SD+HD'),'Color',cols(4,:),'LineStyle','-','HandleVisibility','on'); + plot(xax,netRateDB.HD.NetRate.*1e-9,'DisplayName',sprintf('PAM %d | DB+MLSE HD',M),'Color',cols(M/2,:),'LineStyle',':','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + plot(xax,netRateDB.O_FEC.NetRate.*1e-9,'DisplayName',sprintf('PAM %d | DB+MLSE O-FEC',M),'Color',cols(M/2,:),'LineStyle','-.','HandleVisibility','on','Marker','x','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + % plot(xax,netRateDB.KP4_hamming.NetRate.*1e-9,'DisplayName',sprintf('PAM %d | DB+MLSE KP4-FEC',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','*','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + plot(ief.(sprintf('M%d', M)).baudr_new,netRateIEF_new.SDHD.NetRate,'DisplayName',sprintf('PAM %d | SD+HD',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','o','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + plot(ief.(sprintf('M%d', M)).baudr_new,ief.(sprintf('M%d', M)).ngmi_new.*ief.(sprintf('M%d', M)).baudr_new.*log2M ,'DisplayName',sprintf('PAM %d | AIR',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','o','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + + elseif M == 6 + + log2M = 2.5; + netRateMLSE = TransmissionPerformance().calculateNetRate(log2(M)*baudrate,'BER',best_mlse); + netRateVNLE = TransmissionPerformance().calculateNetRate(log2(M)*baudrate,'NGMI',best_gmi_vnle./log2(M), 'BER',best_vnle_2); + + plot(xax,netRateVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('MLSE - HD'),'Color',cols(M/2,:),'LineStyle',':','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + plot(xax,netRateVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('VNLE SD+HD'),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','o','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + % plot(ief.(sprintf('M%d', M)).baudr_new,netRateIEF_new.HD.NetRate,'DisplayName',sprintf('PAM %d | HD IEF',M),'Color',cols(M/2,:),'LineStyle',':','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + % plot(ief.(sprintf('M%d', M)).baudr_new,netRateIEF_new.SDHD.NetRate,'DisplayName',sprintf('PAM %d | SD+HD IEF',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','o','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + + else + + log2M = log2(M); + + netRateMLSE = TransmissionPerformance().calculateNetRate(log2(M)*baudrate,'BER',best_mlse); + netRateVNLE = TransmissionPerformance().calculateNetRate(log2(M)*baudrate,'NGMI',best_gmi_vnle./log2(M), 'BER',best_vnle_2); + + % plot(xax,netRateVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('MLSE - HD'),'Color',cols(M/2,:),'LineStyle',':','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + % plot(xax,netRateVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('VNLE SD+HD'),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','o','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + plot(ief.(sprintf('M%d', M)).baudr_new,netRateIEF_new.HD.NetRate,'DisplayName',sprintf('PAM %d | HD IEF',M),'Color',cols(M/2,:),'LineStyle',':','HandleVisibility','on','Marker','square','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + plot(ief.(sprintf('M%d', M)).baudr_new,netRateIEF_new.SDHD.NetRate,'DisplayName',sprintf('PAM %d | SD+HD IEF',M),'Color',cols(M/2,:),'LineStyle','-','HandleVisibility','on','Marker','o','MarkerFaceColor',[1,1,1],'MarkerEdgeColor',cols(M/2,:),'MarkerSize',4,'LineWidth',1.4); + + end + set(gca, 'YScale', 'log'); + ylim([150 450]); + xlim([min(xax), max(xax) ]); + xlim([92, 270 ]); + xticks([0:16:280]); + yline([0.8],'HandleVisibility','off'); + % legend + + xlabel('Baudrate in GBd'); + ylabel('Net Rate in Gbps'); + set(findall(gca, '-property', 'Interpreter'), 'Interpreter', 'latex'); + + set(gcf, 'Color', 'w'); + set(gca, 'Box', 'on', 'LineWidth', 0.8); % Thicker border + grid on; + set(gca, 'FontSize', 10, 'FontName', 'Times New Roman'); + +end \ No newline at end of file diff --git a/projects/Messung_Zürich/dsp_ief_file.m b/projects/Messung_Zürich/dsp_ief_file.m new file mode 100644 index 0000000..3fd0754 --- /dev/null +++ b/projects/Messung_Zürich/dsp_ief_file.m @@ -0,0 +1,203 @@ + +function [results] = dsp_ief_file(varargin) + + mu_ffe1 = 0.0001; + mu_ffe2 = 0.0008; + mu_ffe3 = 0.001; + mu_dc = 0.0003; %0.0003; + + mu_dfe_training = 0.0004; + vnle_order1 = 50; + vnle_order2 = 3; + vnle_order3 = 3; + + dfe_mu = 0.0005; + + + tcorrect = 0; + %%% change specific parameter if given in varargin + % Parse optional input arguments + if ~isempty(varargin) + var_s = varargin{1}; + if isstruct(var_s) + fields = fieldnames(var_s); + for i = 1:numel(fields) + if isnumeric(fields{i}) + eval([fields{i}, ' = ', num2str( var_s.(fields{i}) ), ';']); + fprintf("%s <-- %.2f \n", fields{i}, var_s.(fields{i})); + else + eval([fields{i}, ' = ', 'var_s.(fields{',num2str(i),'})' , ';']); + end + + end + else + error('Optional variables should be passed as a struct.'); + end + end + + pathToMeasurement = "C:\Users\Silas\Documents\MATLAB\Datensätze\IEF_Polariton_2025\36_IMDD_Kiel\Data\20250221"; + filename = findFileByCode(pathToMeasurement, filecode{1}); + + pathToTimingRecov = "C:\Users\Silas\Documents\MATLAB\Datensätze\IEF_Polariton_2025\36_IMDD_Kiel\Data\TR_ZIP"; + filename_2 = findFileByCode(pathToTimingRecov, filecode{1}); + + % Extract parameters from the filename using an updated regex + tokens = regexp(filename, 'Pmod_([-0-9p]+)dBm_P_PD_([-0-9p]+)dBm_.*?_(\d+)GBd_(\d+)PAM__\d+T\d+', 'tokens'); + if isempty(tokens) + error('Filename format not recognized.'); + end + + tokens = tokens{1}; + + % Convert values + config.P_laser = str2double(strrep(tokens{1}, 'p', '.')); + config.P_pd = str2double(strrep(tokens{2}, 'p', '.')); + config.fsym = str2double(tokens{3}) * 1e9; % Convert GBd to Hz + config.M = str2double(tokens{4}); + + % Display results + fprintf('Loaded file: %s\n', filename); + fprintf('P_laser: %.3f dBm\n', config.P_laser); + fprintf('P_pd: %.3f dBm\n', config.P_pd); + fprintf('fsym: %.1f GBd\n', config.fsym*1e-9); + fprintf('M: %d\n', config.M); + + %%% Load Data + filepath = fullfile(filename); + ief_ = h5info(filepath); + + config.fs_rx = h5readatt(filepath,'/','fs'); %sampling frequency at Rx + config.fs_tx = h5readatt(filepath, '/','fs_Tx'); %sampling frequency at Tx + config.fsym = h5readatt(filepath, '/','R'); %Baudrate + config.M = h5readatt(filepath, '/','M'); % PAM- 'M' + config.ROF = h5readatt(filepath, '/','ROF'); + config.PulseShape =h5readatt(filepath, '/','PulseShape'); + + yOrg = h5readatt(filepath, ief_.Groups(4).Groups(1).Name, 'YOrg'); + yInc = h5readatt(filepath, ief_.Groups(4).Groups(1).Name, 'YInc'); + + dataRx = double(h5read(filepath, '/Waveforms/Channel 2/Channel 2Data')); % rohdaten des CH4 + bitsTx = h5read(filepath, '/Settings/dataTx'); %Binär + + if config.M ~= 6 + bitsTx = reshape(bitsTx,log2(config.M),[])'; + else + %bitsTx = reshape(bitsTx,5,[])'; + end + + %%% Build Tx Signal (Bits, Pam Map, Symbols) + Tx_bits = Informationsignal(bitsTx); + Tx_symbols = PAMmapper(config.M,0,"eth_style",1).map(Tx_bits); + Tx_symbols.fs = config.fsym; + + if 0 + Rx_bits = PAMmapper(config.M,0,"eth_style",1).demap(Tx_symbols); + [~,~,ber_bw,~] = calc_ber(Rx_bits.signal,Tx_bits.signal(1:length(Rx_bits.signal)),"skip_front",100,"skip_end",150,"returnErrorLocation",1); + end + + %%% Build Rx Signal (Rx, normalize,remove mean) + loadAfterTR = 1; + if loadAfterTR + if 1 + rx_sig = load(filename_2); + rx_sig=rx_sig.signal_TR; + if config.M == 6 + rx_sig_pam6 = zeros(length(rx_sig)*2,1); + rx_sig_pam6(1:2:end) = real(rx_sig); + rx_sig_pam6(2:2:end) = imag(rx_sig); + rx_sig = rx_sig_pam6; + Rx_Sig_resamp = Informationsignal(rx_sig,"fs",config.fsym); + % Tx_symbols.signal = Tx_symbols.signal(1:end/2); + else + Rx_Sig_resamp = Informationsignal(rx_sig,"fs",config.fsym*2); + end + + else + rx_sig = load(string(['testSilas_',char(filecode{1}),'.mat']),"signal_TR3"); + rx_sig=rx_sig.signal_TR3; + Rx_Sig_resamp = Informationsignal(rx_sig,"fs",config.fsym*2); + end + + else + + dataRx = dataRx*yInc+yOrg; + Rx_Sig = Informationsignal(dataRx,"fs",config.fs_rx); + + Rx_Sig.signal = Rx_Sig.signal - mean(Rx_Sig.signal); + + Rx_Sig = Rx_Sig.normalize("mode","rms"); + + %%%%%% Sample to 2x fsym %%%%%% + Rx_Sig_resamp = Rx_Sig.resample("fs_out",config.fsym); + + end + + + %%%%%% Sync Rx signal with reference (S is a cell array with all occurences) %%%%%% + + [~,S,isFlipped] = Rx_Sig_resamp.tsynch("reference",Tx_symbols,"fs_ref",config.fsym,"debug_plots",1); + + output = struct(); + vnle_package = {}; + vnle_pf_package = {}; + dbtgt_package = {}; + + for s = 1%:length(S) + Rx_Sig_sync = S{s}; + Rx_Sig_sync = Rx_Sig_sync.normalize("mode","rms"); + Rx_Sig_sync = Rx_Sig_sync.resample("fs_out",2*config.fsym); + + mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; + vnle_order=[vnle_order1,vnle_order2,vnle_order3]; + eq_ = EQ("Ne",vnle_order,"Nb",[0,0,0],"training_length",4096*4,"training_loops",4,"dd_loops",3,"K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe_training],"DFEmu",dfe_mu,"FFEmu",0,"plotfinal",0,"ideal_dfe",0,"plottrain",0); + + + %%%%% VNLE only (or DFE) %%%% + if 0 + + eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",1001,"sps",1,"decide",0); + [result] = vnle(eq_,config.M,Rx_Sig_sync,Tx_symbols,Tx_bits,"precode_mode",db_mode.no_db,"showAnalysis",1,'eth_style',1,'postFFE',eq_post); + + netRate = TransmissionPerformance().calculateNetRate(log2(config.M)*config.fsym,'NGMI',result.inf_rate_vnle, 'BER',result.ber_vnle); + vnle_package{s} = result; + + end + + %%%%% VNLE + PF + MLSE %%%% + if 1 + + eq_post = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",1001,"sps",1,"decide",0); + pf_ = Postfilter("ncoeff",1,"useBurg",1); + mlse_ = MLSE("duobinary_output",0,'M',config.M,'trellis_states',PAMmapper(config.M,0).levels); + mlse_ = MLSE_viterbi("duobinary_output",0,'M',config.M,'trellis_states',PAMmapper(config.M,0).levels); + doub_mode = db_mode.no_db; + + [result] = vnle_postfilter_mlse(eq_,pf_,mlse_,config.M,Rx_Sig_sync,Tx_symbols,Tx_bits,"precode_mode",doub_mode,'showAnalysis',1,'eth_style_symbol_mapping',1,'postFFE',eq_post); + netRate = TransmissionPerformance().calculateNetRate(log2(config.M)*config.fsym,'NGMI',result.gmi./log2(config.M), 'BER',result.ber_mlse); + fprintf('VNLE SD: %.1f GBd \n',netRate.SDHD.NetRate.*1e-9); + fprintf('MLSE HD: %.1f GBd \n',netRate.HD.NetRate.*1e-9); + + vnle_pf_package{s} = result; + + end + + %%%%% Duobinary Targeting %%%% + if 1 + + mlse_db = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",config.M,"trellis_states",PAMmapper(config.M,0).levels); + doub_mode = db_mode.db_emulate; + + [result] = duobinary_target(eq_, mlse_db, config.M, Rx_Sig_sync, Tx_symbols, Tx_bits, "precode_mode", doub_mode,'showAnalysis',1,'eth_style_symbol_mapping',1); + dbtgt_package{s} = result; + + end + + end + + results.vnle_package = vnle_package; + results.vnle_pf_package = vnle_pf_package; + results.dbtgt_package = dbtgt_package; + results.config = config; + + +end diff --git a/projects/Messung_Zürich/findFileByCode.m b/projects/Messung_Zürich/findFileByCode.m new file mode 100644 index 0000000..7fd9eb4 --- /dev/null +++ b/projects/Messung_Zürich/findFileByCode.m @@ -0,0 +1,40 @@ +function fileName = findFileByCode(folderPath, code) +% FINDFILEBYCODE Searches for a file in a folder structure by a given code. +% fileName = findFileByCode(folderPath, code) searches recursively in +% folderPath for a file containing 'code' in its name and returns the +% full file name if found. +% +% Inputs: +% folderPath - The root directory to search in +% code - The unique code to search for in file names +% +% Output: +% fileName - The full file name if found, empty if not found + + % Initialize output + fileName = ''; + + % Get list of all files and folders in the folderPath + files = dir(folderPath); + + % Iterate through the list + for i = 1:length(files) + % Skip '.' and '..' + if files(i).isdir + if ~startsWith(files(i).name, '.') % Avoid hidden folders + % Recursive search in subdirectories + subFolder = fullfile(folderPath, files(i).name); + fileName = findFileByCode(subFolder, code); + if ~isempty(fileName) + return; % Stop searching if found + end + end + else + % Check if the file name contains the code + if contains(files(i).name, code) + fileName = fullfile(folderPath, files(i).name); + return; + end + end + end +end diff --git a/projects/Messung_Zürich/minimal_pam4_mlse.m b/projects/Messung_Zürich/minimal_pam4_mlse.m new file mode 100644 index 0000000..1a2bb54 --- /dev/null +++ b/projects/Messung_Zürich/minimal_pam4_mlse.m @@ -0,0 +1,19 @@ + + +mlse_sig_sd=load("imdd_simulation\projects\Messung_Zürich\mlse_sig_sd.mat","mlse_sig_sd"); +mlse_sig_sd = mlse_sig_sd.mlse_sig_sd; +tx_symbols=load("imdd_simulation\projects\Messung_Zürich\tx_symbols.mat","tx_symbols"); +tx_symbols = tx_symbols.tx_symbols; + +mlse_ = MLSE_viterbi("duobinary_output",0,'M',4,'trellis_states',PAMmapper(4,0).levels); +mlse_.DIR = [1.0000 0.5452]; +mlse_sig_sd = mlse_.process(mlse_sig_sd); + + +% +% mlse_ = MLSE_viterbi("duobinary_output",0,'M',4,'trellis_states',PAMmapper(4,0).levels); +% mlse_.DIR = [1.0000 0.5452]; +% mlse_sig_sd = mlse_.process(mlse_sig_sd,tx_symbols); + +h = [1.0000 0.5452]; +chatgpt_answer(mlse_sig_sd.signal,tx_symbols.signal,h) diff --git a/projects/Messung_Zürich/pam4.mat b/projects/Messung_Zürich/pam4.mat new file mode 100644 index 0000000..0b9e57d Binary files /dev/null and b/projects/Messung_Zürich/pam4.mat differ diff --git a/projects/Messung_Zürich/pam4_eq.fig b/projects/Messung_Zürich/pam4_eq.fig new file mode 100644 index 0000000..68e2e79 Binary files /dev/null and b/projects/Messung_Zürich/pam4_eq.fig differ diff --git a/projects/Messung_Zürich/submit_handle.m b/projects/Messung_Zürich/submit_handle.m new file mode 100644 index 0000000..2c5e710 --- /dev/null +++ b/projects/Messung_Zürich/submit_handle.m @@ -0,0 +1,69 @@ +function wh = submit_handle(funcHandle, wh, options) + +arguments + funcHandle + wh + options.parallel = 1; +end + +%%% 2) SUBMIT SIMULATION +% Initialize job results +if options.parallel + curpool = gcp('nocreate'); + if isempty(curpool) + parpool; + else + % stop all forgotten or unfetched jobs from queue + if ~isempty(curpool.FevalQueue.QueuedFutures) || ~isempty(curpool.FevalQueue.RunningFutures) + oldq = length(curpool.FevalQueue.QueuedFutures) + length(curpool.FevalQueue.RunningFutures); + curpool.FevalQueue.cancelAll + fprintf('Canceled %d unfetched jobs from old queue.', oldq); + end + end + results = parallel.FevalFuture.empty(); +else + results = []; +end + +fprintf('Requested %d loops\n', wh.getLastLinIndice); + +for lin_idx = 1:wh.getLastLinIndice + optionalVars = struct(); + + if ~isempty(wh.getDimension) + % Build the optionalVars struct + [parametervalues, parameternames] = wh.getPhysIndicesByLinIndex(lin_idx); + for pidx = 1:numel(parameternames) + optionalVars.(parameternames{pidx}) = parametervalues{pidx}; + end + end + + %%% SIMULATION HERE + if options.parallel + numOutputs = 1; + results(lin_idx) = parfeval(funcHandle, numOutputs, optionalVars); + else + finalresults{lin_idx} = feval(funcHandle, optionalVars); + wh.addValueToStorageByLinIdx(finalresults{lin_idx}, 'ber', lin_idx); + end +end + +if options.parallel + %%% 4) Setup waitbar + h = waitbar(0, 'Processing Simulations...'); + updateWaitbar = @(~) waitbar(mean(arrayfun(@(f) strcmp(f.State, 'finished'), results)), h); + + fprintf('Fetching results... \n'); + + updateWaitbarFutures = afterEach(results, updateWaitbar, 0); + afterAll(updateWaitbarFutures, @(~) delete(h), 0); + + %%% 7) Fetch final results after all computations + fetchOutputs(results); + + for ridx = 1:length(results) + wh.addValueToStorageByLinIdx(results(ridx).OutputArguments{1}, 'ber', ridx); + end +end + +end diff --git a/test/duobinary_minimal_example.m b/test/duobinary_minimal_example.m index f817713..8866e69 100644 --- a/test/duobinary_minimal_example.m +++ b/test/duobinary_minimal_example.m @@ -1,11 +1,11 @@ useprbs = 1; -M = 6; +M = 2; randkey = 1; datarate = 448e9; fsym = round(datarate / log2(M)) ; %%%%% PRBS Generation in correct shape for Modulation Format %%%%%% -O = 15; %O of prbs +O = 16; %O of prbs N = 2^(O); %length of prbs [~,seed] = prbs(O,1); %initialize first seed of prbs bitpattern=[]; @@ -47,28 +47,79 @@ Tx_bits = Informationsignal(bitpattern); %%%%% Duobinary %%%%%% +precode = 1; +db_encode = 0; + close all Symbols_tx = PAMmapper(M,0).map(Tx_bits); Symbols_tx.fs = fsym; -Symbols1 = Duobinary().precode(Symbols_tx); -figure;histogram(Symbols1.signal); +%%% precode +if precode + Symbols0 = Duobinary().precode(Symbols_tx); +else + Symbols0 = Symbols_tx; +end -Symbols2 = Duobinary().encode(Symbols1); -figure;histogram(Symbols2.signal); +figure;histogram(Symbols0.signal); -Symbols3 = Duobinary().decode(Symbols2); -figure;histogram(Symbols3.signal); -autoArrangeFigures; +for n = 0:200 -Rx_bits = PAMmapper(M,0).demap(Symbols3); + if db_encode + Symbols1 = Duobinary().encode(Symbols0); + else + Symbols1 = Symbols0; + end -%%%%% Check BER of Bit Sequence %%%%%% + Symbols2 = Symbols1; + pos = 1; + if n~=0 + for pos = 1:n + po = randi(100); + a = Symbols2.signal(100+pos) == Symbols1.signal(100+po); + while a == 1 + po = po+1; + po = randi(100); + a = Symbols2.signal(100+pos) == Symbols1.signal(100+po); + end + Symbols2.signal(100+pos) = Symbols1.signal(100+po); + end + end -[~,error_num,ber,error_pos] = calc_ber(Tx_bits.signal,Rx_bits.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); + % disp(Symbols2.signal(100:100+pos)==Symbols1.signal(100:100+pos)) + + %%% encode + + if db_encode + + % figure;histogram(Symbols2.signal); + + Symbols3 = Duobinary().decode(Symbols2); + % figure;histogram(Symbols3.signal); + % autoArrangeFigures; + elseif precode + Symbols3 = Duobinary().encode(Symbols2); + Symbols3 = Duobinary().decode(Symbols3); + % figure;histogram(Symbols3.signal); + else + Symbols3 = Symbols2; + end + + Rx_bits = PAMmapper(M,0).demap(Symbols3); + + %%%%% Check BER of Bit Sequence %%%%%% + + [~,error_num(n+1),ber,error_pos] = calc_ber(Tx_bits.signal,Rx_bits.signal,"skip_front",10,"skip_end",10,"returnErrorLocation",1); + + % disp(['BER: ',sprintf('%.1E',ber),sprintf(' - Num. Err: %.1d',error_num(n+1)-2),' - - PAM-',num2str(M)]); + fprintf('n: %d - Num. Err: %.1d \n',n,error_num(n+1)); +end -disp(['BER: ',sprintf('%.1E',ber),' - - PAM-',num2str(M)]); -% + +figure() +hold on +scatter(1:length(Symbols3),Symbols3.signal,1,'.'); +scatter(error_pos,Symbols3.signal(error_pos),14,'o'); figure(3); diff --git a/test/matched_filter_minimal.m b/test/matched_filter_minimal.m new file mode 100644 index 0000000..b298fac --- /dev/null +++ b/test/matched_filter_minimal.m @@ -0,0 +1,61 @@ +%%% Run parameters +% TX +M = 4; +fsym = 32e9; + +apply_pulsef = 1; +fdac = 256e9; +fadc = 256e9; +random_key = 1; + +precomp = 0; +db_precode = 0; + +db_encode = 0; + +rcalpha = 0.05; +kover = 16; +vbias_rel = 0.5; +u_pi = 2.9; +vbias = -vbias_rel*u_pi; +laser_wavelength = 1293; +laser_linewidth = 0; +tx_bw_nyquist = 0.8; + +% 1) PRBS Generation +O = 18; %order of prbs +N = 2^(O-1); %length of prbs + +%%%%% MOVE-IT PRMS %%%% +Mi_prms = Moveit_wrapper("prms"); +if M == 6 + Mi_prms.para.bl = 2^(O-2); + Mi_prms.para.dimension = 5; +else + Mi_prms.para.bl = 2^(O-1); + Mi_prms.para.dimension = log2(M); %2.5bits/sym -> 2 bit/sym +end +Mi_prms.para.rand = 0; +Mi_prms.para.order = floor(O / log2(M)); +Mi_prms.para.skip =0; +Mi_prms.para.bruijn = 0; +Mi_prms.para.reset_prms = 0; +Mi_prms.para.method = 1; +bitpattern = Mi_prms.process([]); +if M == 6 + bitpattern = reshape(bitpattern',[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); +end + +bits = Informationsignal(bitpattern); + +symbols = PAMmapper(M,0).map(bits); +symbols.fs = fsym; + +Pform = Pulseformer("fsym",fsym,"fs",fsym,"alpha",0.6,"pulse","rrc","pulselength",64,"matched_sps",4,"output_sps",2); +Pform.process(symbols); + +MF = Pulseformer("fsym",fsym,"fdac",fdac,"pulse","rrc","pulselength",1024,"alpha",rcalpha,"matched",0); +Digi_sig = MF.process(Digi_sig); + +Digi_sig.spectrum("displayname",'Signal after shaping','fignum',1); diff --git a/test/test_mapping_eth.m b/test/test_mapping_eth.m new file mode 100644 index 0000000..592d9db --- /dev/null +++ b/test/test_mapping_eth.m @@ -0,0 +1,39 @@ + +% Setup PRBS parameters +O = 6; +M = 6; +N = 2^(O-1); % Length of PRBS +randkey = 1; % Random key for random stream +use_eth_mapping =1; + +if M ~= 6 + dimension = log2(M); +else + dimension = 5; +end + +[~, seed] = prbs(O, 1); % Initialize first seed of PRBS +bitpattern = []; + +s = RandStream('twister', 'Seed', randkey); +for i = 1:dimension + bitpattern(:, i) = randi(s, [0 1], N, 1); +end +if M == 6 + bitpattern = reshape(bitpattern',[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); +end + +Tx_bits = Informationsignal(bitpattern); + +Digi_Mod = PAMmapper(M, 0,"eth_style",use_eth_mapping); + +% Map bits to symbols +Symbols = Digi_Mod.map(Tx_bits); + +% Demap symbols back to bits +Rx_bits = Digi_Mod.demap(Symbols); + +[~, error_num, ber, ~] = calc_ber(Tx_bits.signal(1:length(Rx_bits.signal)), Rx_bits.signal,"skip_front", 0, "skip_end", 0, "returnErrorLocation", 1); + +fprintf('BER: %.1E \n',ber); \ No newline at end of file