clc;
clear;
close all;

% Defining variables
filename = "enel420_grp_1.txt";
fs = 1024; %Hz
diff1 = 0;
diff2 = 0;
random_mutation = 0;
%i = 0;
z = 0;
SNR_Min = 0;

%load the ECG data
data = load(filename);

while (z < 30)
    z = z + 1;
    i = 0;

    b = [1, -1.59158, 1];
    X = [1, -1.5, 1.0];
    Y = [1, -1.6, 0.9];
    Z = [1, -1.5, 0.95];

    % Applying each filter
    IIR_filter_signal_X = filter(b, X, data);
    IIR_filter_signal_Y = filter(b, Y, data);
    IIR_filter_signal_Z = filter(b, Z, data);

    % Determining the SNR for each filter
    snr_db_X = calculate_snr_after_filter(data, IIR_filter_signal_X);
    snr_db_Y = calculate_snr_after_filter(data, IIR_filter_signal_Y);
    snr_db_Z = calculate_snr_after_filter(data, IIR_filter_signal_Z);

    SNR_Max = max([snr_db_X, snr_db_Y, snr_db_Z]);

    while (SNR_Max < 15) && (i < 30)

    % Incremeting the counter
    i = i + 1;

    % Removing weakest filter
    SNR_Min = min([snr_db_X, snr_db_Y, snr_db_Z]);

    % Shift around so that Z always has the worse values
    if (snr_db_X == SNR_Min)
        X = Z;
    end

    if (snr_db_Y == SNR_Min)
        Y = Z;
    end

    % Determining the new variables
    diff1 = (X(2) - Y(2));
    diff2 = (X(3) - Y(3));
    if diff1 > 0.002
        diff1 = 0.002;
    end
    if diff2 > 0.002
        diff2 = 0.002;
    end
    random_mutation = (0.2 * rand(1,1)) - 0.1;

    X = [1, X(2) + diff1, X(3) + diff2];
    Y = [1, Y(2) - diff1, Y(3) - diff2];
    Z = [1, (X(2) + random_mutation + diff1), (random_mutation + X(3) + diff2)];

    % Applying each filter
    IIR_filter_signal_X = filter(b, X, data);
    IIR_filter_signal_Y = filter(b, Y, data);
    IIR_filter_signal_Z = filter(b, Z, data);

    % Determining the SNR for each filter
    snr_db_X = calculate_snr_after_filter(data, IIR_filter_signal_X);
    snr_db_Y = calculate_snr_after_filter(data, IIR_filter_signal_Y);
    snr_db_Z = calculate_snr_after_filter(data, IIR_filter_signal_Z);

    if (isnan(snr_db_X) == 1) || (snr_db_X < 0)
        X = [1, -1.5, 1.0];
        IIR_filter_signal_X = filter(b, X, data);
        snr_db_X = calculate_snr_after_filter(data, IIR_filter_signal_X);
    end
    if (isnan(snr_db_Y) == 1) || (snr_db_Y < 0)
        Y = [1, -1.6, 0.9];
        IIR_filter_signal_X = filter(b, Y, data);
        snr_db_Y = calculate_snr_after_filter(data, IIR_filter_signal_Y);
    end
    if (isnan(snr_db_Z) == 1) || (snr_db_Z < 0)
        Z = [1, -1.5, 0.95];
        IIR_filter_signal_X = filter(b, Z, data);
        snr_db_Z = calculate_snr_after_filter(data, IIR_filter_signal_Z);
    end

    SNR_Max = max([snr_db_X, snr_db_Y, snr_db_Z]);

    end
    fprintf("%.2f, %.2f \n", SNR_Max, i);
end