00001 /**********************************************************/ 00002 /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/ 00003 /* 00004 Copyright: 1994 - Grupo de Voz (DAET) ETSII/IT-Bilbao 00005 00006 Nombre fuente................ SPL7.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 00012 00013 Codificacion................. Borja Etxebarria 00014 00015 Version dd/mm/aa Autor Proposito de la edicion 00016 ------- -------- -------- ----------------------- 00017 1.1.2 27/05/96 Borja correccion en documentacion 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.0 06/07/93 Borja Codificacion inicial. 00021 00022 ======================== Contenido ======================== 00023 Filtros de fase lineal, por el metodo de las ventanas (ventana 00024 de Kaiser). 00025 La respuesta impulsional que producen es simetrica respecto 00026 al origen, son _siempre_ de longitud _impar_, y las funciones solo 00027 calculan la mitad positiva mas el valor del origen. Por tanto, 00028 cuando se tiene un filtro de longitud {L} (impar), 00029 realmente se almacenan ({L}+1)/2 valores (funcion nel_fil(L)). 00030 00031 . -(L-1)/2 0 (L-1)/2 00032 . --------|-----------------|-----------------|------------>n 00033 . ^^^^^^^^^^^^^^^^^^^ 00034 . L impar (L+1)/2 almacenados 00035 00036 Los filtros se definen de la siguiente forma: 00037 00038 . A| 00039 . | 00040 . 1 ************ } +/- deltaA 00041 . | * 00042 . | * f 00043 . ------------------|-*********************|--------> 00044 . 0 fc fm/2 00045 . \--/ 00046 . deltaF 00047 00048 - {deltaA}: Amplitud del rizado: De 1+{deltaA} a 1-{deltaA}. 00049 . 0<{deltaA}<1. 00050 - {fc}: Frecuencia de corte (A(fc)=0.5) de la banda de transicion 00051 - {deltaF}: Anchura total de la banda de transicion. 00052 00053 Los valores frecuenciales ({deltaF} y {fc}) vienen normalizados 00054 (es decir, en tanto por uno) a la frecuencia de muestreo {fm}. 00055 Por tanto 0<{fc}<0.5 00056 00057 Fijados un rizado {deltaA}, y una banda de transicion {deltaF}, 00058 la funcion fil_get_L(deltaA,deltaF) {devuelve} la longitud {L} minima 00059 que debe tener el filtro para cumplir las especificaciones. Las 00060 funciones fil_get_deltaA(L,deltaF) y fil_get_deltaF(L,deltaA) son 00061 otras conversiones posibles. 00062 00063 Las funciones fil_??? con ???=(lpf, hpf, bpf, rbf) generan la 00064 respuesta impulsional del filtro. Reciben un vector {vh} de longitud 00065 nel_fil(L), el rizado {deltaA}, la frecuencia de corte {fc} (paso bajo 00066 y paso alto) o {fc1} y {fc2} (paso banda y banda eliminada), la longitud 00067 del filtro {L} (impar) y un factor de ganancia {g} (amplificacion nula para 00068 g=1.0). 00069 00070 . fil_lpf() : Filtro paso bajo 00071 . fil_hpf() : Filtro paso alto 00072 . fil_bpf() : Filtro paso banda 00073 . fil_rbf() : Filtro banda eliminada 00074 00075 El vector {vh} se rellena con la respuesta impulsional h(n): 00076 . vh[0]=h(0) vh[1]=h(1)=h(-1) vh[2]=h(2)=h(-2) ..... vh[(L-1)/2] 00077 00078 Para filtrar una muestra s[n] con el filtro {vh} el proceso es: 00079 00080 . i=(L-1)/2 00081 . y[n] = s[n]*vh[0] + SUM ( (s[n-i]+s[n+i])*vh[i] ) 00082 . i=1 00083 00084 La funcion fil_fil() hace este proceso. Recibe: 00085 . - SPL_pFLOAT {vx} : vector de {L} muestras, vx[0] a vx[L-1], del que 00086 . se va a filtrar la muestra vx[(L-1)/2]. 00087 . - SPL_pFLOAT {vh} : Respuesta impulsional de {L} puntos, que realmente 00088 . es un vector de ({L}+1)/2 valores. 00089 . - SPL_INT {L} : Longitud del filtro completo. 00090 La funcion {devuelve} el valor de vx[(L-1)/2] filtrado. 00091 00092 Ref: 00093 . Digital Processing of Speech Signals 00094 . L.R. Rabiner / R.W Schafer 00095 . Prentice Hall 00096 00097 Definir NDEBUG para desconectar la validacion de parametros 00098 con assert(). No definir este simbolo mientras se depuren 00099 aplicaciones, ya que aunque las funciones son algo mas lentas, 00100 el chequeo de validacion de parametros resuelve muchos problemas. 00101 =========================================================== 00102 */ 00103 /*/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\*/ 00104 /**********************************************************/ 00105 00106 #include "spli.h" 00107 00108 /**********************************************************/ 00109 /* a partir de la amplitud de rizado {deltaA} y de la banda 00110 de transicion {deltaF} {devuelve} la longitud L del filtro. 00111 L es siempre impar */ 00112 00113 PUBLIC SPL_INT XAPI fil_get_L( SPL_FLOAT deltaA, SPL_FLOAT deltaF ) 00114 { 00115 SPL_INT L; 00116 assert(deltaA>0.0); 00117 assert(deltaF>0.0); 00118 00119 L=(SPL_INT)(ceil((-a2db(deltaA)-7.95)/(14.36*deltaF)))+1; 00120 if (L&1) /* L debe ser impar */ 00121 return L; 00122 else 00123 return L+1; 00124 } 00125 00126 /**********************************************************/ 00127 /* a partir de la longitud del filtro {L} y de la banda de 00128 transicion {deltaF} {devuelve} la amplitud del rizado deltaA */ 00129 00130 PUBLIC SPL_FLOAT XAPI fil_get_deltaA( SPL_INT L, SPL_FLOAT deltaF ) 00131 { 00132 assert(L>0); 00133 assert(L&1); /* L impar */ 00134 assert(deltaF>0.0); 00135 00136 return db2a(-((L-1)*(14.36*deltaF)+7.95)); 00137 } 00138 00139 /**********************************************************/ 00140 /* a partir de la longitud del filtro {L} y de la amplitud del 00141 rizado {deltaA} {devuelve} la banda de transicion deltaF */ 00142 00143 PUBLIC SPL_FLOAT XAPI fil_get_deltaF( SPL_INT L, SPL_FLOAT deltaA ) 00144 { 00145 assert(L>0); 00146 assert(L&1); /* L impar */ 00147 assert(deltaA>0.0); 00148 00149 return (-a2db(deltaA)-7.95)/((L-1)*14.36); 00150 } 00151 00152 /**********************************************************/ 00153 /* Funcion de uso interno. 00154 a partir de la amplitud de rizado {deltaA} {devuelve} la constante 00155 alpha necesaria para el prototipo del filtro */ 00156 00157 PRIVATE SPL_FLOAT fil_get_alpha( SPL_FLOAT deltaA ) 00158 { 00159 assert(deltaA>0.0); 00160 00161 deltaA = -a2db(deltaA); 00162 if (deltaA>50.0) 00163 return 0.1102*(deltaA-8.7); 00164 else if (deltaA>21.0) 00165 return 0.5842*pow(deltaA-21.0,0.4)+0.07886*(deltaA-21); 00166 else 00167 return 0; 00168 } 00169 00170 /**********************************************************/ 00171 /* Rellena el vector {vh} de ({L}+1)/2 puntos con la mitad positiva (de 00172 0 a ({L}-1)/2) de la respuesta impulsional de un filtro PASO BAJO, 00173 con amplitud de rizado {deltaA}, frecuencia de corte {fc}, y longitud 00174 del filtro de {L} puntos, y con una ganancia {g} (factor multiplicativo 00175 de escalado). 00176 {devuelve} el puntero {vh}. 00177 {vh} tiene ({L}+1)/2 puntos. 00178 0<{deltaA}<1, 0<{fc}<0.5, {L}>0, {L} impar. */ 00179 00180 PUBLIC SPL_pFLOAT XAPI fil_lpf( SPL_pFLOAT vh, SPL_FLOAT deltaA, 00181 SPL_FLOAT fc, SPL_INT L, SPL_FLOAT g ) 00182 { 00183 SPL_INT i, n; 00184 SPL_FLOAT nd, alpha, I0alpha, tmp, wp; 00185 assert(deltaA>0.0); 00186 assert(fc>0.0); 00187 assert(fc<=0.5); 00188 assert(L>0); 00189 assert(L&1); /* L debe ser impar */ 00190 00191 n = (L-1)/2; 00192 nd = n; 00193 alpha = fil_get_alpha(deltaA); 00194 I0alpha = bessel_I0(alpha); 00195 wp = (M_PI*2.0)*fc; 00196 00197 vh[0]= g*2.0*fc; 00198 for (i=1; i<=n; i++) { /* sinc enventanada con Kaiser */ 00199 tmp = i/nd; 00200 vh[i] =g*sin(wp*i)/(M_PI*i)* 00201 bessel_I0(alpha*sqrt(1.0-tmp*tmp))/I0alpha; 00202 } 00203 return vh; 00204 } 00205 00206 /**********************************************************/ 00207 /* Rellena el vector {vh} de ({L}+1)/2 puntos con la mitad positiva (de 00208 0 a ({L}-1)/2) de la respuesta impulsional de un filtro PASO ALTO, 00209 con amplitud de rizado {deltaA}, frecuencia de corte {fc}, y longitud 00210 del filtro de {L} puntos, y con una ganancia {g} (factor multiplicativo 00211 de escalado). 00212 {devuelve} el puntero {vh}. 00213 {vh} tiene ({L}+1)/2 puntos. 00214 0<{deltaA}<1, 0<{fc}<0.5, {L}>0, {L} impar. */ 00215 00216 PUBLIC SPL_pFLOAT XAPI fil_hpf( SPL_pFLOAT vh, SPL_FLOAT deltaA, 00217 SPL_FLOAT fc, SPL_INT L, SPL_FLOAT g ) 00218 { 00219 SPL_INT i, n; 00220 00221 fil_lpf(vh,deltaA,0.5-fc,L,g); 00222 n=(L+1)/2; 00223 00224 for (i=1; i<n; i+=2) 00225 vh[i] = (-vh[i]); /* multiplicar por cos(fm/2*t) */ 00226 00227 return vh; 00228 } 00229 00230 /**********************************************************/ 00231 /* Rellena el vector {vh} de ({L}+1)/2 puntos con la mitad positiva (de 00232 0 a ({L}-1)/2) de la respuesta impulsional de un filtro PASO BANDA, 00233 con amplitud de rizado {deltaA}, frecuencia de corte inferior {fc1} y 00234 superior {fc2}, y longitud del filtro de {L} puntos, y con una 00235 ganancia {g} (factor multiplicativo de escalado). 00236 {devuelve} el puntero {vh}. 00237 {vh} tiene ({L}+1)/2 puntos. 00238 0<{deltaA}<1, 0<{fc1}<0.5, 0<{fc1}<0.5, {fc1}<{fc2}, {L}>0, {L} impar. */ 00239 00240 PUBLIC SPL_pFLOAT XAPI fil_bpf( SPL_pFLOAT vh, SPL_FLOAT deltaA, 00241 SPL_FLOAT fc1, SPL_FLOAT fc2, 00242 SPL_INT L, SPL_FLOAT g ) 00243 { 00244 SPL_INT i, n; 00245 SPL_FLOAT w; 00246 assert(fc2>fc1); 00247 00248 fil_lpf(vh,deltaA,(fc2-fc1)/2.0,L,g); 00249 n=(L+1)/2; 00250 w = (fc2+fc1)*M_PI; 00251 00252 for (i=0; i<n; i++) /* modula cos(w*t) */ 00253 vh[i] *= 2*cos(w*i); 00254 00255 return vh; 00256 } 00257 00258 /**********************************************************/ 00259 /* Rellena el vector {vh} de ({L}+1)/2 puntos con la mitad positiva (de 00260 0 a ({L}-1)/2) de la respuesta impulsional de un filtro BANDA ELIMINADA, 00261 con amplitud de rizado {deltaA}, frecuencia de corte inferior {fc1} y 00262 superior {fc2}, y longitud del filtro de {L} puntos, y con una 00263 ganancia {g} (factor multiplicativo de escalado). 00264 {devuelve} el puntero {vh}. 00265 {vh} tiene ({L}+1)/2 puntos. 00266 0<{deltaA}<1, 0<{fc1}<0.5, 0<{fc1}<0.5, {fc1}<{fc2}, {L}>0, {L} impar. */ 00267 00268 PUBLIC SPL_pFLOAT XAPI fil_rbf( SPL_pFLOAT vh, SPL_FLOAT deltaA, 00269 SPL_FLOAT fc1, SPL_FLOAT fc2, 00270 SPL_INT L, SPL_FLOAT g ) 00271 { 00272 SPL_INT i, n; 00273 SPL_FLOAT w; 00274 assert(fc2>fc1); 00275 00276 fil_lpf(vh,deltaA,(fc2-fc1)/2.0,L,g); 00277 n = (L+1)/2; 00278 w = (fc2+fc1)*M_PI; 00279 00280 for (i=0; i<n; i++) /* modula cos(w*t) */ 00281 vh[i] *= -2*cos(w*i); 00282 vh[0] += g; 00283 00284 return vh; 00285 } 00286 00287 /**********************************************************/ 00288 /* {devuelve} la longitud en puntos que debe tener un vector 00289 para almacenar una respuesta impulsional de longitud {L}, simetrica 00290 respecto al origen */ 00291 00292 PUBLIC SPL_INT XAPI nel_fil( SPL_INT L ) 00293 { 00294 assert(L>0); 00295 assert(L&1); /* L debe ser impar */ 00296 00297 return (L+1)/2; 00298 } 00299 00300 /**********************************************************/ 00301 /* {devuelve} el valor filtrado para la muestra vx[(L-1)/2], 00302 es decir, que por delante de la muestra a filtrar debe haber 00303 (nel_fil(L)-1) muestras. 00304 - {vx} es un vector de longitud L. 00305 - {vh} es la respuesta impulsional del filtro, 00306 . de longitud ({L}+1)/2 (nel_fil(L)). 00307 - {L} es la longitud total del filtro. {L}>0 impar. */ 00308 00309 PUBLIC SPL_FLOAT XAPI fil_fil( SPL_pFLOAT vx, SPL_pFLOAT vh, SPL_INT L ) 00310 { 00311 SPL_INT i,n; 00312 SPL_FLOAT sum; 00313 assert(L>0); 00314 assert(L&1); /* L debe ser impar */ 00315 00316 n=(L-1)/2; 00317 00318 __sum(i,1,n,(vx[n-i]+vx[n+i])*vh[i],sum,vh[0]*vx[n]); 00319 00320 return sum; 00321 } 00322 00323 /**********************************************************/ 00324
1.5.6