function [y,err]=golaycodec(x,enc,ext)
% GOLAYCODEC - encode/decode a binary array using the Golay code with error correction.
%
% golaycodec encodes or decodes a binary array using the Golay code.
% Words of length 12 are encoded as codewords of length 23 (or 24 for the
% extended Golay code). The Golay code is tolerant of up to 3 errors,
% i.e. for any word of length 23 there is a unique codeword that differs
% in at most 3 positions.
%
% The Golay code can be constructed as a polynomial code using one of the
% degree 11 factors of z^23+1 (mod 2). Here we use
% g = 1 + z^2 + z^4 + z^5 + z^6 + z^10 + z^11
% For more details see Wikipedia etc.
%
% Usage:
% y = golaycodec(x)
% [y,err] = golaycodec(x)
% ... = golaycodec(x,enc)
% ... = golaycodec(x,enc,ext)
%
% Inputs:
% x can be an Mx12, Mx23 or Mx24 binary (logical) array or an Mx1 integer array.
% If x is Mx1, the binary representation of each integer is used.
% enc is a logical; if true x is to be encoded, otherwise x is to be decoded.
% If not specified, enc is inferred from the size of x or set to true.
% ext is a logical; if true use the extended Golay code (length 24),
% otherwise use the non-extended Golay code (length 23).
% If not specified, enc is inferred from the size of x or set to false.
% NOTE: the algorithm has so far only been implemented for ext=false
%
% Outputs:
% y is an array consisting of the coded or decoded rows of x
% If x is Mx12 and ext=false then y is Mx23
% If x is Mx12 and ext=true then y is Mx24
% If x is Mx23 or Mx24 then y is Mx12
% If x is Mx1 then y is Mx1
% err (only available if decoding is done, e.g enc=false) is an array
% with the errors that differ each word in x from its closest codeword)
%
% Example: encode some messages, add transmission errors and decode
% n=10;
% x=zeros(n,12);for i=1:n,x(i,ceil(12*rand(1,8)))=1;end; % random message
% y=golaycodec(x);
% err=zeros(n,23);for i=1:n,err(i,ceil(23*rand(1,3)))=1;end; % 3 random errors per message
% y1=xor(y,err); % add transmission error
% [x1,err1]=golaycodec(y1); % decode: should have x1==x and err1==err
%
% Author: Ben Petschel 18/3/2009
%
% Change history:
% 18/3/2009 - created
% 30/4/2009 - updated to implement extended Golay codes
%
encdef = true; % default value of "enc" if not specified for integer input
extdef = false; % default value of "ext" if not specified
% input checking: x
nx=size(x,2);
if any(~isfinite(x(:))) || any(~isreal(x(:)))
error('golaycodec: input array must be finite real values');
end;
if nx==1
if any(x~=round(x))
error('golaycodec: column input must be integers');
end;
if any(x<0)
% check later that integers are small enough
error('golaycodec: column input must be non-negative integers');
end;
intout = true;
elseif nx>1
if any((x(:)~=0)&(x(:)~=1))
error('golaycodec: entries of row or matrix input must be 0 or 1');
end;
if ~any(nx==[12,23,24]),
error('golaycodec: number of columns of input must be 1, 12, 23 or 24');
end;
intout = false; % binary output
else
error('golaycodec: empty input');
end; % if nx==1 elseif ... else ... end;
% input checking: enc, ext must be booleans
if nargin>=2,
if isnumeric(enc),
enc = enc>0;
end;
if ~islogical(enc) || ~isscalar(enc)
error('golaycodec: enc must be 1x1 logical value');
end;
end;
if nargin>=3,
if isnumeric(ext),
ext = ext>0;
end;
if ~islogical(ext) || ~isscalar(ext)
error('golaycodec: enc must be 1x1 logical value');
end;
end;
% input checking: consistency of nx with enc and ext if specified
switch nx
case 1
if nargin == 1,
enc = encdef; % set to default value (usually true)
end;
if nargin <= 2,
ext = extdef; % set to default value (usually false)
end;
case 12
if nargin == 1,
enc = true; % encode 12-digit word
elseif ~enc,
error('golaycodec: cannot decode 12-digit words');
end;
if nargin <=2,
ext = extdef; % set to default value (usually false)
end;
case 23
if nargin == 1,
enc = false; % decode 23-digit word
elseif enc,
error('golaycodec: cannot encode 23-digit words');
end;
if nargin <=2,
ext = false; % not from the extended code
elseif ext,
error('golaycodec: 23 digit words are not in the extended code');
end;
case 24
if nargin == 1,
enc = false; % decode 24-digit word
elseif enc,
error('golaycodec: cannot encode 24-digit words');
end;
if nargin <=2,
ext = true; % use the extended code
elseif ~ext,
error('golaycodec: 24-digit words are not in the non-extended code');
end;
otherwise
error('golaycodec: input must have 1, 12, 23 or 24 columns');
end; % switch nx
% for integer inputs, check that the integers are not too large
if intout
if enc
if any(x(:)>=2^12),
error('golaycodec: integers for coding must be less than 2^12');
end;
elseif any(x(:)>=2^23),
error('golaycodec: integers for decoding must be less than 2^23');
end;
end;
% output checking: nargout>1 only for decoding
if (nargout>2) || ((nargout>1) && enc)
error('golaycodec: can only specify two outputs when decoding');
end;
if nx==1,
% convert integers to the appropriate length
if enc
nx=12;
elseif ext
% decode extended code
nx=24;
else
% decode non-extended code
nx=23;
end;
x = dec2logi(x,nx);
end; % if nx==1
% retain generating polynomial, parity check and error parity table
persistent g h errtab
if isempty(g)
% generating polynomial: g = 1+x^2+x^4+x^5+x^6+x^10+x^11
g=zeros(1,12);g([0,2,4,5,6,10,11]+1)=1;
end;
if isempty(h)
% parity check polynomial is the other factor of x^23-1 (mod 2)
% h = (1+x)*(1+x+x^5+x^6+x^7+x^9+x^11)
h=mod(deconv([1,zeros(1,22),1],g),2);
end;
if isempty(errtab),
% set up the error table for decoding: determine all errors of up to 3 bits
% the error will be determined by simple table lookup of the parity check
% the table is a sparse array indexed by the integer representation of the parity bits
errtab = geterrtab(h,3,23);
end;
% now encode or decode the array
if enc
y = circprod(x,g,23);
if ext
% if using extended code, append a checksum bit (row sums mod 2)
y = [y,mod(sum(y,2),2)];
end;
else
if ext
% drop the checksum bit if using extended code
c = x(:,24);
x = x(:,1:23);
end;
% check parity
p = circprod(x,h,23);
% look up parity in error table
err = dec2logi(errtab(logi2dec(p)+1),23);
y = zeros(size(x,1),12);
for i=1:size(x,1),
y(i,:) = mod(deconv(x(i,:)-err(i,:),g),2);
end;
if ext
% if using extended code, append the checksum bit error
err = [err,mod(c-sum(y,2),2)];
end;
end;
% convert back to integers, if necessary
if intout,
y = logi2dec(y);
if nargout == 2,
err = logi2dec(err);
end;
end;
end % main function golaycodec
function t=geterrtab(h,k,n)
% produces lookup table of parity checks of all errors with weight <=k on n digits
% y is a sparse matrix: index is 1 + (h*err converted to integer)
% and y(i) is (err) converted to binary
%
% h must be row vector of length <= n
% also 0 <= k <= n.
m=length(h);
if (2^n-1)>intmax,
error('golaycodec:geterrtab: too many digits to convert to integer');
end;
if m>n,
error('golaycodec:geterrtab: parity check polynomial is too long');
end;
t = sparse(1,1,0,2^n,1); % initialize with 0*h+1 -> 0
for i=1:k,
ind = nchoosek(1:n,i); % all possible combinations of i error digits
nck = size(ind,1);
errbin = zeros(nck,n); % error in binary format
errbin(sub2ind([nck,n],repmat((1:nck)',i,1),ind(:)))=1; % in row j place 1's in columns ind(j,:)
errint = logi2dec(errbin); % error in integer format
errh = circprod(errbin,h,n); % do parity check using h
t(logi2dec(errh)+1)=errint;
end;
end % helper function geterrtab(...)
function z=circprod(x,h,n)
% does parity check of rows of x using polynomial h
m=size(x,1);
% do circular convolution by fourier transform (fft acts on columns)
% (could also use conv and wrap the tail elements, but conv isn't vectorized)
z = mod(round(ifft(fft(x',n).*repmat(fft(h',n),1,m))'),2);
end % helper function paritycheck
function y=logi2dec(x)
% converts each row of logical array from binary to a decimal number
n=size(x,2);
y=x*(2.^(n-1:-1:0))';
end % helper function logi2dec
function y=dec2logi(x,n)
% converts vector of integers to array of logicals (binary representation)
y=rem(floor(x(:)*2.^(1-n:0)),2); % see dec2bin.m (base matlab) for details
end % helper function dec2logi
