fe_polynomial.h

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

///////////////////////////////////////////////////////////////////////////////
// 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]
///////////////////////////////////////////////////////////////////////////////

#ifndef _FE_POLYNOMIAL_H_
#define _FE_POLYNOMIAL_H_

#include <cstdio>
#include <cstdarg>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <typeinfo>
#include <vector>
using namespace std;
#pragma warning(disable:4786)

#include "FE_complex.h"

template <typename T> class CPolynomial;

// friend template functions
template <typename T> CPolynomial<T> poly_plus(const CPolynomial<T>& A, const CPolynomial<T>& B);
template <typename T> CPolynomial<T> poly_plus(const CPolynomial<T>& A, const double b);
template <typename T> CPolynomial<T> poly_plus(const double b, const CPolynomial<T>& A);

template <typename T> CPolynomial<T> poly_minus(const CPolynomial<T>& A, const CPolynomial<T>& B);
template <typename T> CPolynomial<T> poly_minus(const CPolynomial<T>& A, const double b);
template <typename T> CPolynomial<T> poly_minus(const double b, const CPolynomial<T>& A);

template <typename T> CPolynomial<T> poly_times(const CPolynomial<T>& A, const CPolynomial<T>& B);
template <typename T> CPolynomial<T> poly_times(const CPolynomial<T>& A, const double b);
template <typename T> CPolynomial<T> poly_times(const double b, const CPolynomial<T>& A);

template <typename T> CPolynomial<T> poly_divide(const CPolynomial<T>& A, const CPolynomial<T>& B);
template <typename T> CPolynomial<T> poly_divide(const CPolynomial<T>& A, const double b);
template <typename T> CPolynomial<T> poly_divide(const double b, const CPolynomial<T>& A);

template <typename T> CPolynomial<T> poly_modulo(const CPolynomial<T>& A, const double b);
template <typename T> CPolynomial<T> poly_modulo(const CPolynomial<T>& A, const double b);
template <typename T> CPolynomial<T> poly_modulo(const double b, const CPolynomial<T>& A);

// GNU compiler gcc-2.95.2 cannot handle operator overloading with user-defined class type.
#if 0
template <typename T> CPolynomial<T> poly_plus(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
template <typename T> CPolynomial<T> poly_minus(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
template <typename T> CPolynomial<T> poly_times(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
template <typename T> CPolynomial<T> poly_divide(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
#endif


template<typename T> 
class CPolynomial
{
typedef Complex<T> CComplex;
public:
	CPolynomial();
	CPolynomial(const int v0) { v.resize(1); v[0]=(T)v0; }
	CPolynomial(const double v0) { v.resize(1); v[0]=(T)v0; }
	CPolynomial(int m, const T* a);
	CPolynomial(const char* a);
	CPolynomial(const CPolynomial<T>& a) : v(a.v) {}
	virtual ~CPolynomial();
	void Print() const;
	T Eval(T x, T* pd = NULL, int nd = 0) const;
	void Roots(CComplex* roots) const;
	void zroots(CComplex* a, int m, CComplex* roots, int polish) const;
	void laguer(CComplex* a, int m, CComplex* x, int *its) const;

public:
	vector<T> v;

	CPolynomial<T> operator+();
	CPolynomial<T> operator-();
	CPolynomial<T>& operator=(const CPolynomial<T>& a);
	CPolynomial<T>& operator+=(const CPolynomial<T>& a);
	CPolynomial<T>& operator-=(const CPolynomial<T>& a);
	CPolynomial<T>& operator*=(const CPolynomial<T>& a);
	CPolynomial<T>& operator/=(const CPolynomial<T>& a);

	// GCC-2.95-2 did not allow friend operator overloading definition outside class declaration.
        friend CPolynomial<T> operator+(const CPolynomial<T>& A, const CPolynomial<T>& B) {
                return poly_plus(A,B);
        }
        friend CPolynomial<T> operator+(const CPolynomial<T>& A, const double b) {
                return poly_plus(A,b);
        }
        friend CPolynomial<T> operator+(const double b, const CPolynomial<T>& A) {
                return poly_plus(b,A);
        }
        friend CPolynomial<T> operator-(const CPolynomial<T>& A, const CPolynomial<T>& B) {
                return poly_minus(A,B);
        }
        friend CPolynomial<T> operator-(const CPolynomial<T>& A, const double b) {
                return poly_minus(A,b);
        }
        friend CPolynomial<T> operator-(const double b, const CPolynomial<T>& A) {
                return poly_minus(b,A);
        }
        friend CPolynomial<T> operator*(const CPolynomial<T>& A, const CPolynomial<T>& B) {
                return poly_times(A,B);
        }
        friend CPolynomial<T> operator*(const CPolynomial<T>& A, const double b) {
                return poly_times(A,b);
        }
        friend CPolynomial<T> operator*(const double b, const CPolynomial<T>& A) {
                return poly_times(b,A);
        }
        friend CPolynomial<T> operator/(const CPolynomial<T>& A, const CPolynomial<T>& B) {
                return poly_divide(A,B);
        }
        friend CPolynomial<T> operator/(const CPolynomial<T>& A, const double b) {
                return poly_divide(A,b);
        }
        friend CPolynomial<T> operator/(const double b, const CPolynomial<T>& A) {
                return poly_divide(b,A);
        }
        friend CPolynomial<T> operator%(const CPolynomial<T>& A, const CPolynomial<T>& B) {
                return poly_modulo(A,B);
        }
        friend CPolynomial<T> operator%(const CPolynomial<T>& A, const double b) {
                return poly_modulo(A,b);
        }
        friend CPolynomial<T> operator%(const double b, const CPolynomial<T>& A) {
                return poly_modulo(b,A);
        }

#if 0
        friend CPolynomial<T> operator+(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B) {
                return poly_plus(A,B);
        }
        friend CPolynomial<T> operator+(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B) {
                return poly_minus(A,B);
        }
        friend CPolynomial<T> operator+(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B) {
                return poly_times(A,B);
        }
        friend CPolynomial<T> operator+(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B) {
                return poly_divide(A,B);
        }
#endif

#ifdef __GNUC__
        friend CPolynomial<T> poly_plus<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_plus<T>(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_plus<T>(const double b, const CPolynomial<T>& A);
        friend CPolynomial<T> poly_minus<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_minus<T>(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_minus<T>(const double b, const CPolynomial<T>& A);
        friend CPolynomial<T> poly_times<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_times<T>(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_times<T>(const double b, const CPolynomial<T>& A);
        friend CPolynomial<T> poly_divide<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_divide<T>(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_divide<T>(const double b, const CPolynomial<T>& A);
        friend CPolynomial<T> poly_modulo<T>(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_modulo<T>(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_modulo<T>(const double b, const CPolynomial<T>& A);
#if 0
	friend CPolynomial<T> poly_plus<T>(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
	friend CPolynomial<T> poly_minus<T>(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
	friend CPolynomial<T> poly_times<T>(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
	friend CPolynomial<T> poly_divide<T>(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
#endif
#else
	friend CPolynomial<T> poly_plus(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_plus(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_plus(const double b, const CPolynomial<T>& A);
        friend CPolynomial<T> poly_minus(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_minus(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_minus(const double b, const CPolynomial<T>& A);
        friend CPolynomial<T> poly_times(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_times(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_times(const double b, const CPolynomial<T>& A);
        friend CPolynomial<T> poly_divide(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_divide(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_divide(const double b, const CPolynomial<T>& A);
        friend CPolynomial<T> poly_modulo(const CPolynomial<T>& A, const CPolynomial<T>& B);
        friend CPolynomial<T> poly_modulo(const CPolynomial<T>& A, const double b);
        friend CPolynomial<T> poly_modulo(const double b, const CPolynomial<T>& A);
#if 0
	friend CPolynomial<T> poly_plus(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
	friend CPolynomial<T> poly_minus(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
	friend CPolynomial<T> poly_times(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
	friend CPolynomial<T> poly_divide(const CPolynomial<Complex<T> >& A, const CPolynomial<Complex<T> >& B);
#endif
#endif

#if 0
	friend CPolynomial<T> operator+(const CPolynomial<T>& a, const CPolynomial<T>& b);
	friend CPolynomial<T> operator+(const CPolynomial<T>& a, const double& b);
	friend CPolynomial<T> operator+(const double& a, const CPolynomial<T>& b);
	friend CPolynomial<T> operator+(const CPolynomial<T>& a);
	friend CPolynomial<T> operator-(const CPolynomial<T>& a, const CPolynomial<T>& b);
	friend CPolynomial<T> operator-(const CPolynomial<T>& a, const double& b);
	friend CPolynomial<T> operator-(const double& a, const CPolynomial<T>& b);
	friend CPolynomial<T> operator-(const CPolynomial<T>& a);
	friend CPolynomial<T> operator*(const CPolynomial<T>& a, const CPolynomial<T>& b);
	friend CPolynomial<T> operator*(const CPolynomial<T>& a, const double& b);
	friend CPolynomial<T> operator*(const double& a, const CPolynomial<T>& b);
	friend CPolynomial<T> operator/(const CPolynomial<T>& a, const CPolynomial<T>& b);
	friend CPolynomial<T> operator/(const CPolynomial<T>& a, const double& b);
	friend CPolynomial<T> operator/(const double& a, const CPolynomial<T>& b);
	friend CPolynomial<T> operator%(const CPolynomial<T>& a, const CPolynomial<T>& b);	
	friend CPolynomial<T> operator%(const CPolynomial<T>& a, const double& b);
	friend CPolynomial<T> operator%(const double& a, const CPolynomial<T>& b);
#endif

};

template<typename T>
inline void PolynomialDivision(const CPolynomial<T>& a, const CPolynomial<T>& b, CPolynomial<T>& q, CPolynomial<T>& r)
{
	int an = a.v.size() - 1;
	int bn = a.v.size() - 1;
	int qn = ((an >= bn) ? an-bn : 0);
	int rn = ((bn >= 1) ? ((an <bn-1) ? an : bn-1) : 0);
	q.v.resize(qn+1);
	r.v.resize(rn+1);

	int k,j;
	for(j=0;j<=an;j++){
		r.v[j] = a.v[j];
		q.v[j] = 0;
	}
	for(k=an-bn;k>=0;k--){
		q.v[k] = r.v[bn+k]/b.v[bn];
		for(j=bn+k-1;j>=k;j--) r.v[j] -= q.v[k]*b.v[j-k];
	}
	for(j=bn;j<=an;j++) r.v[j] = 0;
}

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

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

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

template<typename T>
inline CPolynomial<T> operator+(const CPolynomial<T>& a)
{
	CPolynomial<T> c = a;
	return c;
}

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

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

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

template<typename T>
inline CPolynomial<T> operator-(const CPolynomial<T>& a)
{
	CPolynomial<T> c = 0;
	c -= a;
	return c;
}

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

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

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

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

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