diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index b601d40..1ae61d2 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -24,11 +24,11 @@ classdef Signal obj.signal = signal; obj.signal = obj.signal; obj.fs = options.fs; - + [~,obj.gitSHA] = system('git rev-parse HEAD'); [~,obj.gitStatus] = system('git status --porcelain'); % [~,obj.gitPatch] = system('git diff'); - + %%% Stuff for Logbook %%% SignalType = []; TimeStamp = []; @@ -133,6 +133,7 @@ classdef Signal end + %% function plot(obj, options) % signal to plot: obj.signal % fsamp : obj.fs (e.g. 92e9 => 92 GHz) @@ -144,6 +145,7 @@ classdef Signal options.displayname = []; options.timeframe = 0; options.clear = 0; + options.color = []; end figure(options.fignum); % If figure does not exist, create new figure @@ -169,8 +171,11 @@ classdef Signal end hold on; - plot(t, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1); - + if isempty(options.color) + plot(t, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1); + else + plot(t, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color); + end % 2 c) % - xlabel if not already here: time in readable format (1 ms and not 1e-3 s) % - ylabel amplitude @@ -268,7 +273,7 @@ classdef Signal Nase = [0]; SignalCopy = obj.signal; ModifierName = class(CallingModifier); - + cell = {SignalType , TimeStamp , Length , SignalPower(1) , Nase, SignalCopy, ModifierName, ModifierCopy, Description}; @@ -318,17 +323,17 @@ classdef Signal if options.fs_in == options.fs_out desc = ['No need to resample signal from ', num2str(options.fs_in*1e-9), ' GHz to ', num2str(options.fs_out*1e-9), ' GHz' ]; - + obj = obj.logbookentry(desc,obj); else obj.signal = resample(obj.signal,options.fs_out,options.fs_in,options.n,options.beta); - + desc = ['resample signal from ', num2str(options.fs_in*1e-9), ' GHz to ', num2str(options.fs_out*1e-9), ' GHz' ]; - + obj = obj.logbookentry(desc,obj); - + obj.fs = options.fs_out; end @@ -345,62 +350,77 @@ classdef Signal options.color = []; options.normalizeToNyquist = 0; options.normalizeTo0dB = 0; + options.max_num_lines = []; % Leave empty or omit to disable line rotation + options.fft_length = []; end - % spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname); - - N = 2^(nextpow2(length(obj.signal))-9); - - if options.normalizeToNyquist==0 - [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","psd","mean"); - w=w.*1e-9; - else - [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,"centered","psd","mean"); - % p_lin = smooth(p_lin,0.05,'rloess'); + if isempty(options.fft_length) + options.fft_length = 2^(nextpow2(length(obj.signal))-7); 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"); + 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"); + end + if options.normalizeTo0dB p_lin = p_lin./ max(p_lin); - p_dbm = 10*log10(p_lin); %dB to dBm in case of "power" + p_dbm = 10*log10(p_lin); % normalized to 0 dB ylab = "normalized to 0 dB"; else - p_dbm = 10*log10(p_lin); %dB to dBm in case of "power" - ylab = "Power (dBm)"; + p_dbm = 10*log10(p_lin); + ylab = "Power (dB/Hz)"; end - figure(options.fignum); % If figure does not exist, create new figure + figure(options.fignum); ax = gca; hold on if isempty(options.color) - plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1); + hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1); else - plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color); + hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color); end - if options.normalizeToNyquist==0 + % If user wants to limit the number of lines, check and remove old lines + if ~isempty(options.max_num_lines) && options.max_num_lines > 0 + allLines = findall(ax, 'Type', 'Line'); + if length(allLines) > options.max_num_lines + % Sort lines by creation order. Usually, the oldest lines appear first in allLines. + % If needed, you can sort by UserData or other criteria. + numToRemove = length(allLines) - options.max_num_lines; + delete(allLines(1:numToRemove)); + end + end + + if options.normalizeToNyquist == 0 xlabel("Frequency in GHz"); - % xlim([-obj.fs/2 obj.fs/2].*1e-9) edgetick = 2^(nextpow2(obj.fs*1e-9)); - xticks([-edgetick:16:edgetick]); - xlim([100*round( min(w)/100,1)-10,100*round( max(w)/100,1)+10]) + xticks(-edgetick:16:edgetick); + xlim([100*round(min(w)/100,1)-10, 100*round(max(w)/100,1)+10]) + xlim([-128 128]);%256GSa/s else xlabel("Normalized Frequency"); xlim([-pi, pi]); end - - ylabel("Power/frequency (dB/Hz)"); + ylabel(ylab); - - ylim([min(floor( min(p_dbm))-3 , ax.YLim(1)), max(ceil( max(p_dbm) )+(3), ax.YLim(2))]); - yticks([-200:10:10]); - grid on - grid minor - legend('Interpreter','none'); + + 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))]); + end + yticks(-200:10:10); + grid on; grid minor; + legend + % legend('Interpreter','none'); end + function move_it_spectrum(obj,options) arguments obj @@ -624,6 +644,7 @@ classdef Signal if options.mode == delay_mode.samples obj.signal=delayseq(obj.signal,delay); + % obj.signal=circshift(obj.signal,delay); elseif options.mode == delay_mode.time @@ -658,50 +679,63 @@ 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); + [pks,pkpos] = findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend'); shifts = lags(pkpos); - - %Cut occurences of ref signal from signal (only positive shifts) + + shifts = shifts(shifts>=0); + S = {}; - for c = shifts(shifts>=0) - sig = obj.delay(-c,'mode','samples'); - sig.signal = sig.signal(1:length(b)); - S{end+1,1} = sig; - end - - % isFlipped=0; - if all(sign(co(pkpos))) - isFlipped = 1; + + if numel(shifts) > 0 + + %Cut occurences of ref signal from signal (only positive shifts) + + for c = shifts(shifts>=0) + sig = obj.delay(-c,'mode','samples'); + sig.signal = sig.signal(1:length(b)); + S{end+1,1} = sig; + end + + % + + 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)); + obj.signal = S{idx}.signal; + %put signal with highest corr. to first index in S array + swap = S{1}; + S{1} = S{idx}; + S{idx} = swap; + + for c = 1:numel(shifts(shifts>0)) + S{c}.logbook = []; + end + + else + + %do nothing when shifts are negative or there are none... + end - %return/keep the sinal with the highest correlation (only within positive shifts) - [~,idx]=max(pks(shifts>=0)); - obj.signal = S{idx}.signal; - %put signal with highest corr. to first index in S array - swap = S{1}; - S{1} = S{idx}; - S{idx} = swap; - - for c = 1:numel(shifts(shifts>0)) - S{c}.logbook = []; - end - %plot all synced signals and the ref signal debug = 0; - if debug + if debug 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]); plot(b(1000:1100),'LineWidth',1); end end - + end - + %% function obj = filter(obj,a,b) lbdesc = ['Filtering signal with H = a: ',num2str(a),' / b: ',num2str(b)]; @@ -710,6 +744,7 @@ classdef Signal obj.signal = filter(a,b,obj.signal); end + %% function er = extinctionratio(obj,fsym,M) histpoints = 1024; %% verticale resolution histpoints = floor(histpoints/2)*2+1; %% to have the eye digram centered around one point make the vertical resolution uneven @@ -789,7 +824,7 @@ classdef Signal end - + %% function eye(obj,fsym,M,options) arguments @@ -920,85 +955,85 @@ classdef Signal try - hist_interest = plot_data(:,posxall); - hist_interest_smoth = smooth(hist_interest,20); - a = scatter(hist_interest_smoth+posxall,1:length(hist_interest_smoth),4,'.','MarkerEdgeColor','red'); + hist_interest = plot_data(:,posxall); + hist_interest_smoth = smooth(hist_interest,20); + a = scatter(hist_interest_smoth+posxall,1:length(hist_interest_smoth),4,'.','MarkerEdgeColor','red'); - [pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",10,"NPeaks",M,"MinPeakHeight",30,"MinPeakProminence",10); - - scatter(posxall,loc,'red','Marker','x','LineWidth',2); - - yline(loc,'Color','red','LineWidth',1,'LineStyle',':'); + [pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",10,"NPeaks",M,"MinPeakHeight",30,"MinPeakProminence",10); - for i = 1:numel(loc) - ppeak(i) = maxA - (difference/histpoints*loc(i)); - end + scatter(posxall,loc,'red','Marker','x','LineWidth',2); - oma = false; - if isa(obj,'Opticalsignal') - er=10*log10(ppeak(1)/ppeak(end)); - elseif isa(obj,'Electricalsignal') - if mean([ppeak(1),ppeak(end)]) < 1e-2 - oma = true; - er=max(ppeak)-min(ppeak); + yline(loc,'Color','red','LineWidth',1,'LineStyle',':'); + + for i = 1:numel(loc) + ppeak(i) = maxA - (difference/histpoints*loc(i)); + end + + oma = false; + if isa(obj,'Opticalsignal') + er=10*log10(ppeak(1)/ppeak(end)); + elseif isa(obj,'Electricalsignal') + if mean([ppeak(1),ppeak(end)]) < 1e-2 + oma = true; + er=max(ppeak)-min(ppeak); + else + er=10*log10(ppeak(1)/ppeak(end)); + end else er=10*log10(ppeak(1)/ppeak(end)); end - else - er=10*log10(ppeak(1)/ppeak(end)); - end - % Adjust position for the third box (slightly to the right) - boxPosition = [0.59 0.86 0.2 0.05]; % Adjusted position + % Adjust position for the third box (slightly to the right) + boxPosition = [0.59 0.86 0.2 0.05]; % Adjusted position - % Create third annotation box for Vmax - if ~oma - thirdboxstring = ['ER (db):',num2str(er),' dB']; - else - thirdboxstring = ['OMA outer:',num2str(er),' V']; - end + % Create third annotation box for Vmax + if ~oma + thirdboxstring = ['ER (db):',num2str(er),' dB']; + else + thirdboxstring = ['OMA outer:',num2str(er),' V']; + end - plot_infos = 0; - if plot_infos + plot_infos = 0; + if plot_infos - % Define properties - boxPosition = [0.15 0.86 0.2 0.05]; % Position for the first box [x y width height] - boxColor = [0.9 0.9 0.9]; % Light grey background color - boxEdgeColor = 'k'; % Black edge color - boxLineStyle = '--'; % Dashed line style - boxFontWeight = 'bold'; % Bold font + % Define properties + boxPosition = [0.15 0.86 0.2 0.05]; % Position for the first box [x y width height] + boxColor = [0.9 0.9 0.9]; % Light grey background color + boxEdgeColor = 'k'; % Black edge color + boxLineStyle = '--'; % Dashed line style + boxFontWeight = 'bold'; % Bold font - % Create first annotation box for Power - annotation('textbox', boxPosition, ... - 'String', ['Power: ',num2str(pwr_dbm),' dBm'], ... - 'BackgroundColor', boxColor, ... - 'EdgeColor', boxEdgeColor, ... - 'LineStyle', boxLineStyle, ... - 'FontWeight', boxFontWeight, ... - 'HorizontalAlignment', 'center'); + % Create first annotation box for Power + annotation('textbox', boxPosition, ... + 'String', ['Power: ',num2str(pwr_dbm),' dBm'], ... + 'BackgroundColor', boxColor, ... + 'EdgeColor', boxEdgeColor, ... + 'LineStyle', boxLineStyle, ... + 'FontWeight', boxFontWeight, ... + 'HorizontalAlignment', 'center'); - % Adjust position for the second box (slightly to the right) - boxPosition = [0.37 0.86 0.2 0.05]; % Adjusted position + % Adjust position for the second box (slightly to the right) + boxPosition = [0.37 0.86 0.2 0.05]; % Adjusted position - % Create second annotation box for PAPR - annotation('textbox', boxPosition, ... - 'String', ['PAPR(lin):',num2str(papr_),''], ... - 'BackgroundColor', boxColor, ... - 'EdgeColor', boxEdgeColor, ... - 'LineStyle', boxLineStyle, ... - 'FontWeight', boxFontWeight, ... - 'HorizontalAlignment', 'center'); + % Create second annotation box for PAPR + annotation('textbox', boxPosition, ... + 'String', ['PAPR(lin):',num2str(papr_),''], ... + 'BackgroundColor', boxColor, ... + 'EdgeColor', boxEdgeColor, ... + 'LineStyle', boxLineStyle, ... + 'FontWeight', boxFontWeight, ... + 'HorizontalAlignment', 'center'); - - annotation('textbox', boxPosition, ... - 'String',thirdboxstring , ... - 'BackgroundColor', boxColor, ... - 'EdgeColor', boxEdgeColor, ... - 'LineStyle', boxLineStyle, ... - 'FontWeight', boxFontWeight, ... - 'HorizontalAlignment', 'center'); - end + + annotation('textbox', boxPosition, ... + 'String',thirdboxstring , ... + 'BackgroundColor', boxColor, ... + 'EdgeColor', boxEdgeColor, ... + 'LineStyle', boxLineStyle, ... + 'FontWeight', boxFontWeight, ... + 'HorizontalAlignment', 'center'); + end end diff --git a/Classes/01_transmit/PAMmapper.m b/Classes/01_transmit/PAMmapper.m index c233a13..00ed0ca 100644 --- a/Classes/01_transmit/PAMmapper.m +++ b/Classes/01_transmit/PAMmapper.m @@ -23,7 +23,6 @@ classdef PAMmapper obj.scaling = rms(obj.get_levels()); - end function out = map(obj,signal_in) @@ -41,11 +40,25 @@ classdef PAMmapper end - function signalclass_out = demap(obj,signalclass_in) - signalclass_in.signal = obj.demap_(signalclass_in.signal); - lbdesc = ['Demap PAM ',num2str(obj.M),' symbols to bit stream']; - signalclass_in = signalclass_in.logbookentry(lbdesc,obj); - signalclass_out = signalclass_in; + function signal_out = demap(obj,signal_in) + issignalclass = 0; + + if isa(signal_in,'Signal') + signalclass = signal_in; + signal_in = signal_in.signal; + issignalclass = 1; + end + + signal_out = obj.demap_(signal_in); + + if issignalclass + lbdesc = ['Demap PAM ',num2str(obj.M),' symbols to bit stream']; + signal_in = signalclass; + signal_in = signal_in.logbookentry(lbdesc,obj); + signal_in.signal = signal_out; + signal_out = signal_in; + end + end function pam_sig = map_(obj,bitpattern) @@ -156,6 +169,8 @@ classdef PAMmapper 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 @@ -194,7 +209,6 @@ classdef PAMmapper end end - function [data_out] = demap_(obj,data_in) data_in= data_in'; if obj.M ~= 6 @@ -309,14 +323,35 @@ classdef PAMmapper end function [Signal_out] = quantize(obj,Signal_in) + constellation = obj.get_levels(); constellation = constellation ./ rms(constellation); - Signal_out = Signal_in; - dist = abs(Signal_in.signal - constellation); - [~,symbol_idx] = min(dist,[],2); % decision for closest constellation point - Signal_out.signal = constellation(symbol_idx); - Signal_out.signal = reshape(Signal_out.signal,size(Signal_in.signal)); + issignalclass = 0; + if isa(Signal_in,'Signal') + issignalclass = 1; + Sig_class = Signal_in; + Signal_in = Signal_in.signal; + end + + [~,high_dim_sig] = max(size(Signal_in)); + [~,high_dim_const] = max(size(constellation)); + + if high_dim_sig == high_dim_const + Signal_in = Signal_in'; + end + + dist = abs(Signal_in - constellation); + [~,symbol_idx] = min(dist,[],2); % decision for closest constellation point + Signal_out = constellation(symbol_idx); + Signal_out = reshape(Signal_out,size(Signal_in)); + + if issignalclass + Sig_class.signal = Signal_out; + Signal_out = Sig_class; + end + + end diff --git a/Classes/01_transmit/Pulseformer.m b/Classes/01_transmit/Pulseformer.m index a306ef8..c8e8537 100644 --- a/Classes/01_transmit/Pulseformer.m +++ b/Classes/01_transmit/Pulseformer.m @@ -4,6 +4,8 @@ classdef Pulseformer properties(Access=public) fdac + end + properties(Access=private) fsym pulse pulselength diff --git a/Classes/04_DSP/Coding/Duobinary.m b/Classes/04_DSP/Coding/Duobinary.m index 00f8e0c..e93c594 100644 --- a/Classes/04_DSP/Coding/Duobinary.m +++ b/Classes/04_DSP/Coding/Duobinary.m @@ -78,11 +78,12 @@ classdef Duobinary end - function signal = encode(~,signal) + function signal = encode(~,signal,options) arguments ~ signal + options.M = []; end if isa(signal,'Signal') @@ -93,17 +94,19 @@ classdef Duobinary data = data./rms(data); - u = unique(data); - M = numel(u); + if isempty(options.M) + u = unique(data); + options.M = numel(u); + end %make unipolar - if M == 4 + if options.M == 4 data = data .* sqrt(5); - elseif M == 6 + elseif options.M == 6 data = data .* sqrt(10); - elseif M == 8 + elseif options.M == 8 data = data .* sqrt(21); - elseif M == 16 + elseif options.M == 16 data = data .* sqrt(85); warning('Check if PAM16 implementation, mapping and scaling is correct!') end @@ -113,7 +116,7 @@ classdef Duobinary data = data - b; data = data ./ 2; - % assert(isequal((0:M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar + % assert(isequal((0:options.M-1)',unique(data)),'Check Duobinary Precoding'); %seems the signal is not unipolar % duobinary coding (1+D) % coeff = [1,1]; @@ -137,14 +140,14 @@ classdef Duobinary mean_power = sum((unique_points .^ 2) .* probabilities); scaling_factor = sqrt(mean_power); - if M == 4 + if 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 M == 6 + elseif options.M == 6 data = data ./ sqrt(5.8); - elseif M == 8 + elseif options.M == 8 data = data ./ sqrt(10.5); % 15-level constellation weighted with probability after DB code i.e. else - errormsg("Error in: Duobinary encode > scale unipolar to bipolar > The data is not PAM4, PAM6 or PAM 8? ") + error("Error in: Duobinary encode > scale unipolar to bipolar > The data is not PAM4, PAM6 or PAM 8? ") end if isa(signal,'Signal') @@ -165,15 +168,25 @@ classdef Duobinary end - function signalclass = decode(~,signalclass) + function signalclass = decode(~,signalclass,options) + + arguments + ~ + signalclass + options.M = []; + end data = signalclass.signal; u = unique(data); I = numel(u); %number of duobinary coded const. points - M = (I+1)/2; %PAM-M order + if isempty(options.M) + M = (I+1)/2; %PAM-M order + else + M = options.M; %PAM-M order + end %make unipolar - if I == 7 + if I == 7 || I == 6 data = data .* sqrt(2.5); elseif I == 11 data = data .* sqrt(5.8); diff --git a/Classes/04_DSP/Equalizer/Postfilter.m b/Classes/04_DSP/Equalizer/Postfilter.m index 8ec5280..1fed009 100644 --- a/Classes/04_DSP/Equalizer/Postfilter.m +++ b/Classes/04_DSP/Equalizer/Postfilter.m @@ -4,8 +4,8 @@ classdef Postfilter < handle properties(Access=public) ncoeff = 2; - burg_coeff = []; - + coefficients = []; + useBurg = NaN end methods (Access=public) @@ -14,7 +14,9 @@ classdef Postfilter < handle % Detailed explanation goes here arguments - options.ncoeff double = 2 + options.useBurg = NaN + options.coefficients = [] + options.ncoeff double = 2 end % @@ -24,15 +26,37 @@ classdef Postfilter < handle obj.(fn{n}) = options.(fn{n}); end end - + % do more stuff end - function signalclass_out = process(obj,signalclass_in,noiseclass_in) + function signalclass_out = process(obj,signalclass_in,noiseclass_in,options) - obj.burg_coeff = arburg(noiseclass_in.signal,obj.ncoeff); - signalclass_in = signalclass_in.filter(obj.burg_coeff,1); + arguments + obj + signalclass_in + noiseclass_in + options.useBurg = obj.useBurg; + options.coefficients = obj.coefficients; + end + + if ~isnan(options.useBurg) && options.useBurg + + disp('using burg alg') + obj.coefficients = arburg(noiseclass_in.signal,obj.ncoeff); + + elseif ~isempty(options.coefficients) + + disp('using given taps') + obj.coefficients = options.coefficients; + obj.useBurg = 0; + + else + + end + + signalclass_in = signalclass_in.filter(obj.coefficients,1); % append to logbook lbdesc = ['Postfilter']; @@ -43,15 +67,29 @@ classdef Postfilter < handle end - function showFilter(obj,noiseclass_in) + function showFilter(obj,noiseclass_in,options) - noiseclass_in.spectrum('displayname','Noise PSD shifted to 0dBm','fignum',121,'normalizeTo0dB',1); + arguments + obj + noiseclass_in + options.fignum = 121 + options.color = [] + end - [h,w] = freqz(1,obj.burg_coeff,length(noiseclass_in),"whole",noiseclass_in.fs); + % noiseclass_in.spectrum('displayname','Noise PSD shifted to 0dBm','fignum',options.fignum,'normalizeTo0dB',1); + + figure(options.fignum) + [h,w] = freqz(1,obj.coefficients,length(noiseclass_in),"whole",noiseclass_in.fs); h = h/max(abs(h)); hold on w_ = (w - noiseclass_in.fs/2); - plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(obj.burg_coeff), ' ']); + + if isempty(options.color) + plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(obj.coefficients,2)), ' '],'LineWidth',2); + else + plot(w_.*1e-9,20*log10(fftshift(abs(h))),'DisplayName',['Burg Coeffs: ', num2str(round(obj.coefficients,2)), ' '],'Color',options.color,'LineWidth',2); + end + ylim([-30,2]); end diff --git a/Classes/04_DSP/Equalizer/VNLE.m b/Classes/04_DSP/Equalizer/VNLE.m index 397d58e..d95bb98 100644 --- a/Classes/04_DSP/Equalizer/VNLE.m +++ b/Classes/04_DSP/Equalizer/VNLE.m @@ -57,7 +57,7 @@ classdef VNLE < handle end - function [X] = process(obj, X, D) + function [X,N] = process(obj, X, D) % actual processing of the signal (steps 1. - 3.) % 1 normalize RMS @@ -91,6 +91,9 @@ classdef VNLE < handle lbdesc = [num2str(obj.order),' tap FFE']; X = X.logbookentry(lbdesc); % append to logbook + N = X; + N = X - D; + end @@ -156,7 +159,6 @@ classdef VNLE < handle end if ~all(mu==0,'all') %mu has not only zeros - % obj.e = obj.e - (mu * err * x_in) ; % Weight update rule of LMS obj.e = obj.e - ( (mu * x_in) * err ) ; % Weight update rule of LMS else normalizationfactor = (x_in.' * x_in); diff --git a/Classes/04_DSP/Sequence Detection/MLSE.m b/Classes/04_DSP/Sequence Detection/MLSE.m index 3a780da..59b9c33 100644 --- a/Classes/04_DSP/Sequence Detection/MLSE.m +++ b/Classes/04_DSP/Sequence Detection/MLSE.m @@ -34,17 +34,22 @@ classdef MLSE < handle end - function signalclass = process(obj,signalclass) + function [signalclass_hd,signalclass_sd] = process(obj,signalclass,ref_symbolclass) data_in = signalclass.signal; + data_ref = ref_symbolclass.signal; - data_out = obj.process_(data_in); + [data_out_hd,data_out_sd] = obj.process_(data_in,data_ref); - signalclass.signal = data_out; + signalclass_hd = signalclass; + signalclass_hd.signal = data_out_hd; + + signalclass_sd = signalclass; + signalclass_sd.signal = data_out_sd; end - function data_out = process_(obj,data_in) + function [VITERBI_ESTIMATION_SYMBOLS,soft_decisions] = process_(obj,data_in,data_ref) % remove unnecessary zeros at start of impulse response to keep @@ -58,6 +63,17 @@ classdef MLSE < handle obj.DIR = [0 obj.DIR]; end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%% PREPARATIONS %%%%%%%% + + %%%% 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); + + + % impulse respnse i.e. [0.5, 1.0000] obj.DIR = flip(obj.DIR); % RMS normalization of input data @@ -77,6 +93,7 @@ classdef MLSE < handle % Calculate all possible input symbols for the desired impulse % response. Row number is the index of the previous state, % column number is the index of the next state + % noise free received == branch metrics noise_free_received = zeros(length(states),length(states)); count_row = 1; count_col = 1; @@ -101,127 +118,302 @@ classdef MLSE < handle end end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% FORWARD PASS %%%%% - - % - % %% Optimized - % % Preallocate and initialize variables - % data_out = NaN(size(data_in)); % output vector - % sum_path_metrics_res = zeros(length(states), length(data_in)); - % path_idx = zeros(length(states), length(data_in)); - % - % % Precompute repmat size - % num_states = length(states); - % num_signals = numel(noise_free_received); - % - % % First trellis path (initialize) - % sum_path_metrics = zeros(num_states, num_states); - % path_metrics = (abs(data_in(1) - noise_free_received)).^2; % Use broadcasting instead of repmat - % sum_path_metrics = sum_path_metrics + path_metrics; % Compute initial path metrics - % - % [sum_path_metrics_res(:,1), path_idx(:,1)] = min(sum_path_metrics, [], 2); % Best path for first step - % - % % Preallocate path_metrics and sum_path_metrics to avoid reallocating in each loop - % path_metrics = zeros(num_states, num_signals); - % - % % Loop over remaining trellis paths - % for n = 2:length(data_in) - % % Avoid reallocation of sum_path_metrics, reuse the same matrix and update - % previous_sum_path_metrics = sum_path_metrics_res(:,n-1).'; % Transpose once for broadcasting - % sum_path_metrics = repmat(previous_sum_path_metrics, num_states, 1); % Avoid dynamic resizing - % - % % Calculate path metrics using broadcasting - % path_metrics = (abs(data_in(n) - noise_free_received)).^2; % Avoid repmat - % - % % Update sum path metrics - % sum_path_metrics = sum_path_metrics + path_metrics; - % - % % Find the best path for each state and store results - % [sum_path_metrics_res(:,n), path_idx(:,n)] = min(sum_path_metrics, [], 2); - % end - % - % %% Traceback - % ideal_path = NaN(1, length(data_in)+1); % Preallocate ideal path - % [~, ideal_path(length(data_in)+1)] = min(sum_path_metrics_res(:,length(data_in))); % Start from final state - % - % % Trace back through trellis - % for h = length(data_in):-1:1 - % ideal_path(h) = path_idx(ideal_path(h+1),h); % Follow the best path back - % end - % - % % Extract the output indices - % idx_out = ideal_path(2:end); + % Initialize the output vector + pm = zeros(length(states),length(states)); + bm_fw = zeros(length(states),length(states),length(data_in)); + + % Forward Recursion (FSM Computation) + for n = 1:length(data_in) + bm = abs(data_in(n) - noise_free_received).^2; + pm = pm + bm; + [pm_survivor_fw(:,n),pm_survivor_fw_idx(:,n)] = min(pm,[],2); % choose lowest path metric as new state + pm = repmat(pm_survivor_fw(:,n).',length(states),1); % update pm (chosen state to 2nd dimension -> FROM state) + bm_fw(:,:,n) = bm; - % OLD - % initilaize the output vector - data_out = NaN(size(data_in)); - - sum_path_metrics = zeros(length(states),length(states)); - - % first trellis path - % euclidian distance as path metric - path_metrics = (abs(repmat(data_in(1),size(noise_free_received)) - noise_free_received)).^2; - sum_path_metrics = sum_path_metrics + path_metrics; % calculation of all possible sum path metrics - [sum_path_metrics_res(:,1),path_idx(:,1)] = min(sum_path_metrics,[],2); % find the best path to each state, store sum path metric and predecessor for each state - - % remaining trellis paths - for n = 2:length(data_in) - sum_path_metrics = repmat(sum_path_metrics_res(:,n-1).',length(states),1); - % path_metrics = (repmat(data_in(n),size(noise_free_received)) - noise_free_received).^2;%.*prob_mat; - path_metrics = (abs(repmat(data_in(n),size(noise_free_received)) - noise_free_received)).^2; - sum_path_metrics = sum_path_metrics + path_metrics; - [sum_path_metrics_res(:,n),path_idx(:,n)] = min(sum_path_metrics,[],2); end - %% trace back - ideal_path = NaN(1,length(data_in)+1); - % find ideal trellis path by going through the trellis - % backwards - [~,ideal_path(length(data_in)+1)] = min(sum_path_metrics_res(:,length(data_in))); + % we can now get the best path as min + best_fw_path = NaN(1,length(data_in)+1); + + % find ideal trellis path by going through the trellis backwards + [~,best_fw_path(length(data_in)+1)] = min(pm_survivor_fw(:,length(data_in))); + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% BACKWARD PASS %%%%% + + % Initialize the output vector + pm = zeros(length(states),length(states)); + pm_survivor_bw = zeros(length(states),length(data_in)+1); + pm_survivor_bw_idx = zeros(length(states),length(data_in)+1); + bm_bw = zeros(length(states),length(states),length(data_in)+1); % starting with the state that has the lowest sum path % metric, follow the stored information about the % predecessor for h = length(data_in):-1:1 - ideal_path(h) = path_idx(ideal_path(h+1),h); + + best_fw_path(h) = pm_survivor_fw_idx(best_fw_path(h+1),h); + + bm = abs(data_in(h) - noise_free_received).^2; + pm = pm + bm.'; + [pm_survivor_bw(:,h),pm_survivor_bw_idx(:,h)] = min(pm,[],2); % choose lowest path metric as new state + pm = repmat(pm_survivor_bw(:,h).',length(states),1); % update pm (chosen state to 2nd dimension -> FROM state) + + bm_bw(:,:,h) = bm; + end - idx_out = ideal_path(2:length(data_in)+1); + VITERBI_ESTIMATION_IDX(1:length(data_in)) = first_sym(best_fw_path(2:end)); + VITERBI_ESTIMATION_SYMBOLS(1:length(data_in)) = constellation(best_fw_path(2:end)); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%% FORWARD PASS %%%%% + + %calc the log probabilities (llp's) + + for k = 2:length(data_in) + + fsm = repmat(pm_survivor_fw(:,k-1)',[length(states),1]); %pm_survivor_fw size: length(states)xlength(sequence) + bm = bm_fw(:,:,k); %bm_fw size: length(states)xlength(states)xlength(sequence) + bsm = repmat(pm_survivor_bw(:,k)',[length(states),1]); %pm_survivor_bw size: length(states)xlength(sequence)+1 + + llp(:,k) = min(fsm+bm+bsm,[],2); + + end + + % 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'); + [~,~,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); + % fprintf('LLP BER +1: %.2e \n',ber_llp); + % [~,~,ber_llp,~] = calc_ber(circshift(rx_bits,-1),tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + % fprintf('LLP BER -1: %.2e \n',ber_llp); + + %%%% DECIDE based on Viterbi traceback + rx_bits = PAMmapper(numel(constellation),0).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); + % fprintf('Viterbi BER: %.2e \n',ber_viterbi); + + % 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'); + [~,~,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'); + [~,~,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); + % fprintf('BW BER: %.2e \n',ber_viterbi); + % [~,~,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); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + tx_symbolpos = zeros(numel(constellation),length(data_ref)); + est_symbolpos = zeros(numel(constellation),length(data_ref)); + for lvl = 1:numel(constellation) + tx_symbolpos(lvl,data_ref==constellation(lvl)) = 1; + est_symbolpos(lvl,VITERBI_ESTIMATION_IDX==first_sym(lvl)) = 1; + est_symbolpos_llm(lvl,llp_based_state_seq==lvl) = 1; + end + est_symbolpos= logical(est_symbolpos); + tx_symbolpos = logical(tx_symbolpos); + est_symbolpos_llm = logical(est_symbolpos_llm); + figure(299);clf;hold on; + llp_filt = NaN(length(states),length(data_in)); + llp_ = llp - min(llp,[],1); + for c = 2%:numel(constellation) + for c2 = 1:3%numel(constellation) - if obj.duobinary_output - %use duobinary encoder, output is already scaled inside this - %one - data_out = Duobinary().encode(first_sym(idx_out)); + idx = est_symbolpos_llm(c,:); + llp_filt(c2,idx) = (llp(c2,idx)); + + plotidx = 2:200000; + scatter(plotidx,llp_filt(c2,plotidx),1,'o','LineWidth',1,'DisplayName',sprintf('LLP of Symbol %d | %d was transmitted',c2,c)); - else - % - data_out(1:length(data_in)) = first_sym(idx_out); - % scale to rms = 1 using the standard sqrt() expressions - if obj.M == 4 - data_out = data_out./sqrt(5); - elseif obj.M == 6 - data_out = data_out./sqrt(10); - elseif obj.M == 8 - data_out = data_out./sqrt(21); end + end + + % Assume llp is a 4 x N matrix (4 PAM-4 symbols, N time steps) + [numSymbols, numTimeSteps] = size(llp); + P_symbols = zeros(numSymbols, numTimeSteps); % To hold the soft output probabilities + + for k = 1:numTimeSteps + + % Compute the exponentials. (Use -llp_shifted if llp's are costs.) + exponents = -llp_(:, k); + + % Normalize to form a probability vector. + P_symbols(:, k) = exponents / sum(exponents); + end + + % Optionally, if you prefer a single soft output value per symbol (an expectation), + % define your PAM-4 constellation levels, e.g.: + soft_output = sum(P_symbols .* constellation, 1); % 1 x N vector of soft outputs + + + %%%% BITWISE LLR's ??? %%%%% + figure(100); + clf + hold on + title('LLR between inner and outer bits') + bitmapping = PAMmapper(numel(constellation),0).demap(constellation); + for bp = 1:size(bitmapping,2) + pos_bitone = find(bitmapping(:,bp)==1); + pos_bitzero = find(bitmapping(:,bp)~=1); + + llr_bits(bp,:) = min(llp(pos_bitone,:),[],1) - min(llp(pos_bitzero,:),[],1); + + subplot(size(bitmapping,2),1,bp) + scatter(1:length(data_ref),llr_bits(bp,:),1,'.'); + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % From: Log-Likelihood Probabilities LLP --> To: Log-Likelihood Ratios LLR + for s = 1:length(states)-1 + llr(s,:) = llp_(s,:) - llp_(s+1,:); % subtract llp's of successive const. points to get llr's + end + + % Build Sets of LLR's that contain only those values at + % timepoint k, where the symbols: + % a) were actually transmitted: set "S" + % b) were decoded by viterbi: set "SC" + S = NaN(length(states)-1,length(data_in)); + SC = NaN(length(states)-1,length(data_in)); + + llp_inf = llp_; + llp_inf(llp_==0) = Inf; + [~,scnd_idx] = min(llp_inf,[],1) ; + for s = 1:length(states)-1 + + %%% SET "S" + % find all time indices k, where const point s was actually transmitted + indice = tx_symbolpos(s,:)==1; + S(s,indice) = llr(s,indice); + %calc mean of all llr's where symbol s was transmitted + K(s,1) = mean(S(s,indice),'omitnan'); + + % find all time indices k, where const. point s+1 was actually transmitted + indice_plusone = tx_symbolpos(s+1,:)==1; + S(s,indice_plusone) = llr(s,indice_plusone); + K(s,2) = mean(S(s,indice_plusone),'omitnan'); + + + %%% SET "SC" + % find all time indices k, where symbol s was decoded + % idx: find positions in time where a symbol s was decoded + idx = est_symbolpos_llm(s,:); + % idx 2: find positions where the strongest competitor is + % s+1, i.e. the second best llp is at s+1 + idx2 = scnd_idx == s+1; + SC(s,idx&idx2) = llr(s,idx&idx2); + KC(s,1) = mean(SC(s,idx&idx2),'omitnan'); + + % find all time indices k, where symbol s+1 was decoded + % idx: find positions in time where a symbol s+1 was decoded + idx = est_symbolpos_llm(s+1,:); + idx2 = scnd_idx == s; + % idx 2: find positions where the strongest competitor is + % s, i.e. the second best llp is at s + SC(s,idx&idx2) = llr(s,idx&idx2); + KC(s,2) = mean(SC(s,idx&idx2),'omitnan'); + + % scale the sets, using the average (?) of the + llrcn(s,:) = SC(s,:) * ( constellation(s+1)-constellation(s) ) ./ ( K(s,2)-K(s,1)) + decisionLevels(s); + + end + + figure(111) + title("Bla") + clf + for s = 1:length(states)-1 + figure(1111) + hold on + title(sprintf('llrcn = llp %d - llp %d',s, s+1,s, s+1)); + scatter(1:length(llrcn),llrcn(s,:),1,'.'); + + + % STUFENARTIGE LLR'S UND LLP'S %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + figure(111) + + subplot(1,length(states),s) + hold on + title(sprintf('llr = llp %d - llp %d \n SC=llr(tx sym = %d or %d)',s, s+1,s, s+1)); + scatter(1:length(llr),llr(s,:),1,'.'); + scatter(1:length(SC),SC(s,:),1,'.'); + yline([K(s,1),K(s,2)]); + low = min(min(llr)); + hi = max(max(llr)); + ylim([low,hi]); + + subplot(1,length(states),length(states)) + hold on + scatter(1:length(llr),llr(s,:),1,'.'); + ylim([low,hi]); + + + % HISTOGRAMME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + figure(112) + subplot(length(states),1,s) + hold on + title(sprintf('histogram of llr; between const point %d - %d',s, s+1)); + histogram(llrcn(s,:),10000,'EdgeAlpha',0) + % histogram(SC(s,:),10000,'EdgeAlpha',0) + xlim([-20, 20]); + subplot(length(states),1,length(states)) + hold on + histogram(llrcn(s,:),10000,'EdgeAlpha',0) + % histogram(SC(s,:),10000,'EdgeAlpha',0) + % xlim([-20, 20]); end + softdecisions = mean(llrcn,1,'omitnan'); + + distance = abs(VITERBI_ESTIMATION_SYMBOLS-llrcn); + [win_cost,win_idx] = min(distance,[],1); + + for i = 1:length(data_in) + soft_decisions(i) = llrcn(win_idx(i),i); + end + + soft_decisions = max(llrcn,[],1); + + showLevelHistogram(soft_decisions,data_ref) + + figure() + % scatter(1:length(data_in),VITERBI_ESTIMATION_SYMBOLS,1,'.'); + scatter(1:length(data_in),soft_decisions,1,'.'); + + + MLM_ESTIMATION(1:length(data_in)) = PAMmapper(numel(constellation),0).quantize(soft_decisions)'; + rx_bits = PAMmapper(numel(constellation),0).demap(MLM_ESTIMATION'); + [~,~,ber_mlm,~] = calc_ber(rx_bits,tx_bits,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + fprintf('MLM BER: %.2e \n',ber_mlm); end end methods (Access=private) - % Cant be seen from outside! So put all your functions here that can/ - % shall not be called from outside - end end diff --git a/Classes/DataBaseHandler/DBHandler.m b/Classes/DataBaseHandler/DBHandler.m index 6db91a3..1495209 100644 --- a/Classes/DataBaseHandler/DBHandler.m +++ b/Classes/DataBaseHandler/DBHandler.m @@ -430,9 +430,10 @@ classdef DBHandler < handle baseQuery = [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 ' ... - 'LEFT JOIN Equalizer ON BERs.eq_id = Equalizer.eq_id ' ... 'WHERE ']; + % 'LEFT JOIN BERs ON Runs.run_id = BERs.run_id ' ... + % 'LEFT JOIN Equalizer ON BERs.eq_id = Equalizer.eq_id ' ... + % Loop through each table in filterParams filterClauses = []; diff --git a/Datatypes/db_mode.m b/Datatypes/db_mode.m index 728f38c..93421b9 100644 --- a/Datatypes/db_mode.m +++ b/Datatypes/db_mode.m @@ -1,10 +1,11 @@ classdef db_mode < int32 enumeration - no_db (0) - db_precoded (1) - db_encoded (2) - db_emulate (3) + no_db (0) + db_precoded (1) %sequence is precoded + db_encoded (2) %sequence is precoded and encoded + db_emulate (3) %for eq'ing: emulate precode on sequence that was not precoded at tx + db_discard (4) %for eq'ing: discard precode on sequence that was precoded at tx 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 3f2e8be..73f131c 100644 --- a/Functions/EQ_structures/duobinary_signaling.m +++ b/Functions/EQ_structures/duobinary_signaling.m @@ -1,10 +1,11 @@ -function [eq_signal,eq_noise,ber,numErrors] = duobinary_signaling(EQ, MLSE,M ,rx_signal, tx_symbols, tx_bits) +function [eq_package] = duobinary_signaling(eq_, mlse_,M ,rx_signal, tx_symbols, tx_bits) %Duobinary Signaling - [eq_signal, eq_noise] = EQ.process(rx_signal,tx_symbols); + [eq_signal, eq_noise] = eq_.process(rx_signal,tx_symbols); - eq_signal = MLSE.process(eq_signal); + eq_signal = mlse_.process(eq_signal); + eq_signal = Duobinary().encode(eq_signal); eq_signal = Duobinary().decode(eq_signal); % M = numel(unique(eq_signal.signal)); @@ -12,4 +13,6 @@ function [eq_signal,eq_noise,ber,numErrors] = duobinary_signaling(EQ, MLSE,M ,rx [~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + 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 93ad6bf..c2d760e 100644 --- a/Functions/EQ_structures/duobinary_target.m +++ b/Functions/EQ_structures/duobinary_target.m @@ -1,17 +1,77 @@ -function [eq_signal,eq_noise,ber,numErrors] = duobinary_target(EQ, MLSE,M, rx_signal, tx_symbols, tx_bits) +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 + options.precode_mode db_mode + options.showAnalysis = 0; +end %Duobinary Targeting - [eq_signal, eq_noise] = EQ.process(rx_signal,Duobinary().encode(tx_symbols)); + [eq_signal, eq_noise] = eq_.process(rx_signal,Duobinary().encode(tx_symbols)); % dir = [1,1]; - eq_signal = MLSE.process(eq_signal); + mlse_sig_sd = mlse_.process(eq_signal); - eq_signal = Duobinary().decode(eq_signal); + 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) + + % takes: + % -> eq_signal_hd: hard decision signal after eq + % -> tx_symbols: that where used as reference for eq + + switch options.precode_mode + case db_mode.db_emulate + + mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"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); + + case db_mode.db_discard + + % normal dsp for precoded sequence == discard/omit/ignore precode + tx_bits = PAMmapper(M,0).demap(tx_symbols); + + case db_mode.db_encoded + + % normal DB encoded data (only for 10KM) + + case db_mode.db_precoded + + mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M); + + end % M = numel(unique(tx_symbols.signal)); - rx_bits = PAMmapper(M,0).demap(eq_signal); + rx_bits = PAMmapper(M,0).demap(mlse_sig_hd); [~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + + eq_package.ber = ber; + + if options.showAnalysis + eq_noise = eq_noise - mean(eq_noise.signal); + +rx_signal.spectrum("normalizeTo0dB",1,"fignum",250); + + showEQNoisePSD(eq_noise,"fignum",250,"displayname",'Duobinary Target Noise'); + + Duobinary().encode(tx_symbols).spectrum("normalizeTo0dB",1,"fignum",250); + end + + end \ No newline at end of file diff --git a/Functions/EQ_structures/vnle.m b/Functions/EQ_structures/vnle.m index c4e266c..ba705d7 100644 --- a/Functions/EQ_structures/vnle.m +++ b/Functions/EQ_structures/vnle.m @@ -1,4 +1,4 @@ -function [eq_signal,eq_noise,ber,numErrors] = vnle(EQ,M,rx_signal,tx_symbols,tx_bits,emulate_precode) +function [eq_package] = vnle(eq_,M,rx_signal,tx_symbols,tx_bits,options) %VNLE Apply an equalization algorithm to the received signal and calculate BER % This function takes an equalizer object, a received signal, and the % transmitted symbols to apply equalization, map the received signal back to bits, @@ -14,23 +14,73 @@ function [eq_signal,eq_noise,ber,numErrors] = vnle(EQ,M,rx_signal,tx_symbols,tx_ % ber - Bit error rate after equalization % numErrors - Number of bit errors detected - %FFE or VNLE - [eq_signal,eq_noise] = EQ.process(rx_signal,tx_symbols); - - eq_signal = PAMmapper(M,0).quantize(eq_signal); - - if emulate_precode - eq_signal = Duobinary().encode(eq_signal); - eq_signal = Duobinary().decode(eq_signal); - - tx_symbols= Duobinary().encode(tx_symbols); - tx_symbols = Duobinary().decode(tx_symbols); - tx_bits = PAMmapper(M,0).demap(tx_symbols); + arguments + eq_ + M + rx_signal + tx_symbols + tx_bits + options.precode_mode db_mode + options.showAnalysis = 0 end - - % M = numel(unique(tx_symbols.signal)); - rx_bits = PAMmapper(M,0).demap(eq_signal); + + %FFE or VNLE + [eq_signal_sd,eq_noise] = eq_.process(rx_signal,tx_symbols); + + eq_signal_hd = PAMmapper(M,0).quantize(eq_signal_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) + + % takes: + % -> eq_signal_hd: hard decision signal after eq + % -> tx_symbols: that where used as reference for eq + + 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); + + tx_symbols_precoded = Duobinary().encode(tx_symbols); + tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded); + + tx_bits = PAMmapper(M,0).demap(tx_symbols_precoded); + + case db_mode.db_discard + + % normal dsp for precoded sequence == discard/omit/ignore precode + tx_bits = PAMmapper(M,0).demap(tx_symbols); + + case db_mode.db_encoded + + % normal DB encoded data (only for 10KM) + + 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); + + end + + rx_bits = PAMmapper(M,0).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; + + if options.showAnalysis + + fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl); + + fprintf('VNLE BER: %.2e \n',ber); + + 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 83a8b83..ac2493c 100644 --- a/Functions/EQ_structures/vnle_postfilter_mlse.m +++ b/Functions/EQ_structures/vnle_postfilter_mlse.m @@ -1,21 +1,125 @@ -function [eq_signal,eq_noise,ber,numErrors] = vnle_postfilter_mlse(eq_,pf_,mlse_,M,rx_signal,tx_symbols,tx_bits) - +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 + options.precode_mode db_mode + options.showAnalysis = 0; +end + %FFE or VNLE - [eq_signal,eq_noise] = eq_.process(rx_signal,tx_symbols); + [eq_signal_sd,eq_noise] = eq_.process(rx_signal,tx_symbols); - % eq_noise.signal = eq_noise.signal - movmean(eq_noise.signal,[5000,0]); + eq_signal_hd = PAMmapper(M,0).quantize(eq_signal_sd); - eq_signal = pf_.process(eq_signal,eq_noise); + mlse_sig_sd = pf_.process(eq_signal_sd,eq_noise); - %M = numel(unique(tx_symbols.signal)); - mlse_.DIR = pf_.burg_coeff; - mlse_.trellis_states = PAMmapper(M,0).levels; - mlse_.M = M; + 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); - eq_signal = mlse_.process(eq_signal); + 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) + + % takes: + % -> M + % -> eq_signal_hd: hard decision signal after eq + % -> tx_symbols: that where used as reference for eq + + 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); + + mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"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); + + case db_mode.db_discard + + % normal dsp for precoded sequence == discard/omit/ignore precode + tx_bits = PAMmapper(M,0).demap(tx_symbols); + + case db_mode.db_encoded + + % normal DB encoded data (only for 10KM) + + 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); + + mlse_sig_hd = Duobinary().encode(mlse_sig_hd,"M",M); + mlse_sig_hd = Duobinary().decode(mlse_sig_hd,"M",M); + + end + + % 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); - rx_bits = PAMmapper(M,0).demap(eq_signal); + % 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); + + % 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); + + + 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_; + + + if options.showAnalysis + + % fprintf(['VNLE EVM lvl: ',repmat('%.3f ',1,numel(evm_lvl)),' \n'],evm_lvl); + + fprintf('VNLE BER: %.2e \n',ber_vnle); + + fprintf('MLSE BER: %.2e \n',ber_mlse); + + showEQNoisePSD(eq_noise,"fignum",336,"displayname",'VNLE+DFE','postfilter_taps',pf_.coefficients); + + + rx_signal.spectrum("normalizeTo0dB",1,"fignum",337,"displayname",'Rx Signal'); + + tx_symbols.spectrum("normalizeTo0dB",1,"fignum",337,'displayname','Tx Signal'); + + showLevelHistogram(eq_signal_sd,tx_symbols) + % showLevelHistogram(mlse_sig_sd,tx_symbols) + + showEQcoefficients(eq_.e,eq_.e2,eq_.e3,"displayname",'Coefficients'); + + showEQNoiseSNR(tx_symbols,eq_noise,"displayname",'vnle snr','fignum',101); + + %%% EQ SNR Spectrum %230 + %snr + snr_vnle = snr(tx_symbols.signal,eq_noise.signal); + + % showErrorBurstCount(eq_signal_sd,tx_symbols) + + + end - [~,numErrors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - end \ No newline at end of file diff --git a/Functions/beautifyBERplot.m b/Functions/beautifyBERplot.m index a50f626..690c4ba 100644 --- a/Functions/beautifyBERplot.m +++ b/Functions/beautifyBERplot.m @@ -1,9 +1,5 @@ function beautifyBERplot() % BEAUTIFYBERPLOT Enhances a BER plot for publication-quality figures. - % - % This function automatically adjusts the aesthetics of an existing plot. - % It sets limits based on the data already plotted and modifies the - % appearance for publication-ready figures. % Set line properties for all current plot lines lines = findall(gca, 'Type', 'Line'); @@ -11,10 +7,10 @@ function beautifyBERplot() num_markers = length(markers); for i = 1:length(lines) - lines(i).LineWidth = 1.5; % Set a thicker line width - lines(i).LineStyle = '-'; % Solid lines for simplicity + 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 - lines(i).MarkerSize = 6; % Set marker size + lines(i).MarkerSize = 4; % Marker size lines(i).MarkerFaceColor = 'auto'; % Use line color for marker face end @@ -24,34 +20,16 @@ function beautifyBERplot() % Set figure background to white set(gcf, 'Color', 'w'); - % Set axis labels with LaTeX - xlabel('Bit Rate in Gbps', 'Interpreter', 'latex'); - ylabel('Bit Error Rate', 'Interpreter', 'latex'); - - % Determine x and y limits from the data in the plot - x_data = []; - y_data = []; - for i = 1:length(lines) - x_data = [x_data; lines(i).XData(:)]; %#ok % Append x data - y_data = [y_data; lines(i).YData(:)]; %#ok % Append y data - end - x_min = min(x_data); - x_max = max(x_data); - y_min = 1e-4; % Ensure no negative or zero values for log scale - y_max = 0.5; - - % Adjust axis ranges - xlim([x_min, x_max]); - ylim([y_min, y_max]); - - % Set logarithmic scale for y-axis - set(gca, 'YScale', 'log'); + + % 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'); % Customize grid and box appearance - set(gca, 'Box', 'on', 'LineWidth', 1.2); % Add a thicker border + set(gca, 'Box', 'on', 'LineWidth', 0.8); % Thicker border grid on; - grid minor; + % grid minor; % Adjust font size and style for better readability - set(gca, 'FontSize', 12, 'FontName', 'Times New Roman'); + set(gca, 'FontSize', 10, 'FontName', 'Times New Roman'); end diff --git a/Functions/calc_ber.m b/Functions/calc_ber.m deleted file mode 100644 index 142ac1a..0000000 --- a/Functions/calc_ber.m +++ /dev/null @@ -1,62 +0,0 @@ -function [bits,errors,ber,errorIndice] = calc_ber(data_in,data_ref,options) - -arguments(Input) - data_in; - data_ref; - 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)0; - - % Parse optional input arguments - if ~isempty(varargin) - var_s = varargin{1}; - if isstruct(var_s) - fields = fieldnames(var_s); - for i = 1:numel(fields) - eval([fields{i}, ' = ', num2str( var_s.(fields{i}) ), ';']); - fprintf("%s <-- %.2f \n", fields{i}, var_s.(fields{i})); - end - else - error('Optional variables should be passed as a struct.'); - end - end - - %%%% TX Signal - Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rcalpha); - - [Digi_sig,Symbols,Bits] = PAMsource(... - "fsym",fsym,"M",M,"order",19,"useprbs",1,... - "fs_out",fdac,... - "applyclipping",0,"clipfactor",1.5,... - "applypulseform",apply_pulsef,"pulseformer",Pform,... - "randkey",random_key,... - "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); - - %%%%% AWG - % El_sig = M8199A("kover",kover).process(Digi_sig); - El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); - % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); - % El_sig = El_sig.setPower(0,"dBm"); - - %%%%% Low-pass el. components %%%%%% - El_sig = Filter('filtdegree',4,"f_cutoff",tx_bw_nyquist.*f_nyquist,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); - - %%%%% Electrical Driver Amplifier %%%%%% - El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); - El_sig = El_sig.normalize("mode","oneone"); - - %%%%% MODULATE E/O CONVERSION %%%%%% - [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); - - %%%%%% Fiber %%%%%% - Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); - - %%%%%% ROP %%%%%% - Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); - - - - %%%%%% PD Square Law %%%%%% - Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); - - - %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% - Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bw_nyquist.*f_nyquist,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); - - % %%%%%% Low-pass Scope %%%%%% - %dactivated in scope module! - Lp_scpe = Filter('filtdegree',4,"f_cutoff",10e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); - - % Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); - - %%%%%% Scope %%%%%% - Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... - "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... - "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... - "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',0,'H_lpf',Lp_scpe).process(Rx_sig); - - % Scpe_sig.spectrum("displayname",'Digital (256 GSa/s) Rx Spectrum','fignum',100,'normalizeTo0dB',1); - - %%%%%% Sample to 2x fsym %%%%%% - Scpe_sig = Scpe_sig.resample("fs_in",fadc,"fs_out",2*fsym); - - %%%%%% Sync Rx signal with reference %%%%%% - [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); - - Scpe_sig.spectrum("displayname",'Prior DSP (2x fsym) Spectrum','fignum',100,'normalizeTo0dB',1); - - %%% EQUALIZING - ber = struct(); - switch eq_mode - - case equalizer_structure.ffe - - %FFE - if db_precode - Bits_ = PAMmapper(M,0).demap(Symbols); - else - Bits_ = Bits; - end - - eq_ffe = EQ("Ne",ffe_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_ffe = 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_ffe = 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); - - [eq_sig,eq_noise,ber.ber_ffe,totalErrors] = vnle( eq_ffe,M,Scpe_sig ,Symbols, Bits_); - eq_noise.spectrum("displayname",'Noise Spectrum after FFE','fignum',41,'normalizeTo0dB',0); - - case equalizer_structure.vnle - - if db_precode - Bits_ = PAMmapper(M,0).demap(Symbols); - else - Bits_ = Bits; - end - - - %VNLE - 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_sig,eq_noise,ber.ber_vnle,totalErrors] = vnle(eq_vnle,M,Scpe_sig ,Symbols, Bits_); - eq_noise.spectrum("displayname",'Noise Spectrum after VNLE','fignum',41,'normalizeTo0dB',0); - case equalizer_structure.vnle_pf_mlse - - if db_precode - Bits_ = PAMmapper(M,0).demap(Symbols); - else - Bits_ = Bits; - end - - %VNLE + PF + MLSE - eq_mlse = 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_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); - 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); - - pf_ = Postfilter("ncoeff",1); - mlse_ = MLSE("DIR",[0,0],"duobinary_output",0,"M",[],"trellis_states",[]); - [eq_sig,eq_noise,ber.ber_mlse,totalErrors] = vnle_postfilter_mlse(eq_mlse , pf_, mlse_,M, Scpe_sig ,Symbols, Bits_); - pf_.showFilter(eq_noise); - eq_noise.spectrum("displayname",'Noise Spectrum after VNLE+PF','fignum',41,'normalizeTo0dB',0); - case equalizer_structure.db_precoded - - %EQ targets DB => less precompensation; pre-coded - mlse_db_pre = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels); - eq_db_pre = 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_sig,eq_noise,ber.ber_db,totalErrors] = duobinary_target(eq_db_pre, mlse_db_pre,M, Scpe_sig ,Symbols, Bits); - eq_noise.spectrum("displayname",'Noise Spectrum after DB','fignum',41,'normalizeTo0dB',0); - %->append BER to DB - - case equalizer_structure.db_encoded - - %db signaling => db encoded - mlse_db_enc = MLSE("DIR",[1,1],"duobinary_output",1,"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); - [eq_sig,eq_noise,ber.ber_db_enc,totalErrors] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Bits); - %->append BER to DB - end - - % fprintf('BER FFE: %.2e \n',ber.ber_mlse); - % % El_sig.spectrum("displayname",'Tx Spectrum','fignum',10,'normalizeTo0dB',1); - % Scpe_sig.spectrum("displayname",'Rx Spectrum','fignum',100,'normalizeTo0dB',1); - - - - - if ~isempty(curFolder) - cd(curFolder); - end - -end \ No newline at end of file diff --git a/projects/IMDD_base_system/imdd_it.m b/projects/IMDD_base_system/imdd_it.m index 86e8430..b69a8c0 100644 --- a/projects/IMDD_base_system/imdd_it.m +++ b/projects/IMDD_base_system/imdd_it.m @@ -1,26 +1,367 @@ % basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; % db = DBHandler("pathToDB",[basePath,'silas_labor.db']); - -uloops = struct; -uloops.bitrate = [300,330,360,390,420,450,480].*1e9; -uloops.M = 4; - -wh = DataStorage(uloops); -wh.addStorage("ber"); - -wh = submit_simulations(wh,db,"parallel",0,"simulation_mode",0); - -a = wh.getStoValue('ber', uloops.bitrate,uloops.M); -for i = 1:numel(a) - bers(i) = a{i}.ber_mlse; +if 1 + uloops = struct; + uloops.precomp = [0,1]; + uloops.db_precode = [0,1]; + uloops.bitrate = [300,330,360,390,420,450,480].*1e9; %[300,330,360,390,420,450,480] + % uloops.bitrate = 390e9; + % uloops.laser_wavelength = [1293,1297.5,1302,1306.5,1310,1313.4,1318,1322.7,1327.4]; + uloops.laser_wavelength = [1293]; + uloops.M = [4,6,8]; + uloops.link_length = [2]; % 1,2,3,5,6,8,10 + wh = DataStorage(uloops); + wh.addStorage("ber"); + + wh = submit_simulations(wh,"parallel",1,"simulation_mode",0); end -figure(2024) +wh_ana = wh; + +ber_mlse = {}; +ber_vnle = {}; +inf_rate_vnle ={}; + +ber_dbtgt ={}; + +ber_dbenc ={}; +alpha = {}; +ber_dfe = {}; +ngmi = []; + +wavelength=1293; +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 + + +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 -plot(uloops.bitrate.*1e-9,bers); -yline(3.8e-3,'LineWidth',2,'DisplayName','3.8e-3'); -yline(2e-2,'LineWidth',2,'LineStyle','--','DisplayName','2e-2'); +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]); + +figure() +title(sprintf('%d km | 1310 nm | PAM %d | VNLE',uloops.link_length,uloops.M)); +hold on +line([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],[min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)],'Color',[.7,.7,.7],'Marker','none','Handlevisibility','off'); +plot(uloops.bitrate.*1e-9,cellfun(@min, inf_rate_vnle).*uloops.bitrate./log2(uloops.M).*1e-9,'DisplayName',sprintf('AIR'),'Color',cols(1,:),'LineStyle',':'); +plot(uloops.bitrate.*1e-9,netRatesVNLE.SDHD.NetRate.*1e-9,'DisplayName',sprintf('SD+HD'),'Color',cols(2,:),'LineStyle',':'); +plot(uloops.bitrate.*1e-9,netRatesVNLE.HD.NetRate.*1e-9,'DisplayName',sprintf('HD'),'Color',cols(3,:),'LineStyle',':'); +plot(uloops.bitrate.*1e-9,netRatesVNLE.KP4_hamming.NetRate.*1e-9,'DisplayName',sprintf('KP4+Hamming'),'Color',cols(4,:),'LineStyle',':'); +beautifyBERplot() +xlim([min(uloops.bitrate.*1e-9) max(uloops.bitrate.*1e-9)]) +xlabel('Gross Bitrate in Gbps'); +ylabel('Net Bitrate in Gbps'); legend + +% plot(uloops.bitrate.*1e-9,ber_vnle,'DisplayName',sprintf('NGMI MLSE; %d km',len),'Color',cols(1,:),'LineStyle','-'); +% 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'); + + + + + +% 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)); +% ber_vnle_dfe(i) = mean((a{i}.ber_vnle_dfe)); +% 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'); +% xlabel('VNLE 2nd order taps'); +% yline([3.8e-3;1e-2],'LineWidth',1,'HandleVisibility','off'); +% beautifyBERplot() + +% r = 1; +% for rate = uloops.bitrate +% +% i = 1; +% for alpha = uloops.alpha +% cols = linspecer(7);%cbrewer2('Set2',10); +% cnt = 1; +% a = wh.getStoValue('ber',alpha,uloops.vnle_order2,uloops.vnle_order3,precomp, precode, rate,uloops.M); +% ber_mlse(i) = mean((a{1}.ber_mlse)); +% 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) +% all2nd = wh.getStoValue('ber',uloops.vnle_order3(j), uloops.vnle_order2, uloops.bitrate,uloops.M); +% for i = 1:numel(a) +% 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); +% for i = 1:numI +% 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_labdata_example.m b/projects/IMDD_base_system/imdd_labdata_example.m deleted file mode 100644 index a36db10..0000000 --- a/projects/IMDD_base_system/imdd_labdata_example.m +++ /dev/null @@ -1,172 +0,0 @@ -% function [ber] = imdd_labdata_example(varargin) - -% TX - M = 4; - fsym = 224e9; - f_nyquist = fsym/2; - apply_pulsef = 0; - fdac = 2*fsym;%256e9; - fadc = 2*fsym;%256e9; - fdac = 256e9; - fadc = 256e9; - random_key = 1; - db_precode = 1; - 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 = 1.5; - -% Channel - link_length = 0; - -% RX - rop = -7.5; - rx_bw_nyquist = 1.5; - -% EQ - eq_mode = equalizer_structure.vnle_pf_mlse; - ffe_order=[50,0,0]; - vnle_order=[50,7,7]; - dfe_order = [0 0 0]; - - len_tr = 4096*2; - mu_ffe = [0.0004 0.0004 0.0004]; - mu_dfe = 0.0004; - mu_dc = 0.05; - - dfe_ = sum(dfe_order)>0; - - % - % % Parse optional input arguments - % if ~isempty(varargin) - % var_s = varargin{1}; - % if isstruct(var_s) - % fields = fieldnames(var_s); - % for i = 1:numel(fields) - % eval([fields{i}, ' = ', num2str( var_s.(fields{i}) ), ';']); - % fprintf("%s <-- %.2f \n", fields{i}, var_s.(fields{i})); - % end - % else - % error('Optional variables should be passed as a struct.'); - % end - % end - - basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; - useGui = 0; - db = DBHandler("pathToDB",[basePath,'silas_labor.db']); - filterParams = db.tables; - % filterParams.Runs.run_id = 2958; % no db - % filterParams.Runs.run_id = 2937; % no db - filterParams.Configurations = struct( ... - 'bitrate', 300e9, ... - 'db_mode', 0, ... - 'fiber_length', 1, ... - 'interference_attenuation', [], ... - 'interference_path_length', [], ... - 'is_mpi', 0, ... - 'pam_level', 4, ... - 'precomp_amp', [], ... - 'rop_attenuation', 0, ... - 'symbolrate', [], ... - 'v_awg', [], ... - 'v_bias', [], ... - 'wavelength', 1310 ... - ); - 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.fiber_length','Configurations.wavelength','Configurations.precomp_amp','BERs.ber'}; - - [dataTable,sql_query] = db.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)))); - - tx_bits = load([basePath, char(dataTable.tx_bits_path(1))]); - tx_bits = tx_bits.Bits; - tx_symbols = load([basePath, char(dataTable.tx_symbols_path(1))]); - tx_symbols = tx_symbols.Symbols; - - rx_sync = load([basePath, char(dataTable.rx_sync_path(1))]); - rx_sync = rx_sync.S; - - fsym = tx_symbols.fs; - - %%%%%% Sample to 2x fsym %%%%%% - Scpe_sig = rx_sync{1}.resample("fs_in",rx_sync{1}.fs,"fs_out",2*fsym); - - %%%%%% Sync Rx signal with reference %%%%%% - [Scpe_sig,S] = Scpe_sig.tsynch("reference",tx_symbols,"fs_ref",fsym); - - Scpe_sig.spectrum("displayname",'Rx (Scpe+Sync+Resample)','fignum',100,'normalizeTo0dB',0); - - Scpe_sig = Filter('filtdegree',4,"f_cutoff",tx_symbols.fs.*0.6,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); - - Scpe_sig.spectrum("displayname",'Rx (Scpe+Sync+Resample+LPF)','fignum',100,'normalizeTo0dB',0); - - %%% EQUALIZING - ber = struct(); - switch eq_mode - - case equalizer_structure.ffe - - %FFE - if db_precode - Bits_ = PAMmapper(M,0).demap(Symbols); - else - Bits_ = Bits; - end - - eq_ffe = EQ("Ne",ffe_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_ffe = 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_ffe = 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); - - [eq_sig,eq_noise,ber.ber_ffe,totalErrors] = vnle( eq_ffe,M,Scpe_sig ,Symbols, Bits_); - eq_noise.spectrum("displayname",'Noise Spectrum after FFE','fignum',41,'normalizeTo0dB',0); - - case equalizer_structure.vnle - %VNLE - 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_sig,eq_noise,ber.ber_vnle,totalErrors] = vnle(eq_vnle,M,Scpe_sig ,Symbols, Bits_); - eq_noise.spectrum("displayname",'Noise Spectrum after VNLE','fignum',41,'normalizeTo0dB',0); - case equalizer_structure.vnle_pf_mlse - - %VNLE + PF + MLSE - eq_mlse = 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",1); - mlse_ = MLSE("DIR",[0,0],"duobinary_output",0,"M",[],"trellis_states",[]); - [eq_sig,eq_noise,ber.ber_mlse,totalErrors] = vnle_postfilter_mlse(eq_mlse , pf_, mlse_,M, Scpe_sig ,Symbols, Bits_); - pf_.showFilter(eq_noise); - - eq_noise.spectrum("displayname",'Noise Spectrum after VNLE+PF','fignum',41,'normalizeTo0dB',0); - - case equalizer_structure.db_precoded - - %EQ targets DB => less precompensation; pre-coded - mlse_db_pre = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels); - eq_db_pre = 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_sig,eq_noise,ber.ber_db,totalErrors] = duobinary_target(eq_db_pre, mlse_db_pre,M, Scpe_sig ,Symbols, Bits); - eq_noise.spectrum("displayname",'Noise Spectrum after DB','fignum',41,'normalizeTo0dB',0); - %->append BER to DB - - case equalizer_structure.db_encoded - - %db signaling => db encoded - mlse_db_enc = MLSE("DIR",[1,1],"duobinary_output",1,"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); - [eq_sig,eq_noise,ber.ber_db_enc,totalErrors] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Bits); - %->append BER to DB - end - - fprintf('BER FFE: %.2e \n',ber.ber_mlse); - % % El_sig.spectrum("displayname",'Tx Spectrum','fignum',10,'normalizeTo0dB',1); - % Scpe_sig.spectrum("displayname",'Rx Spectrum','fignum',100,'normalizeTo0dB',1); - - - - - -% end \ 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 183f940..49e6723 100644 --- a/projects/IMDD_base_system/imdd_model.m +++ b/projects/IMDD_base_system/imdd_model.m @@ -1,7 +1,4 @@ -function [output] = imdd_model(simulation_mode,database,varargin) - - - +function [output] = imdd_model(simulation_mode,varargin) %%% Change folder curFolder = pwd; @@ -14,48 +11,58 @@ end % TX M = 4; fsym = 180e9; -f_nyquist = fsym/2; + apply_pulsef = 1; fdac = 256e9; fadc = 256e9; random_key = 1; -db_precode = 1; -emulate_precode = 0; -discard_precode = 1; +precomp = 0; +db_precode = 0; + db_encode = 0; -% duob_mode = db_mode.db_emulate; -emulate_db = 1; rcalpha = 0.05; kover = 16; vbias_rel = 0.5; u_pi = 2.9; vbias = -vbias_rel*u_pi; -laser_wavelength = 1310; +laser_wavelength = 1293; laser_linewidth = 0; -tx_bw_nyquist = 0.9; +tx_bw_nyquist = 0.8; % Channel link_length = 1; % RX rop = -8; -rx_bw_nyquist = 0.7; +rx_bw_nyquist = 0.8; -% EQ -eq_mode = equalizer_structure.vnle_pf_mlse; -ffe_order=[50,0,0]; -vnle_order=[50,7,7]; -dfe_order = [2 0 0]; +vnle_order1 = 50; +vnle_order2 = 7; +vnle_order3 = 7; -len_tr = 4096*2; -mu_ffe = [0.0004 0.0004 0.0004]; +vnle_order=[vnle_order1,vnle_order2,vnle_order3]; +dfe_order = [0 0 0]; + + +alpha = 0; + +len_tr = 4096*2; + +mu_ffe1 = 0.0001; +mu_ffe2 = 0.0008; +mu_ffe3 = 0.001; +mu_dc = 0.005; + +mu_ffe = [mu_ffe1 mu_ffe3 mu_ffe3]; mu_dfe = 0.0004; -mu_dc = 0.00; + dfe_ = sum(dfe_order)>0; +doub_mode = db_mode.db_precoded; + %%% change specific parameter if given in varargin % Parse optional input arguments if ~isempty(varargin) @@ -76,30 +83,118 @@ if ~isempty(varargin) 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; + %%% run the simulation or measurement or ... - - if simulation_mode - fsym_ = floor( bitrate*1e-9./log2(M) ).*1e9; + Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"rrcalpha",rcalpha); - 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 + [Digi_sig,Symbols,Tx_bits] = PAMsource(... + "fsym",fsym,"M",M,"order",17,"useprbs",1,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key,... + "db_precode",db_precode,"db_encode",db_encode,... + "mrds_code",0,"mrds_blocklength",512).process(); - tx_simulation; + % Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); - + %%%%% AWG + % El_sig = M8199A("kover",kover).process(Digi_sig); + El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); + % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); + % El_sig = El_sig.setPower(0,"dBm"); + + %%%%% Low-pass el. components %%%%%% + tx_bwl = tx_bw_nyquist.*f_nyquist; + % tx_bwl = 80e9; + El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); + % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1); + + %%%%% Electrical Driver Amplifier %%%%%% + El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); + El_sig = El_sig.normalize("mode","oneone"); + + %%%%% MODULATE E/O CONVERSION %%%%%% + [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); - rx_simulation; + %%%%%% ROP %%%%%% + Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); + + %%%%%% PD Square Law %%%%%% + Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); + + %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% + rx_bwl = rx_bw_nyquist.*f_nyquist; + % rx_bwl = 80e9; + Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); + + % %%%%%% Low-pass Scope %%%%%% + Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); + + % Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); + + %%%%%% Scope %%%%%% + Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc,... + "delay",0,"fixed_delay",0,"filtertype",filtertypes.butterworth,... + "samplingdelay",0,"rand_samplingdelay",0,"freq_offset",0,"samp_jitter",0,... + "adcresolution",8,"quantbuffer",0.1,'block_dc',1,'lpf_active',1,'H_lpf',Lp_scpe).process(Rx_sig); Scpe_cell{1} = Scpe_sig; else - + profile on + basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; + database = DBHandler("pathToDB",[basePath,'silas_labor.db']); + profile off basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; useGui = 0; % db = DBHandler("pathToDB",[basePath,'silas_labor.db']); @@ -119,36 +214,53 @@ else 'symbolrate', [], ... 'v_awg', [], ... 'v_bias', [], ... - 'wavelength', 1310 ... + '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','BERs.ber'}; + 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','Measurements.power_rop','Configurations.v_bias'}; [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)))); - Tx_bits = load([basePath, char(dataTable.tx_bits_path(2))]); + Tx_bits = load([basePath, char(dataTable.tx_bits_path(end))]); Tx_bits = Tx_bits.Bits; - Symbols = load([basePath, char(dataTable.tx_symbols_path(2))]); + Symbols = load([basePath, char(dataTable.tx_symbols_path(end))]); Symbols = Symbols.Symbols; - Scpe_load = load([basePath, char(dataTable.rx_sync_path(2))]); + Scpe_load = load([basePath, char(dataTable.rx_sync_path(end))]); Scpe_cell = Scpe_load.S; + % Raw_signal = load([basePath, char(dataTable.rx_raw_path(1))]); % 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); + % + % Raw_signal.spectrum("normalizeTo0dB",0,"fignum",336,"fft_length",2^12); + % Raw_signal.move_it_spectrum("fignum",334); fsym = Symbols.fs; end -output = struct(); +if db_precode + Symbols_precoded = Symbols; +end -for occ = 1:1%length(Scpe_cell) +output = struct(); +vnle_pf_package = {}; +vnle_dfe_package = {}; +dbtgt_package = {}; + + +proc_occ = min(8,length(Scpe_cell)); +for occ = 1:proc_occ Scpe_sig = Scpe_cell{occ}; @@ -156,195 +268,76 @@ for occ = 1:1%length(Scpe_cell) Scpe_sig = Scpe_sig.resample("fs_out",2*fsym); %%%%%% Sync Rx signal with reference %%%%%% - [Scpe_sig,S] = Scpe_sig.tsynch("reference",Symbols,"fs_ref",fsym); - - % Scpe_sig.spectrum("displayname",'Rx (Scpe+Sync+Resample)','fignum',100,'normalizeTo0dB',0); - - Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.6,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); - - Scpe_sig.spectrum("displayname",'Rx (Scpe+Sync+Resample+LPF)','fignum',110,'normalizeTo0dB',0); - - Scpe_sig.plot("displayname",'Filtered Scope Signal','fignum',111,'clear',1); - - Scpe_sig.eye(fsym,M); + % [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); %%% EQUALIZING - - switch eq_mode - - case equalizer_structure.ffe - - %FFE - if db_precode - Bits_ = PAMmapper(M,0).demap(Symbols); - else - Bits_ = Bits; - end - - eq_ffe = EQ("Ne",ffe_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_ffe = 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_ffe = 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); - - [eq_sig,Eq_noise,ber_ffe(occ),totalErrors] = vnle( eq_ffe,M,Scpe_sig ,Symbols, Bits_); - Eq_noise.spectrum("displayname",'Noise Spectrum after FFE','fignum',41,'normalizeTo0dB',0); - - fprintf('BER FFE: %.2e \n',ber_ffe(occ)); - - case equalizer_structure.vnle - - - - - - %VNLE - 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",0); - - [Eq_sig,Eq_noise] = eq_vnle.process(Scpe_sig,Symbols); - - Eq_sig = PAMmapper(M,0).quantize(Eq_sig); - - if emulate_precode && db_precode == 0 - % emulation - Eq_sig = Duobinary().encode(Eq_sig); - Eq_sig = Duobinary().decode(Eq_sig); - - Symbols= Duobinary().encode(Symbols); - Symbols = Duobinary().decode(Symbols); - Tx_bits = PAMmapper(M,0).demap(Symbols); - elseif db_precode == 1 && db_encode == 0 && discard_precode == 1 - % normal dsp for precoded sequence - Tx_bits = PAMmapper(M,0).demap(Symbols); - - elseif db_precode == 1 && db_encode == 1 - error('not implemented') - - elseif db_precode == 1 && db_encode == 0 - Eq_sig = Duobinary().encode(Eq_sig); - Eq_sig = Duobinary().decode(Eq_sig); - end - - % M = numel(unique(tx_symbols.signal)); - Rx_bits = PAMmapper(M,0).demap(Eq_sig); - - [~,numErrors,ber,~] = calc_ber(Rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); - - Eq_noise.spectrum("displayname",'Noise Spectrum after VNLE','fignum',41,'normalizeTo0dB',0); - - output.ber_vnle(occ) = ber; - - if 1c - figure(51); - clf - title(sprintf('DB coded PAM after EQ ; BER: %1.2e',M, ber )); - constellation = unique(Symbols.signal); - received = NaN(numel(constellation),length(Symbols)); - for lvl = 1:numel(constellation) - %Separate the equalized signal into the - %respective levels based on the actually - %transmitted level! - received(lvl,Symbols.signal==constellation(lvl)) = Eq_sig.signal(Symbols.signal==constellation(lvl)); - intermediate = received(lvl,:); - cnt(lvl) = numel(intermediate(~isnan(intermediate))); - hold on - histogram(received(lvl,:),1000,"EdgeAlpha",0,'DisplayName',['Lvl ',num2str(lvl),' | ',num2str(cnt(lvl)),' entries']); - end - legend - end - - Eq_sig.eye(fsym,M,"displayname",'Eye after EQ','fignum',52); - Eq_sig.spectrum("displayname",'Spectrum after EQ','fignum',53,'normalizeTo0dB',1); - Scpe_sig.spectrum("displayname",'Spectrum before EQ','fignum',53,'normalizeTo0dB',1); - fprintf('BER FFE: %.2e \n',ber); - - - - - - case equalizer_structure.vnle_pf_mlse - - - %VNLE + PF + MLSE - eq_mlse = 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_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]; - pf_ = Postfilter("ncoeff",2); - mlse_ = MLSE("DIR",[0,0],"duobinary_output",0,"M",[],"trellis_states",[]); + % %%%%% 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",1,"ideal_dfe",0); - %FFE or VNLE - [Eq_sig,Eq_noise] = eq_mlse.process(Scpe_sig,Symbols); + [result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",0); + vnle_dfe_package{occ} = result; + - Eq_sig = pf_.process(Eq_sig,Eq_noise); + end + %%%%% VNLE + PF + MLSE %%%% + if 1 - %M = numel(unique(tx_symbols.signal)); - mlse_.DIR = pf_.burg_coeff; - mlse_.trellis_states = PAMmapper(M,0).levels; - mlse_.M = M; - - Eq_sig = mlse_.process(Eq_sig); - pf_.showFilter(Eq_noise); - % Eq_noise.spectrum("displayname",'Noise Spectrum after VNLE+PF','fignum',41,'normalizeTo0dB',0); - - if emulate_precode && db_precode == 0 - % emulation - Eq_sig = Duobinary().encode(Eq_sig); - Eq_sig = Duobinary().decode(Eq_sig); - - Symbols= Duobinary().encode(Symbols); - Symbols = Duobinary().decode(Symbols); - Tx_bits = PAMmapper(M,0).demap(Symbols); - elseif db_precode == 1 && db_encode == 0 && discard_precode == 1 - % normal dsp for precoded sequence - Tx_bits = PAMmapper(M,0).demap(Symbols); - - elseif db_precode == 1 && db_encode == 1 - error('not implemented') - - elseif db_precode == 1 && db_encode == 0 - Eq_sig = Duobinary().encode(Eq_sig); - Eq_sig = Duobinary().decode(Eq_sig); + % 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); + 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); + vnle_pf_package{occ} = result; + end - Rx_bits = PAMmapper(M,0).demap(Eq_sig); - [~,numErrors,ber,~] = calc_ber(Rx_bits.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); + %%%%% Duobinary Targeting %%%% + if 1 - output.ber_mlse(occ) = ber; - output.pf_taps(occ,:) = pf_.burg_coeff; + 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); + dbtgt_package{occ} = result; + + + 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{occ} = result; + end - fprintf('BER MLSE: %.2e \n',ber); - - case equalizer_structure.db_precoded - - %EQ targets DB => less precompensation; pre-coded - mlse_db_pre = MLSE("DIR",[1,1],"duobinary_output",1,"M",M,"trellis_states",PAMmapper(M,0).levels); - eq_db_pre = 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_sig,Eq_noise,ber_db(occ),totalErrors] = duobinary_target(eq_db_pre, mlse_db_pre,M, Scpe_sig ,Symbols, Bits); - Eq_noise.spectrum("displayname",'Noise Spectrum after DB','fignum',41,'normalizeTo0dB',0); - %->append BER to DB - - fprintf('BER VNLE+DB: %.2e \n',ber_db(occ)); - - case equalizer_structure.db_encoded - - %db signaling => db encoded - mlse_db_enc = MLSE("DIR",[1,1],"duobinary_output",1,"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); - [eq_sig,Eq_noise,ber_db_enc,totalErrors] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Bits); - %->append BER to DB - - fprintf('BER DB: %.2e \n',ber_db_enc(occ)); - end - autoArrangeFigures; - + % autoArrangeFigures; + disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') + fprintf('\n') + end - +output.vnle_dfe_package = vnle_dfe_package; +output.vnle_pf_package = vnle_pf_package; +output.dbtgt_package = dbtgt_package; if ~isempty(curFolder) cd(curFolder); diff --git a/projects/IMDD_base_system/submit_simulations.m b/projects/IMDD_base_system/submit_simulations.m index bfc68f1..d8e1d24 100644 --- a/projects/IMDD_base_system/submit_simulations.m +++ b/projects/IMDD_base_system/submit_simulations.m @@ -1,8 +1,7 @@ -function wh = submit_simulations(wh,db,options) +function wh = submit_simulations(wh,options) arguments wh - db options.parallel = 1; options.simulation_mode = 1; end @@ -10,14 +9,24 @@ end %%% 2) SUBMIT SIMULATION % Initialize job results if options.parallel - if isempty(gcp('nocreate')) + 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',wh.getLastLinIndice); lin_idx = 1; for lin_idx = 1:wh.getLastLinIndice @@ -34,11 +43,11 @@ for lin_idx = 1:wh.getLastLinIndice if options.parallel numOutputs = 1; - results(lin_idx) = parfeval(@imdd_model, numOutputs, options.simulation_mode, db, optionalVars); + results(lin_idx) = parfeval(@imdd_model, numOutputs, options.simulation_mode, optionalVars); else - finalresults{lin_idx} = imdd_model(options.simulation_mode,db,optionalVars); + finalresults{lin_idx} = imdd_model(options.simulation_mode,optionalVars); wh.addValueToStorageByLinIdx(finalresults{lin_idx}, 'ber', lin_idx); end diff --git a/projects/IMDD_base_system/tx_simulation.m b/projects/IMDD_base_system/tx_simulation.m index 926c60c..77a1679 100644 --- a/projects/IMDD_base_system/tx_simulation.m +++ b/projects/IMDD_base_system/tx_simulation.m @@ -11,7 +11,7 @@ "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); diff --git a/test/duobinary_emulation_minimal_example.m b/test/duobinary_emulation_minimal_example.m index c390d3f..3f808b3 100644 --- a/test/duobinary_emulation_minimal_example.m +++ b/test/duobinary_emulation_minimal_example.m @@ -1,5 +1,5 @@ useprbs = 1; -M = 4; +M = 6; randkey = 1; datarate = 224e9; fsym = round(datarate / log2(M)) ; @@ -63,6 +63,7 @@ Tx_bits = Bits; disp(['BER: ',sprintf('%.1E',ber),' - - PAM-',num2str(M)]); + figure(200) clf hold on diff --git a/test/duobinary_minimal_example.m b/test/duobinary_minimal_example.m index 79e9f86..f817713 100644 --- a/test/duobinary_minimal_example.m +++ b/test/duobinary_minimal_example.m @@ -1,5 +1,5 @@ useprbs = 1; -M = 4; +M = 6; randkey = 1; datarate = 448e9; fsym = round(datarate / log2(M)) ; @@ -38,19 +38,33 @@ loop = 1; [data_out,state_out] = prms(data_in, state_, para); bitpattern = data_out'; +if M == 6 + bitpattern = reshape(bitpattern',[],1); + bitpattern = bitpattern(1:end-mod(length(bitpattern),5)); +end + Tx_bits = Informationsignal(bitpattern); +%%%%% Duobinary %%%%%% + +close all Symbols_tx = PAMmapper(M,0).map(Tx_bits); Symbols_tx.fs = fsym; Symbols1 = Duobinary().precode(Symbols_tx); +figure;histogram(Symbols1.signal); Symbols2 = Duobinary().encode(Symbols1); +figure;histogram(Symbols2.signal); Symbols3 = Duobinary().decode(Symbols2); +figure;histogram(Symbols3.signal); +autoArrangeFigures; Rx_bits = PAMmapper(M,0).demap(Symbols3); +%%%%% Check BER of Bit Sequence %%%%%% + [~,error_num,ber,error_pos] = calc_ber(Tx_bits.signal,Rx_bits.signal,"skip_front",0,"skip_end",0,"returnErrorLocation",1); disp(['BER: ',sprintf('%.1E',ber),' - - PAM-',num2str(M)]); diff --git a/test/mapping_minimal_example.m b/test/mapping_minimal_example.m index fdc86ac..08e1b3a 100644 --- a/test/mapping_minimal_example.m +++ b/test/mapping_minimal_example.m @@ -7,6 +7,55 @@ randkey = 1; fsym = 112e9; viewresults = 0; +%%%%% Mapping %%%%% +M = 8; +data = 0:M-1; +bitpersymbol = log2(M); + +%Bits +bits = int2bit(data,bitpersymbol, true); + +%Integers +ints = bit2int(bits,bitpersymbol, true); + +Tx_bits = Informationsignal(bits'); +Digi_Mod = PAMmapper(M,0); +Symbols = Digi_Mod.map(Tx_bits); +moveitgray = Symbols.signal; + +%Gray Symbol Mapping +matlabgray = pammod(ints,M,0,'gray'); +scaling_factor = rms(unique(matlabgray)); +matlabgray = matlabgray ./ scaling_factor; + + +%Demod +moveitgray = moveitgray.* scaling_factor; +matlabgray = matlabgray .* scaling_factor; +ints_demap = pamdemod(matlabgray,M,0,'gray'); + +bits_demap = int2bit(ints_demap,bitpersymbol, true); + +% for i = 1:M +% fprintf('%d , %d , %d --> %d \n',bits(1,i),bits(2,i),bits(3,i),matlabgray(i)); +% end +% +% for i = 1:M +% fprintf('%d , %d , %d --> %d \n',bits(1,i),bits(2,i),bits(3,i),moveitgray(i)); +% end + +scatterplot(matlabgray,1,0,'b*'); +for k = 1:M + text(real(matlabgray(k)),imag(matlabgray(k))+0.6,num2str(ints_demap(k)),"Color",[1 1 1]); + + text(real(matlabgray(k)),imag(matlabgray(k))-1.6,num2str(bits(:,k)),"Color",'blue'); + text(real(moveitgray(k)),imag(moveitgray(k))-3,num2str(bits(:,k)),"Color",'green'); +end +axis([-M M -3 2]) + + +symbols = bit2int(bitGroups',bitpersymbol, true); + %%%%% PRBS Generation in correct shape for Modulation Format %%%%%% O = 10; %order of prbs @@ -34,6 +83,9 @@ Tx_bits = Informationsignal(bitpattern); Digi_Mod = PAMmapper(M,0); +symbols = bitMapper(bitpattern, M, 'PAM'); + + Symbols = Digi_Mod.map(Tx_bits); Rx_bits = PAMmapper(M,0).demap(Symbols);