64b/66b encoder

function [ enc ] = enc64b66b( bits )
%ENC64B66B 64b/66b encoder
%   Input: stream of bits to be encoded
%     Should be a multiple of 64 bits
%   Output: stream of 64b/66b encoded bits
%     Will be a multiple of 66 bits
% 
% Alex Forencich <alex@alexforencich.com>
 
% Sources:
% http://cp.literature.agilent.com/litweb/pdf/ads2008/numeric/ads2008/64b66bCoder.html
% http://grouper.ieee.org/groups/802/3/10G_study/public/jan00/walker_1_0100.pdf
 
bits = bits(:)' > 0;
 
rc = size(bits);
 
m = mod(rc(2), 64);
 
if m > 0
    bits = [bits zeros(1,m)];
end
 
rc = size(bits);
 
bc = rc(2)/64;
 
bytes = reshape(bits, 64, bc)';
 
enc = zeros(bc, 66);
 
sreg = ones(1,58);
 
for i = 1:bc
    val = bytes(i,:);
 
    for j = 1:64
        v = xor(sreg(58),sreg(39));
        v = xor(v, val(j));
        sreg = circshift(sreg, [0 1]);
        sreg(1) = v;
        val(j) = v;
    end
 
    enc(i,:) = [0 1 val];
end
 
enc = reshape(enc', 1, bc*66);
 
end