fe_complex.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_COMPLEX_H_
#define _FE_COMPLEX_H_

#include <cstdio>
#include <cmath>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

template <typename T> class Complex;

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

template <typename T> Complex<T> complex_minus(const Complex<T>& A, const Complex<T>& B);
template <typename T> Complex<T> complex_minus(const Complex<T>& A, const double b);
template <typename T> Complex<T> complex_minus(const double b, const Complex<T>& A);

template <typename T> Complex<T> complex_times(const Complex<T>& A, const Complex<T>& B);
template <typename T> Complex<T> complex_times(const Complex<T>& A, const double b);
template <typename T> Complex<T> complex_times(const double b, const Complex<T>& A);

template <typename T> Complex<T> complex_divide(const Complex<T>& A, const Complex<T>& B);
template <typename T> Complex<T> complex_divide(const Complex<T>& A, const double b);
template <typename T> Complex<T> complex_divide(const double b, const Complex<T>& A);


template <typename T>
class Complex
{
public:
	Complex(const T a=(T)0, const T b=(T)0) : r(a), i(b) {}
	Complex(const Complex<T>& z) : r(z.r), i(z.i) {}
	~Complex();

public:
	T r;
	T i;

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

	friend Complex<T> operator+(const Complex<T>& a, const Complex<T>& b) {
		return complex_plus(a,b);
	}
	friend Complex<T> operator+(const Complex<T>& a, const T& b) {
		return complex_plus(a,b);
	}
	friend Complex<T> operator+(const T& a, const Complex<T>& b) {
		return complex_plus(a,b);
	}
	friend Complex<T> operator-(const Complex<T>& a, const Complex<T>& b) {
		return complex_minus(a,b);
	}
	friend Complex<T> operator-(const Complex<T>& a, const T& b) {
		return complex_minus(a,b);
	}
	friend Complex<T> operator-(const T& a, const Complex<T>& b) {
		return complex_minus(a,b);
	}
	friend Complex<T> operator*(const Complex<T>& a, const Complex<T>& b) {
		return complex_times(a,b);
	}
	friend Complex<T> operator*(const Complex<T>& a, const T& b) {
		return complex_times(a,b);
	}
	friend Complex<T> operator*(const T& a, const Complex<T>& b) {
		return complex_times(a,b);
	}
	friend Complex<T> operator/(const Complex<T>& a, const Complex<T>& b) {
		return complex_divide(a,b);
	}
	friend Complex<T> operator/(const Complex<T>& a, const T& b) {
		return complex_divide(a,b);
	}
	friend Complex<T> operator/(const T& a, const Complex<T>& b) {
		return complex_divide(a,b);
	}

#ifdef __GNUC__
        friend Complex<T> complex_plus<T>(const Complex<T>& A, const Complex<T>& B);
        friend Complex<T> complex_plus<T>(const Complex<T>& A, const double b);
        friend Complex<T> complex_plus<T>(const double b, const Complex<T>& A);
        friend Complex<T> complex_minus<T>(const Complex<T>& A, const Complex<T>& B);
        friend Complex<T> complex_minus<T>(const Complex<T>& A, const double b);
        friend Complex<T> complex_minus<T>(const double b, const Complex<T>& A);
        friend Complex<T> complex_times<T>(const Complex<T>& A, const Complex<T>& B);
        friend Complex<T> complex_times<T>(const Complex<T>& A, const double b);
        friend Complex<T> complex_times<T>(const double b, const Complex<T>& A);
        friend Complex<T> complex_divide<T>(const Complex<T>& A, const Complex<T>& B);
        friend Complex<T> complex_divide<T>(const Complex<T>& A, const double b);
        friend Complex<T> complex_divide<T>(const double b, const Complex<T>& A);
#endif

};

template <typename T> 
inline Complex<T> Complex<T>::operator+() {
        return (*this);
}

template <typename T> 
inline Complex<T> Complex<T>::operator-() {
        Complex<T> C=(*this);
        C*=(-1);
        return C;
}

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

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

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

template <typename T>
inline Complex<T> complex_plus(const Complex<T>& a)
{
	Complex<T> c = a;
	return c;
}

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

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

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

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

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

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

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

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

template <typename T>
inline Complex<T> complex_divide(const Complex<T>& a, const T& b)
{
	Complex<T> c = a;
	c /= Complex<T>((T)b,(T)0);
	return c;
}

template <typename T>
inline Complex<T> complex_divide(const T& a, const Complex<T>& b)
{
	Complex<T> c((T)a, (T)0);
	c /= b;
	return c;
}

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

template <typename T>
T real(const Complex<T>& z)
{
	return z.r;
}

template <typename T>
T imag(const Complex<T>& z)
{
	return z.i;
}

template <typename T>
T abs(const Complex<T>& z)
{
	T x,y,ans,temp;
	x=(T)fabs(z.r);
	y=(T)fabs(z.i);
	if(x == (T)0)
		ans=y;
	else if(y == (T)0)
		ans=x;
	else if(x>y){
		temp=y/x;
		ans=(T)(x*sqrt((T)1.0+temp*temp));
	} else {
		temp = x/y;
		ans=(T)(y*sqrt((T)1.0+temp*temp));
	}
	return ans;
}

template <typename T>
T arg(const Complex<T>& z)
{
	if(z.r == (T)0 && z.i == (T)0) return (T)0;
	else if(z.i >= (T)0) return (T)acos(z.r/abs(z));
	else return (T)(2*M_PI-acos(z.r/abs(z)));
}

template <typename T>
Complex<T> sqrt(const Complex<T>& z)
{
	Complex<T> c;
	T x,y,w,r;
	if((z.r == (T)0) && (z.i == (T)0)){
		c.r = (T)0;
		c.i = (T)0;
		return c;
	} else {
		x=(T)fabs(z.r);
		y=(T)fabs(z.i);
		if(x>=y){
			r=y/x;
			w=(T)(sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))));
		} else {
			r=x/y;
			w=(T)(sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))));
		}
		if(z.r >= (T)0){
			c.r=w;
			c.i=z.i/(T)(2*w);
		} else {
			c.i=(z.i >= (T)0) ? w : -w;
			c.r=z.i/(T)(2.0*c.i);
		}
		return c;
	}
}

template <typename T>
Complex<T> conjg(const Complex<T>& z)
{
	Complex<T> c;
	c.r = z.r;
	c.i = -z.i;
	return c;
}

template <typename T>
Complex<T>& Complex<T>::operator=(const Complex<T>& a)
{
	r = a.r;
	i = a.i;
	return *this;
}

template <typename T>
Complex<T>& Complex<T>::operator+=(const Complex<T>& a)
{
	r += a.r;
	i += a.i;
	return *this;
}

template <typename T>
Complex<T>& Complex<T>::operator-=(const Complex<T>& a)
{
	r -= a.r;
	i -= a.i;
	return *this;
}

template <typename T>
Complex<T>& Complex<T>::operator*=(const Complex<T>& a)
{
	T t;
	r = (t=r)*a.r-i*a.i;
	i = i*a.r+t*a.i;
	return *this;
}

template <typename T>
Complex<T>& Complex<T>::operator/=(const Complex<T>& a)
{
	T ratio,den,t;
	if((T)fabs(a.r) >= (T)fabs(a.i)){
		ratio = a.i/a.r;
		den = a.r+ratio*a.i;
		r = ((t=r)+ratio*i)/den;
		i = (i-ratio*t)/den;
	}
	else{
		ratio = a.r/a.i;
		den = a.i+ratio*a.r;
		r = ((t=r)*ratio+i)/den;
		i = (i*ratio-t)/den;
	}

	return *this;
}

#endif