clear all;

global N = 2048;

function rejection = rejection_criteria(log_data, fc)
    N = length(log_data);

    % Index of the first point in the tail fir
    index_tail = round((fc + 0.1) * N) + 1;

    % Index of th last point on the band
    index_band = round((fc - 0.1) * N);

    % Get the worst rejection in stopband
    worst_rejection = -max(log_data(index_tail:end));

    % Get the total of deviation in passband
    worst_band = mean(-1 * abs(log_data(1:index_band)));

    % Compute the rejection
    passband_malus = 10; % weighted value to penalize the deviation in passband
    rejection = worst_band * passband_malus + worst_rejection;
endfunction

function [h, log_curve, rejection] = compute_freqz(filename)
    global N;

    b = load(filename);
    [h, w] = freqz(b, 1, N/2);
    mag = abs(h);
    mag = mag ./ mag(1);
    log_curve = 20 * log10(mag);

    rejection = rejection_criteria(log_curve, 0.5);
endfunction

# Stages
hTotal = ones(N/2, 1);
% [h, curve1, c1] = compute_freqz("filters/fir1/fir1_033_int08"); % 1) -8dB
% [h, curve1, c1] = compute_freqz("filters/fir1/fir1_037_int08"); % 2) -9dB
[h, curve1, c1] = compute_freqz("filters/fir1/fir1_037_int08");
hTotal = hTotal .* h;
% [h, curve2, c2] = compute_freqz("filters/fir1/fir1_033_int10"); % 1) -8dB
% [h, curve2, c2] = compute_freqz("filters/fir1/fir1_033_int10"); % 2) -9dB
[h, curve2, c2] = compute_freqz("filters/fir1/fir1_033_int10");
hTotal = hTotal .* h;
% [h, curve3, c3] = compute_freqz("filters/fir1/fir1_033_int08");
% hTotal = hTotal .* h;
% [h, curve4, c4] = compute_freqz("filters/fir1/fir1_033_int10");
% hTotal = hTotal .* h;
% [h, curve5, c5] = compute_freqz("filters/fir1/fir1_015_int11");
% hTotal = hTotal .* h;

# Log total
mag = abs(hTotal);
mag = mag ./ mag(1);
log_freqz = 20 * log10(mag);
cTotal = rejection_criteria(log_freqz, 0.5);

[ c1+c2 cTotal ]
% [ c1+c2+c3+c4+c5 cTotal ]

clf;
f_axe = [1:N/2] * 2/N;
hold on;
color = [0/255 114/255 189/255];
plot(f_axe, curve1, "linewidth", 1.5, "color", color);
plot([0 1], [-c1 -c1], "--", "linewidth", 1.5, "color", color);

color = [217/255 83/255 25/255];
plot(f_axe, curve2, "linewidth", 1.5, "color", color);
plot([0 1], [-c2 -c2], "--", "linewidth", 1.5, "color", color);
% plot(f_axe, curve3, "linewidth", 1.5);
% plot(f_axe, curve4, "linewidth", 1.5);
% plot(f_axe, curve5, "linewidth", 1.5);

color = [237/255 177/255 32/255];
plot(f_axe, log_freqz, "linewidth", 1.5, "color", color);
plot([0 1], [-cTotal -cTotal], "--", "linewidth", 1.5, "color", color);
plot([0 1], [-(c1 + c2) -(c1 + c2)], ":", "linewidth", 1.5, "color", color);

plot([0.4 0.4], [-500 50], "k:")
plot([0.6 0.6], [-500 50], "k:")
ylim([-200 10])
hold off;

xlabel("Normalized Frequency (a.u.)")
ylabel("Rejection (dB)")
legend("Reponse of 1st filter", "Rejection of 1st filter", "Reponse of 2nd filter", "Rejection of 2nd filter", "Reponse Total", "Actual Rejection", "Expected Rejection", "location", "southwest")