00001 /**********************************************************/ 00002 /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/ 00003 /* 00004 Copyright: 1994 - Grupo de Voz (DAET) ETSII/IT-Bilbao 00005 00006 Nombre fuente................ SPL5e.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.2.0 24/06/97 Borja incorporar lp_k2ta() 00018 1.1.1 30/07/95 Borja scope funciones explicito 00019 1.1.0 08/12/94 Borja revision general (tipos,idx,nel,tnel...) 00020 1.0.1 07/12/93 Borja soporte LP_NEGSUM_NEG() y LP_POSSUM_NEG() 00021 1.0.0 06/07/93 Borja Codificacion inicial. 00022 00023 ======================== Contenido ======================== 00024 Algoritmos de prediccion lineal (LP). Ver SPL5A.C para mas detalles. 00025 00026 Definir NDEBUG para desconectar la validacion de parametros 00027 con assert(). No definir este simbolo mientras se depuren 00028 aplicaciones, ya que aunque las funciones son algo mas lentas, 00029 el chequeo de validacion de parametros resuelve muchos problemas. 00030 =========================================================== 00031 */ 00032 /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/ 00033 /**********************************************************/ 00034 00035 #include "spli.h" 00036 00037 /**********************************************************/ 00038 /* convierte los coeficientes de reflexion {vki} en coeficientes 00039 LPC {vai}. El orden del analisis es {p}. Los vectores {vki} y {vai} tienen 00040 {p} elementos cada uno. El usuario debe crear, pasar y luego destruir 00041 un vector temporal de tnel_lp_k2a(p) elementos en {vtmp}. 00042 La funcion devuelve el propio vector {vai}. 00043 {p}>0 */ 00044 00045 PUBLIC SPL_pFLOAT XAPI lp_k2a( SPL_pFLOAT vki, SPL_pFLOAT vai, 00046 SPL_INT p, SPL_pFLOAT vtmp ) 00047 { 00048 SPL_INT i,j; 00049 SPL_FLOAT ki; 00050 SPL_pFLOAT vo, vn, vx; /* 2 vectores: old, new , y uno auxiliar */ 00051 assert(p>0); 00052 00053 /* inicializa vectores old(vo) y new(vn) en funcion del order p */ 00054 if (p&1) { /* si el orden es impar */ 00055 vn = vai; 00056 vo = vtmp; 00057 } /* al final, vn encaja con vai */ 00058 else { 00059 vo = vai; 00060 vn = vtmp; 00061 } 00062 00063 vn[0] = vki[0]; 00064 for ( i=1; i<p; i++ ) { 00065 __swap(vn,vo,vx); /* a(i) ---> a(i-1) */ 00066 00067 vn[i] = ki = vki[i]; 00068 /* a(i)[j]=a(i-1)[j]+k[i]*a(i-1)[i-j] */ 00069 for ( j=0; j<i; j++ ) 00070 vn[j] = vo[j] + LP_NEGSUM_NEG(ki*vo[i-j-1]); 00071 } 00072 return vai; 00073 } 00074 00075 /**********************************************************/ 00076 00077 PUBLIC SPL_INT XAPI tnel_lp_k2a( SPL_INT p ) 00078 { 00079 assert(p>0); 00080 00081 return (p-1); 00082 } 00083 00084 /**********************************************************/ 00085 /* convierte los coeficientes LPC {vai} en coeficientes de 00086 reflexion {vki}. El orden del analisis es {p}. Los vectores {vki} y {vai} 00087 tienen {p} elementos cada uno. El usuario debe crear, pasar y luego 00088 destruir un vector temporal de tnel_lp_a2k(p) elementos en {vtmp}. 00089 La funcion devuelve el propio vector {vki}. 00090 {p}>0 */ 00091 00092 PUBLIC SPL_pFLOAT XAPI lp_a2k( SPL_pFLOAT vai, SPL_pFLOAT vki, 00093 SPL_INT p, SPL_pFLOAT vtmp ) 00094 { 00095 SPL_INT i,j; 00096 SPL_FLOAT f, ai; 00097 #define vn vtmp 00098 SPL_pFLOAT vo, vx; /* 2 vectores: old, new , y uno auxiliar */ 00099 assert(p>0); 00100 00101 i=p-1; 00102 vki[i] = ai = vai[i]; 00103 f = 1-ai*ai; 00104 i--; 00105 for (j=0; j<=i; j++) 00106 vki[j] = (vai[j]+LP_POSSUM_NEG(ai*vai[i-j]))/f; 00107 00108 /* vn=vtmp en el define */ 00109 vo=vki; 00110 00111 for ( ;i>0; ) { 00112 ai=vo[i]; 00113 f = 1-ai*ai; 00114 i--; 00115 for (j=0; j<=i; j++) 00116 vn[j] = (vo[j]+LP_POSSUM_NEG(ai*vo[i-j]))/f; 00117 00118 vki[i]=vn[i]; 00119 __swap(vn,vo,vx); /* a(i) ---> a(i-1) */ 00120 } 00121 return vki; 00122 #undef vn 00123 } 00124 00125 /**********************************************************/ 00126 00127 PUBLIC SPL_INT XAPI tnel_lp_a2k( SPL_INT p ) 00128 { 00129 assert(p>0); 00130 00131 return (p-1); 00132 } 00133 00134 /**********************************************************/ 00135 /* convierte los coeficientes de reflexion {vki} en areas del 00136 tubo sin perdidas (vti). El orden del analisis es {p}. Los 00137 vectores {vki} y {vti} tienen {p} elementos cada uno. 00138 La funcion devuelve el propio vector {vti}. 00139 {p}>0. {vki}!={vti}. 00140 Notese que son areas, no area-ratios, y se han normalizado 00141 para area=1 (area del primer "tubito")!! 00142 No se puede decir que las areas obtenidas encajen con las 00143 areas reales del tracto vocal, sin embargo, aplicando 00144 preenfasis para eliminar el efecto del pulso glotal y 00145 la radiacion, mas o menos encajara. 00146 Se puede pasar del area al radio de un tubo circular haciendo 00147 radio=sqrt(area/PI). */ 00148 00149 PUBLIC SPL_pFLOAT XAPI lp_k2ta( SPL_pFLOAT vki, SPL_pFLOAT vti, SPL_INT p ) 00150 { 00151 SPL_INT i; 00152 SPL_FLOAT k; 00153 00154 assert(p>0); 00155 p--; 00156 vti[0]=1; 00157 for (i=0; i<p; i++) { 00158 k = LP_POSSUM_NEG(vki[i]); 00159 if (k==-1) 00160 vti[i+1]=vti[i]; 00161 else 00162 vti[i+1] = vti[i]*(1-k)/(1+k); 00163 } 00164 return vti; 00165 } 00166 00167 /**********************************************************/ 00168 /* a partir de los coeficientes LPC de orden {p} en {vai} ({p} elementos) 00169 calcula las {p} primeras autocorrelaciones R(1) a R(p) de la respuesta 00170 impulsional del filtro reconstructor H(z)=1/(1+-sum(ak*z^-k)), 00171 normalizadas a la autocorrelacion en el origen en {vrh} ({p} elementos). 00172 La funcion devuelve el valor de la autocorrelacion en el origen 00173 (valor de normalizacion). 00174 Las autocorrelaciones de la respuesta impulsional coinciden con 00175 las del segmento de analisis utilizadas para el calculo de los 00176 LPC por cualquiera de los metodos de autocorrelacion (Durbin...) 00177 normalizadas a la del origen. Por tanto, a partir de estas 00178 autocorrelaciones (con un 1 por delante) se pueden recuperar los 00179 coeficientes LPC ( es decir, la rutina lp_rh2a() realmente es 00180 cualqiera de los metodos para la autocorrelacion: lpa_cor_dur_a()...) 00181 El vector temporal {vtmp} tiene tnel_lp_a2rh(p) elementos. 00182 {p}>0 */ 00183 00184 PUBLIC SPL_FLOAT XAPI lp_a2rh( SPL_pFLOAT vai, SPL_pFLOAT vrh, 00185 SPL_INT p, SPL_pFLOAT vtmp ) 00186 { 00187 #define vn vtmp 00188 SPL_INT i,j; 00189 SPL_FLOAT ai,ai12,sum; 00190 SPL_pFLOAT vo; 00191 assert(p>0); 00192 00193 vo=vai; 00194 /*vn=vtmp en el define */ 00195 for (i=(--p); i>=1; ) { 00196 ai=vo[i--]; 00197 ai12=1-ai*ai; 00198 for (j=0; j<=i; j++) 00199 vn[j] = (vo[j]+LP_POSSUM_NEG(ai*vo[i-j]))/ai12; 00200 vo=vn; 00201 vn += (i+1); 00202 } 00203 00204 vn--; 00205 for (i=0; i<p; ) { 00206 #ifdef LP_NEGSUM 00207 __xsum(j,0,i,vn[j]*vrh[i-j-1],sum,vn[i]); 00208 #else 00209 __xsub(j,0,i,vn[j]*vrh[i-j-1],sum,-vn[i]); 00210 #endif 00211 vrh[i]=sum; 00212 i++; 00213 vn -= (i+1); 00214 } 00215 #ifdef LP_NEGSUM 00216 __xsum(j,0,p,vai[j]*vrh[p-j-1],sum,vai[p]); 00217 #else 00218 __xsub(j,0,p,vai[j]*vrh[p-j-1],sum,-vai[p]); 00219 #endif 00220 vrh[p]=sum; 00221 00222 #ifdef LP_NEGSUM 00223 __sub(j,0,p,vai[j]*vrh[j],sum,1.0); 00224 #else 00225 __sum(j,0,p,vai[j]*vrh[j],sum,1.0); 00226 #endif 00227 return 1.0/sum; 00228 #undef vn 00229 } 00230 00231 /**********************************************************/ 00232 00233 PUBLIC SPL_INT XAPI tnel_lp_a2rh( SPL_INT p ) 00234 { 00235 assert(p>0); 00236 00237 return ((p-1)*p)/2; 00238 } 00239 00240 /**********************************************************/ 00241 /* a partir de los coeficientes lpc {vai} de orden {p} calcula las 00242 correlaciones de la respuesta impulsional del filtro inverso A(z). 00243 A(z) es un filtro FIR de {p}+1 elementos (la respuesta impulsional 00244 consiste en un 1 seguido de los coeficientes lpc {vai}) por lo que 00245 solo hay {p}+1 autocorrelaciones distintas de 0. 00246 El vector {vai} tiene {p} elementos (coeficientes LPC de a(1) a a(p). 00247 El vector {vra} que se rellena, tiene {p}+1 elementos (correlaciones 00248 desde 0 hasta {p}) 00249 La funcion {devuelve} el vector {vra}. 00250 */ 00251 00252 PUBLIC SPL_pFLOAT XAPI lp_a2ra( SPL_pFLOAT vai, 00253 SPL_pFLOAT vra, SPL_INT p ) 00254 { 00255 SPL_INT i,k; 00256 SPL_FLOAT sum; 00257 assert(p>=0); 00258 00259 __xsum(k,0,p,vai[k]*vai[k],sum,1.0); 00260 vra[0]=sum; 00261 00262 for (i=1; i<=p; i++) { 00263 #ifdef __LP_NEGSUM 00264 __xsum(k,i,p,vai[k]*vai[k-i],sum,-vai[i-1]); 00265 #else 00266 __xsum(k,i,p,vai[k]*vai[k-i],sum,vai[i-1]); 00267 #endif 00268 vra[i]=sum; 00269 } 00270 00271 return vra; 00272 } 00273 00274 /**********************************************************/ 00275 /* calcula los cepstrum LPC de H(z). {vai} el el vector de 00276 coeficientes lpc de orden {p} ({p} elementos) y en {vci} se devuelven 00277 los cepstrum c(1) a c(p) ({p} elementos). Para obtener c(0), llamar 00278 a la funcion lp_c0(e) siendo {e} la potencia del error de prediccion. 00279 La funcion {devuelve} el vector {vci} */ 00280 00281 PUBLIC SPL_pFLOAT XAPI lp_a2c( SPL_pFLOAT vai, SPL_pFLOAT vci, SPL_INT p ) 00282 { 00283 SPL_INT i,k; 00284 SPL_FLOAT sum; 00285 assert(p>=0); 00286 00287 for (i=0; i<p; i++) { 00288 #ifdef LP_NEGSUM 00289 __sum(k,1,i,k/(i+1)*vci[k-1]*vai[i-k],sum,vai[i]); 00290 #else 00291 __sub(k,1,i,k/(i+1)*vci[k-1]*vai[i-k],sum,-vai[i]); 00292 #endif 00293 vci[i]=sum; 00294 } 00295 00296 return vci; 00297 } 00298 00299 /**********************************************************/ 00300 /* {devuelve} el cepstrum c(0). {e} es la potencia del error de 00301 prediccion. 00302 {e}>0 */ 00303 00304 PUBLIC SPL_FLOAT XAPI lp_c0( SPL_FLOAT e ) 00305 { 00306 assert(e>0.0); 00307 00308 return 0.5*log(e); 00309 } 00310 00311 /**********************************************************/ 00312 /* genera el filtro de anchos de banda expandidos {vbwea} a partir 00313 del filtro lpc {vai}, de orden {p}. 00314 00315 . P(z) ---------> P(z/gamma) 00316 . a[k] ---------> a[k]*gamma^k 00317 00318 {vai} y {vbwea} tienen {p} elementos. 00319 {gamma} es el factor de expandido. 1 no hace nada, 0 tiene efecto 00320 maximo. 00321 Los vectores {vai} y {vbwea} SI pueden ser el mismo. 00322 0<={gamma}<=1. 00323 La funcion {devuelve} el vector {vbwea} */ 00324 00325 PUBLIC SPL_pFLOAT XAPI lp_a2bwea( SPL_pFLOAT vai, SPL_pFLOAT vbwea, 00326 SPL_INT p, SPL_FLOAT gamma ) 00327 { 00328 SPL_INT i; 00329 SPL_FLOAT agm; 00330 assert(gamma>=0.0); 00331 assert(gamma<=1.0); 00332 assert(p>=0); 00333 00334 agm=1.0; 00335 for (i=0; i<p; i++) { 00336 agm*=gamma; 00337 vbwea[i]=vai[i]*agm; 00338 } 00339 return vbwea; 00340 } 00341 00342 /**********************************************************/ 00343 /* Predictor P(z). 00344 efectua la prediccion de orden {p} para la muestra actual en 00345 funcion de {p} muestras anteriores(contiguas o no). En {vs} se 00346 envian las {p} muestras utilizadas (p elementos): 00347 . vs[0]=s[n-p] hasta vs[p-1]=s[n-1] (corto termino) 00348 . vs[0]=s[n-p-M+1] hasta vs[p-1]=s[n-M] (largo termino: retardo M) 00349 y predice s[n] (valor que {devuelve} la funcion) utilizando los 00350 coeficientes lpc {vai} ({p} elementos). */ 00351 00352 PUBLIC SPL_FLOAT XAPI lps_P_a( SPL_pFLOAT vs, SPL_pFLOAT vai, SPL_INT p ) 00353 { 00354 SPL_INT i; 00355 SPL_FLOAT sum; 00356 assert(p>=0); 00357 00358 __xsum(i,0,p,vai[i]*vs[p-i-1],sum,0.0); 00359 00360 return LP_POSSUM_NEG(sum); 00361 } 00362 00363 /**********************************************************/ 00364 /* Filtro inverso A(z). 00365 efectua la prediccion de orden {p} para la muestra actual en 00366 funcion de {p} muestras anteriores(contiguas o no). En {vs} se envian las 00367 {p} muestras utilizadas({p} elementos): 00368 . vs[0]=s[n-p] hasta vs[p-1]=s[n-1] (corto termino) 00369 . vs[0]=s[n-p-M+1] hasta vs[p-1]=s[n-M] (largo termino: retardo M) 00370 En {sn} se envia la muestra en el instante actual s[n]. 00371 Se predice la muestra actual s[n] utilizando los coeficientes lpc {vai} 00372 ({p} elementos). 00373 La funcion {devuelve} el error de prediccion e[n] del filtro, restando 00374 al valor de la muestra actual el valor de prediccion. */ 00375 00376 PUBLIC SPL_FLOAT XAPI lps_A_a( SPL_pFLOAT vs, SPL_pFLOAT vai, 00377 SPL_INT p, SPL_FLOAT sn ) 00378 { 00379 SPL_INT i; 00380 SPL_FLOAT sum; 00381 assert(p>=0); 00382 00383 __xsum(i,0,p,vai[i]*vs[p-i-1],sum,0.0); 00384 00385 return sn+LP_NEGSUM_NEG(sum); 00386 } 00387 00388 /**********************************************************/ 00389 /* Filtro reconstructor H(z). 00390 efectua la prediccion de orden {p} para la muestra actual en 00391 funcion de {p} muestras anteriores(contiguas o no). En {vs} se envian las 00392 {p} muestras utilizadas({p} elementos): 00393 . vs[0]=s[n-p] hasta vs[p-1]=s[n-1] (corto termino) 00394 . vs[0]=s[n-p-M+1] hasta vs[p-1]=s[n-M] (largo termino: retardo M) 00395 En {en} se envia el error en el instante actual e[n]. 00396 Se predice la muestra actual s[n] utilizando los coeficientes lpc {vai} 00397 ({p} elementos). 00398 La funcion {devuelve} el valor de prediccion corregido (sumado) con 00399 el error e[n]. */ 00400 00401 PUBLIC SPL_FLOAT XAPI lps_H_a( SPL_pFLOAT vs, SPL_pFLOAT vai, 00402 SPL_INT p, SPL_FLOAT en ) 00403 { 00404 SPL_INT i; 00405 SPL_FLOAT sum; 00406 assert(p>=0); 00407 00408 __xsum(i,0,p,vai[i]*vs[p-i-1],sum,0.0); 00409 00410 return en+LP_POSSUM_NEG(sum); 00411 } 00412 00413 /**********************************************************/ 00414 /* a partir de los coeficientes LPC {vai} de orden {p}, devuelve 00415 en {van} las N primeras muestras de la respuesta impulsional a(n) 00416 del filtro A(z), desde a(0) hasta a(N-1). 00417 La respuesta impulsional es finita, con a(0)=1, a(1)= +-a1 00418 hasta a(p)= +-ap. 00419 +- depende de si LP_NEGSUM esta o no definido. Si esta 00420 definido, es -, si no esta definido, es +. 00421 A partir de a(p+1), la respuesta impulsional es nula. 00422 00423 La funcion {devuelve} {van}. */ 00424 00425 PUBLIC SPL_pFLOAT XAPI lps_A_an_a( SPL_pFLOAT vai, SPL_INT p, 00426 SPL_pFLOAT van, SPL_INT N ) 00427 { 00428 SPL_INT i; 00429 assert(p>=0); 00430 assert(N>=0); 00431 00432 if (N>0) 00433 van[0]=1.0; 00434 00435 if (p>=N) 00436 p=N-1; 00437 00438 for (i=1; i<=p; i++) 00439 van[i]=LP_NEGSUM_NEG(vai[i-1]); 00440 00441 for (i=p+1; i<N; i++) 00442 van[i]=0.0; 00443 00444 return van; 00445 } 00446 00447 /**********************************************************/ 00448 /* a partir de los coeficientes LPC {vai} de orden {p}, devuelve 00449 en {vhn} las N primeras muestras de la respuesta impulsional h(n) 00450 del filtro H(z), desde h(0) hasta h(N-1). 00451 La respuesta impulsional es infinita, con h(0)=1. 00452 00453 La funcion {devuelve} {vhn}. */ 00454 00455 PUBLIC SPL_pFLOAT XAPI lps_H_hn_a( SPL_pFLOAT vai, SPL_INT p, 00456 SPL_pFLOAT vhn, SPL_INT N ) 00457 { 00458 SPL_INT i; 00459 assert(p>=0); 00460 assert(N>=0); 00461 00462 if (N>0) 00463 vhn[0]=1.0; 00464 00465 if (N>p) { 00466 for (i=1; i<=p; i++) 00467 vhn[i]=lps_H_a(vhn+(i-1),vai,i,0.0); 00468 } 00469 00470 for (i=p+1; i<N; i++) 00471 vhn[i]=lps_H_a(vhn+(i-1),vai,p,0.0); 00472 00473 return vhn; 00474 } 00475 00476 /**********************************************************/ 00477 /* a partir de los coeficientes LPC {vai} de orden {p}, devuelve 00478 en {vpn} las {N} primeras muestras de la respuesta impulsional p(n) 00479 del filtro P(z), desde p(0) hasta p(N-1). 00480 La respuesta impulsional es finita, con p(0)=0, p(1)= -+a1 00481 hasta p(p)= -+ap. 00482 -+ depende de si LP_NEGSUM esta o no definido. Si esta 00483 definido, es +, si no esta definido, es -. 00484 A partir de p(p+1), la respuesta impulsional es nula. 00485 00486 La funcion {devuelve} {vpn}. */ 00487 00488 PUBLIC SPL_pFLOAT XAPI lps_P_pn_a( SPL_pFLOAT vai, SPL_INT p, 00489 SPL_pFLOAT vpn, SPL_INT N ) 00490 { 00491 SPL_INT i; 00492 assert(p>=0); 00493 assert(N>=0); 00494 00495 if (N>0) 00496 vpn[0]=0.0; 00497 00498 if (p>=N) 00499 p=N-1; 00500 00501 for (i=1; i<=p; i++) 00502 vpn[i] = LP_POSSUM_NEG(vai[i-1]); 00503 00504 for (i=p+1; i<N; i++) 00505 vpn[i]=0.0; 00506 00507 return vpn; 00508 } 00509 00510 /**********************************************************/ 00511
1.5.6