fe_plcc.cpp

这是一个语音特征提取的程序源码

			
			fsq = f0 * f0;
			r_2 = (float)(fsq + 1.6e5);
			rsss = (float)(fsq * fsq * (fsq + 1.44e6) / (r_2 * r_2 * (fsq + 9.61e6)));
			
			cb[icount] = rsss * cb[icount];
			/* add the qual-loundness curve to the weighting    */
			++icount;
		}
    }
    return (0);
} 


/*
* int cosf_(m, nfilt, wcos)
* computes the cosine weightings for the Inverse Discrete Fourier Transfrom
* m    (in): order of plp auto regressive model
* nfilt(in):
* wcos (in): cosine weightings of inverse dft
*/

int Fe::cosf_ ( int m, int nfilt, float *wcos )
{
	int ii, jj;                       /* Loop counter variables     */
	
	wcos -= 24;                              /* Parameter adjustments      */
	/* [23][16] --> [j=1+24]      */
	
	for (ii = 1; ii <= (m+1); ++ii) {
		for (jj = 2; jj <= (nfilt-1); ++jj) {
			wcos[jj + ii * 23] = (float)(2.0 * cos(2.0 * M_PI * (ii - 1) * (jj - 1) / (2.0 * ( nfilt - 1 ))));
		}
		wcos[nfilt + ii * 23] = (float)(cos(2.0 * M_PI * (ii - 1) * (jj - 1) / (2.0 * ( nfilt - 1 ))));
	}
	
	return (0);
} 


/* 
* int a2gexp_ (a, gexp, i, nc, expon)
* a    (in): pointer to autoregressive coefficients array
* gexp(out):
* i    (in): number of elements in autoregressive coefficient array
* nc   (in): number of elements in output array
* expon(in):  
*/

int Fe::a2gexp_ ( float *a, float *gexp, int i, int nc, float expon)
{
    int i_2;      
	
    float c[257];                       /* Local variables          */
    int j, l, jb, ii, lp;
    float sum;
	
    --a;                                       /* Parameter adjustments    */
    --gexp;
	
    c[0] = (float)log(a[1]);
    c[1] = -a[2] / a[1];
    for (l = 2; l <= nc; ++l){
		lp = l + 1;
		if (l <= i + 1)
			sum = l * a[lp] / a[1];
		else
			sum = 0.0;
		i_2 = l;
		if(l < i + 1) i_2 = l;
		else i_2 = i + 1;
		for (j = 2; j <= i_2; ++j){
			jb = l - j + 2;
			sum += a[j] * c[jb - 1] * (jb - 1) / a[1];
		}
		c[lp - 1] = -sum / l;
	}
    gexp[1] = c[0];
	
    if (expon != 0.0)
		for (ii = 2; ii <= nc; ++ii)
			gexp[ii] = (float) pow((ii-1),expon ) * c[ii - 1];
		else
			for (ii = 2; ii <= nc; ++ii)
				gexp[ii] = c[ii - 1];
			
			return (0);
} 


/*
 * int TrigRecombFFT( cx, y, m )
 * Computers the FFT of a complex signal, using a trigonometric recombination
 * principle in order to speed up matters somewhat.
 * cx (in): pointer to input complex signal of type FeComplex
 * y (out): pointer to output FFT signal of type FeComplex
 * m  (in): the #point fft to computer (2^x format)
 * On exit returns (+1) for success and (-1) for failure
 * Usage: success = TrigRecombFFT( InSignal, FFTout, 8)
 */

int Fe::TrigRecombFFT (FeComplex *cx, FeComplex *y, int m )
{
	FeComplex *ck, *xk, *xnk;             /* array pointer variables      */
	float Realsum, Realdif;             /* real part of trig recombine  */
	float Imagsum, Imagdif;             /* imag part of trig recombine  */
	int num, k;                   /* loop counter variables       */

	/* Do FFT on combined data                  */
	num = 1 << m;
	PlpFFT( cx, m );

	/* Calculate Recombination Factors          */
	if( m != m_plpMofCF ) {
		if( CalculateCF( m ) < 0 ) {
			return( -1 );
		}
	}

	/* DC component, no multiplies              */
	y[0].m_re = cx[0].m_re + cx[0].m_im;
	y[0].m_im = 0.0;

	/* other frequencies by trig recombination  */
	ck = &m_plpCF[0];
	xk = cx + 1;
	xnk = cx + num - 1;

	for( k = 1; k < num; k++ ) {
		Realsum = ( xk->m_re + xnk->m_re ) / (float)2.0;
		Imagsum = ( xk->m_im + xnk->m_im ) / (float)2.0;
		Realdif = ( xk->m_re - xnk->m_re ) / (float)2.0;
		Imagdif = ( xk->m_im - xnk->m_im ) / (float)2.0;

		y[k].m_re = Realsum + ck->m_re * Imagsum -
					ck->m_im * Realdif;

		y[k].m_im = Imagdif - ck->m_im * Imagsum -
					ck->m_re * Realdif;

		ck++;
		xk++;
		xnk--;
	}

	return 1;
}


/*
 * FFT( x, m )
 * In-place radix 2 decimation in time FFT
 * Computes the Fast Fourier transform of a signal
 * x (in): pointer to complex array, power of 2 size of FFT
 * m (in): size fo fft = 2^m
 * Places FFT output on top of input FeComplex array.
 */

int Fe::PlpFFT( FeComplex *x, int m )
{
	FeComplex u, temp, tm;
	FeComplex *xi, *xip, *xj, *wptr;
	int i, j, k, l, le, windex;
	int n;

	/* Calculate Twiddle factors for FFT */
	if( m != m_plpMofW ) {
		if( CalculateW( m ) < 0 ) {
			return( -1 );
		}
	}

	/* n = 2^m = fft length   */
	n = 1 << m;

	/* start fft              */
	le = n;
	windex = 1;
	for( l = 0; l < m; l++ ) {
		le = le / 2;
		
		/* first iteration with no multiplies */
		for( i = 0; i < n; i = i + 2*le ) {
			xi = x + i;
			xip = xi + le;
			temp.m_re = xi->m_re + xip->m_re;
			temp.m_im = xi->m_im + xip->m_im;
			xip->m_re = xi->m_re - xip->m_re;
			xip->m_im = xi->m_im - xip->m_im;
			*xi = temp;
		}

		/* remaining iterations use to store w */
		wptr = &m_plpW[0] + windex - 1;
		for( j = 1; j < le; j++ ) {
			u = *wptr;
			for( i = j; i < n; i = i + 2*le ) {
				xi  = x + i;
				xip = xi + le;
				temp.m_re = xi->m_re + xip->m_re;
				temp.m_im = xi->m_im + xip->m_im;
				tm.m_re = xi->m_re - xip->m_re;
				tm.m_im = xi->m_im - xip->m_im;
				xip->m_re = tm.m_re * u.m_re - 
							tm.m_im * u.m_im;
				xip->m_im = tm.m_re * u.m_im +
							tm.m_im * u.m_re;
				*xi = temp;
			}
			wptr = wptr + windex;
		}
		windex = 2 * windex;
	}

	/* rearrage data by bit reversing */
	j = 0;
	for( i = 1; i < ( n - 1 ); i++ ) {
		k = n / 2;
		while( k <= j ) {
			j = j - k;
			k = k / 2;
		}
		j = j + k;
		if( i < j ) {
			xi = x + i;
			xj = x + j;
			temp = *xj;
			*xj = *xi;
			*xi = temp;
		}
	}

	return 1;
}


/* 
 * int CalculateW( m )
 * This procedure calculates the W twiddle factors for the 
 * FFT routine, thus needed only once. Then check to see if already
 * called, and do not call again.
 * m (in): number of points in FFT = (2^m)
 */

int Fe::CalculateW ( int m )
{
	double w_real, w_imag;
	m_plpMofW = m;
	int n = 1 << m;
	int le = n / 2;	
	m_plpW.resize(le - 1);
	double arg = 4.0 * atan( 1.0 ) / (float)le;
	double wrecur_real = w_real = cos( arg );
	double wrecur_imag = w_imag = -sin( arg );
	FeComplex *xj = &m_plpW[0];
	for(int j = 1; j < le; j++ ) {
		xj->m_re = (float)wrecur_real;
		xj->m_im = (float)wrecur_imag;
		xj++;
		double wtmp_real = wrecur_real * w_real - wrecur_imag * w_imag;
		wrecur_imag = wrecur_real * w_imag + wrecur_imag * w_real;
		wrecur_real = wtmp_real;
	}
	return 1;
}


/*
 * int CalculateCF( m )
 * Calculates the Trigonometric Recombination Coefficients.
 * Used only once, and never called again. With multiple FFT calls
 * saves computation time.
 * m (in): number of points in FFT transform = (2^m)
 */

int Fe::CalculateCF( int m )
{
	m_plpMofCF = m;
	int num = 1 << m;
	m_plpCF.resize(num - 1);
	double factor = 4.0 * atan( 1.0 ) / ((double)num);
	for(int k = 1; k < num; k++ ) {
		double arg = factor * ((double)k);
		m_plpCF[k-1].m_re = (float)cos( arg );
		m_plpCF[k-1].m_im = (float)sin( arg );
	}

	return 1;
}


/*****************************************************************************
*   rasta.c : RASTA processing by using Hermansky or SRI filter constants
*
*   11/30/1992 by Aki Ohshima ([e-mail])
*      Separated from mfcc_lib.c.
*   10/20/1992 by Aki Ohshima ([e-mail])
*      Created. RASTA processing by Nobu Hanai, modified by Aki.
****************************************************************************/

/* --------------------------------------------------------------------------
RastaFilter(mfsc, r_filter, r_f_param); 
 Where 'r_filter' is an integer and supposed to specify a kind of filter, 
 such as
 0 : no rasta filtering(just quit)
 1 : Hermansky's bandpass filetr
 2 : SRI's high pass filter
 'r_f_param' is a float and supposed to specify a filter parameter.
 0.94 : latest parameter for Hermansky's
 0.98 : previous parameter for Hermansky's
 0.97 : SRI's
 rasta filetering reduces the number of frames, and is taken care of in rasta 
 subroutine.
* ------------------------------------------------------------------------- */

/* 	History
* 	Coded by Nobutoshi Hanai      Sep 18, 1992
*    	Rasta filterling
*/

int Fe::RastaFilter(FeSpectrum *mfsc, int r_filter, float r_f_param)
{
	int		i,j;
	float		**in_dat;
	
	if (r_filter == RASTA_NO) return mfsc->m_ntime;

	/* initialize */
	in_dat = mfsc->m_specData;
	vector<float> rasta(mfsc->m_nfreq);
	vector<float> prerasta(mfsc->m_nfreq);
	
	for (j = 0; j < mfsc->m_nfreq; j++) {
		rasta[j] = 0.0;
		prerasta[j] = 0.0;
	}
	
	/* RASTA filtering */
	if (r_filter == RASTA_BY_HERMANSKY_BPF) {
	/*
	printf("RASTA filtering  -- Hermansky's\n");
		*/
		mfsc->m_ntime -= 4;
		for (i = 0; i < mfsc->m_ntime; i++){
			for (j = 0; j < mfsc->m_nfreq; ++j) {
				rasta[j] = (float)(r_f_param * prerasta[j] + 0.2 * *(in_dat[i+4] + j) + 0.1 * *(in_dat[i+3] + j) 
					- 0.1 * *(in_dat[i+1] + j) - 0.2 * *(in_dat[i] + j));
				*(in_dat[i] + j) = rasta[j];
				prerasta[j] = rasta[j];
			}
		}
	}
	else if (r_filter == RASTA_BY_SRI_HPF) {
	/*
	printf("RASTA filtering  -- SRI's\n");
		*/
		mfsc->m_ntime -= 1;
		for (i = 0; i < mfsc->m_ntime; i++){
			for (j = 0; j < mfsc->m_nfreq; ++j) {
				rasta[j] = r_f_param * prerasta[j] + *(in_dat[i+1] + j) - *(in_dat[i] + j);
				*(in_dat[i] + j) = rasta[j];
				prerasta[j] = rasta[j];
			}
		}
	}
	else {
		err_quit("ERROR : Invalid RASTA argument !!\n");
	}	
	return(mfsc->m_ntime);
}