00001 #include <math.h> 00002 #include "c_lpc10.h" 00003 00004 /**********************************************************/ 00005 /* Invert a covariance matrix using Choleski decomposition method 00006 Inputs: 00007 . psi[ORDER] - Column vector to be predicted 00008 In/Outputs: 00009 . phi[ORDER,ORDER] - Covariance matrix 00010 Outputs: 00011 . rc[ORDER] - Pseudo reflection coefficients 00012 Constants: 00013 . ORDER - Analysis order */ 00014 00015 00016 VOID invert( FLOAT phi[ORDER][ORDER], FLOAT psi[ORDER], FLOAT rc[ORDER] ) 00017 { 00018 #define EPS ((FLOAT)1.0e-10) 00019 INDEX i, j, k; 00020 FLOAT save; 00021 00022 /* Decompose PHI into V * D * V' where V is a triangular matrix whose 00023 main diagonal elements are all 1, V' is the transpose of V, and 00024 D is a vector. Here D(n) is stored in location V(n,n). */ 00025 for (j = 0; j < ORDER; j++) { 00026 for (k = 0; k < j; k++) { 00027 save = phi[j][k] * phi[k][k]; 00028 for (i = j; i < ORDER; i++) 00029 phi[i][j] -= phi[i][k] * save; 00030 } 00031 00032 /* Compute intermediate results, which are similar to RC's */ 00033 if ((FLOAT)fabs(phi[j][j]) < EPS) { /* if algorithm terminate early...*/ 00034 for (i = j; i < ORDER; i++) 00035 rc[i] = (FLOAT)0.0; /* ...zero out higher order RC's */ 00036 return; 00037 } 00038 00039 rc[j] = psi[j]; 00040 for (k = 0; k < j; k++) 00041 rc[j] -= rc[k] * phi[j][k]; 00042 00043 phi[j][j] = (FLOAT)1. / phi[j][j]; 00044 rc[j] *= phi[j][j]; 00045 00046 /* Clipping is not necessary. Later we check for stability. */ 00047 /* 00048 if (rc[j] > (FLOAT)0.999) 00049 rc[j] = (FLOAT)0.999; 00050 else if (rc[j] < (FLOAT)-0.999) 00051 rc[j] = (FLOAT)-0.999; 00052 */ 00053 } 00054 } 00055 00056 /**********************************************************/
1.5.6