function code = gilbertcode(bits, words, corrects)
%gilbertcode( bits, words, corrects ):
% [0] Let N = 1<<bits, let radius = 2*corrects
% [1] For n=0 to N-1, allocate tag[n] with tag[n]=LIVE
% [2] Allocate c[0],...,c[words-1] and let c[0] = 0
% [3] For j=0 to words-1, do [4] to [7]
% [4]   For n=c[j]+1 to N-1 do [5]
% [5]     If dist(n,c[j])<=radius, then let tag[n] = DEAD
% [6]   For n=c[j]+1 to N-1 do [7]
% [7]     If tag[n]==LIVE, then let c[j+1]=n and goto [3]
% [8] For j=0 to words-1, print c[j]
  N=2^bits; % Number of available bitstrings
  radius =  2*corrects; % use this distance between codewords
  tag=ones(1,N); % array to indicate codewords LIVE=1
  c=zeros(1,words); % array to hold codewords; c(1)=0 always
  for j=1:words-1 % current codeword 
    for n=(1+c(j)):(N-1) % check remaining bitstrings
    % Hamming distance.  Add N=10...0 (base 2) so that dec2bin()
    % retuns equal-length (=1+bits) bitstrings.
      hdist = sum((dec2bin(N+n))!= dec2bin(N+c(j))); 
      if hdist <= radius
	tag(n+1)=0; % indicate bitstring n is DEAD=0
      end
    end % bitstring removing loop
    for n=(1+c(j)):(N-1) % search remaining bitstrings for a codeword
      if tag(n+1)==1 % Found a LIVE=1 one!
	c(j+1)=n; % store bitstring n as a codeword
	break; % Exit this search loop
      end
    end % surviving codeword search loop
  end
  code=dec2bin(c); % convert codewords to bitstrings
  return;
end

