% generating error patterns using gilbert-elliot model

p_gg = 0.99;
p_bb = 0.1;
p_bg = 1 - p_bb;

% freely chosen bit error probability
e_g = 10^-4;
e_b = 10^-1;

% generate random bitstream
bits_array_len = 10^6;
bits = randi([0, 1], 1, bits_array_len);

% start with a clean error free pipe
channel_state = true;

% preallocated array/vector for the resulting bits
result_stream = zeros(1, bits_array_len);

bit_error_rate_counter = 0;
gap_begin = 1; % stores the last dirty bit idx, meaning the begin of a new gap
gap_tracker = []; % holds a 3xn matrix of begin and end gap values and the gap size between them

state_sequence = strings(0); % stores the state timeline
state_visit_counter = 0;

for it = 1 : bits_array_len
    state_visit_counter = state_visit_counter + 1;
    current_bit = bits(it);

    % two different rands for channel and bit logic
    % to prevent unwanted coupling of these two mechanics
    channel_change_probability = rand();
    bit_error_probability = rand();

    % error and channel
    if channel_state == true
        % channel selection from good state with probability p_gg
        if channel_change_probability < p_gg
            channel_state = 1;
        else
            channel_state = 0;
            state_sequence(end + 1) = string(state_visit_counter) + 'x Good';
            state_visit_counter = 0;
        end
        
        % flip the bit with probability e_g
        if bit_error_probability < e_g
            received_bit = 1 - current_bit;
        else 
            received_bit = current_bit;
        end
    else
        % change channel to good with probability of p_bg
        if channel_change_probability < p_bg
            channel_state = 1;
            state_sequence(end + 1) = string(state_visit_counter) + 'x Bad';
            state_visit_counter = 0;
        else
            channel_state = 0;
        end

        if bit_error_probability < e_b
            received_bit = 1 - current_bit;
        else 
            received_bit = current_bit;
        end
    end

    result_stream(it) = received_bit;

    % BER and Gap
    if current_bit ~= received_bit
        bit_error_rate_counter = bit_error_rate_counter + 1;
        
        gap_tracker(:, end + 1) = [gap_begin; it; it - gap_begin - 1];
        gap_begin = it;
    end
end

gap_tracker(:, 1) = []; % trim first gap, this is a fault entry by design
bit_error_rate = bit_error_rate_counter / bits_array_len;

gap_sizes = gap_tracker(3, :);
gap_sizes = gap_sizes(gap_sizes <= 500); % limit gap size to 500 for histogram

histogram(gap_sizes, 'BinEdges', -0.5:1:500.5);
xlabel('Gap size')
ylabel('Frequency')
title('Histogram of gap sizes between consecutive errors')

disp(state_sequence)
disp(bit_error_rate)