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