fc6c4acd79
Makefiles modified to generate a Linux SICServer. Other files modified to get rid of irritating compiler warning messages.
133 lines
2.9 KiB
C
133 lines
2.9 KiB
C
/*
|
|
*-----------------------------------------------------------------------------
|
|
* file: matdurbn.c
|
|
* desc: Levinson-Durbin algorithm
|
|
* by: ko shu pui, patrick
|
|
* date:
|
|
* revi: 21 may 92 v0.3
|
|
* ref:
|
|
*
|
|
* [1] "Fundementals of Speech Signal Processing," Shuzo Saito,
|
|
* Kazuo Nakata, Academic Press, New York, 1985.
|
|
*
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
|
|
#include "matrix.h"
|
|
|
|
/*
|
|
*-----------------------------------------------------------------------------
|
|
* funct: mat_durbin
|
|
* desct: Levinson-Durbin algorithm
|
|
*
|
|
* This function solve the linear eqns Ax = B:
|
|
*
|
|
* | v0 v1 v2 .. vn-1 | | a1 | | v1 |
|
|
* | v1 v0 v1 .. vn-2 | | a2 | | v2 |
|
|
* | v2 v1 v0 .. vn-3 | | a3 | = | .. |
|
|
* | ... | | .. | | .. |
|
|
* | vn-1 vn-2 .. .. v0 | | an | | vn |
|
|
*
|
|
* where A is a symmetric Toeplitz matrix and B
|
|
* in the above format (related to A)
|
|
*
|
|
* given: R = autocorrelated matrix (v0, v1, ... vn) (dim (n+1) x 1)
|
|
* retrn: x (of Ax = B)
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
MATRIX mat_durbin( R )
|
|
MATRIX R;
|
|
{
|
|
int i, i1, j, ji, p, n;
|
|
MATRIX W, E, K, A, X;
|
|
|
|
p = MatRow(R) - 1;
|
|
W = mat_creat( p+2, 1, UNDEFINED );
|
|
E = mat_creat( p+2, 1, UNDEFINED );
|
|
K = mat_creat( p+2, 1, UNDEFINED );
|
|
A = mat_creat( p+2, p+2, UNDEFINED );
|
|
|
|
W[0][0] = R[1][0];
|
|
E[0][0] = R[0][0];
|
|
|
|
for (i=1; i<=p; i++)
|
|
{
|
|
K[i][0] = W[i-1][0] / E[i-1][0];
|
|
E[i][0] = E[i-1][0] * (1.0 - K[i][0] * K[i][0]);
|
|
|
|
A[i][i] = -K[i][0];
|
|
|
|
i1 = i-1;
|
|
if (i1 >= 1)
|
|
{
|
|
for (j=1; j<=i1; j++)
|
|
{
|
|
ji = i - j;
|
|
A[j][i] = A[j][i1] - K[i][0] * A[ji][i1];
|
|
}
|
|
}
|
|
|
|
if (i != p)
|
|
{
|
|
W[i][0] = R[i+1][0];
|
|
for (j=1; j<=i; j++)
|
|
W[i][0] += A[j][i] * R[i-j+1][0];
|
|
}
|
|
}
|
|
|
|
X = mat_creat( p, 1, UNDEFINED );
|
|
for (i=0; i<p; i++)
|
|
{
|
|
X[i][0] = -A[i+1][p];
|
|
}
|
|
|
|
mat_free( A );
|
|
mat_free( W );
|
|
mat_free( K );
|
|
mat_free( E );
|
|
return (X);
|
|
}
|
|
|
|
/*
|
|
*-----------------------------------------------------------------------------
|
|
* funct: mat_lsolve_durbin
|
|
* desct: Solve simultaneous linear eqns using
|
|
* Levinson-Durbin algorithm
|
|
*
|
|
* This function solve the linear eqns Ax = B:
|
|
*
|
|
* | v0 v1 v2 .. vn-1 | | a1 | | v1 |
|
|
* | v1 v0 v1 .. vn-2 | | a2 | | v2 |
|
|
* | v2 v1 v0 .. vn-3 | | a3 | = | .. |
|
|
* | ... | | .. | | .. |
|
|
* | vn-1 vn-2 .. .. v0 | | an | | vn |
|
|
*
|
|
* domain: where A is a symmetric Toeplitz matrix and B
|
|
* in the above format (related to A)
|
|
*
|
|
* given: A, B
|
|
* retrn: x (of Ax = B)
|
|
*
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
MATRIX mat_lsolve_durbin( A, B )
|
|
MATRIX A, B;
|
|
{
|
|
MATRIX R, X;
|
|
int i, n;
|
|
|
|
n = MatRow(A);
|
|
R = mat_creat(n+1, 1, UNDEFINED);
|
|
for (i=0; i<n; i++)
|
|
{
|
|
R[i][0] = A[i][0];
|
|
}
|
|
R[n][0] = B[n-1][0];
|
|
|
|
X = mat_durbin( R );
|
|
mat_free( R );
|
|
return (X);
|
|
}
|