On 9/10/13 3:36 PM, Randy Yates wrote:
> ...
> I could just write it by the time you and Robert go through your old
> floppies.
sorry Randy. you need to be more explicit (that you want something) for
me to get off my bum. ya gotta unwrap the lines.
have phun.
--
r b-j rbj@audioimagination.com
"Imagination is more important than knowledge."
PAGE 255
IF 0
/**********************************************************************
* *
* Durbin Algorithm for LPC Coefficients *
* *
* input: *
* n = LPC order *
* r -> r[0] to r[n] autocorrelation values *
* *
* output: *
* a -> a[0] to a[n] LPC coefficients, a[0] = 1.0 *
* k -> k[0] to k[n] reflection coefficients, k[0] = unused *
* *
* *
**********************************************************************/
#define N 17
durbin(const int n, const double* r, double* a, double* k)
{
double a_temp[N], alpha, epsilon; /* n <= N = constant */
int i, j;
k[0] = 0.0; /* unused */
a[0] = 1.0;
a_new[0] = 1.0; /* unnecessary but
consistent */
alpha = r[0];
for(i=1; i<=n; i++)
{
epsilon = r[i]; /* epsilon = a[0]*r[i]; */
for(j=1; j<i; j++)
{
epsilon += a[j]*r[i-j];
}
a[i] = k[i] = -epsilon/alpha;
alpha = alpha*(1.0 - k[i]*k[i]);
for(j=1; j<i; j++)
{
a_temp[j] = a[j] + k[i]*a[i-j]; /* update a[] array into
temporary array */
}
for(j=1; j<i; j++)
{
a[j] = a_temp[j]; /* update a[] array */
}
}
}
/*
Less concise version with inside for() loops executed
at least once in the same manner as Fortran version.
*/
durbin(const int n, const double* r, double* a, double* k)
{
double a_temp[N], alpha, epsilon; /* n <= N = constant */
int i, j;
k[0] = 0.0; /* unused */
a[0] = 1.0;
a_new[0] = 1.0; /* unnecessary but
consistent */
a[1] = k[1] = -r[1]/r[0];
alpha = r[0]*(1.0 - k[1]*k[1]);
for(i=2; i<=n; i++)
{
epsilon = 0.0;
for(j=0; j<i; j++)
{
epsilon += a[j]*r[i-j];
}
a[i] = k[i] = -epsilon/alpha;
alpha = alpha*(1.0 - k[i]*k[i]);
for(j=1; j<i; j++)
{
a_temp[j] = a[j] + k[i]*a[i-j]; /* update a[..] array
first into temporary array */
}
for(j=i-1; j>=1; j--)
{
a[j] = a_temp[j]; /* update a[..] array */
}
}
}
(Nearly) original FORTRAN source program: [1]
SUBROUTINE DURBIN(N,R,K,A)
C
C THIS ROUTINE USES THE DURBIN ALGORITHM TO
C TRANSFORM THE AUTOCORRELATION COEFFICIENTS R(0) TO R(N)
C TO THE REFLECTION COEFFICIENTS K(1) TO K(N) AND TO THE FILTER
C COEFFICIENTS A(0) TO A(N) (A(0)=1).
C
REAL K(1..N),R(0..N),A(0..N)
C
C
REAL A_TEMP(0..N),EPSILON
C
A(0) = 1.0
A_TEMP(0) = 1.0
K(1) = -R(1)/R(0)
A(1) = K(1)
ALPHA = R(0)*(1.0 - K(1)*K(1))
DO 400 I = 2,N
EPSILON = 0.0
DO 100 J = 0,I-1
EPSILON = EPSILON + A(J)*R(I-J)
100 CONTINUE
K(I) = - EPSILON / ALPHA
A(I) = K(I)
ALPHA = ALPHA*(1.0 - K(I)*K(I))
DO 200 J = 1,I-1
A_TEMP(J) = A(J) + K(I)*A(I-J)
200 CONTINUE
DO 300 J = 1,I-1
A(J) = A_TEMP(J)
300 CONTINUE
400 CONTINUE
RETURN
END
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.
ENDIF
;
; i/o file definitions
;
file_data_in EQU y:$FFFF
file_data_out EQU y:$FFFE
LPC_order EQU 10 ; 10th order LPC analysis
ORG x:
r DS LPC_order+1 ; autocorrelation coefficients input
ORG x:
k DC 0.0 ; define k[0] = 0.0 (unused)
DS LPC_order ; reflection coefficients output
ORG y:
a_coef DC 0.99999988 ; define a[0] ~= 1.0
DS LPC_order ; LPC filter coef array
ORG x:
a_temp DC 0.99999988 ; define a_temp[0] ~= 1.0
DS LPC_order ; temporary LPC filter coef array
;
; The following registers correspond to the C variables:
;
; r0 - r[i] r4 - a_coef[i]
; r1 - k[i-1] r5 - a_coef[j-1]
; r2 - loop counter r6 - a_coef[i-j+1]
; r3 - r[i-j+1], a_new[j-1]
;
;
move #r,r0 ; r0 -> r[0] input array to algorithm
move #k+1,r1 ; r1 -> k[1] output array of algorithm
move #a_coef+1,n4 ; n4 -> a[1]
move #a_temp,n3 ; n3 -> a_temp[0]
move #LPC_order-1,n2
jsr durbin
durbin
move n4,r4 ; r4 -> a[1]
move n4,r5 ; r5 -> a[1]
move x:(r0)+,y1 ; r[0]
move x:(r0)+,a ; r[1]
move #2,n6 ; n6 = 2
; signed division
abs a a,b ; a = abs(r[1]), b = r[1]
and #$FE,ccr
do #24,_div1
div y1,a ; compute abs(r[1])/r[0]
_div1
neg a a0,x0
eor y1,b a0,a ; sign bit, N = 1 if y1 and b are different signs
tmi x0,a ; negate quotient if r[1] and r[0] have same sign
; a = -r[1]/r[0]
move #<-1.0,b
move a,x0 a,y:(r4)+ ; x0 = k[1], a[1] = k[1] = -r[1]/r[0]
macr x0,x0,b ; b = -1.0 + k[1]^2
move #0,r2 ; r2 = 0
move b,x1 ; x1 = -(1.0 - k[1]^2)
mpyr -x1,y1,b x0,x:(r1)+ ; b = r[0]*(1.0 - k[1]^2), write k[1]
move (r5)- ; r5 -> a[0]
move b,y1 ; alpha = r[0]*(1.0 - k[1]^2)
do n2,_outer_loop ; for (i=2; i<=LPC_order; i++) {
move r0,r3 ; r3 -> r[i]
move (r0)+ ; r0 -> r[i+1]
move (r2)+ ; r2 = i-1
nop
clr a x:(r3)-,x0 y:(r5)+,y0 ; epsilon = 0;
do r2,_loop1 ; for (j=0; j<i-1; j++) { // i-1 iterations
mac x0,y0,a x:(r3)-,x0 y:(r5)+,y0 ; epsilon += a[j]*r[i-j];
_loop1 ; }
mac x0,y0,a ; epsilon += a[i-1]*r[1];
; signed division
abs a a,b ; a = abs(epsilon), b = epsilon
and #$FE,ccr
do #24,_div2
div y1,a ; compute abs(epsilon)/alpha
_div2
neg a a0,x0
eor y1,b a0,a ; sign bit, N = 1 if y1 and b are different signs
tmi x0,a ; negate quotient if epsilon and alpha have same sign
; a = k[i] = -epsilon/alpha
move #<-1.0,b
move a,x0 a,y:(r4)+ ; x0 = k[i], a[i] = k[i]
macr x0,x0,b r4,r6 ; b = -1.0 + k[i]^2, r6 -> a[i+1]
move n4,r5 ; r5 -> a[1]
move b,x1 ; x1 = -(1.0 - k[i]^2)
mpyr -x1,y1,b x0,x:(r1)+ ; b = alpha*(1.0 - k[i]^2), store k[i]
move (r6)-n6 ; r6 -> a[i-1]
move b,y1 ; alpha <- alpha*(1.0 - k[i]^2)
move y:(r6)-,y0 ; y0 = a[i-1]
move y:(r5)+,a ; a = a[1]
move n3,r3 ; r3 -> a_temp[0]
do r2,_loop2 ; for (j=1; j<=i-1; j++) { // i-1 iterations
macr x0,y0,a y:(r6)-,y0 ; // a[j] + k[i]*a[i-j], y0 =
a[i-(j+1)]
move (r3)+ ; // r3 -> a_temp[j]
move a,x:(r3) y:(r5)+,a ; a_temp[j] = a[j] + k[i]*a[i-j] // a =
a[j+1]
_loop2 ; }
move (r5)- ; r5 -> a[i-1]
do r2,_loop3 ; for (j=i-1; j>=1; j--) { // i-1 iterations
move x:(r3)-,y0
move y0,y:(r5)- ; a[j] = a_temp[j]
_loop3 ; }
nop ; r5 -> a[0]
outer_loop ; }
rts