clear;
clf;

function [num, den] = Pbande2ordNum (Q, Fo, Fe)
    a = 1/tan(%pi*Fo/Fe);
    num = [1, 0, -1];
    den = round([1 + Q*a + Q/a, 2*Q/a - 2*Q*a, Q*a - 1 + Q/a]);
endfunction

function [num, den] = Pbas1ordNum (Fc, Fe)
    a = 1/tan(%pi*Fc/Fe);
    num = [1 1];
    den = round([1 + a, 1 - a]);
endfunction

Q  = 7:0.25:9;
F1 = 1000;      F2 = 1200;
Fe = 9.725e3;   Te = 1/Fe;

fc = 110;

ecart_f1 = 999;
ecart_f2 = 999;
bestQ1 = 999;
bestQ2 = 999;
bestDen1 = 0;
bestDen2 =0;
bestf1 = 0;
bestf2 = 0;

for i = 1 : length(Q)
     res = 0.5;
     N = 10000

    [num1, den1] = Pbande2ordNum (Q(i), F1, Fe);
    [num2, den2] = Pbande2ordNum (Q(i), F2, Fe);

    [H1, fr1] = frmag(num1, den1, N);
    [H2, fr2] = frmag(num2, den2, N);

    [val, ind1] = max(H1);
    f1_reel = fr1(ind1)*Fe;
    [val, ind2] = max(H2);
    f2_reel = fr2(ind2)*Fe;

    err_f1 = abs(f1_reel - F1);
    err_f2 = abs(f2_reel - F2);

    if err_f1 < ecart_f1 then
        ecart_f1 = err_f1;
        bestQ1 = Q(i);
        bestDen1 = den1;
        bestf1 = f1_reel;
    end;

    if err_f2 < ecart_f2 then
        ecart_f2 = err_f2;
        bestQ2 = Q(i);
        bestDen2 = den2;
        bestf2 = f2_reel;
    end;

end

[num3, den3] = Pbas1ordNum (fc, Fe);
[H3, fr3] = frmag(num3, den3, 10000);
plot(fr3.*Fe, H3, 'b');

disp(den3, "den3: ");


disp(bestf1, "f1 = ");
disp(bestQ1, "bestQ1 = ");
disp(bestDen1, "den1 = ");
disp(bestf2, "f2 = ");
disp(bestQ2, "bestQ2 = ");
disp(bestDen2, "den2 = ");