fe_polynomial.h

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

}

template<typename T>
inline CPolynomial<T> poly_divide(const double a, const CPolynomial<T>& b)
{
	CPolynomial<T> q,r;
	PolynomialDivision(CPolynomial<T>((T)a),b,q,r);
	return q;
}

template<typename T>
inline CPolynomial<T> poly_modulo(const CPolynomial<T>& a, const CPolynomial<T>& b)
{
	CPolynomial<T> q,r;
	PolynomialDivision(a,b,q,r);
	return r;
}

template<typename T>
inline CPolynomial<T> poly_modulo(const CPolynomial<T>& a, const double b)
{
	CPolynomial<T> q,r;
	PolynomialDivision(a,CPolynomial<T>((T)b),q,r);
	return r;
}

template<typename T>
inline CPolynomial<T> poly_modulo(const double a, const CPolynomial<T>& b)
{
	CPolynomial<T> q,r;
	PolynomialDivision(CPolynomial<T>((T)a),b,q,r);
	return r;
}


#if 0
template<typename T> inline CPolynomial<T> poly_plus(const CPolynomial<Complex<T> >& a, const CPolynomial<Complex<T> >& b)
{
	CPolynomial<T> c = a; c += b; return c;
}

template<typename T> inline CPolynomial<T> poly_minus(const CPolynomial<Complex<T> >& a, const CPolynomial<Complex<T> >& b)
{
	CPolynomial<T> c = a; c -= b; return c;
}

template<typename T> inline CPolynomial<T> poly_times(const CPolynomial<Complex<T> >& a, const CPolynomial<Complex<T> >& b)
{
	CPolynomial<T> c = a; c *= b; return c;
}

template<typename T> inline CPolynomial<T> poly_divide(const CPolynomial<Complex<T> >& a, const CPolynomial<Complex<T> >& b)
{
	CPolynomial<T> c = a; c /= b; return c;
}
#endif

template<typename T>
CPolynomial<T>::CPolynomial(int m, const T* a)
{
	for(int i=0; i<=m; i++){
		v.push_back(a[i]);
	}
}

#if 0
// do not use <cstdarg> and <typeinfo> because they are very fragile
// and depend on compilers
template<typename T>
CPolynomial<T>::CPolynomial(int m, const T v0, ...)
{
	va_list args;
	va_start (args, m); // Although GNU compiler warns here, we can ignore it.
	const type_info& a = typeid(v0);
	if(strcmp(a.name(),"float")==0 || strcmp(a.name(),"double")==0){
		for(int i=1; i<=m; i++){
			v.push_back((T)va_arg(args, double));
		}
	}
	else{
		printf("Unsupported Polynomial type (%s)!\n", a.name());
		exit(1);
	}
    va_end (args);
}
#endif


template<typename T>
CPolynomial<T>::CPolynomial(const char* a)
{
	char *t = strdup(a);
	char *p = strtok(t, " \t\r\n");
	v.push_back((T)atof(p));
	while((p = strtok(NULL, " \t\r\n"))){
		v.push_back((T)atof(p));
	}
	free(t);
}

template<typename T>
CPolynomial<T>::CPolynomial()
{
	v.push_back((T)0);
}

template<typename T>
CPolynomial<T>::~CPolynomial()
{
}

#define EPS 2.0e-6
#define MAXM 100

/* Roots() yields n roots into roots[1]...roots[n] */
template<typename T>
void CPolynomial<T>::Roots(CComplex* roots) const
{
	int n = v.size()-1;
	vector<CComplex> a(n+1);
	for(int i=0;i<=n;i++) a[i] = CComplex(v[i],0);
	zroots(&a[0], n, roots, 1);
}

#define EPSS 1.0e-7
#define MR 8
#define MT 10
#define MAXIT (MT*MR)

#ifndef my_max
#define my_max(a, b) (((a) > (b)) ? (a) : (b)) 
#endif
#ifndef my_min
#define my_min(a, b) (((a) < (b)) ? (a) : (b)) 
#endif


template<typename T>
void CPolynomial<T>::zroots(CComplex* a, int m, CComplex* roots, int polish) const
{
	int i,its,j,jj;
	CComplex x,b,c,ad[MAXM];

	for(j=0;j<=m;j++) ad[j] = a[j];
	for(j=m;j>=1;j--){
		//x = 0;
		x = roots[j];
		laguer(ad,j,&x,&its);
		if(fabs(imag(x)) <= 2.0*EPS*fabs(real(x))) x=CComplex(real(x), 0);
		roots[j] = x;
		b = ad[j];
		for(jj=j-1;jj>=0;jj--){
			c = ad[jj];
			ad[jj] = b;
			b = x*b+c;
		}
	}
	if(polish)
		for(j=1;j<=m;j++)
			laguer(a,m,&roots[j],&its);
	for(j=2;j<=m;j++){
		x = roots[j];
		for(i=j-1;i>=1;i--){
			if(real(roots[i]) <= real(x)) break;
			roots[i+1] = roots[i];
		}
		roots[i+1] = x;
	}
}

template<typename T>
void CPolynomial<T>::laguer(CComplex* a, int m, CComplex* x, int *its) const
{
	int iter,j;
	T abx,abp,abm,err;
	CComplex dx,x1,b,d,f,g,h,sq,gp,gm,g2;
	static T frac[MR+1]={
		(T)0.00,(T)0.50,(T)0.25,
		(T)0.75,(T)0.13,(T)0.38,
		(T)0.62,(T)0.88,(T)1.00
	};

	for(iter=1;iter<=MAXIT;iter++){
		*its = iter;
		b = a[m];
		err = (T)abs(b);
		d = f = 0;
		abx = (T)abs(*x);
		for(j=m-1;j>=0;j--){
			f = (*x)*f+d;
			d = (*x)*d+b;
			b = (*x)*b+a[j];
			err = (T)abs(b)+abx*err;
		}
		err *= (T)EPSS;
		if(abs(b) <= err) return;
		g = d/b;
		g2 = g*g;
		h = g2 - (T)2*(f/b);
		sq = sqrt((T)(m-1)*(((T)m*h)-g2));
		gp = g+sq;
		gm = g-sq;
		abp = (T)abs(gp);
		abm = (T)abs(gm);
		if(abp < abm) gp = gm;
		dx = ((my_max(abp,abm) > 0.0 ? CComplex(m,0)/gp
			: ((T)exp(log(1+abx))*CComplex(cos(iter),sin(iter)))));
		x1 = (*x)-dx;
		if(real(*x) == real(x1) && imag(*x) == imag(x1)) return;
		if(iter%MT) *x = x1;
		else *x = (*x)-frac[iter/MT]*dx;
	}
	printf("too many iteration in laguer");
	return;
}

template<typename T>
T CPolynomial<T>::Eval(T x, T* pd, int nd) const
{
	int nnd, j, i;
	T cnst = 1;
	T ftmp;

	if(!pd){
		nd = 0;
	}
	int n = v.size()-1;
	ftmp = v[n];
	if(pd) pd[0] = v[n];
	for(j=1;j<=nd;j++) pd[j] = 0;
	for(i=n-1;i>=0;i--){
		nnd = (nd < (n-i) ? nd : n-i);
		for(j=nnd;j>=1;j--)
			pd[j] = (T)(pd[j]*x+pd[j-1]);
		ftmp = (T)(ftmp*x+v[i]);
		if(pd) pd[0] = (T)(pd[0]*x+v[i]);
	}
	for(i=2;i<=nd;i++){
		cnst *= i;
		pd[i] *= cnst;
	}
	return ftmp;
}

template<typename T>
void CPolynomial<T>::Print() const
{
	int n=v.size()-1;
	int bAllZero = 1;
	for(int i=n;i>=0;i--){
		if(fabs(v[i]) < 1e-37){
			if(i==0 && bAllZero) printf("0");
			continue;
		}
		bAllZero = 0;
		if(v[i] > 0 && i != n) printf(" + ");
		else if(v[i] < 0) printf(" - ");

		if(fabs(fabs(v[i])-1) > 1e-37){
			printf("%gx(%d)",fabs(v[i]),i);
		}
		else{
			printf("x(%d)",i);	
		}
	}
	printf("\n");
}

template<typename T>
CPolynomial<T>& CPolynomial<T>::operator=(const CPolynomial<T>& a)
{
	int an = a.v.size()-1;
	int n = an;
	v.resize(n+1);
	int m = 0;
	for(int i=0; i<=n; i++){
		v[i] = a.v[i];
		if(v[i] != 0.0) 
			m = i;
	}
	if(m != n) v.resize(m+1);
	return *this;
}

template<typename T>
CPolynomial<T>& CPolynomial<T>::operator+=(const CPolynomial<T>& a)
{
	int m = v.size()-1;
	int an = a.v.size()-1;	
	int n = my_max(m,an);
	v.resize(n+1);
	for(int i=0; i<=n; i++)
		v[i] = ((i<=m)? v[i]: 0) + ((i<=an)? a.v[i]: 0);
	return *this;
}

template<typename T>
CPolynomial<T>& CPolynomial<T>::operator-=(const CPolynomial<T>& a)
{
	int m = v.size()-1;
	int an = a.v.size()-1;	
	int n = my_max(m,an);
	v.resize(n+1);
	for(int i=0; i<=n; i++)
		v[i] = ((i<=m)? v[i]: 0) - ((i<=an)? a.v[i]: 0);
	return *this;
}

template<typename T>
CPolynomial<T>& CPolynomial<T>::operator*=(const CPolynomial<T>& a)
{
	int p,i;
	int m = v.size()-1;
	int an = a.v.size()-1;
	int n = m + an;
	v.resize(n+1);
	for(p=n;p>m;p--) v[p] = 0;
	for(p=n; p>=0; p--){
		T sum = 0;
		for(i=my_min(p,m); i>=0 && p-i<=an; i--){
			sum += v[i] * a.v[p-i];
		}
		v[p] = sum;
	}
	return *this;
}

template<typename T>
CPolynomial<T>& CPolynomial<T>::operator/=(const CPolynomial<T>& a)
{
	CPolynomial<T> q,r;
	PolynomialDivision(*this,a,q,r);
	return q;
}

#endif