Name: DURBIN1
         Type: Assembler Program
      Version: 1.2
 Date Entered: 21-May-87
  Last Change: 3-June-87

Description: Durbin Solution for PARCOR (LPC) Coefficients

This program uses the Durbin algorithm to find a set of PARCOR coefficients
given an autocorrelation vector. These PARCOR coefficients can be used in
speech analysis/synthesis systems such as vocoding, stochastic speech
coding, multipulse speech coding, speech recognition, etc.

This program will solve any order system by equating the order "NK" to the
number of K coefficients desired. The program is written as an illustrative
example and the user may change the routine as needed. The program will:

  1. Read a data file to determine how many sets of PARCOR coefficients
     to calculate.
  2. Read the autocorrelations from a data file.
  3. Perform the Durbin algorithm to find the K's.
  4. Output the K coefficients to a data file.
  5. Continue steps 2,3,4 for the number of specified sets of coefficients.

The explanation of the algorithm is beyond the scope of this help file but
the coding of the algorithm can be easily understood by examining the C
routine below (DURBIN). The C coding of the Durbin algorithm is a modified
version of the original FORTRAN source code given on page 122 of the book
entitled Practical Approaches to Speech Coding by Panos E. Papamichalis,
(Prentice - Hall, 1987).

Approximate execution times for DURBIN1 on a 20.5 MHz DSP56000/1 are:

    NK = 8  (8th order) : 63.5  uS
    NK = 10 (10th order): 85.6  uS
    NK = 16 (16th order): 165.7 uS
  
Original FORTRAN source program: [1]

       SUBROUTINE DURBIN(R,N,K,A)
C
C      THIS ROUTINE USES THE DURBIN ALGORITHM TO
C      TRANSFORM THE AUTOCORRELATION COEFFICIENTS R(1)-R(N+1)
C      TO THE REFLECTION COEFFICIENTS K(1)-K(N) AND TO THE FILTER
C      COEFFICIENTS A(1)-A(N+1) (A(1)=1).
C
       REAL    K(1),R(1),A(1)
C
       A(1) = 1
       K(1) = -R(2) / R(1)
       ALPHA = R(1) * (1 - K(1)**2)
       A(2) = K(1)
       DO 200 I = 2,N
            ERROR = 0.
            DO 100 J = 1,I
                 ERROR = ERROR + A(J) * R(I+2-J)
100         CONTINUE
            K(I) = - ERROR / ALPHA
            ALPHA = ALPHA * (1 - K(I)**2)
            A(I+1) = K(I)
            DO 200 J = 2,I
                 A(J) = A(J) + K(I) * A(I+2-J)
200    CONTINUE
       RETURN
       END

Modified C source for direct DSP56000/1 conversion:

/********************************************************
*                                                       *
* Durbin Algorithm for LPC Coefficients                 *
*                                                       *
********************************************************/
durbin(r,n,k,a)
int n;
double r[], k[], a[];
{
	int i, j;
	double anew[NK+1], alpha, error;

	a[0] = 1.0;
	anew[0] = 1.0;
	k[0] = -r[1]/r[0];
	alpha = r[0] * (1.0 - (k[0] * k[0]));
	a[1] = k[0];
	for( i = 2; i <= n; i++)
		{
		error = 0.0;
		for( j = 1; j <= i; j++)
			error = error + (a[j-1] * r[i-j+1]);
		k[i-1] = -(error / alpha);
		alpha = alpha * (1.0 - (k[i-1] * k[i-1]));
		a[i] = k[i-1];
		for( j = 2; j <= i; j++)
			anew[j-1] = a[j-1] + (k[i-1] * a[i-j+1]);
		for( j = 2; j <= i; j++)
			a[j-1] = anew[j-1];
		}
}

References

[1] Papamichalis, Panos E., Practical Approaches to Speech Coding, Englewood
    Cliffs, NJ: Prentice - Hall, 1987.
[2] Rabiner, L.R. and R.W. Schafer, Digital Processing of Speech Signals,
    Englewood Cliffs, NJ: Prentice - Hall, 1978.
[3] Makhoul, John, "Linear Prediction: A Tutorial Review," Proceedings of
    the IEEE, Volume 63 (April,1975): pp. 561 - 580.