trace1 = load("traces/2412mhz.mat", "iq");
trace2 = load("traces/2432mhz.mat", "iq");
trace3 = load("traces/2452mhz.mat", "iq");
trace4 = load("traces/2472mhz.mat", "iq");

window = 512;
overlap = 64;
fft_precision = 2048;

sample_rate = 20e6; % 20 MS/s -> 20e6 S/s

% 1.
duration1 = length(trace1.iq) / sample_rate;
duration2 = length(trace2.iq) / sample_rate;
duration3 = length(trace3.iq) / sample_rate;
duration4 = length(trace4.iq) / sample_rate;

% 2.
% get all four plot values
[S1,F1,T1] = spectrogram(trace1.iq, window, overlap, fft_precision, sample_rate, "centered");
[S2,F2,T2] = spectrogram(trace2.iq, window, overlap, fft_precision, sample_rate, "centered");
[S3,F3,T3] = spectrogram(trace3.iq, window, overlap, fft_precision, sample_rate, "centered");
[S4,F4,T4] = spectrogram(trace4.iq, window, overlap, fft_precision, sample_rate, "centered");

mag1 = abs(S1).^2;
mag2 = abs(S2).^2;
mag3 = abs(S3).^2;
mag4 = abs(S4).^2;

% begin calibrated color scale
P1_dB = 10 * log10(mag1);
P2_dB = 10 * log10(mag2);
P3_dB = 10 * log10(mag3);
P4_dB = 10 * log10(mag4);

all_P_dB = [P1_dB(:); P2_dB(:); P3_dB(:); P4_dB(:)];

p_dB_min = min(all_P_dB);
p_dB_max = max(all_P_dB);
% end power calibrated color scale

if ~calculations_only
    % plot all figures like example shows
    figure(1);
    imagesc(T1, F1, P1_dB);
    axis ij;
    xlabel("Time [s]");
    ylabel("Freq. [Hz]");
    title("Spectrogram Trace 1");
    cb = colorbar;
    ylabel(cb, "Power [dB]");
    clim([p_dB_min p_dB_max]);
    xlim([0 0.05]);
    
    figure(2);
    imagesc(T2, F2, P2_dB);
    axis ij;
    xlabel("Time [s]");
    ylabel("Freq. [Hz]");
    title("Spectrogram Trace 2");
    cb = colorbar;
    ylabel(cb, "Power [dB]");
    clim([p_dB_min p_dB_max]);
    xlim([0 0.05]);
    
    figure(3);
    imagesc(T3, F3, P3_dB);
    axis ij;
    xlabel("Time [s]");
    ylabel("Freq. [Hz]");
    title("Spectrogram Trace 3");
    cb = colorbar;
    ylabel(cb, "Power [dB]");
    clim([p_dB_min p_dB_max]);
    xlim([0 0.05]);
    
    figure(4);
    imagesc(T4, F4, P4_dB);
    axis ij; % flip so its like the example fig
    xlabel("Time [s]");
    ylabel("Freq. [Hz]");
    title("Spectrogram Trace 4");
    cb = colorbar;
    ylabel(cb, "Power [dB]");
    clim([p_dB_min p_dB_max]);
    xlim([0 0.05]); % reduce view so its like the example fig
end

% 3. 
% taking the 10th percentile on the linear values or the logarithmic ones
% should not make any difference

noise_floor1_linear = prctile(mag1(:), 10);
noise_floor2_linear = prctile(mag2(:), 10);
noise_floor3_linear = prctile(mag3(:), 10);
noise_floor4_linear = prctile(mag4(:), 10);

noise_floor1 = 10 * log10(noise_floor1_linear);
noise_floor2 = prctile(P2_dB(:), 10);
noise_floor3 = prctile(P3_dB(:), 10);
noise_floor4 = prctile(P4_dB(:), 10);

% this is false -> linear avg on logarithmic values
% noise_floor_avg = sum([noise_floor1 noise_floor2, noise_floor3, noise_floor4]) / 4;

noise_floor_avg = 10 * log10(abs(sum([prctile(mag1(:), 10) prctile(mag2(:), 10) prctile(mag3(:), 10) prctile(mag4(:), 10)]) / 4));

% 4.
occupancy1 = sum(P1_dB(:) > noise_floor1 + 10) / numel(P1_dB) * 100;
occupancy2 = sum(P2_dB(:) > noise_floor2 + 10) / numel(P2_dB) * 100;
occupancy3 = sum(P3_dB(:) > noise_floor3 + 10) / numel(P3_dB) * 100;
occupancy4 = sum(P4_dB(:) > noise_floor4 + 10) / numel(P4_dB) * 100;

occupancies = [occupancy1 occupancy2 occupancy3 occupancy4];
channels = [1 5 9 13];

[occupancy_max, max_idx] = max(occupancies);
busiest_channel = channels(max_idx);

if ~calculations_only
    disp("Results Task 1:" + newline)

    disp("-- Trace duration ---")
    disp("Channel  1 duration: " + duration1 + "s")
    disp("Channel  5 duration: " + duration2 + "s")
    disp("Channel  9 duration: " + duration3 + "s")
    disp("Channel 13 duration: " + duration4 + "s" + newline)

    disp("--- Noise floor ---")
    disp("Distinct:")
    disp("Channel  1 noise floor: " + noise_floor1 + "dB")
    disp("Channel  5 noise floor: " + noise_floor2 + "dB")
    disp("Channel  9 noise floor: " + noise_floor3 + "dB")
    disp("Channel 13 noise floor: " + noise_floor4 + "dB" + newline)

    disp("Combined:")
    disp("Channel all noise floor: " + noise_floor_avg + "dB" + newline)

    disp("-- Occupancy ---")
    disp("Channel  1 occupancy: " + occupancy1 + "%")
    disp("Channel  5 occupancy: " + occupancy2 + "%")
    disp("Channel  9 occupancy: " + occupancy3 + "%")
    disp("Channel 13 occupancy: " + occupancy4 + "%" + newline)

    disp("Busiest channel: " + busiest_channel)
end