function [next_state, output]=trellis(G)

% Trellis diagram generator function
% [next_state, output]=trellis(G)
% G: Generator polynomial matrix (N_IN by N_OUT real matrix)
%   G(i,j) is the polynomial associated with j th output bit and i th input bit (decimal unit)
% next_state : next state matrix (N_STATE by 2^N_IN)
%              state is numbered from 0 to N_STATE-1
% output : output codeword matrix (N_STATE by 2^N_IN)
%              output(i,j) means output codeword correponding to input value of (j-1) at (i-1) state
%
% This version of trellis function assumes the feedforward implementation 
% (i.e G should always be the interger matrix)

[N_IN N_OUT]=size(G);   %N_IN : Number of input bit of the convolutional encoder
                        %N_OUT : Number of output bit of the convolutional encoder

for i=1:N_IN
    n_tap(i)=size(dec2bin(G(i,:)),2)-1; % Maximum number of delay elements for i th input bit
end

N_TOTAL_TAP=sum(n_tap);    % Total number of delay elements
N_STATE=2^sum(n_tap)  ;    % Number of states

current_state=0:N_STATE-1;

%Next state and output codeword initialization
next_state=zeros(N_STATE,2^N_IN);
output=zeros(N_STATE,2^N_IN);

%State update and output codeword calculation

%State loop
for i=1:N_STATE
    
    %Input loop
    for j=1:2^(N_IN)
        
        %Output bit vector initialization (bit_out is a binary representation of output(i,j) in the vector form)
        bit_out=zeros(1,N_OUT);
        
        %Specific update for each input bit
        for k=1:N_IN
            
            %start_index and end_index tell the beginning and the end of state bits associated with delay elements for k-th input bit
            %   start_index corresponds to the unit-delayed element
            %   end_index corresponds to the most delayed (i.e n_tap(k)) element
            
            if k==1
                start_index=1;
            else
                start_index=(k-1)*n_tap(k-1)+1;
            end
            
            end_index=start_index+n_tap(k)-1;
            left_shift=N_TOTAL_TAP-sum(n_tap(1:k));
            
            %Contents in delay elements associated with k-th input bit
            small_state=bitand(bitshift(current_state(i),-(start_index-1),N_STATE),2^(end_index+1)-1); 
            
            %Output calcuation
            bit_in=bitget(j-1,k); %k-th input bit for j-th input
            state_and_input=2*small_state+bit_in; %delayed and the current bit of k-th input bit
        
            for l=1:N_OUT
               %Contribution of k-th input bit to l th output bit
                contr=bitand(G(k,l),state_and_input);
                %output bit vector is the sum of each contribution from k-th input bit
                bit_out(l)=mod(bit_out(l)+sum(dec2bin(contr)),2);
            end
            
            %State update
            next_small_state(k)=bitshift(small_state,1,n_tap(k))+bit_in;
            next_state(i,j)=next_state(i,j)+next_small_state(k)*2^(start_index-1);
        end

        output(i,j)=sum(2.^(0:N_OUT-1).*bit_out);
    end
end