fe_basic.cpp

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

///////////////////////////////////////////////////////////////////////////////
// This is a part of the Feature program.
// Version: 1.0
// Date: February 22, 2003
// Programmer: Oh-Wook Kwon
// Copyright(c) 2003 Oh-Wook Kwon. All rights reserved. [e-mail]
///////////////////////////////////////////////////////////////////////////////

#include "StdAfx.h"
#include "FE_feature.h"


typedef CPolynomial<float> Polynomial;
typedef Complex<float> CComplex;


int Fe::fftcep_spectrum_basic(short *sample, int frameSize, float *spectrum, int fftSize, int cepFilterLen)
{
	_fft_spectrum_basic(sample, frameSize, spectrum, fftSize, 1, cepFilterLen);
	return 1;
}


int Fe::lpc_basic(short *sample, int frameSize, float *acf, int norder, float *G)
{
	vector<float> frameA(frameSize);
	preprocessing(sample, frameSize, &frameA[0]);
	m_window.Windowing(&frameA[0], frameSize, WIN_HAMMING);
	_lpc_basic(&frameA[0], frameSize, acf, norder, G);
	return(1);
}


int Fe::lpc_cov_basic(short *sample, int frameSize, float *acf, int norder, float *G)
{
	vector<float> frameA(frameSize);
	for(int i=0;i<frameSize;i++) frameA[i]=sample[i];
	FeMatrix<float> cov(norder+1, norder+1);
	lpc_cov(&frameA[0], frameSize, cov, acf, norder, G);
	return(1);
}


int Fe::parcor_basic(short *sample, int frameSize, float *kcf, int norder)
{
	float G;
	vector<float> acf(norder+1);
	vector<float> frameA(frameSize);
	preprocessing(sample, frameSize, &frameA[0]);
	m_window.Windowing(&frameA[0], frameSize, WIN_HAMMING);
	_lpc_parcor_basic(&frameA[0], frameSize, &acf[0], kcf, norder, &G);
	return(1);
}


int	Fe::lar_basic(short *sample, int frameSize, float *lar, int norder)
{
	float G;
	vector<float> acf(norder+1);
	vector<float> kcf(norder+1);
	vector<float> frameA(frameSize);
	preprocessing(sample, frameSize, &frameA[0]);
	m_window.Windowing(&frameA[0], frameSize, WIN_HAMMING);
	_lpc_parcor_basic(&frameA[0], frameSize, &acf[0], &kcf[0], norder, &G);
	lar[0] = 0;
	for(int i=1; i<= norder; i++){
		lar[i] = (float)log((1-kcf[i])/(1+kcf[i]));
	}
	return(1);
}


int	Fe::lpc_spectrum_basic(short *sample, int frameSize, int norder, float *spectrum, int fftSize)
{
	int i;
	float G;
	vector<float> acf(fftSize);
	lpc_basic(sample, frameSize, &acf[0], norder, &G);
	int log2fftsize = (int)(log(fftSize)/log(2)+0.5);
	float logGdb = (G>0 ? 20*LOG10(G) : FE_MIN_LOG_ENERGY);	/* power spectrum */
	acf[0] = 1;
	/* Note that the definition of a[] is y(n) = a(1)y(n-1) + a(2)y(n-2) + ... */
	for(i=1;i<=norder;i++) acf[i] = -acf[i];
	
	compute_spectrum(&acf[0], spectrum, norder+1, log2fftsize);
	
	/* Note that spectrum is in the power domain. */
	for(i=0;i<fftSize;i++){
		float db=(spectrum[i]>FE_MIN_ENERGY ? 10*LOG10(spectrum[i]) : 10*LOG10(FE_MIN_ENERGY));
		spectrum[i] = (float)(logGdb-db);
	}
	
	return(1);
}


int	Fe::lpccep_spectrum_basic(short *sample, int frameSize, int ceporder, int norder, float *spectrum, int fftSize)
{
	int i;
	float G;
	vector<float> acf(fftSize);
	lpc_basic(sample, frameSize, &acf[0], norder, &G);
	vector<float> lpc_cep(fftSize);

	lpc_to_cepstrum(&acf[0], norder, &lpc_cep[0], ceporder, G);

	lpc_cep[0] = 0.5*lpc_cep[0]; /* we need log magnitude because fft_pow will give power spectrum */
	int log2fftsize = (int)(log(fftSize)/log(2)+0.5);
	compute_spectrum(&lpc_cep[0], spectrum, ceporder+1, log2fftsize);

	/* need 10/log(10) to get dB scale because lpc_cep is in the natural log scale. */
	float log10db=(float)(10/log(10));
	for(i=0;i<fftSize;i++){
		spectrum[i] = log10db*spectrum[i];
	}
	return(1);
}


int	Fe::melcep_spectrum_basic(short *sample, int frameSize, int ceporder, float *spectrum, int fftSize)
{
	int i,k;
	vector<float> mel_cep(ceporder+1);
	mel_cepstrum_basic(sample, frameSize, &mel_cep[0], ceporder, fftSize);
	MfccIDCT (&mel_cep[0], &m_idctMatrix[0], ceporder, m_fbOrder, spectrum);

	/* scale up to be consistent with other kinds of cepstrum coefficients */
	float g=2*fftSize/(float)m_fbOrder;
	for(i=0;i<m_fbOrder;i++) spectrum[i] *= g;

	float log10db=(float)(10/log(10));
	for(i=m_fbOrder;i>=1;i--){
		int kmin=(m_MelCenterIdx[i-1]+m_MelCenterIdx[i])/2;
		int kmax=(m_MelCenterIdx[i]+m_MelCenterIdx[i+1])/2;
		for(k=kmin;k<kmax;k++) spectrum[k]=log10db*spectrum[i-1];
	}
	{
		int kmax=(m_MelCenterIdx[0]+m_MelCenterIdx[1])/2;
		for(k=0;k<kmax;k++) spectrum[k]=spectrum[m_MelCenterIdx[1]];
		int kmin=(m_MelCenterIdx[m_fbOrder]+m_MelCenterIdx[m_fbOrder+1])/2;
		for(k=kmin;k<=fftSize/2;k++) spectrum[k]=spectrum[m_MelCenterIdx[m_fbOrder]];
	}
	return(1);
}


int	Fe::lpc_error_basic(short *sample, int frameSize, float *acf, int norder, float* residual)
{
	int i,n;
	for(n=0;n<norder && n<frameSize;n++){
		short si;
		residual[n] = sample[n];	
		for(i=1;i<=norder;i++){
			si = (n-i>=0) ? sample[n-i] : 0;
			residual[n] -= acf[i]*si;
		}
	}
	for(n=norder;n<frameSize;n++){
		residual[n] = sample[n];
		for(i=1;i<=norder;i++){
			residual[n] -= acf[i]*sample[n-i];
		}
	}
	return(1);
}


int	Fe::formant_basic(short *sample, int frameSize, float *formant, int formantN, int norder)
{
	float G;
	vector<float> acf(norder+1);
	lpc_basic(sample, frameSize, &acf[0], norder, &G);
	vector<CComplex> rootsA;
	int nf = lpc_to_formant(&acf[0], norder, formant, formantN, rootsA);
	for(int i=0; i<nf; i++)
		formant[i] *= (m_sampleRate/(float)(2*M_PI));
	return(1);
}


/*
sigma2 = r[0] - sum_{k=1}^N {a[k]*r[k]}
return value = G = sqrt(sigma2)
*/
float Fe::calc_lpc_gain_basic(float *r, float *acf, int norder)
{
	int i;
	float sigma2 = r[0];
	for(i=1; i<=norder; i++){
		sigma2 -= acf[i]*r[i];
	}
	return (float)sqrt(sigma2);
}


int	Fe::_lpc_parcor_basic(float *sample, int frameSize, float *acf, float *kcf, int norder, float *G)
{
	int		i, winSize;
	vector<float> r(norder+1);
	vector<float> rorg(norder+1);
	winSize = frameSize;
	{
		auto_correlation(sample, &r[0], frameSize, norder);
		
		for(i=0; i<=norder; i++) rorg[i] = r[i];
		
		/* solving the autocorrelation matrix */
		levins(&r[0], kcf, acf, norder);
		
		/* just for completeness */
		kcf[0] = 0;
		acf[0] = 1;
		
		/* G denotes LPC gain G */
		if(G) (*G) = calc_lpc_gain_basic(&rorg[0], acf, norder);
	}
	
	return(1);
}


/*
H(z) = G / (1 - sum_{i=1}^N {a[i]*z^{-i}} )
acf[0] = G, acf[i] = a[i], i=1,...,N
*/
int	Fe::_lpc_basic(float *sample, int frameSize, float *acf, int norder, float *G)
{
	int		i, winSize;
	vector<float> r(norder+1);
	vector<float> rorg(norder+1);
	vector<float> kcf(norder+1);
	winSize = frameSize;
	{
		auto_correlation(sample, &r[0], frameSize, norder);
		
		for(i=0; i<=norder; i++) rorg[i] = r[i];
		
		/* solving the autocorrelation matrix */
		levins(&r[0], &kcf[0], acf, norder);
		
		/* just for completeness */
		kcf[0] = 0; 
		acf[0] = 1;
		
		/* G denotes LPC gain G */
		if(G) (*G) = calc_lpc_gain_basic(&rorg[0], acf, norder);
	}
	
	return(1);
}


int	Fe::_lpc_error_basic(float *sample, int frameSize, float *acf, int norder, float* residual)
{
	int i,n;
	for(n=0;n<norder && n<frameSize;n++){
		float si;
		residual[n] = sample[n];	
		for(i=1;i<=norder;i++){
			si = (n-i>=0) ? sample[n-i] : 0;
			residual[n] -= acf[i]*si;
		}
	}
	for(n=norder;n<frameSize;n++){
		residual[n] = sample[n];
		for(i=1;i<=norder;i++){
			residual[n] -= acf[i]*sample[n-i];
		}
	}
	return(1);
}


/* preprocessing */
int Fe::preprocessing(short *sample, int sampleN, float *out)
{
	for(int i=0;i<sampleN;i++) out[i]=sample[i];
	if (m_dither) Dither(out, sampleN);
	preemphasize(out, sampleN, m_emphFac);
	return 1;
}


/* energy */
float Fe::compute_energy(float *sample, int sampleN, float mean)
{
	float sum = 0;
	for(int i=0; i<sampleN; i++){
		sum += (sample[i]-mean)*(sample[i]-mean);
	}
	float v=((sampleN>0) ? sum/sampleN : (float)0);
	v = v+1; /* to make happy feature extraction and analysis routines */
	return v;
}


/* zero crossing rate */
int Fe::compute_zero_cross_rate(float *sample, int sampleN, int level, int dc_bias)
{
	if(sampleN<1) return 0;
	int zcr=0;	
	float prev = sample[0]-dc_bias-level;
	for(int i=1; i<sampleN; i++) {
		float val  = sample[i]-dc_bias-level;
		float ztmp=val*prev;
		if(ztmp<0) zcr++;
		prev=val;
	}
	return zcr;
}


/* cepstral weighting */
/* cepstral coefficients are ordered as [c0, c1, c2, ..., c12] */
int	Fe::cepstral_window(float *cep, int ceporder, FeLifter lifter)
{
	int i;
	switch(lifter){
	case LIFT_NO:
		break;
	case LIFT_SIN:
		for (i=0; i<=ceporder; i++)
			cep[i] *= (float)(1. + (float)ceporder/2. * sin((M_PI*(i+1))/(float)ceporder));
		break;
	case LIFT_LINEAR:
		for (i=0; i<=ceporder; i++)
			cep[i] *= (i+1);
		break;
	case LIFT_SQRT:
		for (i=0; i<=ceporder; i++)
			cep[i] *= (float)sqrt(i+1);
		break;
	case LIFT_CUBE_ROOT:
		for (i=0; i<=ceporder; i++)
			cep[i] *= (float)pow(i+1,1/3.0);
		break;
	}
	return(1);
}


/* in-place preemphasis */
int Fe::preemphasize(float *sample, int sampleN, float emphFac)
{
	/* Setting emphFac=0 turns off preemphasis. */
	int	i;
	for (i = sampleN-1; i > 0; i--) {
		sample[i] = sample[i] - emphFac * sample[i-1];
	}
	sample[0] = (float)(1.0 - emphFac) * sample[0];
	return(1);
}