Files
imdd_silas/Functions/EQ_visuals/analyzeEQperformance.m
Silas Oettinghaus d099efea03 Add new functions
2025-02-17 21:31:12 +01:00

150 lines
4.7 KiB
Matlab

function analyzeEQperformance(ref_bits,ref_symbols,rx_signal,eq_signal,eq_decisions,fsym,M,options)
arguments
ref_bits
ref_symbols
rx_signal
eq_signal
eq_decisions
fsym
M
options.postfilterclass
options.eqclass
options.mlseclass
options.db_precoded
options.displayname
end
% toolkit to visualize stuff related to DSP of IM/DD
%%% Preps
if isempty(eq_decisions)
eq_decisions = PAMmapper(M,0).quantize(eq_signal);
end
%%% Demap eqlzd signal to determine BER
if options.db_precoded
eq_hd = PAMmapper(M,0).quantize(eq_signal);
eq_hd = Duobinary().encode(eq_hd);
eq_hd = Duobinary().decode(eq_hd,"M",M);
rx_bits = PAMmapper(M,0).demap(eq_hd);
else
rx_bits = PAMmapper(M,0).demap(eq_signal);
end
[~,numerr,ber_sd,errpos] = calc_ber(rx_bits.signal,ref_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% fprintf('SD BER: %.2e \n',ber_sd);
%%% Demap provided decisions to determine BER
rx_bits = PAMmapper(M,0).demap(eq_decisions);
[~,numerr,ber_mlse,errpos] = calc_ber(rx_bits.signal,ref_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1);
% fprintf('MLSE BER: %.2e \n',ber_mlse);
%%% Noise prior to DSP (is thius accurate with resampling?)
rx_resampled = rx_signal.normalize("mode","rms").resample("fs_out",fsym);
rx_noise = rx_resampled-ref_symbols;
%%% Noise after to soft-decision DSP
eq_noise = eq_signal-ref_symbols;
col = cbrewer2('paired',12);
lines = numel(findall(figure(200), 'Type', 'Line'))+1;
darkcoloridx = max(mod(2*lines,12),2);
lightcoloridx = max(mod(2*lines-1,12),1);
%%% Separate Classes
constellation = unique(ref_symbols.signal);
received_sd = NaN(numel(constellation),length(ref_symbols));
received_hd = NaN(numel(constellation),length(ref_symbols));
lvlcol = cbrewer2('Set1',numel(constellation));
for lvl = 1:numel(constellation)
%Separate the equalized signal into the
%respective levels based on the actually
%transmitted level!
received_sd(lvl,ref_symbols.signal==constellation(lvl)) = eq_signal.signal(ref_symbols.signal==constellation(lvl));
received_hd(lvl,ref_symbols.signal==constellation(lvl)) = eq_decisions.signal(ref_symbols.signal==constellation(lvl));
end
%%% bursts
% Find differences between consecutive elements
diff_indices = diff(errpos);
% Identify the start of new sequences (when the difference is not 1)
sequence_starts = [1, find(diff_indices ~= 1) + 1]; % Include the first index
sequence_ends = [sequence_starts(2:end) - 1, length(errpos)]; % Calculate end indices
% Initialize burst count and print bursts longer than 10
burst_len = 1:10;
burst_count = zeros(length(burst_len),1);
for t = 1:numel(burst_len)
for i = 1:length(sequence_starts)
% Extract current sequence
current_burst = errpos(sequence_starts(i):sequence_ends(i));
% Check if the sequence length matches criterion
if length(current_burst) == burst_len(t)
burst_count(t) = burst_count(t) + 1;
end
end
end
burst_symbols = burst_count .* burst_len';
burst_rate = burst_symbols ;%./ length(rx_signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Rx Spectrum %200
rx_signal.spectrum("displayname",sprintf('Rx %d GBd PAM%d',fsym.*1e-9,M),'fignum',200,'normalizeTo0dB',1,'color',col(darkcoloridx,:));
xline([-fsym/2,fsym/2].*1e-9,'Color',col(mod(lines,12)+2,:),'HandleVisibility','off');
ylim([-20,3]);
%%% EQ Spectrum
%
%%% EQ Time Series %210
showEQTimeSignal(eq_signal,ref_symbols)
%%% EQ Noise Spectrum + inverted Postfilter %220
showEQNoisePSD(eq_noise,options.postfilterclass.burg_coeff,"fignum",220,"color",col(darkcoloridx,:));
%%% EQ SNR Spectrum %230
fft_length = 2^13;
[s_lin,w] = pwelch(eq_signal.signal,hanning(fft_length),fft_length/2,fft_length,eq_signal.fs,"centered","psd","mean");
[n_lin,w] = pwelch(eq_noise.signal,hanning(fft_length),fft_length/2,fft_length,eq_noise.fs,"centered","psd","mean");
w = w.*1e-9;
snr_dbm = 10*log10(s_lin./n_lin);
figure(230)
hold on
plot(w,snr_dbm,'DisplayName','SNR after EQ','LineWidth',1,'Color',col(darkcoloridx,:));
xlabel("Frequency in GHz");
edgetick = 2^(nextpow2(eq_signal.fs*1e-9));
xticks(-edgetick:16:edgetick);
xlim([-128 128]);
ylim([-20,35]);
grid minor
yticks(-200:10:100);
grid on; grid minor;
legend('Interpreter','none');
title('Noise of soft decision signal (not MLSE)');
%%% FFE histogram %240
showLevelHistogram(eq_signal,ref_symbols)
%%% Confusion Matrix
showLevelConfusionMatrix(eq_decisions,ref_symbols,"fignum",250)
%%% Burst Count
figure(260)
hold on
plot(burst_len,burst_rate,'Marker','x','Color',col(darkcoloridx,:),"displayname",sprintf('Rx %d GBd PAM%d; %s',fsym.*1e-9,M,options.displayname));
xlabel('length of error burst');
ylabel('occurences')
grid on
set(gca, 'YScale', 'log');
autoArrangeFigures
end