_smcep.c

语音信号处理工具包,里面提供了很多有用的操作。这是最新版本

        
        for( i = 0, re = f, im = f + fftsz; i <= m2; i++ ){

            for( j = 0; j < fftsz; j++ )
	        *(re++) =  cos( ww[j] * i ) * dw[j];
            for( j = 0; j < fftsz; j++ )
	        *(im++) = -sin( ww[j] * i ) * dw[j];

	    re -= fftsz;
            im -= fftsz;

            ifft( re, im, fftsz );

            for( j = 1; j <= m1; j++ )
	      re[j] += re[fftsz-j];

            re += fftsz2;
            im += fftsz2;
        }

        free( ww );
        free( dw );


        /*-------  copy "f" to "g" ----------*/            

        for( i = 0, next = f, pf = f, pg = g; i <= m2; i++ ){
            for( j = 0; j <= m1; j++ )
                *(pg++) = *(pf++);
            next += fftsz2;
            pf = next;
        }
        free( f );
        flag_g = 1;

        for( j = 1; j <= m1; j++ )
            g[j] *= 0.5;

        for( i = 1; i <= m2; i++ )
            g[i*(m1+1)] *= 2.0;        
    }

    for( i = 0, pg = g; i <= m2; i++ )
        for( j = 0, c2[i] = 0.; j <= m1; j++ )
            c2[i] += *(pg++) * c1[j];
}



/***************************************************************

  No.3  ifreqt2    static : *h, size3

   Inverse Frequency Transformation

	void	ifreqt2(c1, m1, c2, m2, fftsz, a, t)

	double	*c1   : minimum phase sequence
	int	m1    : order of minimum phase sequence
	double	*c2   : warped sequence
	int	m2    : order of warped sequence
        int	fftsz : ifft size
	double	a     : all-pass constant
        double	t     : emphasized frequency t * pi(rad)

***************************************************************/

void ifreqt2(c1, m1, c2, m2, fftsz, a, t)
double *c1, *c2, a, t;
int m1, m2, fftsz;
{
    register int 	i, j;
    double		w, b, *ww, *f,
                        *re, *im, *pl, *pr, *plnxt, *prnxt,
                        *pf, *ph, *next, warp(), derivw();
    int                 size_h, size_f, fftsz2, m12, m11;
    static double	*h = NULL;
    static int		size3, flag_h = 1;

    b = M_2PI / (double)fftsz;

    size_h = ( m2 + 1 ) * ( m1 + 1 );

    if( h == NULL ){
        flag_h = 0;
        size3 = size_h;
        h = dgetmem( size3 );
    }
    else if( size_h != size3 ){
        free(h);
        flag_h = 0;        
        size3 = size_h;
        h = dgetmem( size3 );
    }

    /*-------  if "h" is not defined  ----------*/    

    if( flag_h == 0 ){

        ww = dgetmem( fftsz );

        for( j = 0, w = 0.; j < fftsz; j++, w+=b )
            ww[j] = warp( w, a, t );

        fftsz2 = fftsz + fftsz;       /* size of (re + im) */

	m12 = m1 + m1 + 1;
        size_f = m12 * fftsz2;      /* size of array "f" */
        f = dgetmem( size_f );
        

        for( i = -m1, re = f, im = f + fftsz; i <= m1; i++ ){

	  for( j = 0; j < fftsz; j++ )
	    *(re++) =  cos( ww[j] * i );
	  
	  for( j = 0; j < fftsz; j++ )
	    *(im++) = -sin( ww[j] * i );

	  re -= fftsz;
	  im -= fftsz;

	  ifft( re, im, fftsz );

	  re += fftsz2;
	  im += fftsz2;
        }

        free( ww );

        /*------- b'(n,m)=b(n,m)+b(n,-m) ----------*/            

	pl = f;
	pr = f + ( m12 - 1 ) * fftsz2;

	for( i = 0, plnxt = pl, prnxt = pr; i < m1; i++ ){
	  plnxt += fftsz2;
	  prnxt -= fftsz2;

	  for( j = 0; j <= m2; j++ )
	    *(pr++) += *(pl++);

	  pl = plnxt;
	  pr = prnxt;
	}

        /*-------  copy "f" to "h" ----------*/            

	m11 = m1 + 1;
	pf = f + m1 * fftsz2;

	for( j = 0, next = pf;  j <= m1; j++ ){

	  next += fftsz2;

	  for( i = 0; i <= m2; i++ )
	    h[m11*i+j] = *(pf++);

	  pf = next;
	}
        free( f );
        flag_h = 1;

        for( j = 1; j <= m1; j++ )
            h[j] *= 0.5;

        for( i = 1; i <= m2; i++ )
            h[i*m11] *= 2.0;
    }

    for( i = 0, ph = h; i <= m2; i++ )
      for( j = 0, c2[i] = 0.; j <= m1; j++ )
	c2[i] += *(ph++) * c1[j];
}



/***************************************************************

  No.4  frqtr2    static : *k, size4

    Frequency Transformation for Calculating Coefficients

	void	frqtr2(c1, m1, c2, m2, fftsz, a, t)

	double	*c1   : minimum phase sequence
	int	m1    : order of minimum phase sequence
	double	*c2   : warped sequence
	int	m2    : order of warped sequence
	int	fftsz  : frame length (fft size)
	double	a     : all-pass constant
	double	t     : emphasized frequency

***************************************************************/

void frqtr2(c1, m1, c2, m2, fftsz, a, t)
double *c1, *c2, a, t;
int m1, m2, fftsz;
{
    register int 	i, j;
    double		w, b, *ww, *f, *tc2,
                        *re, *im, *pf, *pk, *next, warp();
    int                 size_k, size_f, fftsz2;
    static double	*k = NULL;
    static int		size4, flag_k = 1;

    b = M_2PI / (double)fftsz;

    size_k = ( m2 + 1 ) * ( m1 + 1 );

    if( k == NULL ){
        flag_k = 0;
        size4 = size_k;
        k = dgetmem( size4 );
    }
    else if( size_k != size4 ){
        free(k);
        flag_k = 0;        
        size4 = size_k;
        k = dgetmem( size4 );
    }

    /*-------  if "k" is not defined  ----------*/    

    if( flag_k == 0 ){

        ww = dgetmem( fftsz );

        for( j = 0, w = 0.; j < fftsz; j++, w+=b )
            ww[j] = warp( w, a, t );

        fftsz2 = fftsz + fftsz;         /* size of (re + im) */
        size_f = ( m2 + 1 ) * fftsz2; /* size of array "f" */
        f = dgetmem( size_f );
        
        for( i = 0, re = f, im = f + fftsz; i <= m2; i++ ){

            for( j = 0; j < fftsz; j++ )
	      *(re++) =  cos( ww[j] * i );
            for( j = 0; j < fftsz; j++ )
	      *(im++) = -sin( ww[j] * i );

	    re -= fftsz;
            im -= fftsz;

            ifft( re, im, fftsz );

            for( j = 1; j <= m1; j++ )
	      re[j] += re[fftsz-j];
	    
            re += fftsz2;
            im += fftsz2;
        }

        free( ww );


        /*-------  copy "f" to "k" ----------*/            

        for( i = 0, next = f, pf = f, pk = k; i <= m2; i++ ){
            for( j = 0; j <= m1; j++ )
                *(pk++) = *(pf++);
            next += fftsz2;
            pf = next;
        }
        free( f );
        flag_k = 1;
    }

    tc2 = dgetmem( m2 + 1 );   /*  tmp of c2  */
    
    for( i = 0, pk = k; i <= m2; i++ )
        for( j = 0, tc2[i] = 0.; j <= m1; j++ )
            tc2[i] += *(pk++) * c1[j];

    movem(tc2, c2, sizeof(*c2), m2+1);

    free(tc2);
}


/***************************************************************

  Warping Function and Its Derivative
  
  double   warp(w, a, t)  &  derivw(w, a, t)

  double  w     : frequency
  double  a     : all-pass constant
  double  t     : emphasized frequency
  
***************************************************************/

double warp(w, a, t)
double w, a, t;
{
    double              ww, x, y;

    x = w-t;
    y = w+t;
    
    ww = w + atan2( (a * sin(x)), (1. - a * cos(x)) )
           + atan2( (a * sin(y)), (1. - a * cos(y)) );
    return(ww);
}


/*============================================================*/

double derivw(w, a, t)
double w, a, t;
{
    double dw, x, y, a2, aa;
 
    x = w-t;
    y = w+t;

    a2 = a+a;
    aa = a*a;

    dw = 1. + ( a * cos(x) - aa )/( 1. - a2 * cos(x) + aa )
            + ( a * cos(y) - aa )/( 1. - a2 * cos(y) + aa );
    return(dw);
}