function DIST = lpcdistance(X1, X2, FRAMESIZE, ORDER)
% Calculate the average symmetric Itakura distance between the signals X1 and X2.
% 
% DIST = lpcdistance(X1, X2, FRAMESIZE, ORDER);
%
% DIST (scalar) = the average Itakura-Saito distance.
% X1 = the first signal.
% X2 = the second signal.
% FRAMESIZE = the frame length to use (default is 160 samples).
% ORDER = the LPC order (default is 10).
%

% Make sure the inputs are column vectors
X1 = X1(:);
X2 = X2(:);

% Figure out the frame length and LPC order
if nargin<4,
  ORDER = 10;
  if nargin<3,
    FRAMESIZE = 160;
  end;
end;
  
% Create a Hamming window called WINDOW
WINDOW = hamming(FRAMESIZE);

% Calculate the number of frames.  Call it NFRAMES.
DURATION = min(length(X1),length(X2));
NFRAMES = floor(DURATION/FRAMESIZE);

% MAIN LOOP
DIST = 0;
NDIST = 0;
for i=1:NFRAMES,

  % Get the current frame of speech from X1 and X2
  x1 = X1((i-1)*FRAMESIZE+1:i*FRAMESIZE) .* WINDOW;
  x2 = X2((i-1)*FRAMESIZE+1:i*FRAMESIZE) .* WINDOW;
  
  % Check to make sure that both signals are nonzero
  if ((sum(x1.^2) > 0) & (sum(x2.^2) > 0)),
    NDIST = NDIST + 1;
    
    % Calculate LPC coefficients
    A1 = lpc(x1, ORDER);
    A2 = lpc(x2, ORDER);
  
    % Calculate the autocorrelation
    r1 = xcorr(x1);
    r1 = r1((FRAMESIZE+1):(FRAMESIZE+1+ORDER));
    R1 = toeplitz(r1);
    r2 = xcorr(x2);
    r2 = r2((FRAMESIZE+1):(FRAMESIZE+1+ORDER));
    R2 = toeplitz(r2);
    
    % Calculate d12 and d21, the two different Itakura distances
    d12 = 10*log10((A2*R1*A2') / (A1*R1*A1'));
    d21 = 10*log10((A1*R2*A1') / (A2*R2*A2'));
  
    % Accumulate a dB average
    DIST = DIST + (d12 + d21)/2;
    
  end;
end;


DIST = DIST / NDIST;


