AhoLib: Fichero Fuente xfft4.c

00001 
00002 
00028 #include "_xfft.h"
00029 
00030 
00031 
00032 
00033 SPL_FLOAT Hz2Log ( SPL_FLOAT f )
00034 { return log(f); }
00035 
00036 SPL_FLOAT Log2Hz ( SPL_FLOAT f )
00037 { return exp(f); }
00038 
00039 SPL_FLOAT Hz2Log10 ( SPL_FLOAT f )
00040 { return log10(f); }
00041 
00042 SPL_FLOAT Log102Hz ( SPL_FLOAT f )
00043 { return pow(10, f); }
00044 
00045 SPL_FLOAT Hz2Mel ( SPL_FLOAT f )
00046 { return 1127.0*log(1.0+f/700.0); }
00047 
00048 SPL_FLOAT Mel2Hz ( SPL_FLOAT f )
00049 { return 700.0*(exp(f/1127.0)-1.0); }
00050 
00051 SPL_FLOAT Hz2Bark ( SPL_FLOAT f )
00052 {
00053     SPL_FLOAT b = (26.81*f) / (1960+f) - 0.53;
00054     if (b < 2) b = 0.85*b + 0.3;
00055     else if (b > 20.1) b = 1.22*b - 4.422;
00056     return b;
00057 }
00058 
00059 SPL_FLOAT Bark2Hz ( SPL_FLOAT f )
00060 {
00061     if (f < 2) f = (f-0.3)/0.85;
00062     else if (f > 20.1) f = (f+4.422)/1.22;
00063     return 1960*(f+0.53)/(26.28-f);
00064 }
00065     
00066 SPL_FLOAT Hz2Bin ( SPL_FLOAT f, SPL_FLOAT fres )
00067 { return floor(f/fres); }
00068 
00069 SPL_FLOAT Bin2Hz ( SPL_FLOAT f, SPL_FLOAT fres )
00070 { return f*fres; }
00071 
00085 SPL_FLOAT FbandScaleConvert ( SPL_FLOAT f, XFBAND_SCALE from, XFBAND_SCALE to, SPL_FLOAT Fres )
00086 {
00087     /*Primero convertimos lo que sea a lineal*/
00088     switch (from)
00089     {
00090         case XFBAND_SCALE_MEL:
00091             f = Mel2Hz(f);
00092             break;
00093         case XFBAND_SCALE_BARK:
00094             f = Bark2Hz(f);
00095             break;
00096         case XFBAND_SCALE_BIN:
00097             f = Bin2Hz(f, Fres);
00098             break;
00099         case XFBAND_SCALE_LOG:
00100             f = Log2Hz(f);
00101             break;
00102         case XFBAND_SCALE_LOG10:
00103             f = Log102Hz(f);
00104             break;
00105         default:
00106             break;
00107     }
00108     
00109     /*Ahora convertimos de lineal a la salida*/
00110     switch (to)
00111     {
00112         case XFBAND_SCALE_MEL:
00113             f = Hz2Mel(f);
00114             break;
00115         case XFBAND_SCALE_BARK:
00116             f = Hz2Bark(f);
00117             break;
00118         case XFBAND_SCALE_BIN:
00119             f = Hz2Bin(f, Fres);
00120             break;
00121         case XFBAND_SCALE_LOG:
00122             f = Hz2Log(f);
00123             break;
00124         case XFBAND_SCALE_LOG10:
00125             f = Hz2Log10(f);
00126             break;
00127         default:
00128             break;
00129     }
00130     
00131     return f;
00132 
00133 }
00134 
00135 
00139 /* Permite recalcular los índices de banda y los pesos para cada bin.
00140    Devuelve SPL_TRUE si todo va bien y SPL_FALSE si ha habido errores.
00141    En caso de error, libera toda la memoria que se haya reservado ya.
00142 */
00143 PRIVATE SPL_BOOL fband_calcfilters( pFBAND r )
00144 {
00145     assert(r!=NULL);
00146     assert(r->_filt == XFBAND_FILTER_RECT || r->_filt == XFBAND_FILTER_TRIANG);
00147     
00148     if (r->_fft == NULL)
00149         return SPL_FALSE;
00150 
00151     SPL_INT fftnsamples = rfft_getfftnp(r->_fft);
00152     
00153     /*Primero liberamos la memoria*/
00154     if (r->_findex) xfft_free(r->_findex);
00155     if (r->_fweight) xfft_free(r->_fweight);
00156     if (r->_pvec) xfft_free(r->_pvec);      //Podría sustituirse en una función externa setnband
00157     
00158     r->_findex = (SPL_pINT)xfft_malloc(sizeof(SPL_INT)*fftnsamples);
00159     r->_fweight = (SPL_pFLOAT)xfft_malloc(sizeof(SPL_FLOAT)*fftnsamples);
00160     r->_pvec = (SPL_pFLOAT)xfft_malloc(sizeof(SPL_FLOAT)*r->_nf);
00161     
00162     if (r->_findex == NULL || r->_fweight == NULL || r->_pvec == NULL) {
00163         if (r->_findex) xfft_free(r->_findex);
00164         if (r->_fweight) xfft_free(r->_fweight);
00165         if (r->_pvec) xfft_free(r->_pvec);
00166         return SPL_FALSE;
00167     }
00168     
00169     /* Necesitamos minfreq y maxfreq convertidos a la escala solicitada
00170        para crear los filtros */
00171     SPL_FLOAT fres = r->_srate/rfft_getfftnp(r->_fft);
00172     SPL_FLOAT outFreqMin = FbandScaleConvert(r->_minfreqidx, XFBAND_SCALE_BIN, r->_scale, fres);
00173     SPL_FLOAT outFreqMax = FbandScaleConvert(r->_maxfreqidx, XFBAND_SCALE_BIN, r->_scale, fres);
00174     
00175     switch (r->_filt)
00176     {
00177         case XFBAND_FILTER_RECT:
00178         {
00179             SPL_FLOAT filtSize = (outFreqMax-outFreqMin)/r->_nf;    /*Anchura de cada filtro*/
00180             for (SPL_INT i = r->_minfreqidx; i <= r->_maxfreqidx; ++i)
00181             {
00182                 /*Calculamos la frecuencia escalada a la que corresponde ese bin*/
00183                 SPL_FLOAT outBin = FbandScaleConvert(i, XFBAND_SCALE_BIN, r->_scale, fres);
00184         
00185                 /*El peso y el filtro al que corresponde*/
00186                 SPL_FLOAT tmp = (outBin-outFreqMin)/filtSize;
00187                 SPL_FLOAT tmp2;
00188                 modf (tmp, &tmp2);
00189                 r->_fweight[i] = 1.0;   /* Peso unidad por ser filtros rectangulares */
00190                 r->_findex[i] = (SPL_INT)tmp2;
00191         
00192                 /*Para Fs nos queda que le corresponde el filtro r->_nf, el cual no existe.
00193                     El problema es que está en la frontera entre dos filtros. Podemos ponerselo al anterior*/
00194                 if (r->_findex[i] == r->_nf)
00195                     r->_findex[i] = r->_nf-1;
00196         
00197                 /*Escalamos los pesos por el número de bins para terminar teniendo potencia independiente
00198                     de la NFFT (que cambiará con la Fs)*/
00199                 r->_fweight[i] /= fftnsamples;
00200             }
00201             break;
00202         }
00203 
00204         case XFBAND_FILTER_TRIANG:
00205         {
00206             SPL_FLOAT filtSize = (outFreqMax-outFreqMin)/(r->_nf+1)*2;  /*Anchura de cada filtro*/
00207             for (SPL_INT i = r->_minfreqidx; i <= r->_maxfreqidx; ++i)
00208             {
00209                 /*Calculamos la frecuencia escalada a la que corresponde ese bin*/
00210                 SPL_FLOAT outBin = FbandScaleConvert(i, XFBAND_SCALE_BIN, r->_scale, fres);
00211         
00212                 /*El peso y el filtro al que corresponde*/
00213                 SPL_FLOAT tmp = (outBin-outFreqMin)*2/filtSize;
00214                 SPL_FLOAT tmp2;
00215                 r->_fweight[i] = modf (tmp, &tmp2);
00216                 r->_findex[i] = (SPL_INT)tmp2;
00217         
00218                 /*Para Fs nos queda que le corresponde el filtro r->_nf+1, el cual no existe.
00219                     El problema es que está en la frontera entre dos filtros. Podemos ponerselo al anterior
00220                     con pesos modificados*/
00221                 if (r->_findex[i] == r->_nf+1) {
00222                     r->_findex[i] = r->_nf;
00223                     r->_fweight[i] = 1-r->_fweight[i];
00224                     /*peso 0 para el filtro r->_nf+1, peso 1 para r->_nf*/
00225                 }
00226         
00227                 /*Escalamos los pesos por el número de bins para terminar teniendo potencia independiente
00228                     de la NFFT (que cambiará con la Fs)*/
00229                 r->_fweight[i] /= fftnsamples;
00230             }
00231         
00232             break;
00233         }
00234     }
00235     
00236     return SPL_TRUE;
00237 }
00238 
00239 
00282 PUBLIC pFBAND XAPI fband_construct( SPL_INT np, SPL_INT nband, XFBAND_FILTER filter,
00283         XFBAND_SCALE scale, SPL_FLOAT minfreq, SPL_FLOAT maxfreq, SPL_FLOAT fs,
00284         XFFT_WIN win, XFFT_WIN_FUNC winfunc )
00285 {
00286     assert(np>0);
00287     assert(nband>0);
00288     assert(minfreq>=0);
00289     assert(maxfreq>=0);
00290     assert(scale!=XFBAND_SCALE_BIN);    //No es una escala válida
00291     
00292     pFBAND r = NULL;
00293     r = (pFBAND)xfft_malloc(sizeof(FBAND));
00294     if (r!=NULL){
00295         
00296         /* Primero creamos el objeto FFT */
00297         /* El número de puntos de la FFT será el de defecto: Siguiente potencia de 2 de np */
00298         r->_fft = rfft_construct( np, 0, SPL_FALSE, 32768, win, winfunc );
00299         if (r->_fft == NULL) {
00300             xfft_free(r);
00301             r=NULL;
00302         } else {
00303             
00304             r->_nf = nband;
00305             r->_filt = filter;
00306             r->_scale = scale;
00307             r->_srate = fs;
00308             
00309             /* Preparamos los punteros para fband_calcfilters */
00310             r->_findex = NULL;
00311             r->_fweight = NULL;
00312             r->_pvec = NULL;
00313             
00314             SPL_FLOAT fres = fs/rfft_getfftnp(r->_fft);
00315             maxfreq = maxfreq == 0 ? fs/2 : maxfreq;    /* Si no hay límite superior, se toma Nyquist */
00316             
00317             /* Índice del bin FFT de la minfreq */
00318             r->_minfreqidx = (SPL_INT) FbandScaleConvert(minfreq, XFBAND_SCALE_LIN, XFBAND_SCALE_BIN, fres);
00319             //  if (r->_minfreqidx < 1) r->_minfreqidx = 1;
00320     
00321             /* Índice del bin FFT de la maxfreq */
00322             r->_maxfreqidx = (SPL_INT) FbandScaleConvert(maxfreq, XFBAND_SCALE_LIN, XFBAND_SCALE_BIN, fres);
00323             //  if (r->_maxfreqidx >=  fftnsamples/2) r->_minfreqidx = fftnsamples/2;
00324             
00325             /*Calculamos los filtros*/
00326             if (fband_calcfilters( r ) == SPL_FALSE){
00327                 fband_destruct(r);
00328                 r = NULL;
00329             }
00330         }
00331     }
00332     
00333     return r;
00334 }
00335 
00346 PUBLIC SPL_VOID XAPI fband_destruct( pFBAND r )
00347 {
00348     assert(r!=NULL);
00349 
00350     if (r->_fft)
00351         rfft_destruct(r->_fft);
00352     if (r->_findex)
00353         xfft_free(r->_findex);
00354     if (r->_fweight)
00355         xfft_free(r->_fweight);
00356     if (r->_pvec)
00357         xfft_free(r->_pvec);
00358     xfft_free(r);
00359 }
00360 
00388 PUBLIC SPL_BOOL XAPI fband_setwin( pFBAND r, XFFT_WIN win,
00389         XFFT_WIN_FUNC winfunc )
00390 {
00391     assert(r!=NULL);
00392     
00393     return rfft_setwin(r->_fft, win, winfunc);
00394 }
00395 
00410 PUBLIC SPL_BOOL XAPI fband_setnp( pFBAND r, SPL_INT np )
00411 {
00412     assert(r!=NULL);
00413     assert(np>0);
00414     
00415     SPL_INT n = r->_fft->_np;
00416     if (np!=n) {
00417         if (rfft_setnp(r->_fft, np) == SPL_FALSE)
00418             return SPL_FALSE;
00419         if (r->_fft->_np != n) {    /*Tal vez no haya cambiado después de todo*/
00420             if (fband_calcfilters (r) == SPL_FALSE)
00421             return SPL_FALSE;
00422         }
00423     }
00424     return SPL_TRUE;
00425 }
00426 
00427 
00440 PUBLIC SPL_BOOL XAPI fband_setnband( pFBAND r, SPL_INT nband )
00441 {
00442     assert(r!=NULL);
00443     assert(nband>0);
00444     
00445     if (nband!=r->_nf) {
00446         r->_nf = nband;
00447         return fband_calcfilters (r);
00448     }
00449     return SPL_TRUE;
00450 }
00451 
00464 PUBLIC SPL_BOOL XAPI fband_setfilters( pFBAND r, XFBAND_FILTER filter, XFBAND_SCALE scale )
00465 {
00466     assert(r!=NULL);
00467     assert(scale!=XFBAND_SCALE_BIN);    //No es una escala válida
00468     
00469     r->_filt = filter;
00470     r->_scale = scale;
00471     return fband_calcfilters (r);
00472 }
00473 
00486 PUBLIC SPL_BOOL XAPI fband_setfreqlimits( pFBAND r, SPL_FLOAT minfreq, SPL_FLOAT maxfreq )
00487 {
00488     assert(r!=NULL);
00489     assert(minfreq>=0);
00490     assert(maxfreq>=0);
00491     
00492     SPL_FLOAT fres = r->_srate/rfft_getfftnp(r->_fft);
00493     maxfreq = maxfreq == 0 ? r->_srate/2 : maxfreq; /* Si no hay límite superior, se toma Nyquist */
00494 
00495     /* Por si acaso los bins no cambian, en cuyo caso nada cambia*/
00496     SPL_INT oldmin = r->_minfreqidx;
00497     SPL_INT oldmax = r->_maxfreqidx;
00498     
00499     /* Índice del bin FFT de la minfreq */
00500     r->_minfreqidx = (SPL_INT) FbandScaleConvert(minfreq, XFBAND_SCALE_LIN, XFBAND_SCALE_BIN, fres);
00501     //  if (r->_minfreqidx < 1) r->_minfreqidx = 1;
00502 
00503     /* Índice del bin FFT de la maxfreq */
00504     r->_maxfreqidx = (SPL_INT) FbandScaleConvert(maxfreq, XFBAND_SCALE_LIN, XFBAND_SCALE_BIN, fres);
00505     //  if (r->_maxfreqidx >=  fftnsamples/2) r->_minfreqidx = fftnsamples/2;
00506 
00507     /*Calculamos los filtros sólo si algo cambia*/
00508     if (oldmin != r->_minfreqidx || oldmax != r->_maxfreqidx)
00509         return fband_calcfilters (r);
00510     else
00511         return SPL_TRUE;
00512 }
00513 
00524 PUBLIC SPL_BOOL XAPI fband_setminfreq( pFBAND r, SPL_FLOAT minfreq )
00525 {
00526     assert(r!=NULL);
00527     assert(minfreq>=0);
00528     
00529     return fband_setfreqlimits (r, minfreq, FbandScaleConvert(r->_maxfreqidx, XFBAND_SCALE_BIN, XFBAND_SCALE_LIN, r->_srate/rfft_getfftnp(r->_fft)));
00530 }
00531 
00542 PUBLIC SPL_BOOL XAPI fband_setmaxfreq( pFBAND r, SPL_FLOAT maxfreq )
00543 {
00544     assert(r!=NULL);
00545     assert(maxfreq>=0);
00546     
00547     return fband_setfreqlimits (r, FbandScaleConvert(r->_minfreqidx, XFBAND_SCALE_BIN, XFBAND_SCALE_LIN, r->_srate/rfft_getfftnp(r->_fft)), maxfreq);
00548 }
00549 
00550 
00561 PUBLIC SPL_BOOL XAPI fband_setsrate( pFBAND r, SPL_FLOAT srate )
00562 {
00563     assert(r!=NULL);
00564     assert(srate>0);
00565 
00566     /*Calculamos los filtros sólo si algo cambia*/
00567     if (srate == r->_srate) 
00568         return SPL_TRUE;
00569     
00570     r->_srate = srate;
00571     return fband_calcfilters (r);
00572 }
00573 
00574 
00585 PUBLIC SPL_INT XAPI fband_getnp( pFBAND r )
00586 {
00587     assert(r!=NULL);
00588 
00589     return rfft_getnp(r->_fft);
00590 }
00591 
00602 PUBLIC SPL_INT XAPI fband_getnband( pFBAND r )
00603 {
00604     assert(r!=NULL);
00605 
00606     return (r->_nf);
00607 }
00608 
00609 
00617 PUBLIC XFFT_WIN XAPI fband_getwin( pFBAND r )
00618 {
00619     assert(r!=NULL);
00620 
00621     return rfft_getwin(r->_fft);
00622 }
00623 
00635 PUBLIC SPL_pFLOAT XAPI fband_getwinvec( pFBAND r )
00636 {
00637     assert(r!=NULL);
00638 
00639     return rfft_getwinvec(r->_fft);
00640 }
00641 
00649 PUBLIC XFFT_WIN_FUNC XAPI fband_getwinfunc( pFBAND r )
00650 {
00651     assert(r!=NULL);
00652 
00653     return rfft_getwinfunc(r->_fft);
00654 }
00655 
00656 
00670 PUBLIC SPL_pFLOAT XAPI fband_getpband( pFBAND r )
00671 {
00672     assert(r!=NULL);
00673 
00674     return r->_pvec;
00675 }
00676 
00692 PUBLIC SPL_FLOAT XAPI fband_getminfreq( pFBAND r )
00693 {
00694     assert(r!=NULL);
00695 
00696     return FbandScaleConvert(r->_minfreqidx, XFBAND_SCALE_BIN, XFBAND_SCALE_LIN, r->_srate/rfft_getfftnp(r->_fft));
00697 }
00698 
00714 PUBLIC SPL_FLOAT XAPI fband_getmaxfreq( pFBAND r )
00715 {
00716     assert(r!=NULL);
00717 
00718     return FbandScaleConvert(r->_maxfreqidx, XFBAND_SCALE_BIN, XFBAND_SCALE_LIN, r->_srate/rfft_getfftnp(r->_fft));
00719 }
00720 
00721 
00729 PUBLIC SPL_FLOAT XAPI fband_getsrate( pFBAND r )
00730 {
00731     assert(r!=NULL);
00732 
00733     return r->_srate;
00734 }
00735 
00736 
00740 /* Realiza el filtrado por bandas. Supone que la FFT ya está calculada
00741    y en formato módulo/argumento
00742 */
00743 PUBLIC SPL_VOID XAPI fband_dofband( pFBAND r )
00744 {
00745     assert(r!=NULL);
00746     
00747     /*Recuperamos el módulo de la FFT*/
00748     SPL_FLOAT* fftMod = rfft_getrevec(r->_fft);
00749     
00750     /*Borramos todo el contenido del vector de salida*/
00751     for (SPL_INT i = 0; i < r->_nf; ++i)
00752         r->_pvec[i] = 0;
00753 
00754     /*Recorremos la FFT, metiendo su aportación dentro de cada filtro
00755       Escalamos al número de bins (para tener potencia independientemente de la Fs)*/
00756     for (SPL_INT i = r->_minfreqidx; i <= r->_maxfreqidx; ++i)
00757     {
00758         /*El primer y último filtro son especiales*/
00759         if (r->_findex[i] != r->_nf)
00760             /*Para filtros rectangulares ya nos hemos asegurado de que no sea mayor que r->nf-1
00761               luego siempre entraremos aquí.*/
00762             r->_pvec[r->_findex[i]] += r->_fweight[i]*fftMod[i];
00763 
00764         if (r->_findex[i] != 0)
00765             /*Gracias a que en los filtros rectangulares el peso siempre es 1/fftNsamples,
00766               el método de calcular la aportación al filtro anterior es consistente para
00767               los filtros triangulares y rectangulares.*/
00768             r->_pvec[r->_findex[i]-1] += (1.0/rfft_getfftnp(r->_fft) - r->_fweight[i])*fftMod[i];
00769     }
00770 
00771     /*Aplicamos el log*/
00772     for (SPL_INT i = 0; i < r->_nf; ++i)
00773         r->_pvec[i] = 20*log10(r->_pvec[i]);    /*20 para que sea potencia*/
00774 }
00775 
00776 
00792 PUBLIC SPL_VOID XAPI fband_fband( pFBAND r, SPL_pFLOAT invec )
00793 {
00794     assert(r!=NULL);
00795     
00796     /*Calculamos su FFT en formato mod-arg*/
00797     rfft_rfft(r->_fft, invec);
00798     rfft_magarg(r->_fft);
00799     
00800     fband_dofband(r);
00801 }
00802 
00803 
00814 PUBLIC SPL_VOID XAPI fband_fband_i( pFBAND r, SPL_pINT invec )
00815 {
00816     assert(r!=NULL);
00817 
00818     /*Calculamos su FFT en formato mod-arg*/
00819     rfft_rfft_i(r->_fft, invec);
00820     rfft_magarg(r->_fft);
00821     
00822     fband_dofband(r);
00823 }
00824 
00835 PUBLIC SPL_VOID XAPI fband_fband_i16( pFBAND r, pINT16 invec )
00836 {
00837     assert(r!=NULL);
00838 
00839     /*Calculamos su FFT en formato mod-arg*/
00840     rfft_rfft_i16(r->_fft, invec);
00841     rfft_magarg(r->_fft);
00842     
00843     fband_dofband(r);
00844 }
00845 
00846 
00857 PUBLIC SPL_VOID XAPI fband_fband_i32( pFBAND r, pINT32 invec )
00858 {
00859     assert(r!=NULL);
00860 
00861     /*Calculamos su FFT en formato mod-arg*/
00862     rfft_rfft_i32(r->_fft, invec);
00863     rfft_magarg(r->_fft);
00864     
00865     fband_dofband(r);
00866 }
00867 
00878 PUBLIC SPL_VOID XAPI fband_fband_u32( pFBAND r, pUINT32 invec )
00879 {
00880     assert(r!=NULL);
00881 
00882     /*Calculamos su FFT en formato mod-arg*/
00883     rfft_rfft_u32(r->_fft, invec);
00884     rfft_magarg(r->_fft);
00885     
00886     fband_dofband(r);
00887 }
00888 
00899 PUBLIC SPL_VOID XAPI fband_fband_f32( pFBAND r, pFLOAT32 invec )
00900 {
00901     assert(r!=NULL);
00902 
00903     /*Calculamos su FFT en formato mod-arg*/
00904     rfft_rfft_f32(r->_fft, invec);
00905     rfft_magarg(r->_fft);
00906     
00907     fband_dofband(r);
00908 }
00909 
00910 
00930 PUBLIC SPL_VOID API fband2mfcc (pFBAND r, SPL_pFLOAT mfcc, SPL_INT nceps, SPL_INT lift)
00931 {
00932     assert(r!=NULL);
00933     assert(nceps>0);
00934 
00935     SPL_FLOAT c = sqrt(2.0/fband_getnband(r));
00936     for (SPL_INT j = 0; j < nceps; ++j)
00937     {
00938         mfcc[j] = 0;
00939         SPL_FLOAT cc = M_PI*(j+1)/fband_getnband(r);
00940         SPL_pFLOAT fb = fband_getpband(r);
00941         for (SPL_INT k = 0; k < fband_getnband(r); ++k)
00942             mfcc[j] += fb[k]*cos(cc*(k+0.5));
00943     
00944             //Aplicamos el liftering, si se pide
00945         if (lift != 0)
00946             mfcc[j] *= c*(1+lift/2*sin(M_PI*(j+1)/lift));
00947         else
00948             mfcc[j] *= c;
00949     }
00950 }
00951 
00952 
00953 /**********************************************************/
00954 
00955