function C = vqcodebook(X, N, I)
% C= vqcodebook(X, N, I)  Design a VQ codebook
%       X is a set of data rows.  Design an N-entry codebook 
%       by iterating I times.
% 2001-03-19 dpwe@ee.columbia.edu

if nargin < 2
  N = 64;
end
if nargin < 3
  I = 20;
end

[ndata, ndims] = size(X);

% Initialize codebook with random vectors
rp = randperm(ndata);
C = X(rp(1:N),:);

% Iteratively update codebook 

for it = 1:I
  
  Dtot = 0;
  
  % Assign each datum to codeword
  cofx = zeros(ndata,1);
  for da = 1:ndata
    dists = sum( (((ones(N,1)*X(da,:)) - C).^2)' );
    best = find(dists == min(dists));
    cofx(da) = best(1);
    Dtot = Dtot + min(dists);
  end
  
  % Report distortion
  Drms = sqrt(Dtot/ndata);
  disp(['RMS error before it ',num2str(it),' = ',num2str(Drms)]);
  
  % Re-estimate each codeword as mean
  for cw = 1:N
    da = find(cofx == cw);
    if length(da) > 0
      C(cw,:) = mean(X(da,:));
    end
  end
  
end