00001 /**********************************************************/ 00002 /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/ 00003 /* 00004 Copyright: 1994 - Grupo de Voz (DAET) ETSII/IT-Bilbao 00005 00006 Nombre fuente................ SPL5a.C 00007 Nombre paquete............... SPL 00008 Lenguaje fuente.............. C (Borland C/C++ 3.1) 00009 Estado....................... Completado 00010 Dependencia Hard/OS.......... - 00011 Codigo condicional........... NDEBUG, LP_NEGSUM 00012 00013 Codificacion................. Borja Etxebarria 00014 00015 Version dd/mm/aa Autor Proposito de la edicion 00016 ------- -------- -------- ----------------------- 00017 1.1.1 30/07/95 Borja scope funciones explicito 00018 1.1.0 08/12/94 Borja revision general (tipos,idx,nel,tnel...) 00019 1.0.1 07/12/93 Borja soporte LP_NEGSUM_NEG() y LP_POSSUM_NEG() 00020 1.0.0 06/07/93 Borja Codificacion inicial. 00021 00022 ======================== Contenido ======================== 00023 Algoritmos de prediccion lineal (LP). 00024 00025 El funcionamiento de todos los algoritmos LP esta determinado 00026 por el #define LP_NEGSUM en SPL.H : 00027 00028 + Si LP_NEGSUM esta definido: 00029 . H(z) = G / ( 1 - sum{ak * z**-k} ) 00030 . Estructuras lattice propagan error con -Ki 00031 . s'(n) = sum{ak*s(n-k)} 00032 00033 + Si por el contrario LP_NEGSUM no esta definido: 00034 . H(z) = G / ( 1 + sum{ak * z**-k} ) 00035 . Estructuras lattice propagan error con +Ki 00036 . s'(n) = - sum{ak*s(n-k)} 00037 00038 (x**y indica x elevado a y) 00039 en cualquiera de los dos casos, en un predictor de orden p se 00040 cumple ap=Kp, siendo ap el ultimo coeficiente del predictor 00041 lineal, y Kp el correspondiente coeficiente de reflexion. 00042 Realmente los coeficientes Ki *son* los coeficientes de reflexion 00043 del modelo del tubo sin perdidas solo si LP_NEGSUN *no* esta 00044 definido. Si LP_NEGSUM *si* esta definido, entonces los Ki son 00045 realmente los coeficientes de reflexion autenticos cambiados de signo. 00046 00047 Se han definido dos macros para gestionar automaticamente los cambios 00048 de signo en funcion de LP_NEGSUM: 00049 - LP_NEGSUM_NEG(n) convierte {n} en -{n} solo si LP_NEGSUM esta definido 00050 - LP_POSSUM_NEG(n) convierte {n} en -{n} solo si LP_NEGSUM NO esta definido 00051 00052 Para un analisis LP de orden {p} (funciones lpa_???), definimos: 00053 00054 - {vac} es el vector de las {p}+1 primeras autocorrelaciones, desde 00055 R(0) hasta R(p) ==> {p}+1 elementos, vac[0]=R(0). 00056 00057 - {mcv} es la 'pseudo matriz' de covarianzas de orden {p}, obtenida 00058 con las rutinas xcovm_???(). Este vector tiene tnel_xcovm(p) 00059 elementos. Tambien valen las rutinas xwcovm_??? (Metodo mejorado 00060 segun Atal), que utiliza covarianzas modificadas aunque son mucho 00061 mas lentas. 00062 00063 - {vki} es un vector de resultado donde se devuelven los coeficientes 00064 de refrexion k(1) hasta k(p) ==> {p} elementos, vki[0]=k(1). 00065 00066 - {vai} es un vector de resultado donde se devuelven los coeficientes 00067 LPC desde a(1) hasta a(p) ==> {p} elementos, vai[0]=a(1). 00068 00069 - {vti} es un vector de resultado donde se devuelven las areas 00070 del tubo sin perdidas A(1) hasta A(p) ==> {p} elementos, vti[0]=A(1). 00071 00072 - {p} es el orden del analisis. Aunque algunas funciones soportan {p}==0, 00073 la mayor parte no funciona correctamente, por tanto {p}>0. 00074 00075 - {vtmp} es un vector temporal que debe pasarse a la rutina para 00076 realizar operaciones temporales, y que debe tener, como minimo, 00077 espacio para el numero de ELEMENTOS indicado en cada caso. 00078 Se indican elementos, no bytes! Cada rutina de analisis tiene 00079 asociada otra del mismo nombre y prefijo 'tnel_' que devuelve 00080 el numero de elementos SPL_FLOAT que debe tener este vector para un 00081 orden {p} concreto. 00082 00083 Las funciones lpa_??? {devuelven} la potencia del error de prediccion 00084 E=G**2 (si H(z)=G/A(z) es el filtro reconstructor). 00085 00086 Hay varias funciones para cada metodo de calculo: ???_ka, ???_k, 00087 ???_a, ???_e, ???_a_ne... que pueden devolver los coeficientes de 00088 reflexion(k) los coeficientes lpc(a), solo la potencia del error (e), 00089 sin la potencia del error (ne), o los lpc y los coef. reflexion a 00090 la vez (ka), etc. Las combinaciones en algunos casos no implementadas 00091 se pueden obtener con llamadas a otras funciones sin perdida de 00092 eficiencia, y por eso no se han implementado. 00093 00094 Ref prediccion lineal, LPC, LSF/LSP: 00095 . Digital Processing of Speech Signals 00096 . L.R. Rabiner / R.W Schafer 00097 . Prentice Hall 00098 00099 . J. Makhoul 00100 . Linear Prediction: A Tutorial Review 00101 . Proc. IEEE vol 63, pp. 561-580, Apr. 1075 00102 00103 . PUBLIC. Delsarte, Y. V. Genin 00104 . The Split Levinson Algorithm 00105 . IEEE trans. ASSP, vol 34, n. 3, Jun. 1986 00106 00107 . Furui 00108 . ?????? desaparecido en combate 00109 00110 . K.K Paliwal, B.S. Atal 00111 . Efficient Vector Quantization of LPC Parameters at 24 Bits/Frame 00112 . IEEE. Trans. Speech, and Audio Processing, Vol 1, N.1 Jan 93 00113 00114 . F.K. Soong B.H. Juang 00115 . Optimal Quantization of LSP Parameters 00116 . IEEE. Trans. Speech, and Audio Processing, Vol 1, N.1 Jan 93 00117 00118 Ref. covarianzas modificadas, metodo Cholesky estabilizado: 00119 . S. Singhal and B. S. Atal 00120 . Improving performance of multi-pulse LPC coders at low bit rates. 00121 . Proc. Int. Conf. ASSP vol 1, Paper 1.3 00122 00123 Ref. metodo Cholesky estabilizado, noise shapping: 00124 . B. S. Atal, M. R. Schroeder 00125 . Predictive Coding of Speech Signals and Subjective error criteria 00126 . IEEE Trans. ASSP, vol 27, N.3, Jun. 1979 00127 00128 . B. S. Atal 00129 . Predictive Coding of Speech at Low bit Rates 00130 . IEEE Trans. Comm. Vol COM-30, N.4, April 1982 00131 00132 00133 Definir NDEBUG para desconectar la validacion de parametros 00134 con assert(). No definir este simbolo mientras se depuren 00135 aplicaciones, ya que aunque las funciones son algo mas lentas, 00136 el chequeo de validacion de parametros resuelve muchos problemas. 00137 =========================================================== 00138 */ 00139 /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/ 00140 /**********************************************************/ 00141 00142 #include "spli.h" 00143 00144 /**********************************************************/ 00145 /* Metodo de Durbin para las correlaciones. 00146 A partir de las {p}+1 primeras autocorrelaciones {vac} 00147 calcula los coeficientes de reflexion {vki} y los coeficientes 00148 lpc {vai} para analisis de orden {p}. 00149 La funcion {devuelve} la potencia del error de prediccion E=G**2 */ 00150 00151 PUBLIC SPL_FLOAT XAPI lpa_cor_dur_ka( SPL_pFLOAT vac, SPL_pFLOAT vki, 00152 SPL_pFLOAT vai, SPL_INT p, SPL_pFLOAT vtmp ) 00153 { 00154 SPL_INT i,j; 00155 SPL_FLOAT E; /*error prediccion E(i)*/ 00156 SPL_FLOAT sum, ki; 00157 SPL_pFLOAT vo, vn, vx; /* 2 vectores: old, new , y uno auxiliar */ 00158 assert(p>0); 00159 00160 /* inicializa vectores old(vo) y new(vn) en funcion del order p */ 00161 if (p&1) { /* si el orden es impar */ 00162 vn = vai; 00163 vo = vtmp; 00164 } /* al final, vn encaja con vai */ 00165 else { 00166 vo = vai; 00167 vn = vtmp; 00168 } 00169 00170 E = (vac[0]+vac[1])* 00171 (1.0+LP_NEGSUM_NEG(vn[0]=vki[0]=LP_POSSUM_NEG(vac[1]/vac[0]))); 00172 00173 for ( i=1; i<p; i++ ) { 00174 __swap(vn,vo,vx); /* a(i) ---> a(i-1) */ 00175 00176 /* Sumatorio[j=1,j=i-1]( R[i-j] * a(i-1)[j] ) */ 00177 __xsum(j,0,i,vac[i-j]*vo[j],sum,0.0); 00178 00179 /* a(i)[i] = k[i] = -(R[i]+Sum[...])/E(i-1) */ 00180 vn[i] = vki[i] = ki = LP_POSSUM_NEG(vac[i+1]+LP_NEGSUM_NEG(sum))/E; 00181 00182 /* a(i)[j]=a(i-1)[j]+k[i]*a(i-1)[i-j] */ 00183 for ( j=0; j<i; j++ ) 00184 vn[j] = vo[j] + LP_NEGSUM_NEG(ki*vo[i-j-1]); 00185 00186 /* E(i) = (1-k[i]^2)*E(i-1) */ 00187 E *= (1.0-ki*ki); 00188 } 00189 00190 return E; /* devuelve G*G */ 00191 } 00192 00193 /**********************************************************/ 00194 00195 PUBLIC SPL_INT XAPI tnel_lpa_cor_dur_ka( SPL_INT p ) 00196 { 00197 assert(p>0); 00198 00199 return (p-1); 00200 } 00201 00202 /**********************************************************/ 00203 /* Metodo de Durbin para las correlaciones. 00204 A partir de las {p}+1 primeras autocorrelaciones {vac} 00205 calcula los coeficientes de reflexion {vki} para analisis 00206 de orden {p}. 00207 La {funcion} devuelve la potencia del error de prediccion */ 00208 00209 PUBLIC SPL_FLOAT XAPI lpa_cor_dur_k( SPL_pFLOAT vac, SPL_pFLOAT vki, 00210 SPL_INT p, SPL_pFLOAT vtmp ) 00211 { 00212 SPL_INT i,j; 00213 SPL_FLOAT E; /*error prediccion E(i)*/ 00214 SPL_FLOAT sum, ki; 00215 SPL_pFLOAT vo, vn, vx; /* 2 vectores: old, new , y uno auxiliar */ 00216 assert(p>0); 00217 00218 vo = vtmp; 00219 vn = vtmp+(--p); 00220 00221 E = (vac[0]+vac[1])* 00222 (1.0+LP_NEGSUM_NEG(vn[0]=vki[0]=LP_POSSUM_NEG(vac[1]/vac[0]))); 00223 00224 for ( i=1; i<p; i++ ) { 00225 __swap(vn,vo,vx); /* a(i) ---> a(i-1) */ 00226 00227 /* Sumatorio[j=1,j=i-1]( R[i-j] * a(i-1)[j] ) */ 00228 __xsum(j,0,i,vac[i-j]*vo[j],sum,0.0); 00229 00230 /* a(i)[i] = k[i] = -(R[i]+Sum[...])/E(i-1) */ 00231 vn[i] = vki[i] = ki = LP_POSSUM_NEG(vac[i+1]+LP_NEGSUM_NEG(sum))/E; 00232 00233 /* a(i)[j]=a(i-1)[j]+k[i]*a(i-1)[i-j] */ 00234 for ( j=0; j<i; j++ ) 00235 vn[j] = vo[j] + LP_NEGSUM_NEG(ki*vo[i-j-1]); 00236 00237 /* E(i) = (1-k[i]^2)*E(i-1) */ 00238 E *= (1.0-ki*ki); 00239 } 00240 00241 __xsum(j,0,p,vac[p-j]*vn[j],sum,0.0); 00242 vki[p] = LP_POSSUM_NEG(vac[p+1]+LP_NEGSUM_NEG(sum))/E; 00243 00244 return E*(1.0-vki[p]*vki[p]); /* devuelve G*G */ 00245 } 00246 00247 /**********************************************************/ 00248 00249 PUBLIC SPL_INT XAPI tnel_lpa_cor_dur_k( SPL_INT p ) 00250 { 00251 assert(p>0); 00252 00253 return (2*(p-1)); 00254 } 00255 00256 /**********************************************************/ 00257 /* Metodo de Durbin para las correlaciones. 00258 A partir de las {p}+1 primeras autocorrelaciones {vac} 00259 calcula los coeficientes lpc {vai} para analisis de orden {p}. 00260 La funcion {devuelve} la potencia del error de prediccion */ 00261 00262 PUBLIC SPL_FLOAT XAPI lpa_cor_dur_a( SPL_pFLOAT vac, SPL_pFLOAT vai, 00263 SPL_INT p, SPL_pFLOAT vtmp ) 00264 { 00265 SPL_INT i,j; 00266 SPL_FLOAT E; /*error prediccion E(i)*/ 00267 SPL_FLOAT sum, ki; 00268 SPL_pFLOAT vo, vn, vx; /* 2 vectores: old, new , y uno auxiliar */ 00269 assert(p>0); 00270 00271 /* inicializa vectores old(vo) y new(vn) en funcion del order p */ 00272 if (p&1) { /* si el orden es impar */ 00273 vn = vai; 00274 vo = vtmp; 00275 } /* al final, vn encaja con vai */ 00276 else { 00277 vo = vai; 00278 vn = vtmp; 00279 } 00280 00281 E = (vac[0]+vac[1])*(1.0+LP_NEGSUM_NEG(vn[0]=LP_POSSUM_NEG(vac[1]/vac[0]))); 00282 00283 for ( i=1; i<p; i++ ) { 00284 __swap(vn,vo,vx); /* a(i) ---> a(i-1) */ 00285 00286 /* Sumatorio[j=1,j=i-1]( R[i-j] * a(i-1)[j] ) */ 00287 __xsum(j,0,i,vac[i-j]*vo[j],sum,0.0); 00288 00289 /* a(i)[i] = k[i] = -(R[i]+Sum[...])/E(i-1) */ 00290 vn[i] = ki = LP_POSSUM_NEG(vac[i+1]+LP_NEGSUM_NEG(sum))/E; 00291 00292 /* a(i)[j]=a(i-1)[j]+k[i]*a(i-1)[i-j] */ 00293 for ( j=0; j<i; j++ ) 00294 vn[j] = vo[j] + LP_NEGSUM_NEG(ki*vo[i-j-1]); 00295 00296 /* E(i) = (1-k[i]^2)*E(i-1) */ 00297 E *= (1.0-ki*ki); 00298 } 00299 00300 return E; /* devuelve G*G */ 00301 } 00302 00303 /**********************************************************/ 00304 00305 PUBLIC SPL_INT XAPI tnel_lpa_cor_dur_a( SPL_INT p ) 00306 { 00307 assert(p>0); 00308 00309 return (p-1); 00310 } 00311 00312 /**********************************************************/ 00313 /* Metodo de Durbin para las correlaciones. 00314 A partir de las {p}+1 primeras autocorrelaciones {vac} 00315 la funcion devuelve la potencia del error de prediccion para 00316 analisis de orden {p}. */ 00317 00318 PUBLIC SPL_FLOAT XAPI lpa_cor_dur_e( SPL_pFLOAT vac, SPL_INT p, 00319 SPL_pFLOAT vtmp ) 00320 { 00321 SPL_INT i,j; 00322 SPL_FLOAT E; /*error prediccion E(i)*/ 00323 SPL_FLOAT sum, ki; 00324 SPL_pFLOAT vo, vn, vx; /* 2 vectores: old, new , y uno auxiliar */ 00325 assert(p>0); 00326 00327 vo = vtmp; 00328 vn = vtmp+(--p); 00329 00330 E = (vac[0]+vac[1])*(1.0 + 00331 LP_NEGSUM_NEG(vn[0]=LP_POSSUM_NEG(vac[1]/vac[0]))); 00332 00333 for ( i=1; i<p; i++ ) { 00334 __swap(vn,vo,vx); /* a(i) ---> a(i-1) */ 00335 00336 /* Sumatorio[j=1,j=i-1]( R[i-j] * a(i-1)[j] ) */ 00337 __xsum(j,0,i,vac[i-j]*vo[j],sum,0.0); 00338 00339 /* a(i)[i] = k[i] = -(R[i]+Sum[...])/E(i-1) */ 00340 vn[i] = ki = LP_POSSUM_NEG(vac[i+1]+LP_NEGSUM_NEG(sum))/E; 00341 00342 /* a(i)[j]=a(i-1)[j]+k[i]*a(i-1)[i-j] */ 00343 for ( j=0; j<i; j++ ) 00344 vn[j] = vo[j] + LP_NEGSUM_NEG(ki*vo[i-j-1]); 00345 00346 /* E(i) = (1-k[i]^2)*E(i-1) */ 00347 E *= (1.0-ki*ki); 00348 } 00349 00350 __xsum(j,0,p,vac[p-j]*vn[j],sum,0.0); 00351 ki = LP_POSSUM_NEG(vac[p+1]+LP_NEGSUM_NEG(sum))/E; 00352 00353 return E*(1.0-ki*ki); /* devuelve G*G */ 00354 } 00355 00356 /**********************************************************/ 00357 00358 PUBLIC SPL_INT XAPI tnel_lpa_cor_dur_e( SPL_INT p ) 00359 { 00360 assert(p>0); 00361 00362 return (2*(p-1)); 00363 } 00364 00365 /**********************************************************/ 00366 /* Metodo de Atal para las correlaciones. Mas lento que Durbin. 00367 A partir de las {p}+1 primeras autocorrelaciones {vac} 00368 calcula los coeficientes de reflexion {vki} para analisis 00369 de orden {p}. 00370 La funcion {devuelve} la potencia del error de prediccion */ 00371 00372 PUBLIC SPL_FLOAT XAPI lpa_cor_atl_a( SPL_pFLOAT vac, SPL_pFLOAT vai, 00373 SPL_INT p, SPL_pFLOAT vtmp ) 00374 { 00375 SPL_INT i,j; 00376 SPL_FLOAT nsum, dsum; 00377 SPL_pFLOAT vo, vn, vx; /* 2 vectores: old, new , y uno auxiliar */ 00378 assert(p>0); 00379 00380 /* inicializa vectores old(vo) y new(vn) en funcion del order p */ 00381 if (p&1) { /* si el orden es impar */ 00382 vn = vai; 00383 vo = vtmp; 00384 } /* al final, vn encaja con vai */ 00385 else { 00386 vo = vai; 00387 vn = vtmp; 00388 } 00389 00390 nsum = dsum = vac[1]*(vn[0] = LP_POSSUM_NEG(vac[1]/vac[0])); 00391 00392 for ( i=1; i<p; i++ ) { 00393 __swap(vn,vo,vx); /* a(i) ---> a(i-1) */ 00394 00395 vn[i]= LP_POSSUM_NEG(vac[i+1]+LP_NEGSUM_NEG(nsum))/ 00396 (vac[0]+LP_NEGSUM_NEG(dsum)); 00397 nsum = vac[1]*vn[i]; 00398 dsum = vac[i+1]*vn[i]; 00399 00400 for ( j=0; j<i; j++ ) { 00401 vn[j] = vo[j]+LP_NEGSUM_NEG(vn[i]*vo[i-j-1]); 00402 nsum+=vac[i-j+1]*vn[j]; 00403 dsum+=vac[j+1]*vn[j]; 00404 } 00405 } 00406 00407 return vac[0]+LP_NEGSUM_NEG(dsum); 00408 } 00409 00410 /**********************************************************/ 00411 00412 PUBLIC SPL_INT XAPI tnel_lpa_cor_atl_a( SPL_INT p ) 00413 { 00414 assert(p>0); 00415 00416 return (p-1); 00417 } 00418 00419 /**********************************************************/ 00420
1.5.6