%DISTANCE : perform spectral distortion measure.
%Function: y=distance2(a,b,xa,zi,method) returns the distortion measure between
%            two frames of speech data.
%
%INPUT:  a == LPC coefficients for the reference speech
%        b == LPC coefficients for the test speech
%        xa == reference speech
%        zi == initial value of the Linear Prediction filter
%        method==1 calculate the ITAKURA MEASURE.
%              ==2 calculate the WEIGHTED CEPTRAL DISTANCE.
%              ==3 calculate the LOG SPECTRAL DISTANCE.
%
%OUTPUT:  y == Integrate{ sqare[ log(A((w))-log(B(w)) ] };
%      or
%      or y == beta/alpha-1; beta=filter(b,1,xa)^2; alpha=filter(a,1,xa);
%
%REFERENCE: see "Fundamentals of speech recognition" by Rabiner & Juang
%           pg. 195 for details.
       
function y=distance2(a,b,xa,zi,method)


if method==1
%%%%%%%%%%%%%%%% ITAKURA MEASURE

  xa=xa(:)';
  err=filter(b,1,xa,zi);
  beta=err*err';
  err=filter(a,1,xa,zi);
  alpha=err*err';
  y=beta/alpha-1;
  if y<0
     y=abs(y)+1;
  end

elseif method==2
%%%%%%%%%%%%%%%% WEIGHTED CEPTRAL DISTANCE

   % 1. construct the cepstral coefficeint for the reference speech
   order=length(a)-1;
   c1=zeros(1,2*order);
   a=a(:);
   a=[a;zeros(order,1)];
   for i=1:order
       c1(i)=a(i+1)+(1:i-1)/i.*c1(1:i-1)*a(i:-1:2);
   end
   for i=order+1:2*order
       c1(i)=(1:i-1)/i.*c1(1:i-1)*a(i:-1:2);
   end

   % 2. construct the cepstral coefficient for the test speech
   c2=zeros(1,2*order);
   b=b(:);
   b=[b;zeros(order,1)];
   for i=1:order
       c2(i)=b(i+1)+(1:i-1)/i.*c2(1:i-1)*b(i:-1:2);
   end
   for i=order+1:2*order
       c2(i)=(1:i-1)/i.*c2(1:i-1)*b(i:-1:2);
   end

   % 3. construct the weighting function. See function (4.36) at pg. 169
   w=zeros(1,2*order);
   for i=1:order
       w(i)=1+order/2*sin(i*pi/order);
       w(i+order)=0;
   end

   % 4. calculate the LIFTER CEPTRAL DISTANCE

   y=w.*(c1-c2)*(c1-c2)';

elseif method==3

%%%%%%%%%%%% LOG SPECTRAL DISTANCE

   Ha=freqz(1,a,128);
   Hb=freqz(1,b,128);
   HH=log( abs(Ha) )-log( abs(Hb) );
   HH=HH(:)';
   y=(HH*HH')/128;

end